Dynamic dispatch for sparse matrices

March 19, 2014
By

(This article was first published on Rcpp Gallery, and kindly contributed to R-bloggers)

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<S4>(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<arma::sp_mat>(Xr), Y) ;
        } ;
        if(Rf_inherits(Xr, "indMatrix")) {
            return matmult_ind(Xr, Y) ; 
        } ;
        stop("unknown class of Xr") ;
    } else {
        return matmult_dense(as<arma::mat>(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 on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Sponsors

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)