使用 Rcpp 的高效矩阵子集

我正在尝试找到一种有效的方法来使用 Rcpp 对一组非连续的行和列进行矩阵子集：

m <- matrix(1:20000000, nrow=5000)
rows <- sample(1:5000, 100)
cols <- sample(1:4000, 100)

在 R 中，矩阵可以使用rows和cols向量直接子集：

matrix_subsetting <- function(m, rows, cols){
return(m[rows, cols])
}
m[rows, cols]
# or
matrix_subsetting(m, rows, cols)

到目前为止，我能够找到的最快的Rcpp方法是：

Rcpp::cppFunction("

NumericMatrix cpp_matrix_subsetting(NumericMatrix m, NumericVector rows, NumericVector cols){

int rl = rows.length();
int cl = cols.length();
NumericMatrix out(rl, cl);

for (int i=0; i<cl; i++){
NumericMatrix::Column org_c = m(_, cols[i]-1);
NumericMatrix::Column new_c = out(_, i);
for (int j=0; j<rl; j++){
new_c[j] = org_c[rows[j]-1];
}
}
return(out);
}

")

但相比之下，Rcpp 版本要慢得多：

> microbenchmark::microbenchmark(matrix_subsetting(m, rows, cols), cpp_matrix_subsetting(m, rows, cols), times=500)
Unit: microseconds
expr       min        lq       mean    median         uq        max neval
matrix_subsetting(m, rows, cols)    23.269    90.127   107.8273   130.347   135.3285    605.235   500
cpp_matrix_subsetting(m, rows, cols) 69191.784 75254.277 88484.9328 90477.448 95611.9090 178903.973   500

有什么想法，至少要获得与 Rcpp 相当的速度吗？

我已经尝试了RcppArmadilloarma::mat::submat功能，但它比我的版本慢。

使用IntegerMatrix而不是NumericMatrix实现cpp_matrix_subsetting函数。

新基准：

> microbenchmark::microbenchmark(matrix_subsetting(m, rows, cols), cpp_matrix_subsetting(m, rows, cols), times=1e4)
Unit: microseconds
expr    min     lq     mean median      uq      max neval
matrix_subsetting(m, rows, cols) 41.110 60.261 66.88845 61.730 63.8900 14723.52 10000
cpp_matrix_subsetting(m, rows, cols) 43.703 61.936 71.56733 63.362 65.8445 27314.11 10000