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

Question

我处理大型全球降水数据集（很细的空间分辨率）来计算使用的R.

SPEI package

我的标准化降水指数通常我的问题是指使用非常大的数据优化数据处理。我在其他帖子中发现了一些讨论（here,here和here fo实例），但没有一个与我的情况类似。

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

下面的代码示例：

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

：

在另一方面，如果你有机会到多核计算机（因为它是最有可能的情况下现在），你可以通过parallelising计算可知在计算时间的改进您可能需要将更高的值设置为ncpu以获得更高的速度增益，具体取决于您的计算机支持的线程数。 但是,sfApply，不会解决你的记忆问题与非常大的数据集。这是因为函数会将数据集分配给分配的cpus数量。由于系统的总内存不会改变，因此您的内存不足也会相同。

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

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。有一些尝试和错误很容易。