arrays - 矩阵中的累积和

Question

我有一个矩阵

A= [ 1 2 4
     2 3 1
     3 1 2 ]

我想逐行逐列计算它的累积总和，也就是说，我希望结果是

B = [ 1  3  7 
      3  8  13 
      6  12 19 ]

关于如何在 R 中快速实现这一点的任何想法？（可能使用函数 cumsum）（我有巨大的矩阵）

谢谢！

Answer 1

单线：

t(apply(apply(A, 2, cumsum)), 1, cumsum))

基本观察是，您可以首先计算列上的累积和，然后计算行上该矩阵的累积和。

注意：在做行时，您必须转置结果矩阵。

你的例子：

> apply(A, 2, cumsum)
     [,1] [,2] [,3]
[1,]    1    2    4
[2,]    3    5    5
[3,]    6    6    7

> t(apply(apply(A, 2, cumsum), 1, cumsum))
     [,1] [,2] [,3]
[1,]    1    3    7
[2,]    3    8   13
[3,]    6   12   19

关于性能：我现在知道这种方法扩展到大矩阵的效果有多好。在复杂性方面，这应该接近最优。通常，apply性能也没有那么差。

---

编辑

现在我开始好奇 - 哪种方法更好？一个简短的基准：

> A <- matrix(runif(1000*1000, 1, 500), 1000)
> 
> system.time(
+   B <- t(apply(apply(A, 2, cumsum), 1, cumsum))
+ )
       User      System     elapsed 
      0.082       0.011       0.093 
> 
> system.time(
+   C <- lower.tri(diag(nrow(A)), diag = TRUE) %*% A %*% upper.tri(diag(ncol(A)), diag = TRUE)
+ )
       User      System     elapsed 
      1.519       0.016       1.530 

因此：Apply 的性能优于矩阵乘法 15 倍。（只是为了比较：MATLAB 需要 0.10719 秒。）结果并不令人惊讶，因为apply-version 可以在 O(n^2) 中完成，而矩阵乘法将需要大约。O(n^2.7) 次计算。因此，如果 n 足够大，矩阵乘法提供的所有优化都应该丢失。

Answer 2

这是使用 matrixStats 包和更大的示例矩阵的更有效的实现：

library(matrixStats)
A <- matrix(runif(10000*10000, 1, 500), 10000)

# Thilo's answer
system.time(B <- t(apply(apply(A, 2, cumsum), 1, cumsum)))
user  system elapsed 
3.684   0.504   4.201

# using matrixStats
system.time(C <- colCumsums(rowCumsums(A)))
user  system elapsed 
0.164   0.068   0.233 

all.equal(B, C)
[1] TRUE

Answer 3

我的解决方案：函数 cumsum_row()（见下文）采用矩阵 M 并返回 M 行的累积和的矩阵。函数 cumsum_col() 对列做同样的事情。

cumsum_row <- function(M) {
  M2 <- c()
  for (i in 1:nrow(M))
    M2 <- rbind(M2, cumsum(M[i,]))
  return (M2)
}

cumsum_col <- function(M) {
  return (t(cumsum_row(t(M))))
}

例子：

  > M <- matrix(rep(1, 9), nrow=3)
  > M
         [,1] [,2] [,3]
    [1,]    1    1    1
    [2,]    1    1    1
    [3,]    1    1    1

  > cumsum_row(M)
         [,1] [,2] [,3]
    [1,]    1    2    3
    [2,]    1    2    3
    [3,]    1    2    3