R:快速多变量优化包？

Question

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

我想知道,有没有任何好的,快速的优化套件供人们使用,我可能会丢失？理想情况下,这些将是围绕优秀的C++/Fortran库的薄包装形式,因此最小的纯R代码.(虽然这不应该是相关的,但是在尝试将4参数分布拟合到一组值时,通过最小化某个拟合优度度量,我的优化问题出现了).

在过去,我发现R的优化库非常慢,最后编写了一个瘦R包装器来调用商业优化库的C++ API.那么最好的图书馆必然是商业图书馆吗？

UPDATE.这是我正在查看的代码的简化示例:

###########
## 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- 这使用最大似然期望最大化(EM)方法.这很快就邪恶!

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%的时间可归因于目标函数!所以Dirk和Mike在下面的评论中是正确的.我应该更仔细地估计在我的目标函数中花费的时间,打印点不是一个好主意!另外我怀疑MLE(EM)方法非常快,因为它使用分析(闭合形式)作为对数似然函数.

Answer 1

当问题存在时,最大似然估计器总是比任何语言的全局优化器都快.

无论算法如何,全局优化器通常将一些随机跳转与局部最小化例程组合在一起.不同的算法可以在群体(遗传算法),退火,迁移等方面对此进行讨论,但它们在概念上都是相似的.

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

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

通常,您选择的优化程序需要根据您的问题进行选择.在任何编程语言中都没有完美的"最优优化器",并且选择适合您的问题的优化算法可能会比实现的任何轻微低效率产生更大的差异.