Есть ли эффективный способ распараллеливания 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, потому что я не смог ее установить. Однако его пример работает, даже после настройки его функции.), но не заставил его работать. Однако, поскольку он использует 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)

Это может потребовать некоторой настройки, но по сути это разбивает ваш df ровно на столько бит, сколько имеется ядер, и запускает 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,3 с для применения uniroot отдельно для каждого. Это также объединило P и B0 в одно значение раньше времени, поскольку именно так они вводят уравнение.

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

Answer 3

это старая тема, но, к вашему сведению, parallel::mcmapply документ здесь . не забудьте установить 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)