在 R 中快速地对矩阵进行子集操作

Question

我遇到了以下问题：我需要从一个大矩阵中获取许多子集。实际上，我只需要作为另一个函数 f() 的输入视图，因此我不需要更改值。然而，似乎 R 在这个任务上非常慢，或者我做错了什么（后者更可能）。玩具示例说明了选择列所需的时间以及在另一个函数中使用它们的时间（在这个玩具示例中是原始函数 sum()）。作为“基准测试”，我还测试了将整个矩阵相加所需的计算时间，这是出人意料的快。我还尝试使用 ref 包，但无法获得任何性能增益。因此，关键问题是如何对矩阵进行子集处理而不复制它？感谢任何帮助！

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 倍。

- maxatSOflow

---

当我尝试你的示例时，出现错误：r0[, 1:900]中的错误：维度数量不正确。 - 5th

library('ref') 修复了它。 - dvantwisk

好的，那么解决方案的标准是什么？它只需要比你现在做的更快，还是必须比整个矩阵求和更快？ - Hack-R

使用Rcpp或compiler可能会给你带来所需的性能提升。 - Hack-R

@Hack-R 编译器通常并不是很有用。 - F. Privé

@F.Privé 不知道，它仍然提升了50%的性能，但是毋庸置疑，你的方法肯定更好。编译器在这个非常相似的线程中表现得更好：https://dev59.com/OYvda4cB1Zd3GeqPfuqD 好吧，无论如何，编译器比Rcpp更容易使用，但没有什么比你提出的第一个解决方案更好(sum(colSums(m0)[1:900]))。 - Hack-R

Answer 1

你的解决方案存在问题，就是子集分配了另一个矩阵，这会花费时间。

你有两种解决方法：

如果对整个矩阵使用sum仍在可接受范围内，你可以对整个矩阵使用colSums并对结果进行子集分配：

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

或者您可以使用Rcpp来计算具有子集的sum，而无需复制矩阵。

#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版本。

- F. Privé

---

是的，这就是答案。+1。我认为在我们实施之前指出它的评论也应该得到+1 ;) 但真的，如果sum(colSums(m0)[1:900])做得很好，我们甚至不需要更繁琐的Rcpp。 - Hack-R

感谢您的回复。是的，我认为您是完全正确的。主要问题在于为矩阵副本分配内存。在R中是否有可能不复制矩阵而完成此操作？我需要将矩阵子集作为另一个自编函数的输入，sum()函数只是这里的示例。 - maxatSOflow

基本上，我所有的函数都可以在矩阵的子集上运作。如果你不想做一个副本，你只需要这样的函数。 - F. Privé

Answer 2

我使用编译器编写了一个函数，它的速度比你使用的其他方法快2倍（而不是16倍sum(m0)的速度为8倍）：

require(compiler)

compiler_sum <- cmpfun({function(x) {
     tmp <- 0
     for (i in 1:900)
         tmp <- tmp+sum(x[,i])
     tmp
}})

microbenchmark( 
               sum(m0),
               compiler_sum(m0)
               )

Unit: milliseconds
             expr      min       lq     mean   median      uq       max
          sum(m0) 1.016532 1.056030 1.107263 1.084503 1.11173  1.634391
 compiler_sum(m0) 7.655251 7.854135 8.000521 8.021107 8.29850 16.760058
 neval
   100
   100

- Hack-R