为啥 scipy.optimize.minimize （默认）在不使用 Skyfield 的情况下报告成功？
Posted
技术标签:
【中文标题】为啥 scipy.optimize.minimize （默认）在不使用 Skyfield 的情况下报告成功？【英文标题】：Why does scipy.optimize.minimize (default) report success without moving with Skyfield?为什么 scipy.optimize.minimize （默认）在不使用 Skyfield 的情况下报告成功？ 【发布时间】：2016-07-06 18:41:38 【问题描述】：scipy.optimize.minimize 使用默认方法返回初始值作为结果，没有任何错误或警告消息。虽然使用this answer 建议的 Nelder-Mead 方法可以解决问题，但我想了解：
为什么默认方法会返回没有警告的错误答案作为答案的起点 - 有没有一种方法可以防止“没有警告的错误答案” 在这种情况下避免这种行为？
注意，函数separation 使用python 包Skyfield 生成要最小化的值，这不能保证平滑，这可能是Simplex在这里更好的原因。
结果：
测试结果：[2.14159739]'正确'：2.14159265359初始：0.0
默认结果：[ 10000.] '正确'：13054 初始值：10000
Nelder-Mead 结果：[ 13053.81011963] '正确'：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
【问题讨论】：
这是一个related Skyfield 的问题。 抱歉，错过了。我想问题是minimize 默认使用的算法要求你的函数是平滑的。如果你的函数不流畅，你最终会陷入垃圾进垃圾出的境地。
我不明白你在问什么。如果您的函数可能不平滑，则默认算法根本不适合您的问题。那你为什么要使用它？你问如何判断函数是否平滑？
文档指出（默认）“BFGS 已证明即使对于非平滑优化也有良好的性能”。也许应该更新文档以更清楚地建议对非平滑问题使用不同的求解器？
不平滑的种类有很多种，所以对于其他问题可能还可以。但是，当导数不存在时（在关键点），建议基于导数的优化器可能不是一个好的一般建议。
【参考方案1】：
1)
您的函数是分段常数（具有小规模“阶梯”模式）。 它并非处处可微。
初始猜测时函数的梯度为零。
默认的 BFGS 优化器看到零梯度，并根据其标准（基于关于输入函数的假设在这种情况下不正确，例如可微性）确定它是局部最小值。
基本上，完全平坦的区域会轰炸优化器。优化器在初始点周围的小而平坦的区域中探测函数，在那里一切看起来函数只是一个常数，所以它认为你给了它一个常数函数。因为你的函数不是处处可微的，所以几乎所有的初始点都可能在这样的平坦区域内，所以在选择初始点时不会运气不好。
还请注意，Nelder-Mead 并非对此免疫 --- 它只是碰巧它的初始单纯形大于楼梯的大小，因此它会在更大的点周围探测函数。如果初始单纯形小于阶梯大小，则优化器的行为类似于 BFGS。
2)
一般答案：局部优化器返回局部最优值。这些是否与真正的最优值一致取决于函数的属性。
一般来说，要查看您是否陷入局部最优，请尝试不同的初始猜测。
此外，对不可微函数使用基于导数的优化器也不是一个好主意。如果函数在“大”尺度上可微分，则可以调整数值微分的步长。
因为没有廉价/通用的方法来检查一个函数是否处处可微，所以没有进行这样的检查 --- 相反，它是优化方法中的一个假设，必须由输入目标函数并选择的人确保优化方法。
【讨论】：
好的，这正是我需要理解的！感谢您花时间清楚地解释发生了什么。我在“补充”答案中添加了一个情节，表明即使 JulianDate 低于 1E-03 秒也是离散的。这对于包的大多数用途都很好，但当然不适用于衍生产品。 我已更改标题并删除了“失败”一词，因为这正是预期的结果。【参考方案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()
【讨论】：
【参考方案3】：我不能太高估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 位浮点数的精确度，并且它们都被设计为在该点之前停止。但是你隐藏了精确度：你告诉例行程序“给我几秒钟”，这让他们认为他们可以摆弄即使是非常微小的秒值低端数字，而实际上秒数正在与不仅仅是几小时和几天，还有几年，在这个过程中，任何微小的精确到秒的底部都会丢失——尽管最小化器不知道！
所以让我们将真正的浮点时间暴露给算法。在这个过程中，我会做三件事：
让我们避免使用您正在执行的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 已经做到了这一点，而且它在天文学编程中是传统的。
【讨论】：
好的，这太棒了！事实上，我本可以使用print，但是我永远不会像我在这里通过在 SE 中寻求更好的理解来了解这么多“幕后”发生的事情。我很期待1.0——这些新方法真的值得等待。我认为能够搜索日全食或掩星等标准非常棒。时间对象会有自己的简单add method 吗？
哦，我回头一看，发现 float() # seems necessary 实际上只需要 .topos()
以上是关于为啥 scipy.optimize.minimize （默认）在不使用 Skyfield 的情况下报告成功？的主要内容，如果未能解决你的问题，请参考以下文章
为啥使用 glTranslatef？为啥不直接更改渲染坐标？
为啥 DataGridView 上的 DoubleBuffered 属性默认为 false，为啥它受到保护？