6

scipy.optimize.minimize使用默认方法返回初始值作为结果，没有任何错误或警告消息。虽然使用此答案建议的 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
4

3 回答 3

5

1)

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

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

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

基本上，完全平坦的区域会轰炸优化器。优化器在初始点周围的小而平坦的区域中探测函数，在那里一切看起来函数只是一个常数，所以它认为你给了它一个常数函数。因为你的函数不是处处可微的，所以几乎所有的初始点都可能在这样的平坦区域内，因此在选择初始点时可能会发生这种情况而不会运气不好。

另请注意，Nelder-Mead 也不能对此免疫 --- 它只是碰巧它的初始单纯形大于楼梯的大小，因此它会在更大的点周围探测函数。如果初始单纯形小于楼梯大小，则优化器的行为类似于 BFGS。

2)

一般答案：局部优化器返回局部最优值。这些是否与真正的最优值一致取决于函数的属性。

一般来说，要查看您是否陷入局部最优，请尝试不同的初始猜测。

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

因为没有廉价/通用的方法来检查一个函数是否处处可微，所以没有进行这样的检查——相反，它是优化方法中的一个假设，必须由输入目标函数并选择优化方法的人来确保.

于 2016-03-21T22:15:59.693 回答
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()
于 2016-03-22T01:16:02.673 回答
2

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 位浮点数的精确度，并且它们都被设计为在该点之前停止。但是你隐藏了精确度：你告诉例行程序“给我几秒钟”，这让他们认为他们可以摆弄即使是非常微小的秒值低端数字，而实际上秒数正在与不仅仅是几小时和几天，还有几年，在这个过程中，任何微小的精确到秒的底部都会丢失——尽管最小化器不知道！

因此，让我们将真正的浮点时间暴露给算法。在这个过程中，我会做三件事：

  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 已经做到了这一点，它在天文学编程中是传统的。

于 2016-04-03T00:12:28.580 回答