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

1)

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

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

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

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

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

2)

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

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

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

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

Принятый ответ @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
    ...

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

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

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

Вот, для сравнения, весь поиск, выполняемый алгоритмом по умолчанию:

[ 10000.]
[ 10000.00000001]
[ 10000.]

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

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

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

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

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

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

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

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

   Но это легко исправить: теперь, когда Skyfield имеет встроенный 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()

['A.D. 2016-Mar-09 03:37:33.8224 UT']

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

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