Распараллеливание внутренней петли в функции R

Question

Я строю пробоотборник Марковской цепи Монте-Карло в R, идея в том, что у меня есть популяция цепей, которые обновляются параллельно (независимо) и в какой-то момент они взаимодействуют. Что я хочу сделать, это распараллеливать независимое обновление, так как весь код занимает много времени. Однако foreach, похоже, не уместен, так как он должен возвращать значения, мне просто нужно их обновить, у кого-нибудь возникла эта проблема и придумал умное решение?

population_MCMC <- function(niter, burnin,thin ,th0, T_N ,Sig, y0, log_target) 
{ 
    # th0 will be updated at each step, th will contain the output of interest (that is, when T_N = 1) 
    th <- matrix(nrow= ceiling((niter-burnin)/thin), ncol=2) 

    nacp = 0 # number of accepted moves 

    for(i in 1:(niter)) 
    { 
     for(j in 1:length(T_N)){ # <-- THIS IS THE FOR LOOP I WANT TO PARALLELIZE! 

      #this is the local change 
      delta = as.vector(rmvnorm(1, mean = th0[j,], sig = Sig)) 
      lacp <- log_target(th = delta, y_obs = y_obs, y0 = y0, t_n=T_N[j]) 
      lacp <- lacp - log_target(th = th0[j,], y_obs = y_obs, y0 = y0, t_n=T_N[j]) 
      #cat(lacp,"\n") 
      lgu <- log(runif(1)) 
      if(!(is.na(lacp)) & lgu < lacp) 
      { 
      th0[j,] <- delta 
      nacp = nacp + 1 
      } 
     } 

    # Try to exchange theta_l and theta_m where m = l+1 or m= l-1 if l=! 1 and l=! length(T_N) 
    ..... some other markovian stuff ..... 

    if(i>burnin & (i-burnin)%%thin==0){ 
     th[(i-burnin)/thin,] = th0[length(T_N),] 
    } 

    if(i%%1000==0) cat("*** Iteration number ", i,"/", niter, "\n") 
    } 
    cat("Acceptance rate =", nacp/niter, "\n") 
    return(th) 
} 

EDIT: если он может быть полезен для сравнительного анализа, вы можете получить работающую версию моего кода здесь

https://github.com/mariob6/progetto_bayes/blob/master/population_basedMC.R требует этого исходного файла https://github.com/mariob6/progetto_bayes/blob/master/simple_oscillator.R

Answer 1

Существует не какая-либо конкретная проблема с параллелизацией внутреннего цикла, если для ваших работников достаточно задач, и они занимают разумное количество времени вычисления. Параллельность внешнего цикла обычно более эффективна, но если внешний цикл нельзя безопасно распараллелить, вы распараллеливаете внутреннюю.

Реструктуризация кода, позволяющего соответствующим образом обновлять структуры данных, может быть сложной задачей. В этом случае я предлагаю, чтобы тело цикла foreach возвращало список, который указывает, как должна обновляться соответствующая строка «th0», чтобы функция объединения (которая выполняется на главном компьютере) могла выполнять фактическое обновление.

Вот разобранный пример, который демонстрирует эту технику:

example <- function() { 
    th0 <- matrix(0, 4, 4) 
    nacp <- 0 

    updatemaster <- function(ignore, ...) { 
    for (r in list(...)) { 
     if (! is.null(r$delta)) { 
     cat('updating row', r$j, '\n') 
     th0[r$j,] <<- r$delta 
     nacp <<- nacp + 1 
     } else { 
     cat('not updating row', r$j, '\n') 
     } 
    } 
    NULL # combine function called strictly for side effect 
    } 

    for (i in 1:2) { 
    ignore <- 
     foreach(j=1:nrow(th0), .combine='updatemaster', 
       .init=NULL, .multicombine=TRUE) %dopar% { 
     delta <- rnorm(ncol(th0)) 
     if (rnorm(1) >= 0) 
      list(j=j, delta=delta) 
     else 
      list(j=j, delta=NULL) 
     } 

    print(th0) 
    print(nacp) 
    } 
} 

suppressMessages(library(doSNOW)) 
cl <- makeSOCKcluster(2) 
registerDoSNOW(cl) 
example() 

Это может быть улучшено только отправка куска «TH0» к каждому из рабочих один раз, а затем многократно использовать его для каждой итерации внешнего петля. Это может значительно снизить накладные расходы, но также значительно усложняется.