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

Jeffrey Pollock jeffpollock9 at gmail.com
Fri Sep 16 17:31:49 CEST 2011


Hey Dirk,

Thanks for the informative and speedy reply! Having a look through the
archives now....

I should hopefully have time to post up quite a few benchmarks of work I'm
doing at the moment/as I learn more of Rcpp.

Thanks again,

Jeff

On Fri, Sep 16, 2011 at 3:34 PM, Dirk Eddelbuettel <edd at debian.org> wrote:

>
> 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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20110916/2e385060/attachment.htm>


More information about the Rcpp-devel mailing list