R - 加速数组维度上的循环答案

【问题标题】：R - Speeding up loop over array dimensionsR - 加速数组维度上的循环
【发布时间】：2020-04-12 14:27:04
【问题描述】：

我正在处理具有维度的数组

[1] 290 259  55   4

对于最后三个维度的每次重复，我想对第一个维度的 290 个元素执行滚动平均，将元素数量减少到 289 个。最后，我需要创建一个包含更新值的数据框。

下面的代码实现了我所需要的，但是运行时间很长（其实我要在结束前打断它）。

library(zoo)

# Generate random data with same dimensions as mine
my.array <- array(1:16524200, dim=c(290,259,55,4))

# Get dimension sizes
dim2 <- dim(my.array)[2]
dim3 <- dim(my.array)[3]
dim4 <- dim(my.array)[4]

# Pre-allocate data frame to be used within the loop
df2 <- data.frame()

# Loop over dimensions
for (i in 1:dim4) {
  for (j in 1:dim3) {
    for (k in 1:dim2) {

      # Take rolling average
      u <- rollapply(my.array[,k,j,i], 2, mean)

      # Assemble data frame
      df1 <- data.frame(time=i, level=j, lat=k, wind=u)
      df2 <- rbind(df2, df1)

    }
  }
}
# Very slow, and uses only one machine core

我觉得可以通过使用矢量化甚至某种并行性来改善这段代码的处理时间，但我不知道如何。

有什么建议可以让这段代码更高效吗？

【问题讨论】：

永远不要迭代地构建 data.frames。每次调用rbind 时，它都会将整个帧 复制到一个新对象中并覆盖df2。它可能适用于几十个，但是（如您所见）这可怕地。
@r2evans 说得通，但是......还有什么替代方案？
一般来说，something <- lapply(list_of_stuff, somefunc)，然后是do.call(rbind, something)（虽然这个问题需要多一点）。

标签： r performance for-loop multidimensional-array

【解决方案1】：

在使用data.table 计算滚动平均值之前，flatten the multidimensional array first 的另一个选项

library(data.table)
system.time({
    ans <- setDT(as.data.frame.table(my.array))[
        , .(wind=((Freq + shift(Freq)) / 2)[-1L]), 
        .(time=Var4, level=Var3, lat=Var2)]
    cols <- c("time", "level", "lat")
    ans[, (cols) := lapply(.SD, function(x) match(x, unique(x))), .SDcols=cols]
})
ans

输出：

          time level lat       wind
       1:    1     1   1        1.5
       2:    1     1   1        2.5
       3:    1     1   1        3.5
       4:    1     1   1        4.5
       5:    1     1   1        5.5
      ---                          
16467216:    4    55 259 16524195.5
16467217:    4    55 259 16524196.5
16467218:    4    55 259 16524197.5
16467219:    4    55 259 16524198.5
16467220:    4    55 259 16524199.5

时间安排：

   user  system elapsed 
   4.90    1.16    5.66

为了比较：

library(zoo)
system.time({
    as.data.frame.table(apply(my.array, c(2,3,4), rollmean, 2))  
})
#   user  system elapsed 
#  21.89    0.63   22.51

【讨论】：

谢谢，我忘记了data.table 的速度有多快！我将此标记为官方答案，因为我的数据可能比我给出的示例大 10 倍，使用此答案可以为我节省大量时间。
我正在努力将结果数据表转换为维度为 289 259 55 4 的数组（类似于原始数据表）。你碰巧对此有任何提示吗？
我猜它使用数组？也许发布另一个问题？我现在没有电脑

【解决方案2】：

apply() 适用于任意数量的维度，因此您可以使用 as.data.frame.table() 中包含的以下内容更快地获得相同的结果，从而有效地将输出从数组转换为数据框：

library(zoo)
df <- as.data.frame.table(apply(my.array, c(2,3,4), rollmean, 2))

并非绝对必要，但可以整理以匹配您的原始输出：

idx <- sapply(df, is.factor)
df[idx] <- sapply(df[idx], as.integer)

df <- setNames(df[c(4,3,2,5)], c("time", "level", "lat", "wind"))

检查结果是否相同：

identical(df2, df)
[1] TRUE

【讨论】：

我忘记了 as.data.frame.table，它在加快速度方面做得非常好（在初始测试中超过 2 倍）。

【解决方案3】：

在前面，您正在遭受 R 的地狱的第二圈 (https://www.burns-stat.com/pages/Tutor/R_inferno.pdf)：正在生长的物体。每次调用rbind 时，它都会制作帧的完整副本，执行 r 绑定，然后覆盖原始变量名称的完整副本。因此，虽然它可能在前几十个没有明显减速的情况下工作，但它会减速超过 100 次左右......而你正在这样做 56,980 次。

通常最好将事物处理成list，然后在整个列表的末尾执行一次rbind，如do.call(rbind, list_of_frames)。诚然，你仍然可能面临做一些可能很难的事情的计算挑战......幸运的是，zoo 与窗口操作的效率差不多，而且这不是不可能的难。

我将演示一个显着减少的问题集（因为我认为我们查看 16M 或 1.5M 迭代并不重要。

my.array <- array(1:1502200, dim=c(290,259,5,4))
eg <- do.call(expand.grid, lapply(dim(my.array)[-1], seq_len))
dim(eg)
# [1] 5180    3
head(eg)
#   Var1 Var2 Var3
# 1    1    1    1
# 2    2    1    1
# 3    3    1    1
# 4    4    1    1
# 5    5    1    1
# 6    6    1    1

system.time({
  list_of_frames <- Map(function(i,j,k) {
    u <- zoo::rollapply(my.array[,i,j,k], 2, mean)
    data.frame(i, j, k, wind = u)
  }, eg[[1]], eg[[2]], eg[[3]])
})
#    user  system elapsed 
#    5.79    0.00    5.80 
head(list_of_frames[[5]])
#   i j k   wind
# 1 5 1 1 1161.5
# 2 5 1 1 1162.5
# 3 5 1 1 1163.5
# 4 5 1 1 1164.5
# 5 5 1 1 1165.5
# 6 5 1 1 1166.5

system.time({
  out <- do.call(rbind, list_of_frames)
})
#    user  system elapsed 
#    0.50    0.03    0.53 
nrow(out)
# [1] 1497020
rbind(head(out), tail(out))
#           i j k      wind
# 1         1 1 1       1.5
# 2         1 1 1       2.5
# 3         1 1 1       3.5
# 4         1 1 1       4.5
# 5         1 1 1       5.5
# 6         1 1 1       6.5
# 1497015 259 5 4 1502194.5
# 1497016 259 5 4 1502195.5
# 1497017 259 5 4 1502196.5
# 1497018 259 5 4 1502197.5
# 1497019 259 5 4 1502198.5
# 1497020 259 5 4 1502199.5

解释：

do.call(expand.grid, ...) 正在创建一个包含您需要的所有 i,j,k 组合的框架，动态地在您的数组的维度上。
Map(f, is, js, ks) 使用 is、js 和 ks 的第一个参数运行函数 f（此项目符号的概念），因此 Map 看起来像：
```
f(is[1], js[1], ks[1])
f(is[2], js[2], ks[2])
f(is[3], js[3], ks[3])
# ...
```
然后我们使用do.call(rbind, ...) 将它们组合在一个调用中。我们真的必须在这里使用do.call，因为这个调用类似于
```
rbind(list_of_frames[[1]], list_of_frames[[2]], ..., list_of_frames[[5180]])
```
（如果你想写出这个版本，交给你）。

【讨论】：