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

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