R中矩阵的快速子集 [英] Fast subsetting of a matrix in R

问题描述

我面临以下问题:我需要一个大矩阵的许多子集.实际上，我只需要将视图作为另一个函数f()的输入，因此不需要更改值.但是似乎R的执行速度非常慢，或者我做错了(似乎更有可能).玩具示例说明了选择列然后在另一个函数中使用它们需要花费多少时间(在此玩具示例中为原始函数sum()).作为基准"，我还测试了计算时间是否不加总整个矩阵，这出乎意料地更快.我还尝试了ref包，但是无法获得任何性能提升. 所以关键问题是如何在不复制矩阵的情况下对矩阵进行子集化?感谢您的帮助，谢谢！

I face the following problem: I need many subsets of a big matrix. Actually I just need views as input for another function f(), so I don't need to change the values. However it seems, that R is terribly slow for this task, or I'm doing something wrong (which seems more likely). The toy example illustrates how much time it takes just to select the columns, and then use them in another function (in this toy example the primitive function sum()). As 'benchmark' I also test the calculation time against summing up the whole matrix, which is surprisingly faster. I also experimented with the package ref, however could not achieve any performance gain. So the key question is how to subset the matrix without copying it? I appreciate any help, Thanks!

library(microbenchmark)
library(ref)

m0 <- matrix(rnorm(10^6), 10^3, 10^3)
r0 <- refdata(m0)
microbenchmark(m0[, 1:900], sum(m0[, 1:900]), sum(r0[,1:900]), sum(m0))

Unit: milliseconds
             expr       min        lq      mean    median        uq
      m0[, 1:900] 10.087403 12.350751 16.697078 18.307475 19.054157
 sum(m0[, 1:900]) 11.067583 13.341860 17.286514 19.123748 19.990661
 sum(r0[, 1:900]) 11.066164 13.194244 16.869551 19.204434 20.004034
          sum(m0)  1.015247  1.040574  1.059872  1.049513  1.067142
       max neval
 58.238217   100
 25.664729   100
 23.505308   100
  1.233617   100

对整个矩阵求和的基准任务需要1.059872毫秒，比其他函数快16倍.

The benchmark task of summing the whole matrix takes 1.059872 milliseconds and is about 16 times faster than the other functions.

推荐答案

解决方案的问题是子集分配了另一个矩阵，这需要时间.

The problem with your solution is that the subsetting is allocating another matrix, which takes times.

您有两种解决方案:

如果在整个矩阵上使用sum花费的时间还可以，则可以在整个矩阵上使用colSums并将结果子集化:

If the time taken with sum on the whole matrix is okay with you, you could use colSums on the whole matrix and subset the result:

sum(colSums(m0)[1:900])

或者您可以使用Rcpp通过子集计算sum，而无需复制矩阵.

Or you could use Rcpp to compute the sum with subsetting without copying the matrix.

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
double sumSub(const NumericMatrix& x,
              const IntegerVector& colInd) {

  double sum = 0;

  for (IntegerVector::const_iterator it = colInd.begin(); it != colInd.end(); ++it) {
    int j = *it - 1;
    for (int i = 0; i < x.nrow(); i++) {
      sum += x(i, j);
    }
  }

  return sum;
}

    microbenchmark(m0[, 1:900], sum(m0[, 1:900]), sum(r0[,1:900]), sum(m0),
                   sum(colSums(m0)[1:900]),
                   sumSub(m0, 1:900))
Unit: milliseconds
                    expr      min       lq     mean   median       uq       max neval
             m0[, 1:900] 4.831616 5.447749 5.641096 5.675774 5.861052  6.418266   100
        sum(m0[, 1:900]) 6.103985 6.475921 7.052001 6.723035 6.999226 37.085345   100
        sum(r0[, 1:900]) 6.224850 6.449210 6.728681 6.705366 6.943689  7.565842   100
                 sum(m0) 1.110073 1.145906 1.175224 1.168696 1.197889  1.269589   100
 sum(colSums(m0)[1:900]) 1.113834 1.141411 1.178913 1.168312 1.201827  1.408785   100
       sumSub(m0, 1:900) 1.337188 1.368383 1.404744 1.390846 1.415434  2.459361   100

您可以使用展开优化进一步优化Rcpp版本.

You could use unrolling optimization to further optimize the Rcpp version.

这篇关于R中矩阵的快速子集的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！