R関数のinner forループを並列化する

Question

私はMarkov Chain Monte CarloサンプラーをRで構築しています。その背景には、並列に（独立して）更新され、ある時点で相互作用するChainsの集団があります。 私がしたいのは、コード全体が実行に時間がかかるため、独立した更新を並列化することです。 しかし、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"の対応する行がどのように更新されるべきかを示すリストを返し、結合機能（マスター上で実行する）が実際の更新を実行できるようにすることをお勧めします。これは一度だけの労働者のそれぞれに「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() 

：

は、ここでは、この技術を実証ストリップダウン例ですループ。オーバーヘッドを大幅に減らすことができますが、かなり複雑になります。