2013-02-12 51 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