[Rcpp-devel] Easy way to subset a matrix

Jeffrey Pollock jeffpollock9 at gmail.com
Wed Sep 21 15:06:37 CEST 2011


I had a bash at this as I'm currently learning some C++ / Rcpp. I made a
function where you pass in the row start index and row end index, and the
desired matrix column (which is what I assumed you would want).

Here is the code I came up with including some benchmarks I ran,
interestingly, using a bigger matrix meant that the speedup from using
RcppArmadillo was less, I guess this is to do with the question I asked
earlier on the mailing list (
http://lists.r-forge.r-project.org/pipermail/rcpp-devel/2011-September/002876.html)
but I've not yet had time to fully understand the solutions.

Jeff

library(inline)
library(Rcpp)

src <- '
        arma::mat mat = as<arma::mat>(sMat);
        int rowFrom = as<int>(sRowFrom) - 1;
        int rowTo = as<int>(sRowTo) - 1;
        int col = as<int>(sCol) - 1;

        arma::colvec column = mat(arma::span(rowFrom, rowTo),
arma::span(col, col));

        return wrap(arma::var(column));
        '

cppFunction <- cxxfunction(signature(sMat = "NumericMatrix", sRowFrom =
"integer", sRowTo = "integer", sCol = "integer"),
        src, "RcppArmadillo")

rFunction <- function(mat, rowFrom, rowTo, col) {
    return(var(mat[rowFrom:rowTo, col]))
}

# test
mat <- matrix(rnorm(100), 10, 10)
all.equal(cppFunction(mat, 2, 5, 2), rFunction(mat, 2, 5, 2))

# benchmark
library(rbenchmark)

benchmark(cppFunction(mat, 2, 5, 2), rFunction(mat, 2, 5, 2),
        columns = c("test", "replications", "elapsed", "relative",
"user.self", "sys.self"),
        order = "relative", replications = 1e5)

                       test replications elapsed relative user.self sys.self
1 cppFunction(mat, 2, 5, 2)       100000    0.82 1.000000      0.81        0
2   rFunction(mat, 2, 5, 2)       100000    5.41 6.597561      5.41        0

# big matrix benchmark
bigMat <- matrix(rnorm(10000), 100, 100)

benchmark(cppFunction(bigMat, 2, 5, 2), rFunction(bigMat, 2, 5, 2),
        columns = c("test", "replications", "elapsed", "relative",
"user.self", "sys.self"),
        order = "relative", replications = 1e5)

                          test replications elapsed relative user.self
sys.self
1 cppFunction(bigMat, 2, 5, 2)       100000    2.19 1.000000
2.19        0
2   rFunction(bigMat, 2, 5, 2)       100000    5.45 2.488584
5.43        0


On Wed, Sep 21, 2011 at 1:35 PM, Dirk Eddelbuettel <edd at debian.org> wrote:

>
> On 20 September 2011 at 23:00, Noah Silverman wrote:
> | Nice suggestion, but it may not work for what I'm trying to do.
> |
> | I'm building up a matrix of values over a time series as part of a big
> loop.
> |  At certain iterations, I need to calculate some summary statistics on a
> few
> | things.  In R, it is trivial to subset any part of a matrix.  Not sure to
> do
> | that in Rcpp.
>
> RcppArmardillo makes Armadillo and Rcpp (and hence R) interoperable.
>
> Which means Zarrar has just solved your problem for you. I don't quite see
> where that 'may not work' as you claim.  Care to expand?
>
> You also may want to look at things like the Rcpp::SubMatrix class. Here is
> a
> quick cut and paste from one of the unit test:
>
>       NumericMatrix xx(4, 5);
>       xx(0,0) = 3;
>       xx(0,1) = 4;
>       xx(0,2) = 5;
>       xx(1,_) = xx(0,_);
>       xx(_,3) = xx(_,2);
>       SubMatrix<REALSXP> yy = xx( Range(0,2), Range(0,3) ) ;
>       NumericMatrix res = yy ;
>       return res;
>
> Supply a matrix and two ranges, get back a submatrix.  So it's there.  Also
> note the sugar syntax used here accessing entire rows and colums, just like
> in R.
>
> Dirk
>
> |
> |
> | --
> | Noah Silverman
> | UCLA Department of Statistics
> | 8117 Math Sciences Building
> | Los Angeles, CA 90095
> |
> | On Sep 20, 2011, at 10:09 PM, Zarrar Shehzad wrote:
> |
> |
> |     I am not sure if there are native functions in Rcpp but you could use
> |     RcppArmadillo to solve your problem.
> |
> |     So say Xs = x:
> |
> |     // Convert from SEXP => Rcpp => Arma
> |     Rcpp::NumericMatrix Xr(Xs);
> |     arma::mat X(Xr.begin(), Xr.nrow(), Xr.ncol(), false);
> |
> |     // Get subset of matrix and calculate variance
> |     // (i.e., do var(x[3:10,5]))
> |     arma::var( X.submat(3, 10, 5, 5) );
> |
> |     The armadillo docs are really nice so you can see those for more info
> on
> |     the var and submat functions.
> |
> |     Cheers
> |     Zarrar
> |
> |     On Wed, Sep 21, 2011 at 12:16 AM, Noah Silverman <
> noahsilverman at ucla.edu>
> |     wrote:
> |
> |         Hello,
> |
> |         I want to calculate the variance of a subset of a matrix column.
> |
> |         For example, if I wanted the variance of items 3-10 in column 5.
> |
> |         In R, this would be:
> |
> |         x <- matrix(rnorm(100), nrow=10, ncol=10)
> |         varx <- var(x[3:10,5])
> |
> |         In Rcpp, I can construct a matrix object Rcpp::NumericMatrix x
> |         The var function is available thanks to the great Sugar
> implementation
> |
> |         But, how can I easily reference a subset of a column to calculate
> the
> |         variance?
> |
> |         Ideas?
> |
> |
> |         --
> |         Noah Silverman
> |         UCLA Department of Statistics
> |         8117 Math Sciences Building #8208
> |         Los Angeles, CA 90095
> |
> |
> |         _______________________________________________
> |         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
> |
> |
> |
> |
> |
> | ----------------------------------------------------------------------
> | _______________________________________________
> | 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
> _______________________________________________
> 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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20110921/47831e3e/attachment.htm>


More information about the Rcpp-devel mailing list