Dynamic dispatch for sparse matrices

[This article was first published on Rcpp Gallery, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

We want to do matrix multiplication for 3 cases:

  • dense times dense
  • sparse times dense for sparse matrices of class dgCMatrix
  • sparse times dense for sparse matrices of class indMatrix,

using R’s Matrix package for sparse matrices in R and RcppArmadillo for C++ linear algebra:

// [[Rcpp::depends(RcppArmadillo)]]
#include 
using namespace Rcpp ;

arma::mat matmult_sp(const arma::sp_mat X, const arma::mat Y){
    arma::mat ret = X * Y;
    return ret;
};
arma::mat matmult_dense(const arma::mat X, const arma::mat Y){
    arma::mat ret = X * Y;
    return ret;
};
arma::mat matmult_ind(const SEXP Xr, const arma::mat Y){
    // pre-multiplication with index matrix is a permutation of Y's rows: 
    arma::uvec perm =  as(Xr).slot("perm");
    arma::mat ret = Y.rows(perm - 1);
    return ret;
};

//[[Rcpp::export]]
arma::mat matmult_cpp(SEXP Xr, const arma::mat Y) {
    if (Rf_isS4(Xr)) {
        if(Rf_inherits(Xr, "dgCMatrix")) {
            return matmult_sp(as(Xr), Y) ;
        } ;
        if(Rf_inherits(Xr, "indMatrix")) {
            return matmult_ind(Xr, Y) ; 
        } ;
        stop("unknown class of Xr") ;
    } else {
        return matmult_dense(as(Xr), Y) ;
    } 
}

Set up test cases:

library(Matrix)

Loading required package: methods

library(rbenchmark)
set.seed(12211212)
n <- 1000
d <- 50
p <- 30  

X <- matrix(rnorm(n*d), n, d)
X_sp <- as(diag(n)[,1:d], "dgCMatrix")
X_ind <- as(sample(1:d, n, rep=TRUE), "indMatrix")
Y <- matrix(1:(d*p), d, p)

Check exception handling:

matmult_cpp(as(X_ind, "ngTMatrix"), Y)

Error: unknown class of Xr

Dense times dense:

all.equal(X%*%Y, matmult_cpp(X, Y))

[1] TRUE

benchmark(X%*%Y, 
          matmult_cpp(X, Y),
          replications=100)[,1:4]

               test replications elapsed relative
2 matmult_cpp(X, Y)          100   0.086     1.00
1           X %*% Y          100   0.098     1.14

dgCMatrix times dense:

all.equal(as(X_sp%*%Y, "matrix"), matmult_cpp(X_sp, Y),
          check.attributes = FALSE)

[1] TRUE

benchmark(X_sp%*%Y, 
          matmult_cpp(X_sp, Y),
          replications=100)[,1:4]

                  test replications elapsed relative
2 matmult_cpp(X_sp, Y)          100   0.009    1.000
1           X_sp %*% Y          100   0.025    2.778

indMatrix times dense:

all.equal(X_ind%*%Y, matmult_cpp(X_ind, Y))

[1] TRUE

benchmark(X_ind%*%Y, 
          matmult_cpp(X_ind, Y),
          replications=100)[,1:4]

                   test replications elapsed relative
2 matmult_cpp(X_ind, Y)          100   0.013    1.000
1           X_ind %*% Y          100   0.025    1.923

Based on this Q&A on StackOverflow, thanks again to Kevin Ushey for his helpful comment.

To leave a comment for the author, please follow the link and comment on their blog: Rcpp Gallery.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)