使累积和更快答案

【问题标题】：Make cumulative sum faster使累积和更快
【发布时间】：2015-06-12 15:20:50
【问题描述】：

我正在尝试对矩阵的每一列进行累积总和。这是我在 R 中的代码：

testMatrix = matrix(1:65536, ncol=256);
microbenchmark(apply(testMatrix, 2, cumsum), times=100L);

Unit: milliseconds
                         expr      min       lq     mean  median       uq      max neval
 apply(testMatrix, 2, cumsum) 1.599051 1.766112 2.329932 2.15326 2.221538 93.84911 10000

我使用了 Rcpp 进行比较：

cppFunction('NumericMatrix apply_cumsum_col(NumericMatrix m) {
    for (int j = 0; j < m.ncol(); ++j) {
        for (int i = 1; i < m.nrow(); ++i) {
            m(i, j) += m(i - 1, j);
        }
    }
    return m;
}');
microbenchmark(apply_cumsum_col(testMatrix), times=10000L);

Unit: microseconds
                         expr     min      lq     mean  median      uq      max neval
 apply_cumsum_col(testMatrix) 205.833 257.719 309.9949 265.986 276.534 96398.93 10000

所以 C++ 代码的速度是原来的 7.5 倍。在纯 R 中是否有可能比 apply(testMatrix, 2, cumsum) 做得更好？感觉就像我无缘无故地有一个数量级的开销。

【问题讨论】：

您可以尝试通过compile:cmpfun 和该库中的其他工具进行编译。但是，众所周知 R 与 MATLAB 和类似语言一样，由于在命令时“编译”而具有大量开销。这就是为什么人们在需要最大速度时会在 Fortran 或 cpp 中编写函数的核心。
apply(testMatrix, 2, cumsum) 的快速替代方案是 matrixStats::colCumsums(testMatrix)。
@Khashaa，我想你的比 R 更快，因为它也使用 C 代码。我相信作者是在询问严格的 R 实现。
Khashaa 的方法要求用户仅调用 R 函数，我认为这是 OP 正在寻找的。许多“R”函数实际上是 C 或 C++。在我的电脑上，matrixStats 包中的colCumsums 实际上比 Rcpp 函数快大约 35%。
@cdeterman 很高兴知道所有选项！

标签： r rcpp

【解决方案1】：

仅使用 R 代码很难击败 C++。我能想到的最快方法是如果您愿意将矩阵拆分为列表。这样，R 使用原始函数并且不会在每次迭代时复制对象（apply 本质上是一个漂亮的循环）。您可以看到 C++ 仍然胜出，但如果您真的只想使用 R 代码，使用 list 方法可以显着加快速度。

fun1 <- function(){
    apply(testMatrix, 2, cumsum)
}

testList <- split(testMatrix, col(testMatrix))

fun2 <- function(){
    lapply(testList, cumsum)
}

microbenchmark(fun1(),
               fun2(),
               apply_cumsum_col(testMatrix),
               times=100L)


Unit: microseconds
                         expr      min        lq      mean   median        uq      max neval
                       fun1() 3298.534 3411.9910 4376.4544 3477.608 3699.2485 9249.919   100
                       fun2()  558.800  596.0605  766.2377  630.841  659.3015 5153.100   100
 apply_cumsum_col(testMatrix)  219.651  282.8570  576.9958  311.562  339.5680 4915.290   100

编辑请注意，如果包含将矩阵拆分为列表的时间，则此方法比 fun1 慢。

【讨论】：

难道不应该考虑将矩阵拆分为列表所需的时间吗？
@nrussell，那么即使fun1 也会更慢。但这是我能想到的唯一可以提高cumsum 速度的场景。绝对愿意接受其他答案:)

【解决方案2】：

使用字节编译的 for 循环比在我的系统上调用 apply 稍快。我预计它会更快，因为它的工作量比apply 少。正如预期的那样，R 循环仍然比您编写的简单 C++ 函数慢。

colCumsum <- compiler::cmpfun(function(x) {
  for (i in 1:ncol(x))
    x[,i] <- cumsum(x[,i])
  x
})

testMatrix <- matrix(1:65536, ncol=256)
m <- testMatrix
require(microbenchmark)
microbenchmark(colCumsum(m), apply_cumsum_col(m), apply(m, 2, cumsum), times=100L)
# Unit: microseconds
#                 expr      min        lq    median        uq       max neval
#      matrixCumsum(m) 1478.671 1540.5945 1586.1185 2199.9530 37377.114   100
#  apply_cumsum_col(m)  178.214  192.4375  204.3905  234.8245  1616.030   100
#  apply(m, 2, cumsum) 1879.850 1940.1615 1991.3125 2745.8975  4346.802   100
all.equal(colCumsum(m), apply(m, 2, cumsum))
# [1] TRUE

【讨论】：

【解决方案3】：

也许为时已晚，但我会写下我的答案，以便其他人可以看到。

首先，在您的 C++ 代码中，您需要克隆您的矩阵，否则您将写入 R 的内存，并且 CRAN 禁止这样做。所以你的代码变成了：

rcpp_apply<-cppFunction('NumericMatrix apply_cumsum_col(NumericMatrix m) {
    NumericMatrix g=clone(m);
    for (int j = 0; j < m.ncol(); ++j) {
        for (int i = 1; i < m.nrow(); ++i) {
            g(i, j) += g(i - 1, j);
        }
    }
    return g;
}');

由于您的矩阵是 typeof integer，那么您可以将 C++ 的参数更改为 IntegerMatrix。

rcpp_apply_integer<-cppFunction('IntegerMatrix apply_cumsum_col(IntegerMatrix m) {
    NumericMatrix g=clone(m);
    for (int j = 0; j < m.ncol(); ++j) {
        for (int i = 1; i < m.nrow(); ++i) {
            g(i, j) += g(i - 1, j);
        }
    }
    return g;
}');

这将代码改进了大约 2 倍。这是一个基准：

microbenchmark::microbenchmark(R=apply(testMatrix, 2, cumsum),Rcpp=rcpp_apply(testMatrix),Rcpp_integer=rcpp_apply_integer(testMatrix), times=10)

Unit: microseconds
        expr      min       lq      mean    median       uq      max neval
           R 1552.217 1706.165 1770.1264 1740.0345 1897.884 1940.989    10
        Rcpp  502.900  523.838  637.7188  665.0605  699.134  743.471    10
Rcpp_integer  220.455  274.645  274.9327  275.8770  277.930  316.109    10



all.equal(rcpp_apply(testMatrix),rcpp_apply_integer(testMatrix))
[1] TRUE

如果你的矩阵有很大的值，那么你必须使用NumericMatrix。

【讨论】：