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

Question

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

[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

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

有什么建议可以使此代码更有效吗？

Answer 1

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

Answer 2

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

通常最好将事物处理成 alist然后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)运行功能f与每个的第一个参数is，js和ks（名义为这个子弹），使地图看起来是这样的：

  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]])
  

  （如果您想写出此版本，则交给您）。