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

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

1. Попробуйте 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')

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

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

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

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

Среди различных процедур оптимизации, доступных в SciPy, тот, который называется TNC (сопряженный градиент Ньютона с усечением), довольно устойчив к числовым значениям, связанным с проблемой.

симплексный метод Нелдера-Мида (предложенный Cristián Antuña в комментариях выше) хорошо известен как хороший выбор для оптимизации (возможно, плохой) ведет себя) функции без знания производных (см.

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

Одним из приятных аспектов Nelder-Mead является то, что функция выполняет серию шагов, которые настолько интуитивно понятны. Некоторые из этих точек, очевидно, могут выбросить вас из региона, но это легко изменить. Вот реализация 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 из любого псевдокода, и вы можете добавить имитацию отжига В любом случае.

Масштабирование

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

–