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.

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...