Есть ли эффективный способ распараллеливать mapply?

Question

У меня много строк и на каждой строке я вычисляю uniroot нелинейной функции. У меня есть четырехъядерный компьютер Ubuntu, который не прекратил работу моего кода в течение двух дней. Неудивительно, что я ищу способы ускорить работу ;-)

После некоторых исследований я заметил, что в настоящее время используется только одно ядро, и распараллеливание - это то, что нужно сделать. Копая глубже, я пришел к выводу (может быть, неправильно?), Что пакет foreach на самом деле не предназначен для моей проблемы, потому что создается слишком много служебных данных (см., Например, SO). Хорошей альтернативой, по-видимому, является multicore для Unix-машин. В частности, функция pvec представляется наиболее эффективной после проверки страницы справки.

Однако, если я правильно понял, эта функция принимает только один вектор и соответственно разбивает его. Мне нужна функция, которая может быть парализована, но принимает несколько векторов (или вместо data.frame), как и функция mapply. Есть что-то там, что я пропустил?

Вот небольшой пример того, что я хочу сделать: (Обратите внимание, что здесь я использую пример plyr, потому что он может быть альтернативой базовой функции mapply, и у него есть опция параллелизации. Однако он медленнее в моем внедрение и внутренне, он вызывает foreach распараллеливание, так что я думаю, что это не поможет. Разве это правильно?)

library(plyr) 
library(foreach) 
n <- 10000 
df <- data.frame(P = rnorm(n, mean=100, sd=10), 
       B0 = rnorm(n, mean=40, sd=5), 
       CF1 = rnorm(n, mean=30, sd=10), 
       CF2 = rnorm(n, mean=30, sd=5), 
       CF3 = rnorm(n, mean=90, sd=8)) 

get_uniroot <- function(P, B0, CF1, CF2, CF3) { 

    uniroot(function(x) {-P + B0 + CF1/x + CF2/x^2 + CF3/x^3}, 
      lower = 1, 
      upper = 10, 
      tol = 0.00001)$root 

} 

system.time(x1 <- mapply(get_uniroot, df$P, df$B0, df$CF1, df$CF2, df$CF3)) 
    #user system elapsed 
    #0.91 0.00 0.90 
system.time(x2 <- mdply(df, get_uniroot)) 
    #user system elapsed 
    #5.85 0.00 5.85 
system.time(x3 <- foreach(P=df$P, B0=df$B0, CF1=df$CF1, CF2=df$CF2, CF3=df$CF3, .combine = "c") %do% { 
    get_uniroot(P, B0, CF1, CF2, CF3)}) 
    #user system elapsed 
    # 10.30 0.00 10.36 
all.equal(x1, x2$V1) #TRUE 
all.equal(x1, x3) #TRUE 

Кроме того, я пытался реализовать функцию Райан Томпсон chunkapply по ссылке SO выше (только избавилась из doMC, потому что я не смог установить его. Его пример работает даже после настройки его функции.), , но didn не получится. Однако, поскольку он использует foreach, я думал, что те же самые аргументы, упомянутые выше, применяются, поэтому я не пробовал слишком долго.

#chunkapply(get_uniroot, list(P=df$P, B0=df$B0, CF1=df$CF1, CF2=df$CF2, CF3=df$CF3)) 
#Error in { : task 1 failed - "invalid function value in 'zeroin'" 

PS: Я знаю, что я мог бы просто увеличить tol сократить количество шагов, которые необходимо найти uniroot. Тем не менее, я уже установил tol как можно больше.

Answer 1

Я бы использовал пакет parallel, который встроен в R 2.14 и работает с матрицами. Вы могли бы просто использовать mclapply как это:

dfm <- as.matrix(df) 
result <- mclapply(seq_len(nrow(dfm)), 
      function(x) do.call(get_uniroot,as.list(dfm[x,])), 
      mc.cores=4L 
     ) 
unlist(result) 

Это в основном делает то же самое mapply делает, но параллельно.

Но ...

Имейте в виду, что распараллеливание всегда рассчитывает на некоторые накладные расходы, а также. Как я объяснил в вопросе, на который вы ссылаетесь, параллельная работа только окупается, если ваша внутренняя функция вычисляет значительно больше, чем задействованные служебные данные. В вашем случае ваша функция uniroot работает довольно быстро. Затем вы можете рассмотреть возможность сократить кадр данных в больших кусках и объединить как mapply, так и mclapply. Возможный способ сделать это:

ncores <- 4 
id <- floor(
     quantile(0:nrow(df), 
       1-(0:ncores)/ncores 
     ) 
    ) 
idm <- embed(id,2) 

mapply_uniroot <- function(id){ 
    tmp <- df[(id[1]+1):id[2],] 
    mapply(get_uniroot, tmp$P, tmp$B0, tmp$CF1, tmp$CF2, tmp$CF3) 
} 
result <-mclapply(nrow(idm):1, 
        function(x) mapply_uniroot(idm[x,]), 
        mc.cores=ncores) 
final <- unlist(result) 

Это может потребоваться некоторые настройки, но это существенно нарушает ваш ФР в точности так, как много бит, поскольку есть ядра, и запустить mapply на каждом ядре.Для того, чтобы показать, что это работает:

> x1 <- mapply(get_uniroot, df$P, df$B0, df$CF1, df$CF2, df$CF3) 
> all.equal(final,x1) 
[1] TRUE 

Answer 2

Это не совсем наилучшая практика предложение, но значительное ускорение можно получить, идентифицирующий корень для всех параметров в моде «векторизованной». Например,

bisect <- 
    function(f, interval, ..., lower=min(interval), upper=max(interval), 
      f.lower=f(lower, ...), f.upper=f(upper, ...), maxiter=20) 
{ 
    nrow <- length(f.lower) 
    bounds <- matrix(c(lower, upper), nrow, 2, byrow=TRUE) 
    for (i in seq_len(maxiter)) { 
     ## move lower or upper bound to mid-point, preserving opposite signs 
     mid <- rowSums(bounds)/2 
     updt <- ifelse(f(mid, ...) > 0, 0L, nrow) + seq_len(nrow) 
     bounds[updt] <- mid 
    } 
    rowSums(bounds)/2 
} 

, а затем

> system.time(x2 <- with(df, { 
+  f <- function(x, PB0, CF1, CF2, CF3) 
+   PB0 + CF1/x + CF2/x^2 + CF3/x^3 
+  bisect(f, c(1, 10), PB0, CF1, CF2, CF3) 
+ })) 
    user system elapsed 
    0.180 0.000 0.181 
> range(x1 - x2) 
[1] -6.282406e-06 6.658593e-06 

против около 1.3s для применения uniroot отдельно к каждому. Это также комбинировало P и B0 в одно значение заблаговременно, так как именно они входят в уравнение.

Оценки конечного значения: +/- diff(interval) * (.5^maxiter) или около того. Более реалистичная реализация заменила бы биссектризацию линейной или квадратичной интерполяцией (как в ссылке, приведенной в ?uniroot), но тогда было бы сложнее организовать равномерную эффективную конвергенцию (и во всех случаях обработки ошибок).

Answer 3

это старая тема, но у вас есть parallel::mcmapply doc is here. не забудьте установить mc.cores в настройках. Обычно я использую mc.cores=parallel::detectCores()-1, чтобы освободить один процессор для операций ОС.

x4 <- mcmapply(get_uniroot, df$P, df$B0, df$CF1, df$CF2, df$CF3,mc.cores=parallel::detectCores()-1)