优化/并行化 R - 处理在 R 中计算 SPI 的大型数据集

Question

我正在处理一个大型全球降水数据集（具有非常精细的空间分辨率），以使用 R 中的SPEI package 计算标准化降水指数。

我的问题一般是指使用非常大的数据优化数据处理。我在其他帖子中发现了一些讨论（here、here 和 here fo 实例），但似乎没有什么与我的情况相似。

我的输入是一个矩阵，其中包含超过 20 年的每月观测（>20*12 行）降水时间序列 > 1,000,000 个点（列）。 SPI 的计算为每个时间序列执行一系列步骤，并将指数计算为与中位数的标准差。 输出是一个列表，其结果矩阵 ($fitted) 与输入矩阵的大小相同。

这里是代码示例：

require(SPEI)

#generating a random values matrix 
data<-replicate(200, rnorm(240))
# erasing negative values
data[data<=0]=0

spi6 <- spi(data, 6, kernel = list(type = 'rectangular', shift = 0), distribution = 'PearsonIII', fit = 'ub-pwm', na.rm = FALSE, ref.start=NULL, ref.end=NULL, x=FALSE, params=NULL)

#testing the results
plot(spi6$fitted[,67])
#taking my results out
results <- t(spi6$fitted)

此脚本运行良好，但如果我增加点数（在本例中为列），处理时间会成倍增加。直到遇到内存不足问题：

Warning messages:
1: In std[ff, s] <- qnorm(cdfpe3(acu.pred[ff], p3par)) :
  Reached total allocation of 16253Mb: see help(memory.size)
2: In std[ff, s] <- qnorm(cdfpe3(acu.pred[ff], p3par)) :
  Reached total allocation of 16253Mb: see help(memory.size)
3: In NextMethod("[<-") :
  Reached total allocation of 16253Mb: see help(memory.size)
4: In NextMethod("[<-") :
  Reached total allocation of 16253Mb: see help(memory.size)
5: In std[ff, s] <- qnorm(pze + (1 - pze) * pnorm(std[ff, s])) :
  Reached total allocation of 16253Mb: see help(memory.size)
6: In std[ff, s] <- qnorm(pze + (1 - pze) * pnorm(std[ff, s])) :
  Reached total allocation of 16253Mb: see help(memory.size)
7: In NextMethod("[<-") :
  Reached total allocation of 16253Mb: see help(memory.size)
8: In NextMethod("[<-") :
  Reached total allocation of 16253Mb: see help(memory.size)

如何拆分我的输入矩阵（或在输入矩阵上拆分过程）以并行处理列向量组（每个列向量都是特定点的完整时间序列），而不会丢失信息（或弄乱我的数据）？ 谢谢。

Answer 1

我发现处理时间是线性的，而不是指数的：

system.time(spi(data[,1], 6))
system.time(spi(data[,1:10], 6))
system.time(spi(data[,1:100], 6))

如果您看到指数增长，则可能是由于 RAM 分配过多的问题。

一个简单的解决方案是在矩阵上拆分计算：

spi6 <- data*NA
system.time(
  for (i in 1:100) spi6[,i] <- spi(data[,i], 6)$fitted
)

或者，以类似的效率：

system.time(
  spi6 <- apply(data[,1:100], 2, function(x) spi(x, 6)$fitted)
)

如您所见，计算时间与将整个矩阵作为输入提供给spi() 函数的原始选项非常相似。 但是，如果您遇到内存问题，这可能会解决它们。

另一方面，如果您可以使用多核计算机（现在很可能是这种情况），您可能会通过并行计算看到计算时间的改进：

library(snowfall)
sfInit(parallel=TRUE, cpus=2)
sfLibrary(SPEI)

system.time(
  spi6 <- sfApply(data[,1:100], 2, function(x) spi(x, 6)$fitted)
)

sfStop()

您可能希望将更高的值设置为ncpu 以获得更高的速度增益，具体取决于您的计算机支持的线程数。 但是，sfApply 无法解决您对非常大的数据集的内存问题。之所以如此，是因为该函数在分配的 CPU 数量之间拆分数据集。由于系统的总内存不变，所以同样会耗尽内存。

解决方案是合并这两种方法：拆分数据集，然后并行化。像这样的：

data <- replicate(10000, rnorm(240))

sfInit(parallel=TRUE, cpus=10)
sfLibrary(SPEI)

spi6 <- data*NA
for (i in 0:9) {
    chunk <- (i*1000)+(1:1000)
    spi6[,chunk] <- sfApply(data[,chunk], 2, function(x) spi(x, 6)$fitted)
} 

sfStop()

现在，您只需要找到块的最大大小，即您传递给sfApply 的数据原始数量，以免溢出您的可用 RAM。通过一些尝试和错误，这很容易。