[Rcpp-devel] Comparison of R and rcpp for backsolve
Gabor Grothendieck
ggrothendieck at gmail.com
Sat Jun 28 03:46:38 CEST 2014
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>
More information about the Rcpp-devel
mailing list