R中的并行优化

Posted 2023-03-05

【中文标题】R中的并行优化

【英文标题】：Parallelized optimization in R

【发布时间】：2013-03-02 02:24:57

【问题描述】：

我在具有 8 个多核处理器的 linux 机器上运行 R，并且有一个优化问题，我想通过并行化优化例程本身来加快速度。重要的是，这个问题涉及 (1) 多个参数，以及 (2) 固有的缓慢 模型运行。一个相当普遍的问题！

有谁知道用于这种情况的并行优化器？

更具体地说，像 nlm() 这样的求解器会在算法每次在参数空间中迈出一步时运行多个模型评估（每个参数值两个），因此并行化多个模型运行的实例将在这些情况下大大加快速度，当更多比几个参数值适合。

似乎使用包parallel 的代码可以编写成用户必须最少 修改代码以从使用nlm() 或optim() 移动到这个并行优化例程。也就是说，似乎可以基本上不做任何更改地重写这些例程，除了多次调用模型的步骤（这在基于梯度的方法中很常见）将并行完成。

理想情况下，像 nlmPara() 这样的代码会采用看起来像这样的代码

fit <- nlm(MyObjFunc, params0);

并且只需要很小的修改，例如，

fit <- nlmPara(MyObjFunc, params0, ncores=6);

想法/建议？

PS：我已采取措施加速这些模型运行，但由于各种原因它们运行缓慢（即我不需要关于加速模型运行的建议！;-)）。

【问题讨论】：

更多地阅读不同的优化器，看起来这种破解需要重写 C 代码（例如，重写 OPTIF9 例程的 C 端口以使用多个线程）或编写本机R 优化器以利用更高级别的并行化选项，例如 parallel、multicore、snow 等。

optimx/optimplus 包有很多优化算法的原生-R 版本：也许最容易从那里开始......？

感谢 Ben :-) optimx 允许您输入梯度函数。我会尝试一下，看看我是否不能只给它一个并行化的代码块，这应该可以解决问题。

我还有一些想法——可能有一些并行+记忆技巧？一些内置的optim() 优化器也采用可选的gr 参数

rgenoud 包对你有用吗？该包中的 genoud 函数采用 cluster 参数，该参数支持通过 snow 包进行并行计算，尽管您必须使用 cluster=rep('localhost', 6) 而不是 cluster=6。

【参考方案1】：

这是一个粗略的解决方案，至少有一些希望。非常感谢 Ben Bolker 指出许多/大多数优化例程允许用户指定的梯度函数。

具有更多参数值的测试问题可能会显示出更显着的改进，但在 8 核机器上，使用并行梯度函数的运行时间大约是串行版本的 70%。请注意，这里使用的粗略梯度近似似乎会减慢收敛速度，因此会增加一些时间。

## Set up the cluster
require("parallel");
.nlocalcores = NULL; # Default to "Cores available - 1" if NULL.
if(is.null(.nlocalcores))  .nlocalcores = detectCores() - 1; 
if(.nlocalcores < 1)  print("Multiple cores unavailable! See code!!"); return()
print(paste("Using ",.nlocalcores,"cores for parallelized gradient computation."))
.cl=makeCluster(.nlocalcores);
print(.cl)


# Now define a gradient function: both in serial and in parallel
mygr <- function(.params, ...) 
  dp = cbind(rep(0,length(.params)),diag(.params * 1e-8)); # TINY finite difference
  Fout = apply(dp,2, function(x) fn(.params + x,...));     # Serial 
  return((Fout[-1]-Fout[1])/diag(dp[,-1]));                # finite difference 


mypgr <- function(.params, ...)  # Now use the cluster 
  dp = cbind(rep(0,length(.params)),diag(.params * 1e-8));   
  Fout = parCapply(.cl, dp, function(x) fn(.params + x,...)); # Parallel 
  return((Fout[-1]-Fout[1])/diag(dp[,-1]));                  #



## Lets try it out!
fr <- function(x, slow=FALSE)  ## Rosenbrock Banana function from optim() documentation.
  if(slow)  Sys.sleep(0.1);    ## Modified to be a little slow, if needed.
  x1 <- x[1]
  x2 <- x[2]
  100 * (x2 - x1 * x1)^2 + (1 - x1)^2


grr <- function(x, slow=FALSE)  ## Gradient of 'fr'
  if(slow)  Sys.sleep(0.1);    ## Modified to be a little slow, if needed.
  x1 <- x[1]
  x2 <- x[2]
  c(-400 * x1 * (x2 - x1 * x1) - 2 * (1 - x1),
    200 *      (x2 - x1 * x1))


## Make sure the nodes can see these functions & other objects as called by the optimizer
fn <- fr;  # A bit of a hack
clusterExport(cl, "fn");

# First, test our gradient approximation function mypgr
print( mypgr(c(-1.2,1)) - grr(c(-1.2,1)))

## Some test calls, following the examples in the optim() documentation
tic = Sys.time();
fit1 = optim(c(-1.2,1), fr, slow=FALSE);                          toc1=Sys.time()-tic
fit2 = optim(c(-1.2,1), fr, gr=grr, slow=FALSE, method="BFGS");   toc2=Sys.time()-tic-toc1
fit3 = optim(c(-1.2,1), fr, gr=mygr, slow=FALSE, method="BFGS");  toc3=Sys.time()-tic-toc1-toc2
fit4 = optim(c(-1.2,1), fr, gr=mypgr, slow=FALSE, method="BFGS"); toc4=Sys.time()-tic-toc1-toc2-toc3


## Now slow it down a bit
tic = Sys.time();
fit5 = optim(c(-1.2,1), fr, slow=TRUE);                           toc5=Sys.time()-tic
fit6 = optim(c(-1.2,1), fr, gr=grr, slow=TRUE, method="BFGS");    toc6=Sys.time()-tic-toc5
fit7 = optim(c(-1.2,1), fr, gr=mygr, slow=TRUE, method="BFGS");   toc7=Sys.time()-tic-toc5-toc6
fit8 = optim(c(-1.2,1), fr, gr=mypgr, slow=TRUE, method="BFGS");  toc8=Sys.time()-tic-toc5-toc6-toc7

print(cbind(fast=c(default=toc1,exact.gr=toc2,serial.gr=toc3,parallel.gr=toc4),
            slow=c(toc5,toc6,toc7,toc8)))

【参考方案2】：

我是 R 包 optimParallel 的作者。它提供了optim() 的基于梯度的优化方法的并行版本。该包的主要功能是optimParallel()，其用法和输出与optim()相同。如下图所示，使用optimParallel()可以显着减少优化时间（p是参数个数）。

请参阅https://cran.r-project.org/package=optimParallel 和http://arxiv.org/abs/1804.11058 了解更多信息。

【参考方案3】：

由于您尚未接受答案，因此这个想法可能会有所帮助： 对于全局优化，包DEoptim() 具有用于并行优化的内置选项。好消息是，它易于使用且文档编写良好。

参考 http://www.jstatsoft.org/v40/i06/paper（目前关闭）

http://cran.r-project.org/web/packages/DEoptim/index.html

注意：差分 Evolglobal 优化器可能仍会遇到本地优化器。

【讨论】：

顺便说一下，如果您按照lhs 包中所述为您的初始群体使用拉丁超立方体采样，则使用 DE 进行的优化通常收敛得更快。

谢谢！我将不得不使用 DEoptim() 并看看它的比较

【参考方案4】：

我使用包 doSNOW 在 8 个内核上运行代码。 我可以复制和粘贴引用此包的代码部分。 希望对您有所帮助！

    # use multicore libraries
      # specify number of cores to use
    cores<- 8
      cluster <- makeCluster(cores, type="SOCK")
      registerDoSNOW(cluster)

      # check how many cores will be used
      ncores <- getDoParWorkers()
    print(paste("Computing algorithm for ", cores, " cores", sep=""))
      fph <- rep(-100,12)

      # start multicore cicle on 12  subsets
      fph <- foreach(i=1:12, .combine='c') %dopar% 
        PhenoRiceRun(sub=i, mpath=MODIS_LOCAL_DIR, masklocaldir=MASK_LOCAL_DIR, startYear=startYear, tile=tile, evismoothopt=FALSE)
      


  stopCluster(cluster) # check if gives error
  gc(verbose=FALSE)

【讨论】：

感谢 FraNut，但问题不在于如何在 R 中并行运行某些东西，我专门寻找一个基于梯度的优化例程来自动并行化梯度计算。请参阅上面我的问题下方的 cmets - 我认为 Ben 的建议导致了一个可行的解决方案

具体来说，optimx 允许用户自定义梯度函数，这仍然需要用户使用parallel（或doSNOW 等）编写代码，但想必这将是一个相当简单的解决方案解决问题。