Распараллеливание больших симуляций над сеткой в ​​R

Question

Я запускаю серию больших симуляций по сетке. Я выполняю моделирование по строкам, и я обнаружил, что мои функции выборки являются узким местом. Я пытался использовать библиотеки foreach и doMC для ускорения процесса, но обнаружил, что либо параллельный метод медленнее, либо я не смог написать функцию, которая была бы правильно интерпретирована foreach.

Глядя на некоторые другие сообщения, кажется, что мой подход с использованием 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 ) )

Во втором случае я пытаюсь выбрать из набора эмпирических распределений, которые индексируются по столбцу в матрице. Значение элемента grid-row соответствует столбцу для выборки.

#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. Спасибо за вашу помощь.

Answer 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 

Answer 2

Вот небольшое улучшение вашей первой симуляции. Больше 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