为什么scipy.optimize.minimize（默认）不随Skyfield移动就报成功？
Why does scipy.optimize.minimize (default) report success without moving with Skyfield?
scipy.optimize.minimize 使用默认方法返回初始值作为结果，没有任何错误或警告消息。虽然使用 建议的 Nelder-Mead 方法解决了问题，但我想了解：
为什么默认方法returns 没有警告的错误答案作为答案的起点 - 有没有办法我可以防止"wrong answer without warning"在这种情况下避免这种行为？
请注意，函数 separation 使用 python 程序包 Skyfield 来生成要最小化的值，但不能保证平滑，这可能是 Simplex 在这里更好的原因。
结果：
测试结果：[2.14159739] 'correct'：2.14159265359初始：0.0
默认结果：[10000.] 'correct'：13054 初始值：10000
Nelder-Mead 结果：[13053.81011963] 'correct'：13054 初始值：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
1)
您的函数是分段常数（具有小规模 "staircase" 模式）。
它不是处处可微的。
初始猜测时函数的梯度为零。
默认的 BFGS 优化器看到零梯度并根据其标准（基于在这种情况下不正确的输入函数假设，例如可微性）确定它是局部最小值。
基本上，完全平坦的区域会轰炸优化器。优化器在初始点周围的小而平坦的区域中探测函数，其中一切看起来都像函数只是一个常数，所以它认为你给了它一个常数函数。因为你的函数不是处处可微的，所以可能几乎所有的初始点都在这样的平坦区域内，这样就可以在选择初始点时运气不佳。
另请注意，Nelder-Mead 并非 对此免疫——恰好它的初始单纯形大于楼梯的大小，因此它会探测周围的函数一个更大的地方。如果初始单纯形小于阶梯大小，优化器的行为将与 BFGS 类似。
2)
一般答案：局部优化器return局部最优。这些是否与真正的最优值一致取决于函数的属性。
一般来说，要查看您是否陷入局部最优，请尝试不同的初始猜测。
此外，在不可微函数上使用基于导数的优化器也不是一个好主意。如果函数在 "large" 尺度上可微分，您可以调整数值微分的步长。
因为没有 cheap/general 方法来检查函数是否处处可微，所以没有进行此类检查 --- 相反，它是优化方法中的一个假设，必须由输入该函数的任何人确保objective 函数并选择优化方法。
@pv 接受的答案。解释说 Skyfield 有一个 "staircase" 响应，这意味着它 returns 的一些值在局部是平坦的，除了离散跳跃。
我在第一步做了一个小实验 - 将时间转换为 JulianDate 对象，确实看起来每次增量大约 40 微秒，或者大约 5E-10 天。这是合理的，考虑到 JPL 数据库 跨越数千年 。虽然这对于几乎所有一般天文规模的应用程序来说可能都很好，但实际上并不顺利。正如答案所指出的那样——局部平坦度将在一些（可能很多）最小化器中触发 "success"。这是预期和合理的，绝不是该方法的失败。
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.]
就是这样。它尝试稍微离开 1e-8 秒，看不到它得到的答案有任何不同，然后放弃——正如这里正确断言的其他几个答案。
有时您可以通过告诉算法 (a) 开始时采取更大的步长和 (b) 一旦步长变小就停止测试——比如，当步长下降时，您可以解决这种情况到一毫秒。您可以尝试类似的方法：
out_s_def = minimize(separation, sec_init, args=(32.5, 215.1),
tol=1e-3, options={'eps': 500})
在这种情况下，即使有这种帮助，默认的最小化技术似乎仍然太脆弱，无法建设性地找到最小值，所以我们可以做点别的：我们可以告诉最小化函数它真正需要多少位一起玩
你看，这些最小化例程通常是在相当明确的情况下编写的，即在没有更多精度可用之前可以将 64 位浮点数推到多精确，并且它们都设计为在该点之前停止。但是你隐藏了精度：你告诉例程“给我几秒钟”，这让他们认为他们可以 fiddle 使用秒值的非常小的低端数字，而实际上秒不仅与小时和天相结合，而且与年相结合，在这个过程中，秒数底部的任何微小精度都会丢失——尽管最小化器不知道！
所以让我们把真正的浮点时间暴露给算法。在这个过程中，我会做三件事：
让我们避免需要您正在做的 float() 操作。我们的 print 语句显示了问题：即使您提供了一个标量浮点数，最小化器仍然会将其转换为 NumPy 数组：
(array([ 10000.]), 32.5, 215.1)
但这很容易解决：既然 Skyfield 有一个内置的 separation_from() 可以很好地处理数组，我们将使用它：
sepa = mpos.separation_from(spos)
return sepa.degrees
我将切换到创建日期的新语法，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 做到了这一点，在天文编程中是传统的。
scipy.optimize.minimize 使用默认方法返回初始值作为结果，没有任何错误或警告消息。虽然使用
为什么默认方法returns 没有警告的错误答案作为答案的起点 - 有没有办法我可以防止"wrong answer without warning"在这种情况下避免这种行为？
请注意，函数 separation 使用 python 程序包 Skyfield 来生成要最小化的值，但不能保证平滑，这可能是 Simplex 在这里更好的原因。
结果：
测试结果：[2.14159739] 'correct'：2.14159265359初始：0.0
默认结果：[10000.] 'correct'：13054 初始值：10000
Nelder-Mead 结果：[13053.81011963] 'correct'：13054 初始值：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
1)
您的函数是分段常数（具有小规模 "staircase" 模式）。 它不是处处可微的。
初始猜测时函数的梯度为零。
默认的 BFGS 优化器看到零梯度并根据其标准（基于在这种情况下不正确的输入函数假设，例如可微性）确定它是局部最小值。
基本上，完全平坦的区域会轰炸优化器。优化器在初始点周围的小而平坦的区域中探测函数，其中一切看起来都像函数只是一个常数，所以它认为你给了它一个常数函数。因为你的函数不是处处可微的，所以可能几乎所有的初始点都在这样的平坦区域内，这样就可以在选择初始点时运气不佳。
另请注意，Nelder-Mead 并非 对此免疫——恰好它的初始单纯形大于楼梯的大小，因此它会探测周围的函数一个更大的地方。如果初始单纯形小于阶梯大小，优化器的行为将与 BFGS 类似。
2)
一般答案：局部优化器return局部最优。这些是否与真正的最优值一致取决于函数的属性。
一般来说，要查看您是否陷入局部最优，请尝试不同的初始猜测。
此外，在不可微函数上使用基于导数的优化器也不是一个好主意。如果函数在 "large" 尺度上可微分，您可以调整数值微分的步长。
因为没有 cheap/general 方法来检查函数是否处处可微，所以没有进行此类检查 --- 相反，它是优化方法中的一个假设，必须由输入该函数的任何人确保objective 函数并选择优化方法。
@pv 接受的答案。解释说 Skyfield 有一个 "staircase" 响应，这意味着它 returns 的一些值在局部是平坦的，除了离散跳跃。
我在第一步做了一个小实验 - 将时间转换为 JulianDate 对象，确实看起来每次增量大约 40 微秒，或者大约 5E-10 天。这是合理的，考虑到 JPL 数据库 跨越数千年 。虽然这对于几乎所有一般天文规模的应用程序来说可能都很好，但实际上并不顺利。正如答案所指出的那样——局部平坦度将在一些（可能很多）最小化器中触发 "success"。这是预期和合理的，绝不是该方法的失败。
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.]
就是这样。它尝试稍微离开 1e-8 秒，看不到它得到的答案有任何不同，然后放弃——正如这里正确断言的其他几个答案。
有时您可以通过告诉算法 (a) 开始时采取更大的步长和 (b) 一旦步长变小就停止测试——比如，当步长下降时，您可以解决这种情况到一毫秒。您可以尝试类似的方法：
out_s_def = minimize(separation, sec_init, args=(32.5, 215.1),
tol=1e-3, options={'eps': 500})
在这种情况下，即使有这种帮助，默认的最小化技术似乎仍然太脆弱，无法建设性地找到最小值，所以我们可以做点别的：我们可以告诉最小化函数它真正需要多少位一起玩
你看，这些最小化例程通常是在相当明确的情况下编写的，即在没有更多精度可用之前可以将 64 位浮点数推到多精确，并且它们都设计为在该点之前停止。但是你隐藏了精度：你告诉例程“给我几秒钟”，这让他们认为他们可以 fiddle 使用秒值的非常小的低端数字，而实际上秒不仅与小时和天相结合，而且与年相结合，在这个过程中，秒数底部的任何微小精度都会丢失——尽管最小化器不知道！
所以让我们把真正的浮点时间暴露给算法。在这个过程中，我会做三件事：
让我们避免需要您正在做的
float()操作。我们的print语句显示了问题：即使您提供了一个标量浮点数，最小化器仍然会将其转换为 NumPy 数组：(array([ 10000.]), 32.5, 215.1)但这很容易解决：既然 Skyfield 有一个内置的
separation_from()可以很好地处理数组，我们将使用它：sepa = mpos.separation_from(spos) return sepa.degrees我将切换到创建日期的新语法，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 做到了这一点，在天文编程中是传统的。