为什么使用Skyfield时，scipy.optimize.minimize（默认）没有移动但却报告成功？

Question

scipy.optimize.minimize使用默认方法时，返回的结果是初始值，没有任何错误或警告消息。而使用Nelder-Mead方法（如这个答案所建议的）可以解决问题，但我想了解：

为什么默认方法会返回起始点作为答案，有没有办法在这种情况下避免这种行为？

请注意，函数separation使用python包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

- uhoh

---

这里有一个相关的Skyfield问题。 - uhoh

哦，抱歉，我漏掉了。我猜问题在于minimize默认使用需要您的函数平滑的算法。如果您的函数不平滑，您最终会陷入垃圾输入垃圾输出的情况。 - cel

我不明白你的问题。如果你的函数可能是非平滑的，那么默认算法就不适合解决你的问题。那你为什么要使用它呢？你是在问如何判断函数是否平滑吗？ - cel

再次回答已经很明确了：“为什么默认方法会返回错误答案而没有警告 - 是否有办法在这种情况下防止‘无警告的错误答案’行为？” - uhoh

文档说明（默认）“BFGS已经证明在非平滑优化中表现良好”。也许应该更新文档，更清楚地建议使用不同的求解器来解决非平滑问题？ - Ted Pudlik

有不同种类的不光滑性，因此对于其他问题可能会正常工作。但是当关键点处不存在导数时，建议使用基于导数的优化器可能不是一个好的普遍建议。 - pv.

Answer 1

1)

您的函数是分段常数（具有小尺度的“阶梯”模式）。它不是处处可微的。

初始猜测时，函数的梯度为零。

默认的BFGS优化器看到零梯度，并根据其标准（基于关于输入函数的假设，这些假设在这种情况下不成立，例如可微性）决定它是局部最小值。

基本上，完全平坦的区域会炸毁优化器。优化器在围绕初始点的小区域内探测函数，在这个区域内，一切都像函数只是一个常数，所以它认为你给了它一个常数函数。因为您的函数并非处处可微，所以几乎所有的初始点都可能在这样的平坦区域内，因此在选择初始点时，这种情况可能会发生而没有坏运气。

还要注意，Nelder-Mead并不免疫这种问题---它只是它的初始单纯形比楼梯的大小大，所以它在更大的区域内探测函数。如果初始单纯形比楼梯的大小小，则优化器的行为类似于BFGS。

2)

一般来说，局部优化器返回局部最优解。这些是否与真正的最优解相符取决于函数的属性。

通常，要查看是否陷入局部最优解，请尝试不同的初始猜测。

此外，在非可微函数上使用基于导数的优化器不是一个好主意。如果函数在“大”尺度上可微，则可以调整数值微分的步长。

由于没有便宜/通用的方法来数值检查一个函数是否处处可微，因此不进行此类检查---而是将其视为优化方法中的假设，必须由输入目标函数和选择优化方法的人员确保。

- pv.

---

好的，这正是我需要理解的！感谢您花时间清楚地解释发生了什么。我已经在一个“补充”答案中添加了一个图，显示即使在1E-03秒以下，Julian日期也是离散的。对于该软件包的大多数用途来说，这很好，但当然不适用于导数。 - uhoh

我已经更改了标题并删除了“失败”一词，因为这正是预期的结果。 - uhoh

Answer 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()

- uhoh

Answer 3

我非常推崇使用print语句来观察算法在运行过程中的表现价值。如果你在separation()函数顶部添加一个print语句，那么你就可以看到最小化程序逐步接近答案：

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位浮点数推动到多精度，它们都被设计为在那个点之前停止。 但是你隐藏了精度：你告诉程序“给我一些秒”，这让它们认为它们可以玩弄甚至非常微小的秒值底端数字，但实际上秒数不仅与小时和日数相结合，而且与年份相结合，在这个过程中，任何在秒底部的微小精度都会丢失 - 尽管最小化器并不知道！

所以让我们向算法公开真正的浮点时间。 在此过程中，我要做三件事：

1. Let's avoid the need for the float() maneuver you are doing. Our print statement shows the problem: that even though you provided a scalar float, the minimizer turns it into a NumPy array anyway:

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

   But that is easy to fix: now that Skyfield has a separation_from() built in that can handle arrays just fine, we will use it:

   sepa = mpos.separation_from(spos)
   return sepa.degrees
   

2. I will switch to the new syntax for creating dates, that Skyfield has adopted as it heads towards 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实现的相当贫乏的优化器。他们需要注意的主要问题是这里发生的事情：优化器需要给出浮点数字以摆动，这些数字在整个过程中都是显著的。

也许我应该研究一下允许时间对象从两个浮点对象组成其时间，以便可以表示更多的秒数位数。我认为AstroPy已经这样做了，在天文学编程中是传统的。

- Brandon Rhodes

---

好的这很棒！确实我可以使用print，但是通过在SE这里寻求更好的理解，我从这里学到了很多“底层”知识。我期待1.0版本 - 这些新方法真的值得等待。我认为能够搜索像日全食或掩星这样的标准非常棒。时间对象是否会有自己的简单添加方法？ - uhoh

哦，我回头看了一下，发现 float() # 似乎是必要的 实际上只需要用于 .topos()，如果我不小心输入纬度或经度为整数而没有小数点，则无法正常工作。 - uhoh