Сципи не найдет оптимального (простая функция косинуса)

Question

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

Код + демонстрация:

import numpy
import scipy

def solver(data):
    Z=numpy.zeros(len(data))
    a=0.003
    for i in range(len(data)):
        def minimizer(b):
            return numpy.abs(data[i]-numpy.cos(b))
        Z[i]=scipy.optimize.minimize(minimizer,a,bounds=[(0,numpy.pi)],method="L-BFGS-B").x[0]
    return Z

Y=numpy.zeros(100)
for i in range(100):
   Y[i]=numpy.cos(i/25)

solver(Y)

Результат не очень хороший, когда аргумент функции cos достигает значений выше 2, оценка «пропускает» значения и возвращает максимальное значение аргумента.

array([0.        , 0.04      , 0.08      , 0.12      , 0.16      ,
       0.2       , 0.24      , 0.28      , 0.32      , 0.36      ,
       0.4       , 0.44      , 0.48      , 0.52      , 0.56      ,
       0.6       , 0.64      , 0.67999999, 0.72      , 0.75999999,
       0.8       , 0.83999999, 0.88      , 0.92      , 0.95999999,
       1.        , 1.04      , 1.08      , 1.12      , 1.16      ,
       1.2       , 1.24      , 1.28      , 1.32      , 1.36      ,
       1.4       , 1.44      , 1.48      , 1.52      , 1.56      ,
       1.6       , 1.64      , 1.68      , 1.72      , 1.76      ,
       1.8       , 1.84      , 1.88      , 1.91999999, 1.95999999,
       2.        , 2.04      , 3.14159265, 3.14159265, 3.14159265,
       3.14159265, 3.14159265, 3.14159265, 3.14159265, 3.14159265,...

Что вызывает это явление? Существуют ли другие оптимизаторы / настройки, которые могут помочь с этой проблемой?

Answer 1

Помимо общего метода minimize, SciPy имеет minimize_scalar специально для одномерных задач, подобных здесь, и least_squares для минимизации определенного вида функций, которые измеряют разницу между двумя величинами (например, разницу между cos(b) и diff[i] здесь). Последний работает хорошо здесь, даже без тонкой настройки.

for i in range(len(data)):
    Z[i] = scipy.optimize.least_squares(lambda b: data[i] - numpy.cos(b), a, bounds=(0, numpy.pi)).x[0]

Функция, переданная в least_squares - это то, что мы хотели бы, чтобы она была близка к 0 без абсолютного значения. Я добавлю, что a = 0,003 кажется неоптимальным выбором для начальной точки, находящейся так близко к границе; тем не менее это работает.

Кроме того, как уже сообщал a_guest, метод поиска скалярного корня должен делать то же самое, вызывая меньше сюрпризов, учитывая, что у нас уже есть хороший интервал брекетинга [0, pi]. Бисекция надежна, но медленна; Метод Брента - это то, что я, вероятно, буду использовать.

for i in range(len(data)):
    Z[i] = scipy.optimize.brentq(lambda b: data[i] - numpy.cos(b), 0, numpy.pi)

Answer 2

Причина в том, что для функции (например) f = abs(cos(0.75*pi) - cos(z)) градиент f' оказывается равным нулю при z = pi, как видно из следующего графика:

Если вы проверите результат процедуры оптимизации, то увидите, что:

      fun: array([0.29289322])
 hess_inv: <1x1 LbfgsInvHessProduct with dtype=float64>
      jac: array([0.])
  message: b'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
     nfev: 16
      nit: 2
   status: 0
  success: True
        x: array([3.14159265])

Итак, процедура оптимизации достигла одного из критериев сходимости.Более подробную информацию о критерии можно найти в документации L-BFGS-B .Он говорит, что

gtol: float

Итерация остановится, когда max {| proj g_i |i = 1, ..., n} <= gtol, где pg_i - это i-й компонент проецируемого градиента. </p>

Таким образом, в конечном итоге он достигает точки z >= pi, которая затем проецируется обратно вz = pi из-за ограничения и в этот момент градиент функции равен нулю и, следовательно, он останавливается.Вы можете наблюдать это, регистрируя обратный вызов, который печатает текущий вектор параметров:

def new_callback():
    step = 1

    def callback(xk):
        nonlocal step
        print('Step #{}: xk = {}'.format(step, xk))
        step += 1

    return callback

scipy.optimize.minimize(..., callback=new_callback())

, который выводит:

Step #1: xk = [0.006]
Step #2: xk = [3.14159265]

Таким образом, на втором шаге он нажимает z >= pi, который проецируется обратнодо z = pi.

Вы можете обойти эту проблему, уменьшив границы, например, до bounds=[(0, 0.99*np.pi)].Это даст ожидаемый результат, однако метод не будет сходиться;вы увидите что-то вроде:

      fun: array([1.32930966e-09])
 hess_inv: <1x1 LbfgsInvHessProduct with dtype=float64>
      jac: array([0.44124484])
  message: b'ABNORMAL_TERMINATION_IN_LNSRCH'
     nfev: 160
      nit: 6
   status: 2
  success: False
        x: array([2.35619449])

Обратите внимание на сообщение ABNORMAL_TERMINATION_IN_LNSRCH.Это связано с природой abs(x) и тем фактом, что его производная имеет разрыв в x = 0 (вы можете прочитать больше об этом здесь ).

Альтернативный подход (поискroot)

Для всех строк выше мы пытались найти значение z, для которого cos(z) == cos(0.75*pi) (или abs(cos(z) - cos(0.75*pi)) < eps).Эта проблема на самом деле находит корень функции f = cos(z) - cos(0.75*pi), где мы можем использовать тот факт, что cos является непрерывной функцией.Нам нужно установить границы a, b так, чтобы f(a)*f(b) < 0 (т.е. они имели противоположный знак).Например, используя bisect метод :

res = scipy.optimize.bisect(f, 0, np.pi)