使用 RcppGSL 对矩阵进行快速子集化答案

【问题标题】：Fast subsetting of a matrix with RcppGSL使用 RcppGSL 对矩阵进行快速子集化
【发布时间】：2017-09-05 09:07:30
【问题描述】：

查看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 对矩阵进行子集化？

【问题讨论】：

标签： r rcpp gsl

【解决方案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）。

【讨论】：