15

考虑数组a：

> a <- array(c(1:9, 1:9), c(3,3,2))
> a
, , 1

     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9

, , 2

     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9

我们如何有效地计算由第三维索引的矩阵的行和，使得结果是：

     [,1] [,2]
[1,]   12   12
[2,]   15   15
[3,]   18   18

??

通过 的'dims'参数，列总和很容易colSums()：

> colSums(a, dims = 1)

但我找不到在数组上使用来达到预期结果的方法，因为rowSums()它对.'dims'colSums()

使用以下方法计算所需的行总和很简单：

> apply(a, 3, rowSums)
     [,1] [,2]
[1,]   12   12
[2,]   15   15
[3,]   18   18

但这只是隐藏循环。是否有其他有效的、真正矢量化的方法来计算所需的行总和？

4

4 回答 4

13

@Fojtasek 的回答提到拆分数组让我想起了aperm()允许置换数组维度的函数。作为colSums()作品，我们可以使用输出交换前两个维度aperm()并colSums()在输出上运行。

> colSums(aperm(a, c(2,1,3)))
     [,1] [,2]
[1,]   12   12
[2,]   15   15
[3,]   18   18

这个和其他建议的基于 R 的答案的一些比较时间：

> b <- array(c(1:250000, 1:250000),c(5000,5000,2))
> system.time(rs1 <- apply(b, 3, rowSums))
   user  system elapsed 
  1.831   0.394   2.232 
> system.time(rs2 <- rowSums3d(b))
   user  system elapsed 
  1.134   0.183   1.320 
> system.time(rs3 <- sapply(1:dim(b)[3], function(i) rowSums(b[,,i])))
   user  system elapsed 
  1.556   0.073   1.636
> system.time(rs4 <- colSums(aperm(b, c(2,1,3))))
   user  system elapsed 
  0.860   0.103   0.966

因此，在我的系统上，该aperm()解决方案的出现速度略快：

> sessionInfo()
R version 2.12.1 Patched (2011-02-06 r54249)
Platform: x86_64-unknown-linux-gnu (64-bit)

但是，rowSums3d()没有给出与其他解决方案相同的答案：

> all.equal(rs1, rs2)
[1] "Mean relative difference: 0.01999992"
> all.equal(rs1, rs3)
[1] TRUE
> all.equal(rs1, rs4)
[1] TRUE
于 2011-02-27T22:53:55.000 回答
6

您可以将数组分成二维，计算行总和，然后按照您想要的方式将输出重新组合在一起。像这样：

rowSums3d <- function(a){
    m <- matrix(a,ncol=ncol(a))
    rs <- rowSums(m)
    matrix(rs,ncol=2)
}

> a <- array(c(1:250000, 1:250000),c(5000,5000,2))
> system.time(rowSums3d(a))
   user  system elapsed 
   1.73    0.17    1.96 
> system.time(apply(a, 3, rowSums))
   user  system elapsed 
   3.09    0.46    3.74
于 2011-02-27T20:49:50.160 回答
4

我不知道最有效的方法，但sapply似乎做得很好

a <- array(c(1:9, 1:9), c(3,3,2))
x1 <- sapply(1:dim(a)[3], function(i) rowSums(a[,,i]))
x1
     [,1] [,2]
[1,]   12   12
[2,]   15   15
[3,]   18   18

x2 <- apply(a, 3, rowSums)
all.equal(x1, x2)
[1] TRUE

这提供了如下速度改进：

> a <- array(c(1:250000, 1:250000),c(5000,5000,2))

> summary(replicate(10, system.time(rowSums3d(a))[3]))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  2.784   2.799   2.810   2.814   2.821   2.862 

> summary(replicate(10, system.time(apply(a, 3, rowSums))[3]))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  2.730   2.755   2.766   2.776   2.788   2.839 

> summary(replicate(10, system.time( sapply(1:dim(a)[3], function(i) rowSums(a[,,i])) )[3]))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.840   1.852   1.867   1.872   1.893   1.914

时间安排在：

# Ubuntu 10.10
# Kernal Linux 2.6.35-27-generic
> sessionInfo()
R version 2.12.1 (2010-12-16)
Platform: x86_64-pc-linux-gnu (64-bit)
于 2011-02-27T21:14:56.873 回答
1

如果您有一个多核系统，您可以编写一个简单的 C 函数并使用 Open MP 并行线程库。我为我的一个问题做了类似的事情，我在 8 核系统上得到了 8 倍的增长。该代码仍然可以在单处理器系统上运行，甚至可以在没有 OpenMP 的系统上编译，可能会到处乱加#ifdef _OPENMP。

当然，如果您知道这是大部分时间所花费的，那么它唯一值得做的事情。在优化之前对您的代码进行分析。

于 2011-02-27T20:40:22.843 回答