[Rcpp-devel] arma::mat and Rcpp::NumericMatrix

[Dirk Eddelbuettel]

Hi Jeffrey,

And welcome! Thanks for a detailed post with examples and benchmarks!

On 16 September 2011 at 14:48, Jeffrey Pollock wrote:
| Hello all,
| 
| I've recently started looking into the speed gains from using Rcpp and have
| began looking into the arma library too. Ive written a small loop to fill in a
| matrix and done so using the arma::mat class and the Rcpp::NumericMatrix class.
| In this case, the Rcpp::NumericMatrix class was slightly faster. Is this
| usually the case, or have I written quite inefficient code in some parts? It's

This may not be worth worrying about.  See the list archives just two or
three weeks back: Darren Cook noticed that even within Rcpp, the times
between accessing vector (not matrix) elements using () and [] for indexing
differed far more (which may well be indicative of us doing something silly,
there may be a cache lookup gone wrong...)

For all real applications, _actual work_ is bound to dominate.  So as the
difference is small...   Also, check the Arma docs. Conrad is very generous
is his API, there are often several different ways to do similar things

| my understanding that the advantage of the arma::mat is that is has more
| methods associated with it, is this correct?

Armadillo is a full-blown linear algebra library. If you want to do 'math'
with your matrix, it is a good choice to deploy arma matrices

Rcpp::NumericMatrix is a very convenient container to get matrix data back
and from R, and its vector variant also goes easily to STL types etc so that
you can (relatively) easily exchange data with other C++ libraries etc.

So it depends.
 
| I've also added an R function to do this and as per usual, the speed gain of
| simply using Rcpp is incredible!

Yes, compiled loops often win rather handily.
 
| Any input massively appreciated

Hope this helps,  Dirk
 
| Jeff
| 
| library(RcppArmadillo)
| library(inline)
| library(rbenchmark)
| 
| # start with R function....
| 
| fillMiddleTheta <- function(thetaMatrix, thetaStartIndex, thetaEndIndex,
| alphas, alphaCounterStart, timerem) {
|     alphaCounter <- alphaCounterStart
|     for (i in (thetaStartIndex + 1L):(thetaEndIndex - 1L)) {
|         addition <- exp(alphas[alphaCounter] + alphas[alphaCounter + 1L]
| * timerem)
|         thetaMatrix[ , i] <- thetaMatrix[ , i - 1L] + addition
|         alphaCounter <- alphaCounter + 2L
|     }
|     return(thetaMatrix)
| }
| 
| # cpp function....
| 
| src <- '
|         Rcpp::NumericMatrix thetaMatrix(sthetaMatrix);
|        
|         Rcpp::NumericVector alphas(salphas);
|         Rcpp::NumericVector timerem(stimerem);
|        
|         int thetaStartIndex = as<int>(sthetaStartIndex);
|         int thetaEndIndex = as<int>(sthetaEndIndex);
|         int alphaCounterStart = as<int>(salphaCounterStart);
|        
|         int alphaCounter = alphaCounterStart;
|        
|         for (int i = thetaStartIndex + 1; i <= thetaEndIndex - 1; i++) {
|        
|         thetaMatrix(_, i - 1) = thetaMatrix(_, i - 2) + exp(alphas
| (alphaCounter - 1) + alphas(alphaCounter) * timerem);
|         alphaCounter = alphaCounter + 2;
|        
|         }
|        
|         return wrap(thetaMatrix);
|         '
| 
| cppFillMiddleTheta <- cxxfunction(signature(sthetaMatrix = "Matrix",
| sthetaStartIndex = "integer",
|                 sthetaEndIndex = "integer", salphas =
| "numericVector", salphaCounterStart = "integer", stimerem = "numericVector"),
|         src, "Rcpp")
| 
| # cpp function with armadillo....
| 
| armaSrc <- '
|         arma::mat thetaMatrix = Rcpp::as<arma::mat>(sthetaMatrix);
|        
|         arma::colvec alphas = Rcpp::as<arma::colvec>(salphas);
|         arma::colvec timerem = Rcpp::as<arma::colvec>(stimerem);
|        
|         int thetaStartIndex = as<int>(sthetaStartIndex);
|         int thetaEndIndex = as<int>(sthetaEndIndex);
|         int alphaCounterStart = as<int>(salphaCounterStart);
|        
|         int alphaCounter = alphaCounterStart;
|        
|         for (int i = thetaStartIndex + 1; i <= thetaEndIndex - 1; i++) {
|        
|         thetaMatrix.col(i - 1) = thetaMatrix.col(i - 2) + exp(alphas
| (alphaCounter - 1) + alphas(alphaCounter) * timerem);
|         alphaCounter = alphaCounter + 2;
|        
|         }
|        
|         return wrap(thetaMatrix);
|         '
| 
| cppArmaFillMiddleTheta <- cxxfunction(signature(sthetaMatrix = "Matrix",
| sthetaStartIndex = "integer",
|                 sthetaEndIndex = "integer", salphas =
| "numericVector", salphaCounterStart = "integer", stimerem = "numericVector"),
|         armaSrc, "RcppArmadillo")
| 
| # test
| n <- 1e6
| 
| thetaMatrix <- matrix(0, ncol = 10, nrow = 10)
| thetaMatrix[, 1L] <- 1
| thetaMatrix[, 10L] <- 100
| thetaStartIndex <- 1L
| thetaEndIndex <- 10L
| alphas <- seq(0, 0.06, length = 16)
| alphaCounterStart <- 1L
| timerem <- rnorm(10, mean = 40, sd = 2)
| 
| r <- fillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex, alphas,
| alphaCounterStart, timerem)
| cpp <- cppFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex, alphas,
| alphaCounterStart, timerem)
| cppArma <- cppArmaFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex,
| alphas, alphaCounterStart, timerem)
| all.equal(r, cpp, cppArma)
| 
| benchmark(fillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex, alphas,
| alphaCounterStart, timerem),
|         cppFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex,
| alphas, alphaCounterStart, timerem),
|         cppArmaFillMiddleTheta(thetaMatrix, thetaStartIndex,
| thetaEndIndex, alphas, alphaCounterStart, timerem),
|        
|         columns = c("test", "replications", "elapsed", "relative",
| "user.self", "sys.self"),
|         order = "relative",
|         replications = n)
| 
| Here is some of my output;
| 
| > n <- 1e6
| > r <- fillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex, alphas,
| alphaCounterStart, timerem)
| > cpp <- cppFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex,
| alphas, alphaCounterStart, timerem)
| > cppArma <- cppArmaFillMiddleTheta(thetaMatrix, thetaStartIndex,
| thetaEndIndex, alphas, alphaCounterStart, timerem)
| > all.equal(r, cpp, cppArma)
| [1] TRUE
| 
| > benchmark(fillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex,
| alphas, alphaCounterStart, timerem),
| +         cppFillMiddleTheta(thetaMatrix, thetaStartIndex,
| thetaEndIndex, alphas, alphaCounterStart, timerem),
| +         cppArmaFillMiddleTheta(thetaMatrix, thetaStartIndex,
| thetaEndIndex, alphas, alphaCounterStart, timerem),
| +       
| +         columns = c("test", "replications", "elapsed", "relative",
| "user.self", "sys.self"),
| +         order = "relative",
| +         replications = n)
| > benchmark(fillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex,
| alphas, alphaCounterStart, timerem),
| +         cppFillMiddleTheta(thetaMatrix, thetaStartIndex,
| thetaEndIndex, alphas, alphaCounterStart, timerem),
| +         cppArmaFillMiddleTheta(thetaMatrix, thetaStartIndex,
| thetaEndIndex, alphas, alphaCounterStart, timerem),
| +       
| +         columns = c("test", "replications", "elapsed", "relative",
| "user.self", "sys.self"),
| +         order = "relative",
| +         replications = n)
| 
|                                                                                                     
| test replications elapsed relative user.self sys.self
| 2     cppFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex,
| alphas, alphaCounterStart, timerem)      1000000   12.22
| 1.000000     12.17     0.01
| 3 cppArmaFillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex, alphas,
| alphaCounterStart, timerem)      1000000   15.07 1.233224    
| 15.01     0.02
| 1        fillMiddleTheta(thetaMatrix, thetaStartIndex, thetaEndIndex,
| alphas, alphaCounterStart, timerem)      1000000   86.62
| 7.088380     86.02     0.10
| 
| ----------------------------------------------------------------------
| _______________________________________________
| 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
-- 
New Rcpp master class for R and C++ integration is scheduled for 
San Francisco (Oct 8), more details / reg.info available at
http://www.revolutionanalytics.com/products/training/public/rcpp-master-class.php