矩陣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。也許存在矩陣（如Eigen）的地圖視圖來代替&（我不知道GSL）。