R: Пакеты многовариантной оптимизации FAST?

Я ищу локальный минимум скалярной функции от 4 переменных, и у меня есть ограничения диапазона для переменных ("ограничения блока"). Для производной функции нет замкнутой формы, поэтому о методах, требующих аналитической производной функции, не может быть и речи. Я пробовал несколько вариантов и параметров управления с функцией optim , но все они кажутся очень медленными. В частности, кажется, что они проводят много времени между вызовами моей (определенной R) целевой функции, поэтому я знаю, что узким местом является не моя целевая функция, а «размышления» между вызовами моей целевой функции. Я посмотрел на представление задач CRAN для оптимизации и попробовал несколько из этих вариантов ( DEOptim из RcppDE и т. Д.), Но ни один из них не оказался подходящим. Я бы хотел попробовать пакет nloptr (оболочка R для библиотеки NLOPT), но он, похоже, недоступен для 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 - в нем используется метод максимизации ожидания максимального правдоподобия (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% времени приходится на целевую функцию! Итак, Дирк и Майк были правы в своих комментариях ниже. Мне следовало более тщательно оценить время, потраченное на мою целевую функцию, а печатать точки - не лучшая идея! Также я подозреваю, что метод MLE (EM) очень быстр, потому что он использует аналитическую (замкнутую форму) для функции логарифмического правдоподобия.