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

[Jeffrey Pollock]

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 my understanding that the advantage of the arma::mat is that is
has more methods associated with it, is this correct?

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

Any input massively appreciated

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