Почему scipy.optimize.minimize (default) сообщает об успехе, не двигаясь с Skyfield?

scipy.optimize.minimize с использованием метода по умолчанию возвращает исходное значение в качестве результата без каких-либо сообщений об ошибках или предупреждениях. Хотя метод Nelder-Mead, предложенный в этом ответе, решает проблему, я хотел бы понять:

Почему метод по умолчанию возвращает неверный ответ, не предупреждая отправную точку в качестве ответа – и есть ли способ защитить от «неправильного ответа без предупреждения», чтобы избежать такого поведения в этом случае?

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

РЕЗУЛЬТАТЫ:

результат теста: [ 2.14159739 ] 'correct': 2.14159265359 initial: 0.0

результат по умолчанию: [ 10000. ] 'correct': 13054 initial: 10000

Результат Nelder -Mead: [ 13053.81011963 ] 'correct': 13054 initial: 10000

 FULL OUTPUT using DEFAULT METHOD: status: 0 success: True njev: 1 nfev: 3 hess_inv: array([[1]]) fun: 1694.98753895812 x: array([ 10000.]) message: 'Optimization terminated successfully.' jac: array([ 0.]) nit: 0 FULL OUTPUT using Nelder-Mead METHOD: status: 0 nfev: 63 success: True fun: 3.2179306044608054 x: array([ 13053.81011963]) message: 'Optimization terminated successfully.' nit: 28 

Вот полный скрипт:

 def g(x, a, b): return np.cos(a*x + b) def separation(seconds, lat, lon): lat, lon, seconds = float(lat), float(lon), float(seconds) # necessary it seems place = earth.topos(lat, lon) jd = JulianDate(utc=(2016, 3, 9, 0, 0, seconds)) mpos = place.at(jd).observe(moon).apparent().position.km spos = place.at(jd).observe(sun).apparent().position.km mlen = np.sqrt((mpos**2).sum()) slen = np.sqrt((spos**2).sum()) sepa = ((3600.*180./np.pi) * np.arccos(np.dot(mpos, spos)/(mlen*slen))) return sepa from skyfield.api import load, now, JulianDate import numpy as np from scipy.optimize import minimize data = load('de421.bsp') sun = data['sun'] earth = data['earth'] moon = data['moon'] x_init = 0.0 out_g = minimize(g, x_init, args=(1, 1)) print "test result: ", out_g.x, "'correct': ", np.pi-1, "initial: ", x_init # gives right answer sec_init = 10000 out_s_def = minimize(separation, sec_init, args=(32.5, 215.1)) print "default result: ", out_s_def.x, "'correct': ", 13054, "initial: ", sec_init sec_init = 10000 out_s_NM = minimize(separation, sec_init, args=(32.5, 215.1), method = "Nelder-Mead") print "Nelder-Mead result: ", out_s_NM.x, "'correct': ", 13054, "initial: ", sec_init print "" print "FULL OUTPUT using DEFAULT METHOD:" print out_s_def print "" print "FULL OUTPUT using Nelder-Mead METHOD:" print out_s_NM 

    3 Solutions collect form web for “Почему scipy.optimize.minimize (default) сообщает об успехе, не двигаясь с Skyfield?”

    1)

    Ваша функция кусочно-постоянной (имеет мелкомасштабную «лестницу»). Это не везде дифференцируемо.

    Градиент функции в исходном предположении равен нулю.

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

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

    Обратите также внимание на то, что Nelder-Mead не защищен от этого – просто случается, что его первоначальный симплекс больше размера лестницы, поэтому он исследует функцию вокруг большего места. Если начальный симплекс будет меньше размера лестницы, оптимизатор будет вести себя аналогично BFGS.

    2)

    Общий ответ: локальные оптимизаторы возвращают локальные оптимумы. Совместимы ли они с истинным оптимумом, зависит от свойств функции.

    В общем, чтобы убедиться, что вы застряли в локальном оптимуме, попробуйте разные исходные предположения.

    Кроме того, использование оптимизатора на основе производных на недифференцируемой функции не является хорошей идеей. Если функция дифференцируема в «большом» масштабе, вы можете отрегулировать размер шага численного дифференцирования.

    Поскольку нет дешевого / общего способа проверить численно, если функция везде дифференцируема, такая проверка не выполняется — вместо этого это предположение в методе оптимизации, которое должно быть обеспечено тем, кто вводит целевую функцию и выбирает метод оптимизации ,

    Принятый ответ by @pv. объясняет, что Skyfield имеет ответ «лестницы», что означает, что некоторые значения, которые он возвращает, локально плоские, за исключением дискретных прыжков.

    Я сделал небольшой эксперимент на первом этапе – преобразование времени в объекты JulianDate, действительно, оно выглядит примерно 40 микросекунд на каждый шаг или около 5E-10 дней. Это разумно, учитывая, что базы данных JPL охватывают тысячи лет . Хотя это, вероятно, отлично подходит практически для любого приложения с общими астрономическими масштабами, оно на самом деле не является гладким. Как указывает ответ – локальная плоскость вызовет «успех» в некоторых (возможно, многих) минимизаторах. Это ожидается и разумно и никоим образом не является отказом метода.

    дискретное время в небоскребе

     from skyfield.api import load, now, JulianDate import numpy as np import matplotlib.pyplot as plt t = 10000 + np.logspace(-10, 2, 25) # logarithmic spacing jd = JulianDate(utc=(2016, 3, 9, 0, 0, t)) dt = t[1:] - t[:-1] djd = jd.tt[1:] - jd.tt[:-1] t = 10000 + np.linspace(0, 0.001, 1001) # linear spacing jd = JulianDate(utc=(2016, 3, 9, 0, 0, t)) plt.figure() plt.subplot(1,2,1) plt.plot(dt, djd) plt.xscale('log') plt.yscale('log') plt.subplot(1,2,2) plt.plot(t, jd.tt-jd.tt[0]) plt.show() 

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

     def separation(seconds, lat, lon): print seconds ... 

    Добавление этой строки позволит вам увидеть, что метод Nelder-Mead выполняет тщательный поиск диапазона секунд, шаг вперед с шагом в 500 секунд, прежде чем он начнет воспроизводиться:

     [ 10000.] [ 10500.] [ 11000.] [ 11500.] [ 12500.] ... 

    Конечно, он не знает, что это 500-секундные приращения, потому что для такого решателя проблема не имеет единиц. Эти корректировки могут составлять 500 метров, или 500 ангстрем, или 500 лет. Но он натыкается слепо вперед и, в случае с Nelder-Mead, видит достаточно того, как выход изменяется с входом, чтобы отточить ответ на понравившийся вам ответ.

    Здесь для сравнения используется весь поиск по алгоритму по умолчанию:

     [ 10000.] [ 10000.00000001] [ 10000.] 

    Вот и все. Он пытается немного отступить на 1е-8 секунд, не может видеть ничего другого в ответе, который он получает, и отказывается – так как некоторые из других ответов здесь правильно утверждаются.

    Иногда вы можете исправить ситуацию, как это, сказав алгоритму: (a) предпринять больший шаг для начала и (b) прекратить тестирование, как только размер шага, который он создает, становится небольшим – скажем, когда он опускается до миллисекунды. Вы можете попробовать что-то вроде:

     out_s_def = minimize(separation, sec_init, args=(32.5, 215.1), tol=1e-3, options={'eps': 500}) 

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

    Понимаете, эти процедуры минимизации часто пишутся с довольно явным пониманием того, насколько точный 64-битный поплавок может быть нажат до того, как не будет больше точности, и все они предназначены для остановки до этого момента. Но вы скрываете точность: вы говорите рутине «дайте мне несколько секунд», что заставляет их думать, что они могут играть с очень маленькими младшими разрядами значения секунд, когда на самом деле секунды объединяются с а не только часы и дни, но с годами, и при этом теряется всякая крошечная точность вниз в нижней части секунд, хотя минимизатор не знает!

    Итак, давайте разоблачим реальное время с плавающей запятой алгоритма. В этом процессе я сделаю три вещи:

    1. Давайте избежим необходимости маневра float() вы делаете. Наш оператор print показывает проблему: несмотря на то, что вы предоставили скалярный float, минимизатор превращает его в массив NumPy:

       (array([ 10000.]), 32.5, 215.1) 

      Но это легко исправить: теперь, когда Skyfield имеет встроенный separ_from separation_from() который отлично справляется с массивами, мы будем использовать его:

       sepa = mpos.separation_from(spos) return sepa.degrees 
    2. Я переключусь на новый синтаксис для создания дат, который Skyfield принял, когда он направляется к 1.0.

    Это дает нам что-то вроде (но заметьте, что это было бы быстрее, если бы вы только построили topos один раз и передали его, вместо того, чтобы перестраивать его и заставлять его выполнять свою математику каждый раз):

     ts = load.timescale() ... def separation(tt_jd, lat, lon): place = earth.topos(lat, lon) t = ts.tt(jd=tt_jd) mpos = place.at(t).observe(moon).apparent() spos = place.at(t).observe(sun).apparent() return mpos.separation_from(spos).degrees ... sec_init = 10000.0 jd_init = ts.utc(2016, 3, 9, 0, 0, sec_init).tt out_s_def = minimize(separation, jd_init, args=(32.5, 215.1)) 

    Результат – успешное премирование, дающее – я думаю, если бы вы могли проверить меня здесь? – ответ, который вы ищете:

     print ts.tt(jd=out_s_def.x).utc_jpl() ['AD 2016-Mar-09 03:37:33.8224 UT'] 

    Надеюсь, скоро скоро будет добавлено несколько заранее разработанных процедур минификсации в Skyfield – на самом деле, большая причина для написания его, чтобы заменить PyEphem, была нужна возможность развязать мощные оптимизаторы SciPy и быть в состоянии отказаться от довольно анемичных, которые PyEphem реализуется в C. Главное, что им нужно будет быть осторожным, – это то, что произошло здесь: что оптимизатору должны быть даны цифры с плавающей запятой для покачивания, которые значительны до конца.

    Возможно, мне стоит изучить возможность создания объектов Time из двух объектов с плавающей запятой, чтобы можно было представить еще несколько секунд. Я думаю, что AstroPy сделала это, и это традиционно в программировании астрономии.

    Python - лучший язык программирования в мире.