優化/並行化的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。有一些嘗試和錯誤很容易。