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

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

microbenchmark(backsolve_R(rep, R), backsolve_rcpp(rep, R), times=10)