[Rcpp-devel] Comparison of R and rcpp for backsolve
Peter Rossi
perossichi at gmail.com
Sat Jun 28 01:14:49 CEST 2014
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
-------------- next part --------------
A non-text attachment was scrubbed...
Name: backsolve_test.cpp
Type: text/x-c++src
Size: 368 bytes
Desc: not available
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20140627/2651d5ef/attachment-0001.cpp>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: backsolve_test.Rmd
Type: application/octet-stream
Size: 1613 bytes
Desc: not available
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20140627/2651d5ef/attachment-0001.obj>
More information about the Rcpp-devel
mailing list