嵌套最大化与在R中使用全局变量的需求并行

Question

我有一个带有两个嵌套优化的R代码。有一个内部函数。外部函数将某些参数传递给内部函数，内部函数对另一组参数进行优化。然后将这些参数发送到外部函数，该外部函数根据内部函数中估计的参数优化目标函数。然后，将外部函数的估计值传递给内部函数，内部函数会在内部函数中找到新的最佳参数集，然后将其传递给外部函数。重复这些循环，直到最小化外部循环中的目标函数为止。

该代码通过将内部参数设置为全局变量来工作，因此在外部循环中最大化代码后，代码会将这些全局变量传递给内部循环。

我想针对不同的数据集并行运行此过程。 I understand that I cannot use the global variables in parallel，并且我正在考虑在每个循环中保存具有不同文件名的文本文件：我将在外部循环的末尾保存一个带有参数值的文件，并在外部循环的开头重新打开它。但是，有没有更有效的方法可以做到这一点？我认为使用list无效。谢谢。

示例：

require(nloptr)
y = rnorm(100)
x = runif(100)*5

inner <- function(beta){
  mean((y-beta*x)^2)
}
outer <- function(alpha){
  if( !exists("storage") | is.null(storage$solution) ){
    beta <- runif(1)
  }else{
    beta <- storage$solution
  }

  sol.inner <-nloptr(x0 = beta,
                     eval_f = inner,
                     opts = list("algorithm"="NLOPT_LN_BOBYQA",
                                 ftol_rel = 1.e-6, ftol_abs = 1.e-7, xtol_rel = 1.e-6, 
                                 xtol_abs = 0, maxeval = 1000000)) 

  storage <-c()
  storage <<- append(storage,sol.inner)
  beta <- sol.inner$solution

  mean(x^2 - alpha* x + beta)^2
}

alpha0 <- runif(1)
storage <- c()
sol.outer <-nloptr(x0 = alpha0,
                   eval_f = outer,
                   opts = list("algorithm"="NLOPT_LN_BOBYQA",
                               ftol_rel = 1.e-6, ftol_abs = 1.e-7, xtol_rel = 1.e-6, 
                               xtol_abs = 0, maxeval = 1000000)) 
sol.outer

Answer 1

虽然非常简洁，但我一般不建议使用<<-运算符。如果要修改函数中的元素，以便在函数退出后可以使用它们，则建议您改用环境。

具有并行处理的东西是，正如在parallel包中实现的那样，每个线程/子代/子代都在其自己的会话中运行，这意味着它们不会相互交互。在这种情况下，您可以在每个后代过程中完成几乎所有您想做的事情。这是您要尝试执行的操作的示例：

# Simulating 4 random datasets
set.seed(131)
datasets <- replicate(4, {
  list(
    y = rnorm(100),
    x = runif(100)*5
  )
}, simplify = FALSE)


inner <- function(beta, x, y) mean((y-beta*x)^2)

outer <- function(alpha, storage, x, y) {

  if (!length(storage$solution))
    beta <- runif(1)
  else
    # Take the first value, which is the latest to be
    # stored (see below)
    beta <- storage$solution[[1]]

  sol.inner <- nloptr(
    x0     = beta,
    eval_f = inner,
    opts   = list(
      algorithm = "NLOPT_LN_BOBYQA",
      ftol_rel  = 1.e-6,
      ftol_abs  = 1.e-7,
      xtol_rel  = 1.e-6, 
      xtol_abs = 0,
      maxeval = 1000000
    ),
    y = y,
    x = x
  ) 

  # We can append the latest beta as a list
  storage$solution <- c(list(sol.inner$solution), storage$solution)
  beta    <-  sol.inner$solution

  mean(x^2 - alpha* x + beta)^2

}

# Parallel solution with PSOCKcluster --------------------
library(parallel)

# Setting up the cluster object
cl <- makePSOCKcluster(4)

# We need to export the objects we plan to use within
# each session this includes loading the needed packages
clusterExport(cl, c("outer", "inner"))
invisible(clusterEvalQ(cl, library(nloptr)))
invisible({
  clusterEvalQ(cl, {
    # Be careful about random numbers in parallel!
    # This example is not reproducible right now
    alpha0    <- runif(1)

    # This should be an environment, which is easier to handle
    storage   <- new.env() 
  })
})



# You can send data to the offspring sessions and 
# these will be evaluated in separate R sessions 
ans <- parLapply(cl, datasets, function(d) {

  # Making the variables available to the program
  y <- d$y
  x <- d$x

  sol.outer <- nloptr(
    x0     = alpha0,
    eval_f = outer,
    opts   = list(
      algorithm ="NLOPT_LN_BOBYQA",
      ftol_rel  = 1.e-6,
      ftol_abs  = 1.e-7,
      xtol_rel  = 1.e-6,
      xtol_abs  = 0,
      maxeval   = 1000000
    ),
    x = d$x,
    y = d$y,
    # Passing the environment as an extra
    #  argument to the function
    storage = storage
  ) 

  list(
    sol     = sol.outer,
    storage = storage
  )

})

# Stopping the R sessions
stopCluster(cl)

# Checking out the storage vectors
lapply(ans, function(x) unlist(x$storage$solution))
#> [[1]]
#>  [1] -0.04112901 -0.04112901 -0.04112901 -0.04112901 -0.04112901
#>  [6] -0.04112901 -0.04112901 -0.04112901 -0.04112901 -0.04112901
#> [11] -0.04112901 -0.04112901 -0.04112901 -0.04112901 -0.04112901
#> [16] -0.04112901
#> 
#> [[2]]
#>  [1] -0.06877397 -0.06877397 -0.06877397 -0.06877397 -0.06877397
#>  [6] -0.06877397 -0.06877397 -0.06877397 -0.06877397 -0.06877397
#> [11] -0.06877397 -0.06877397 -0.06877397 -0.06877397 -0.06877397
#> [16] -0.06877397 -0.06877397
#> 
#> [[3]]
#>  [1] 0.004505708 0.004505708 0.004505708 0.004505708 0.004505708
#>  [6] 0.004505708 0.004505708 0.004505708 0.004505708 0.004505708
#> [11] 0.004505708 0.004505708 0.004505708 0.004505708 0.004505708
#> [16] 0.004505708 0.004505708
#> 
#> [[4]]
#>  [1] -0.02001445 -0.02001445 -0.02001445 -0.02001445 -0.02001445
#>  [6] -0.02001445 -0.02001445 -0.02001445 -0.02001445 -0.02001445
#> [11] -0.02001445 -0.02001445 -0.02001445 -0.02001445 -0.02001445
#> [16] -0.02001445

^{由reprex package（v0.3.0）在2019-11-20创建}

这里要注意的一件事是，我修改了您的函数，以便显式传递参数，因此在这种情况下，我们将不处理范围界定。这通常更安全，并且R的最新版本足够聪明，可以避免在传递给函数时复制对象。

最后要指出的是，如果您的数据集很大，最好将其实际加载到子代会话中，以避免重复存储（通常，如果使用makeForkCluster，这最后一点不是问题，但是仅适用于基于UNIX的系统）。