I’m extracting a submatrix B from a matrix A using Rcpp. Benchmarking has shown that this is slower than the ordinary R way.
Why is it slow and is there a faster way with Rcpp? (the indices are not contiguous in general). I have also tested with RcppEigen, but the result is the same.
The reason why extracting a submatrix using Rcpp is slower than the ordinary R way is because Rcpp performs the extraction in a row-wise manner, which can be inefficient for non-contiguous indices.
To achieve faster submatrix extraction with Rcpp, you can use the Eigen library in RcppEigen. Eigen provides a more efficient way to extract submatrices using a column-major storage format, which aligns with how matrices are stored in R.
Here’s an example code snippet using RcppEigen to extract a submatrix from a matrix A:
#include <RcppEigen.h>
// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::export]]
Eigen::MatrixXd extractSubmatrix(Eigen::MatrixXd A, Eigen::VectorXi rows, Eigen::VectorXi cols) {
return A(rows, cols);
}
In this code, A is the input matrix, and rows and cols are vectors specifying the row and column indices of the submatrix to extract. The function extractSubmatrix returns the extracted submatrix as a MatrixXd object.
You can then call this function from R to extract the submatrix efficiently:
library(Rcpp)
sourceCpp("path/to/cpp/file.cpp")
A <- matrix(1:9, nrow = 3)
rows <- c(2, 3)
cols <- c(1, 3)
B <- extractSubmatrix(A, rows, cols)
By using RcppEigen and the Eigen library, you should be able to achieve faster submatrix extraction compared to the ordinary R way.