CODEWITHSUNDEEP
rrcppgsl
Stéphane Laurent
Stéphane Laurent

Reputation: 84599

Fast subsetting of a matrix with RcppGSL

After viewing this post, I tried to subset a matrix with Rcpp.

With RcppArmadillo:

// [[Rcpp::depends(RcppArmadillo)]]
#include "RcppArmadillo.h"
// [[Rcpp::export]]
arma::mat submatrix(const arma::mat& m1in, int fromin, int toin){
  arma::mat s1 = m1in.cols(fromin-1,toin-1);
  return(s1);
}

Then submatrix(M, 1, 900) is a bit faster than M[,1:900].

With RcppGSL:

#include <RcppGSL.h>
#include <gsl/gsl_matrix.h>
// [[Rcpp::export]]
gsl_matrix_const_view submatrix(const RcppGSL::Matrix & X, int k1, int k2, int n1, int n2) {
  return gsl_matrix_const_submatrix(X, k1, k2, n1, n2);
}

Here submatrix(M, 0, 0, 1000, 900) is slower than M[,1:900]:

> microbenchmark(M[,1:900], submatrix(M, 0, 0, 1000, 900))
Unit: milliseconds
                          expr       min       lq     mean   median       uq      max neval
                    M[, 1:900]  8.035749 10.20265 13.25657 11.75554 14.27586 117.2533   100
 submatrix(M, 0, 0, 1000, 900) 16.597605 19.55858 23.04454 21.52959 23.98431 141.6158   100

Is there a faster way to subset a matrix with RcppGSL?

Upvotes: 1

Views: 115

Answers (1)

F. Priv&#233;
F. Priv&#233;

Reputation: 11728

I think the reason is that your matrix is not passed by reference (maybe because R matrices and GSL matrices are not compatible).

To prove my point, test this:

// [[Rcpp::depends(RcppGSL)]]
#include <RcppGSL.h>
#include <gsl/gsl_matrix.h>

// [[Rcpp::export]]
gsl_matrix_const_view submatrix(RcppGSL::Matrix & X, int k1, int k2, int n1, int n2) {
  X(0, 0) = 1;
  return gsl_matrix_const_submatrix(X, k1, k2, n1, n2);
}

/*** R
M <- matrix(0, 1000, 1000)
test <- submatrix(M, 0, 0, 1000, 900)
M[1, 1]
*/

If I'm correct, you will have the same problem everytime you use RcppGSL. Maybe there exists a Map view of a matrix (like in Eigen) to use instead of & (I don't know GSL).

Upvotes: 1

Related Questions