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

[Gabor Grothendieck]

On Fri, Jun 27, 2014 at 7: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
>

Try eigen.  On my test run it took roughly half the time of R for 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) 76.62980 85.39488 86.71858 88.01366 89.92510    10
 backsolve_rcpp(rep, R) 41.51224 46.72658 47.10239 47.13911 47.35077    10

See attachment for code.


-- 
Statistics & Software Consulting
GKX Group, GKX Associates Inc.
tel: 1-877-GKX-GROUP
email: ggrothendieck at gmail.com
-------------- next part --------------
A non-text attachment was scrubbed...
Name: backsolve_eigen_test.cpp
Type: text/x-c++src
Size: 883 bytes
Desc: not available
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20140627/2034608b/attachment.cpp>