optimization - R:FAST 多变量优化包？

Question

我正在寻找 4 个变量的标量函数的局部最小值，并且我对变量有范围约束(“框约束”)。函数导数没有封闭形式，因此需要解析导数函数的方法是不可能的。我已经用 optim 尝试了几个选项和控制参数功能，但所有这些似乎都很慢。具体来说，他们似乎在调用我的(R 定义的)目标函数之间花费了大量时间，所以我知道瓶颈不是我的目标函数，而是调用我的目标函数之间的“思考”。我查看了 CRAN 任务 View 以进行优化并尝试了其中的几个选项( DEOptim 来自 RcppDE 等)，但它们似乎都没有任何好处。我很想试试 nloptr包(NLOPT 库的 R 包装器)，但它似乎不适用于 Windows。

我想知道，是否有我可能会遗漏的人们使用的任何好的、快速的优化包？理想情况下，这些将采用围绕良好 C++/Fortran 库的薄包装器的形式，因此只有最少的纯 R 代码。 (虽然这不应该是相关的，但我的优化问题是在尝试通过最小化某个拟合优度来将 4 参数分布拟合到一组值时出现的)。

在过去，我发现 R 的优化库非常慢，最终编写了一个薄的 R 包装器，调用了商业优化库的 C++ API。那么最好的图书馆一定是商业图书馆吗？

更新。 这是我正在查看的代码的简化示例:

###########
## given a set of values x and a cdf, calculate a measure of "misfit":
## smaller value is better fit
## x is assumed sorted in non-decr order;
Misfit <- function(x, cdf) {
  nevals <<- nevals + 1
  thinkSecs <<- thinkSecs + ( Sys.time() - snapTime)
  cat('S')
  if(nevals %% 20 == 0) cat('\n')
  L <- length(x)
  cdf_x <- pmax(0.0001, pmin(0.9999, cdf(x)))
  measure <- -L - (1/L) * sum( (2 * (1:L)-1 )* ( log( cdf_x ) + log( 1 - rev(cdf_x))))
  snapTime <<- Sys.time()
  cat('E')
  return(measure)  
}
## Given 3 parameters nu (degrees of freedom, or shape), 
## sigma (dispersion), gamma (skewness),
## returns the corresponding 4-parameter student-T cdf parametrized by these params
## (we restrict the location parameter mu to be 0).
skewtGen <- function( p ) {
  require(ghyp)
  pars = student.t( nu = p[1], mu = 0, sigma = p[2], gamma = p[3] )
  function(z) pghyp(z, pars)
}

## Fit using optim() and BFGS method
fit_BFGS <- function(x, init = c()) {
  x <- sort(x)
  nevals <<- 0
  objFun <- function(par) Misfit(x, skewtGen(par))
  snapTime <<- Sys.time() ## global time snap shot
  thinkSecs <<- 0 ## secs spent "thinking" between objFun calls
  tUser <- system.time(
              res <- optim(init, objFun,
                           lower = c(2.1, 0.1, -1), upper = c(15, 2, 1),
                           method = 'L-BFGS-B',
                           control = list(trace=2, factr = 1e12, pgtol = .01 )) )[1]
  cat('Total time = ', tUser, 
      ' secs, ObjFun Time Pct = ', 100*(1 - thinkSecs/tUser), '\n')
  cat('results:\n')
  print(res$par)
}

fit_DE <- function(x) {
  x <- sort(x)
  nevals <<- 0
  objFun <- function(par) Misfit(x, skewtGen(par))
  snapTime <<- Sys.time() ## global time snap shot
  thinkSecs <<- 0 ## secs spent "thinking" between objFun calls
  require(RcppDE)
  tUser <- system.time(
              res <- DEoptim(objFun,
                             lower = c(2.1, 0.1, -1),
                             upper = c(15, 2, 1) )) [1]
  cat('Total time = ',             tUser,
      ' secs, ObjFun Time Pct = ', 100*(1 - thinkSecs/tUser), '\n')
  cat('results:\n')
  print(res$par)
}

让我们生成一个随机样本:

set.seed(1)
# generate 1000 standard-student-T points with nu = 4 (degrees of freedom)
x <- rt(1000,4)

第一次使用 fit.tuv (对于“T UniVariate”)函数在 ghyp 中包 - 这使用最大似然期望最大化(E-M)方法。这是邪恶的快!

require(ghyp)
> system.time( print(unlist( pars <- coef( fit.tuv(x, silent = TRUE) ))[c(2,4,5,6)]))
         nu          mu       sigma       gamma 
 3.16658356  0.11008948  1.56794166 -0.04734128 
   user  system elapsed 
   0.27    0.00    0.27 

现在我试图以不同的方式拟合分布:通过使用标准 optim() 最小化上面定义的“失配”度量。函数在基 R 中。请注意，结果通常不会相同。我这样做的原因是将这两个结果在整个类别的情况下进行比较。我将上面的最大似然估计作为此优化的起点。

> fit_BFGS( x, init = c(pars$nu, pars$sigma, pars$gamma) )
N = 3, M = 5 machine precision = 2.22045e-16
....................
....................
.........
iterations 5
function evaluations 7
segments explored during Cauchy searches 7
BFGS updates skipped 0
active bounds at final generalized Cauchy point 0
norm of the final projected gradient 0.0492174
final function value 0.368136

final  value 0.368136 
converged
Total time =  41.02  secs, ObjFun Time Pct =  99.77084 
results:
[1] 3.2389296 1.5483393 0.1161706

我也试图适应 DEoptim()但它跑得太久了，我不得不杀死它。从上面的输出可以看出，99.8% 的时间都归因于目标函数!所以德克和迈克在下面的评论中是正确的。我应该更仔细地估计我的目标函数所花费的时间，打印点不是一个好主意!此外，我怀疑 MLE(E-M) 方法非常快，因为它对对数似然函数使用分析(封闭形式)。

Answer 1

对于您的问题，最大似然估计器在任何语言中总是比全局优化器更快。

无论采用何种算法，全局优化器通常都会将一些随机跳转与局部最小化例程结合起来。不同的算法可能会从种群(遗传算法)、退火、迁移等方面讨论这一点，但它们在概念上都是相似的。

在实践中，这意味着如果你有一个平滑的函数，其他一些优化算法可能是最快的。您的问题函数的特征将决定它是二次、线性、圆锥形还是其他类型的优化问题，其中存在精确(或接近精确)的解析解，或者您是否需要应用全局优化器那必然更慢。

通过使用 ghyp，您是说您的 4 变量函数产生了一个可能适合广义双曲分布的输出，并且您正在使用最大似然估计器来找到与您提供的数据最接近的广义双曲分布。但是如果你这样做，恐怕我不明白你怎么会有一个需要优化的非光滑表面。

一般来说，你选择的优化器需要根据你的问题来选择。在任何编程语言中都没有完美的“最佳优化器”，并且选择适合您的问题的优化算法可能会比任何轻微的低效率实现更重要。