одиночный скаляр scipy оптимизировать минимизировать

У меня проблема с использованием функции minimize() scipy, и я недостаточно понимаю в оптимизации, чтобы понять, что здесь не так.

У меня есть функция, которая вызывает scipy.optimize.minimize(). Он отлично работает и предоставляет мне именно те выходные данные, которые мне нужны, когда x0 представляет собой массив размером > 1, но когда x0 равен ровно 1, он терпит неудачу. В документации говорится, что x0 должен быть np.ndarray размера (n,), но не указывает, что он должен быть> 1, поэтому я предположил, что все будет в порядке. Уменьшенная версия моего кода, вызывающего функцию с оптимальным значением:

import numpy as np
from scipy.optimize import minimize

def to_freq(*arrays):
    # Better version of `convert_to_freq()`
    out = []
    for a in arrays:
        converted = np.array([(x + i / len(a)) / (max(a)+1) for i, x in enumerate(a, start=1)])
        out.append(converted)
    return out

def likelihood(x, x_freq, expected, x_max):
    # Better version, supports vectorisation
    a = 2 * x * np.log(x_freq / expected) 
    b = 2 * (x_max - x) * np.log((1 - x_freq) / (1 - expected))
    return a + b

def objective(x0, labels, a, b):
    R = x0[labels=='R'].item()

    a_c, b_c = np.cumsum(a), np.cumsum(b)
    a_f, b_f = to_freq(a_c, b_c)

    # Get the expected values for signals and noises
    exp_a = ((1 - R) * b_f + R)[:-1]
    exp_b = b_f[:-1]

    # Compute the gsquared using the dual process model parameters
    #   Still getting runtime warnings about division. Function only works with numpy, so can't use math.
    a_lrat = likelihood(x=a_c[:-1], x_freq=a_f[:-1], expected=exp_a, x_max=a_c.max())
    b_lrat = likelihood(x=b_c[:-1], x_freq=b_f[:-1], expected=exp_b, x_max=b_c.max())

    return sum(a_lrat + b_lrat)

# Observations
a = [508,224,172,135,119,63]
b = [102,161,288,472,492,308]
x0 = np.array([0.520274590415736]) # Optimal value for variable
labels = np.array(['R'])

# Gives correct iotimized value of 163.27525607890783
objective(x0, labels, a, b)

А теперь случайная инициализация x0 для случаев, когда оптимальное значение неизвестно:

x0 = np.random.uniform(-.5,0.5, len(labels)) # random initialization

# Without method='nelder-mead' occasionally gives correct value of fun, but frequently fails
opt = minimize(fun=objective, x0=x0, args=(labels, a, b), tol=1e-4)
print(opt)

Неудачный результат оптимизации:

      fun: nan
 hess_inv: array([[1]])
      jac: array([nan])
  message: 'Desired error not necessarily achieved due to precision loss.'
     nfev: 336
      nit: 1
     njev: 112
   status: 2
  success: False
        x: array([1034.74])

Но если я продолжаю запускать это и случайным образом устанавливаю начальное значение, оно иногда выдает хороший результат:

      fun: 163.27525607888913
 hess_inv: array([[4.14149525e-05]])
      jac: array([-1.90734863e-05])
  message: 'Optimization terminated successfully.'
     nfev: 27
      nit: 7
     njev: 9
   status: 0
  success: True
        x: array([0.52027462])

Если я укажу method='nelder-mead' (решение возможно, не связанная с этим проблема) в вызове minimize() в моей более крупной функции, он также на самом деле дает мне ожидаемый результат:

 final_simplex: (array([[0.52026029],
       [0.52031204]]), array([163.27525856, 163.27527298]))
           fun: 163.2752585612531
       message: 'Optimization terminated successfully.'
          nfev: 32
           nit: 16
        status: 0
       success: True
             x: array([0.52026029])

Я действительно не понимаю, как лучше всего это реализовать, поскольку я очень неопытен в оптимизации.

[Сноска]: Алгоритм минимизации иногда пытается использовать значения, несовместимые с моей функцией (например, ‹ 0 или > 1), и вызов np.log() заканчивается предупреждением, но обычно я просто подавляю это, так как кажется, что он работает независимо ...


person fffrost    schedule 05.05.2020    source источник
comment
То, что вы пишете в сноске, вообще не очень хорошая идея. Когда вы используете процедуры оптимизации, вы обычно гарантируете, что оптимизируемые переменные могут принимать любое значение от -inf до +inf. Существуют алгоритмы, которые позволяют вам определять границы, но это должно быть указано явно и работает только с подмножеством алгоритмов оптимизации. Часто легко перепараметризировать, чтобы выполнить условие. Например. если функция f(x) действительна только для значений › 0, вы можете использовать y := exp(x), и f(y) будет действительна для всех y на вещественной оси.   -  person cel    schedule 05.05.2020
comment
Я не совсем уверен, как я мог бы это реализовать - ситуация возникает во время вызова likelihood(), где используется функция np.log()   -  person fffrost    schedule 05.05.2020
comment
Метод Нелдера-Мида не требует вычисления градиентов функций, поэтому он работает в вашем случае. Я нахожу Nelder-Mead хорошим во многих случаях на реальных данных. Если вы можете вычислить градиент своей функции, то алгоритмы, требующие градиентного спуска, будут работать лучше. Здесь есть хорошее объяснение этих алгоритмов: scipy-lectures.org/advanced/mathematical_optimization   -  person Paddy Harrison    schedule 05.05.2020