Минимизация многопараметрической функции с помощью scipy. Производные не известны

Question

У меня есть функция, которая на самом деле является вызовом другой программы (некоторый код Fortran). Когда я вызываю эту функцию (run_moog), я могу проанализировать 4 переменные, и она возвращает 6 значений. Эти значения должны быть близки к 0 (для минимизации). Однако я объединил их вот так: np.sum(results**2). Теперь у меня есть скалярная функция. Я хотел бы свести к минимуму эту функцию, т. Е. Получить np.sum(results**2) как можно ближе к нулю.
Примечание: Когда эта функция (run_moog) принимает 4 входных параметра, она создает входной файл для кода Fortran, который зависит от этих параметров.

Я попытался несколько способов оптимизировать это от the scipy docs. Но никто не работает так, как ожидалось. Минимизация должна иметь ограничения на 4 переменные. Вот попытка:

from scipy.optimize import minimize # Tried others as well from the docs 
x0 = 4435, 3.54, 0.13, 2.4 
bounds = [(4000, 6000), (3.00, 4.50), (-0.1, 0.1), (0.0, None)] 
a = minimize(fun_mmog, x0, bounds=bounds, method='L-BFGS-B') # I've tried several different methods here 
print a 

Это тогда дает мне

status: 0 
success: True 
    nfev: 5 
    fun: 2.3194639999999964 
     x: array([ 4.43500000e+03, 3.54000000e+00, 1.00000000e-01, 
     2.40000000e+00]) 
message: 'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL' 
    jac: array([ 0., 0., -54090399.99999981, 0.]) 
    nit: 0 

Третий параметр изменяется незначительно, в то время как другие точно так же. Также было 5 вызовов функций (nfev), но нет итераций (nit). output from scipy is shown here.

Answer 1

Пара возможностей:

1. Try COBYLA. Он должен быть без производных и поддерживать ограничения неравенства.
2. Вы не можете использовать разные эпсилоны через обычный интерфейс; поэтому попробуйте масштабировать свою первую переменную на 1e4. (Разделите это происходит в умножьте приходит обратно.)
3. Пропустить нормальный автоматический конструктор якобиан, и сделать свой собственный:

Допустим, вы пытаетесь использовать SLSQP, и вы не обеспечивают якобиан функция. Это делает его для вас. Код для него находится в approx_jacobian в slsqp.py. Вот сжатая версия:

def approx_jacobian(x,func,epsilon,*args): 
    x0 = asfarray(x) 
    f0 = atleast_1d(func(*((x0,)+args))) 
    jac = zeros([len(x0),len(f0)]) 
    dx = zeros(len(x0)) 
    for i in range(len(x0)): 
     dx[i] = epsilon 
     jac[i] = (func(*((x0+dx,)+args)) - f0)/epsilon 
     dx[i] = 0.0 

    return jac.transpose() 

Можно попробовать заменить эту петлю:

for (i, e) in zip(range(len(x0)), epsilon): 
     dx[i] = e 
     jac[i] = (func(*((x0+dx,)+args)) - f0)/e 
     dx[i] = 0.0 

Вы не можете обеспечить это как якобиан к minimize, но фиксируя его на который прост:

def construct_jacobian(func,epsilon): 
    def jac(x, *args): 
     x0 = asfarray(x) 
     f0 = atleast_1d(func(*((x0,)+args))) 
     jac = zeros([len(x0),len(f0)]) 
     dx = zeros(len(x0)) 
     for i in range(len(x0)): 
      dx[i] = epsilon 
      jac[i] = (func(*((x0+dx,)+args)) - f0)/epsilon 
      dx[i] = 0.0 

     return jac.transpose() 
    return jac 

Вы можете назвать minimize как:

minimize(fun_mmog, x0, 
     jac=construct_jacobian(fun_mmog, [1e0, 1e-4, 1e-4, 1e-4]), 
     bounds=bounds, method='SLSQP') 

Answer 2

Похоже, что ваша целевая функция не имеет хороших производных. Строка на выходе jac: array([ 0., 0., -54090399.99999981, 0.]) означает, что изменение только третьего значения переменной является значительным. А так как производная w.r.t. эта переменная практически бесконечна, вероятно, что-то не так в функции. Именно поэтому третье значение переменной заканчивается в максимуме.

Я предлагаю вам взглянуть на производные, по крайней мере, в нескольких точках вашего пространства параметров. Вычислите их с использованием конечных разностей и размера шага по умолчанию для SciPy's fmin_l_bfgs_b, 1e-8. Here - пример того, как вы могли бы вычислить производные.

Попробуйте также указать целевую функцию. Например, сохраните два параметра постоянными и позвольте двум другим варьироваться. Если функция имеет несколько локальных оптимумов, вам не следует использовать методы на основе градиента, такие как BFGS.

Answer 3

Насколько сложно получить аналитическое выражение для градиента? Если у вас есть, то вы можете затем аппроксимировать произведение гессиана вектором с использованием конечной разности. Затем вы можете использовать другие подпрограммы оптимизации.

Среди различных подпрограмм оптимизации, доступных в SciPy, тот, который называется TNC (Newton Conjugate Gradient with Truncation), достаточно устойчив к числовым значениям, связанным с проблемой.

Answer 4

Nelder-Mead Simplex Method (предложенный Cristián Antuña в вышеуказанных комментариях) хорошо известен, чтобы быть хорошим выбором для оптимизации (posibly плохо вели себя) функции без знания производных (см Numerical Recipies In C, Chapter 10).

Есть два конкретных аспекта вашего вопроса. Во-первых, это ограничения на входы, а вторая - проблема масштабирования. Следующее предлагает решения для этих точек, но вам может потребоваться вручную перебирать между ними несколько раз, пока все не сработает.

Входные Ограничения

Если предположить, что ваши входные ограничения образуют convex region (как ваши примеры выше показывают, но я хотел бы обобщить его немного), то вы можете написать функцию

is_in_bounds(p): 
    # Return if p is in the bounds 

Используя эту функцию, предположим, что алгоритм хочет перейти от пункта from_ к точке to, где известно, что from_ находится в регионе. Тогда следующая функция будет эффективно найти дальнюю точку на линии между двумя точками, на которых он может продолжаться:

from numpy.linalg import norm 

def progress_within_bounds(from_, to, eps): 
    """ 
    from_ -- source (in region) 
    to -- target point 
    eps -- Eucliedan precision along the line 
    """ 

    if norm(from_, to) < eps: 
     return from_ 
    mid = (from_ + to)/2 
    if is_in_bounds(mid): 
     return progress_within_bounds(mid, to, eps) 
    return progress_within_bounds(from_, mid, eps) 

(Обратите внимание, что эта функция может быть оптимизирована для некоторых регионов, но вряд ли стоит беспокоиться, так как он даже не вызывает вашу первоначальную функцию объекта, которая является дорогой.)

Одним из приятных аспектов Nelder-Mead является то, что функция выполняет ряд шагов, которые являются настолько интуитивными. Некоторые из этих пунктов, очевидно, могут выкинуть вас из региона, но это легко изменить. Вот implementation of Nelder Mead с изменениями, сделанные отмечены между парами строк вида ##################################################################:

import copy 

''' 
    Pure Python/Numpy implementation of the Nelder-Mead algorithm. 
    Reference: https://en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method 
''' 


def nelder_mead(f, x_start, 
     step=0.1, no_improve_thr=10e-6, no_improv_break=10, max_iter=0, 
     alpha = 1., gamma = 2., rho = -0.5, sigma = 0.5): 
    ''' 
     @param f (function): function to optimize, must return a scalar score 
      and operate over a numpy array of the same dimensions as x_start 
     @param x_start (numpy array): initial position 
     @param step (float): look-around radius in initial step 
     @no_improv_thr, no_improv_break (float, int): break after no_improv_break iterations with 
      an improvement lower than no_improv_thr 
     @max_iter (int): always break after this number of iterations. 
      Set it to 0 to loop indefinitely. 
     @alpha, gamma, rho, sigma (floats): parameters of the algorithm 
      (see Wikipedia page for reference) 
    ''' 

    # init 
    dim = len(x_start) 
    prev_best = f(x_start) 
    no_improv = 0 
    res = [[x_start, prev_best]] 

    for i in range(dim): 
     x = copy.copy(x_start) 
     x[i] = x[i] + step 
     score = f(x) 
     res.append([x, score]) 

    # simplex iter 
    iters = 0 
    while 1: 
     # order 
     res.sort(key = lambda x: x[1]) 
     best = res[0][1] 

     # break after max_iter 
     if max_iter and iters >= max_iter: 
      return res[0] 
     iters += 1 

     # break after no_improv_break iterations with no improvement 
     print '...best so far:', best 

     if best < prev_best - no_improve_thr: 
      no_improv = 0 
      prev_best = best 
     else: 
      no_improv += 1 

     if no_improv >= no_improv_break: 
      return res[0] 

     # centroid 
     x0 = [0.] * dim 
     for tup in res[:-1]: 
      for i, c in enumerate(tup[0]): 
       x0[i] += c/(len(res)-1) 

     # reflection 
     xr = x0 + alpha*(x0 - res[-1][0]) 
     ################################################################## 
     ################################################################## 
     xr = progress_within_bounds(x0, x0 + alpha*(x0 - res[-1][0]), prog_eps) 
     ################################################################## 
     ################################################################## 
     rscore = f(xr) 
     if res[0][1] <= rscore < res[-2][1]: 
      del res[-1] 
      res.append([xr, rscore]) 
      continue 

     # expansion 
     if rscore < res[0][1]: 
      xe = x0 + gamma*(x0 - res[-1][0]) 
      ################################################################## 
      ################################################################## 
      xe = progress_within_bounds(x0, x0 + gamma*(x0 - res[-1][0]), prog_eps) 
      ################################################################## 
      ################################################################## 
      escore = f(xe) 
      if escore < rscore: 
       del res[-1] 
       res.append([xe, escore]) 
       continue 
      else: 
       del res[-1] 
       res.append([xr, rscore]) 
       continue 

     # contraction 
     xc = x0 + rho*(x0 - res[-1][0]) 
     ################################################################## 
     ################################################################## 
     xc = progress_within_bounds(x0, x0 + rho*(x0 - res[-1][0]), prog_eps) 
     ################################################################## 
     ################################################################## 
     cscore = f(xc) 
     if cscore < res[-1][1]: 
      del res[-1] 
      res.append([xc, cscore]) 
      continue 

     # reduction 
     x1 = res[0][0] 
     nres = [] 
     for tup in res: 
      redx = x1 + sigma*(tup[0] - x1) 
      score = f(redx) 
      nres.append([redx, score]) 
     res = nres 

Примечание Эта реализация GPL, что либо нормально для вас или нет. Однако очень легко модифицировать NM из любого псевдокода, но в любом случае вы можете захотеть в simulated annealing.

Scaling

Это является хитрой проблемой, но jasaarim сделал интересный момент относительно этого. Как только модифицированный алгоритм NM найдет точку, вы можете запустить matplotlib.contour, фиксируя несколько измерений, чтобы увидеть, как работает функция. На этом этапе вам может потребоваться изменить масштаб одного или нескольких параметров и перезапустить модифицированный NM.