[Rcpp-devel] Comparison of R and rcpp for backsolve

[Matteo Fasiolo]

Hi Hao,

 from what I understand Armadillo doesn't have a forward or back-solver.

In any case it's pretty easy to create one, for example you can have a look
here:

https://github.com/mfasiolo/mvnfast/blob/master/src/mahaCpp.cpp

Have a nice weekend,

Matteo


On Sat, Jun 28, 2014 at 12:25 AM, Hao Ye <hye at ucsd.edu> wrote:

> Hi Peter,
>
> I think Soren had a similar issue in February (about Armadillo not having
> an optimized solver for triangular matrices):
>
> http://lists.r-forge.r-project.org/pipermail/rcpp-devel/2014-February/007147.html
>
> I don't think there was ever a reply though...
>
> In contrast, the notes for backsolve() show:
> > This is a wrapper for the level-3 BLAS routine dtrsm.
>
> which *is* specific for triangular matrices. So I'm not surprised that
> backsolve() is faster than RcppArmadillo's solve() in your situation.
>
> Best,
> --
> Hao Ye
> hye at ucsd.edu
>
> On Jun 27, 2014, at 4:14 PM, Peter Rossi <perossichi at gmail.com> wrote:
>
> > Folks-
> >
> > In my package, bayesm, I use backsolve() to invert upper-triangular
> > arrays.  I am in the process of converting my package to rccp-arma.
> >
> > I thought I would test the analogous operation in arma by declaring
> > the matrix as upper-triangular and inverting using solve().
> >
> > To my surprise, pure R code was considerably faster than
> > rcpp-armadillo for 50 x 50 and larger matrices.
> >
> > I attach an rmd and cpp files necessary to run this benchmark.
> >
> > I assume I'm doing something terribly wrong or naive in my use of
> > rcpp-armadillo and would appreciate any thoughts on what causes these
> > unexpected results.
> >
> > p
> >
> > Here are the results for those who do not want to compile the rmd file:
> >
> > 10 by 10 uppertriangular matrix;
> >
> > # 10 by 10
> > set.seed(777)
> > k <- 10
> > m <- k + 10
> > A <- matrix(rnorm(k*m), m, k)
> > R <- chol(t(A)%*%A)
> >
> > rep <- 500
> > microbenchmark(backsolve_R(rep, R),backsolve_rcpp(rep, R),times=10)
> >
> > ## Unit: microseconds
> > ##                    expr    min     lq median     uq    max neval
> > ##     backsolve_R(rep, R) 3781.3 4002.4 4131.0 4958.8 6146.3    10
> > ##  backsolve_rcpp(rep, R)  809.1  812.5  825.1  851.4  889.7    10
> >
> > 50 by 50 uppertriangular matrix
> >
> > # 50 by 50
> > set.seed(777)
> > k <- 50
> > m <- k + 10
> > A <- matrix(rnorm(k*m), m, k)
> > R <- chol(t(A)%*%A)
> >
> > rep <- 500
> > microbenchmark(backsolve_R(rep, R),backsolve_rcpp(rep, R),times=10)
> >
> > ## Unit: milliseconds
> > ##                    expr   min    lq median    uq   max neval
> > ##     backsolve_R(rep, R) 18.92 19.40  19.97 20.94 44.22    10
> > ##  backsolve_rcpp(rep, R) 36.67 36.88  37.13 37.28 37.41    10
> >
> > 100 by 100 uppertriangular matrix
> >
> > # 100 by 100
> > set.seed(777)
> > k <- 100
> > m <- k+10
> > A <- matrix(rnorm(k*m), m, k)
> > R <- chol(t(A)%*%A)
> >
> > rep <- 500
> > microbenchmark(backsolve_R(rep, R),backsolve_rcpp(rep,R),times=10)
> >
> > ## Unit: milliseconds
> > ##                    expr   min    lq median    uq   max neval
> > ##     backsolve_R(rep, R) 101.5 102.4  103.2 104.4 126.8    10
> > ##  backsolve_rcpp(rep, R) 251.6 251.8  251.8 252.0 254.0    10
> >
> >
> >
> > --
> >  Peter E. Rossi
> >
> <backsolve_test.cpp><backsolve_test.Rmd>_______________________________________________
> > Rcpp-devel mailing list
> > Rcpp-devel at lists.r-forge.r-project.org
> > https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
>
> _______________________________________________
> Rcpp-devel mailing list
> Rcpp-devel at lists.r-forge.r-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20140628/e2338c65/attachment.html>