在R函数中并行化内部for循环答案

【问题标题】：Parallelize inner for loop in R function在R函数中并行化内部for循环
【发布时间】：2017-07-01 06:37:33
【问题描述】：

我正在 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)
}

编辑：如果它对基准测试有用，您可以在此处获取我的代码的运行版本

https://github.com/mariob6/progetto_bayes/blob/master/population_basedMC.R 需要这个源文件 https://github.com/mariob6/progetto_bayes/blob/master/simple_oscillator.R

【问题讨论】：

parallel 包中的parLapply 怎么样。

标签： r for-loop foreach parallel-processing

【解决方案1】：

只要有足够的任务供您的工作人员使用并且它们需要合理的计算时间，并行化内部循环就没有任何特别的问题。并行化外循环通常效率更高，但如果不能安全地并行化外循环，则并行化内循环。

重构代码以允许正确更新数据结构可能会很棘手。在这种情况下，我建议 foreach 循环的主体返回一个列表，该列表指示应如何更新“th0”的相应行，以便组合函数（在 master 上执行）可以执行实际更新。

这是一个演示此技术的精简示例：

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”块然后在外循环的每次迭代中重用它来改进。这可以显着减少开销，但也会变得更加复杂。

【讨论】：