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

[Peter Rossi]

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

// [[Rcpp::depends("RcppArmadillo")]]
#include <RcppArmadillo.h>
#include <Rcpp.h>

using namespace arma;

// [[Rcpp::export]]
mat backsolve_rcpp(int rep, mat R) { 
// solve() and trimatu()
      
  mat RI;
  int k = R.n_rows;
  mat I = eye(k,k);
  mat R_upptri = trimatu(R);
      
  for(int i=0; i<rep; i++){
    RI = solve(R_upptri,I);
  }
  
  return RI;
}

backsolve_test.Rmd
Description: Binary data

_______________________________________________
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel