用RcppGSL实现矩阵的快速子集

Question

在查看了this post之后，我尝试用Rcpp对一个矩阵进行子集。

使用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);
}

那么submatrix(M, 1, 900)比M[,1:900]要快一点。

使用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);
}

这里submatrix(M, 0, 0, 1000, 900)比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

有没有一种更快的方法用RcppGSL对矩阵进行子集

Answer 1

我认为原因是您的矩阵不是通过引用传递的(可能是因为R矩阵和GSL矩阵不兼容)。

为了证明我的观点，测试一下：

// [[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]
*/

如果我没记错的话，你每次使用RcppGSL都会遇到同样的问题。也许有一个矩阵的Map视图(就像在Eigen中)可以用来代替& (我不知道GSL)。