2013-02-12 17 views
1

我正在網格上運行一系列大型仿真。我正在逐行實施模擬,並且發現我的採樣功能是一個瓶頸。我試圖用foreach和doMC庫來加速這個過程,但是我發現並行方法比較慢,或者我一直無法編寫一個可以被foreach正確解釋的函數。在R中的網格中進行大型仿真的並行化

綜觀其他一些帖子,看來使用的foreach我的做法可能會誤導該職位我試圖人數大大超過可用處理器的數量。我想知道在我的情況下,人們是否會對如何最好地實現並行化提出一些建議。我的模擬通常有兩種類型。首先,我計算一個矩陣,該矩陣包含我正在處理的網格行內的每個元素的採樣間隔(行)。然後使用runif進行採樣(在真實模擬中,我的行包含〜9000個單元,並且我正在執行10000次模擬)。

#number of simulations per element 
n = 5 

#Generate an example sampling interval. 
m.int1 <- matrix (seq (1, 20, 1), ncol=10, nrow=2) 

#Define a function to sample over the interval defined in m.int1 
f.rand1 <- function(a) { 
return (runif (n, a[1], a[2])) 
} 

#run the simulation with each columns corresponding to the row element and rows 
#the simultions. 
sim1 <- round(apply (m.int1, 2, f.rand1)) 

在第二種情況下,我試圖從一組矩陣中​​按列索引的經驗分佈中抽樣。網格行元素的值對應於要採樣的列。

#number of simulations per element 
n = 5 

#generate a vector represeting a row of grid values 
v.int2 <- round(runif(10,1,3)) 

#define matrix of data that contains the distributions to be sampled. 
m.samples<-cbind(rep(5,10),rep(4,10),rep(3,10)) 

f.sample <- function(a) { 
return (sample (m.samples [ ,a], n,)) 
} 

#Sample m.samples indexed by column number. 
sim2<- sapply(v.int2,f.sample) 

在第二個例子,我能夠利用的foreach()%dopar%並行運行,但仿真花基本上長於串行代碼。在上面的第一個例子中,我無法寫出一個正確的函數來利用foreach並行化。我將把我在第二種情況下使用的代碼放在一起來展示我的想法 - 但現在我意識到我的方法在開銷方面太昂貴了。

library(foreach) 
library(doMC) 
registerDoMC(2) 

n = 5 

#Sample m.samples indexed by column number using parallel method. 
sim2.par <- foreach (i = 1 : length (v.int2), 
    .combine="cbind") %dopar% sample ( 
    m.samples [ , v.int2 [i] ] , n) 

我會很感激的做法提出了一些建議(和一些代碼!),這將幫助我有效地利用並行。再次,我正在處理的行通常包含約9000個元素,我們正在對每個元素進行10000次模擬。所以我的輸出仿真矩陣一般在10000 X 9000的數量級。感謝您的幫助。

+0

在情況下,個人迭代短,開銷可能會非常昂貴,相對來說。這就是爲什麼在許多內核上運行時沒有看到任何提升。換句話說,工作速度如此之快,以至於溝通比實際工作花費更多時間。 – 2013-02-12 19:52:21

回答

1

這是你第一次模擬略有改善。更大的n可能在運行時產生更大的收益。

> n <- 1000 
> m.int1 <- matrix (seq (1, 20, 1), ncol=10, nrow=2) 
> f.rand1 <- function(a) { 
+ return(runif(n, a[1], a[2])) 
+ } 
> system.time(x1 <- replicate(n, round(apply(m.int1, 2, f.rand1)))) 
    user system elapsed 
    2.84 0.06 2.95 
> system.time(x2 <- replicate(n, matrix(round(runif(n*10, min = m.int1[1, ], max = m.int1[2, ])), ncol = 10, byrow = TRUE))) 
    user system elapsed 
    2.48 0.06 2.61 
> head(x1[,,1]) 
    [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] 
[1,] 1 4 5 7 10 12 13 16 17 20 
[2,] 1 3 6 7 10 11 13 16 17 19 
[3,] 1 3 6 7 10 12 14 16 18 20 
[4,] 2 4 5 7 9 12 14 16 17 19 
[5,] 1 4 5 7 10 12 14 16 17 20 
[6,] 1 4 6 8 9 11 13 15 18 20 
> head(x2[,,1]) 
    [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] 
[1,] 2 4 6 7 9 12 14 16 17 20 
[2,] 1 3 6 8 10 12 14 15 18 20 
[3,] 2 4 5 7 9 11 13 15 17 20 
[4,] 2 3 5 7 9 11 14 15 17 19 
[5,] 2 3 6 7 9 12 13 16 17 20 
[6,] 2 4 6 7 10 12 14 16 17 20 
1

嘗試使用此代替兩步過程。它跳過apply步:

f.rand2 <- function(a) { 
    matrix(runif (n*ncol(a), rep(a[1,], n) , rep(a[2,], n)), nrow=ncol(a)) 
        } 

f.rand2(m.int1) 
      [,1]  [,2]  [,3]  [,4]  [,5] 
[1,] 1.693183 1.404336 1.067888 1.904476 1.161198 
[2,] 3.411118 3.852238 3.621822 3.969399 3.318809 
[3,] 5.966934 5.466153 5.624387 5.646181 5.347473 
[4,] 7.317181 7.106791 7.403022 7.442060 7.161711 
[5,] 9.491231 9.656023 9.518498 9.569379 9.812931 
[6,] 11.843074 11.594308 11.706276 11.744094 11.994256 
[7,] 13.375382 13.599407 13.416135 13.634053 13.539246 
[8,] 15.948597 15.532356 15.692132 15.442519 15.627716 
[9,] 17.856878 17.208313 17.804288 17.875288 17.232867 
[10,] 19.214776 19.689534 19.732680 19.813718 19.866297 

對於我來說,切割時間縮短了一半:

> system.time(x1 <- replicate(n, round(apply(m.int1, 2, f.rand1)))) 
    user system elapsed 
    1.088 0.470 1.550 

> system.time(x1 <- replicate(n, f.rand2(m.int1))) 
    user system elapsed 
    0.559 0.256 0.811