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

[Hao Ye]

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