使用 Rcpp 的高效矩阵子集

Posted 2023-02-22

【中文标题】使用 Rcpp 的高效矩阵子集

【英文标题】：Efficient Matrix subsetting with Rcpp

【发布时间】：2019-12-11 10:54:05

【问题描述】：

我正在尝试找到一种有效的方法来使用 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

【问题讨论】：

您尝试过更大的矩阵吗？我们在这里讨论的是毫秒。

我尝试了更大的矩阵，但差异变得更大了。 （我编辑了问题）

【参考方案1】：

这是因为您有一个 integer 类型的矩阵 m（不是 double，因为 NumericMatrix 所期望的）所以这会复制整个矩阵（这需要很多时间）。

例如，尝试改用m <- matrix(1:20000000 + 0, nrow=5000)。

【讨论】：

我用IntegerMatrix 代替NumericMatrix 重写了Rcpp 函数，现在它和R 版本一样快。

你应该创建一个检测SEXP类型的函数并调用2个不同的函数。