理解scipy的basin hopping优化函数的示例

Question

我在scipy中遇到了一个叫做盆地跳跃算法的东西，于是我创建了一个简单的问题来理解它的用法，但它似乎在这个问题上没有正确工作。可能是我做错了什么。

这是我的代码：

import scipy.optimize as spo
import numpy as np
minimizer_kwargs = {"method":"BFGS"}    
f1=lambda x: (x-4)
def mybounds(**kwargs):
    x = kwargs["x_new"]
    tmax = bool(np.all(x <= 1.0))
    tmin = bool(np.all(x >= 0.0))
    print x
    print tmin and tmax
    return tmax and tmin


def print_fun(x, f, accepted):
      print("at minima %.4f accepted %d" % (f, int(accepted)))
x0=[1.]     
spo.basinhopping(f1,x0,accept_test=mybounds,callback=print_fun,niter=200,minimizer_kwargs=minimizer_kwargs)

它给出的结果是 x: array([ -1.80746874e+08])

Answer 1

Yike Lu已经指出了问题：你的边界限制只在最顶层生效，但局部优化器BFGS对此一无所知。

一般来说，在优化过程中使用“硬性”边界往往不是个好主意，因为大多数算法不会允许任何可能直接到达你允许的空间边界的路径触碰到边界，否则就会被终止。你可以想象在上面的例子中，找到最优解（x=0）是多么困难，如果不尝试x=-0.0000001，根本无法发现自己走得太远了，需要稍微退回来一点。

现在，有些算法可以通过转换输入数据来解决这个问题（在scipy.optimize中，这些算法接受边界作为参数），但一般的解决方案是这样的：

你需要更新你的成本函数，让它在输入超出允许范围时迅速增加：

def f1(x):
    cost_raw = (x-4)
    if   x >= 1.0: cost_overrun = (1000*(x-1))**8
    elif x <= 0.0: cost_overrun = (1000*(-x))**8
    else: cost_overrun = 0.0

    return(cost_raw + cost_overrun)

这样，任何优化器都会看到成本函数在增加，一旦超出边界，就会立刻回到允许的范围。这并不是严格的强制，但优化器本身就是在不断接近最优解，所以根据你的需求，可以调整惩罚函数，让增加的幅度更大或更小。有些优化器喜欢连续的导数（因此使用幂函数），而有些则可以接受固定的步长——在这种情况下，你可以简单地在超出边界时加上10000。

Answer 2

你写的这个“盆地跳跃”算法的工作方式是把一些小的扰动和局部最小化结合在一起。

你的程序一直产生不理想的结果，主要是因为局部优化的部分。简单来说，你用的BFGS算法是完全没有限制的，所以它会沿着梯度一直走到负无穷。这种结果又会反馈到你的检查器里。

所以，不管你的“盆地跳跃”点x1在哪里，BFGS部分总是会得到一个非常大的负值。

你使用的基准函数x - 4并不是一个理想的目标。可以看看比如Rastrigin函数。如果你真的需要优化一个线性函数，还有一整套算法可以做到这一点（可以参考维基百科上的线性规划）。

Answer 3

你正在测试的这个函数使用了一种叫做Metropolis-Hastings的方法，这种方法可以改造成一种叫做模拟退火的程序，能够以随机的方式优化函数。

它的工作原理是这样的。首先，你选择一个点，比如你的点x0。从这个点出发，你会生成一个随机的扰动（这个叫做“提议”）。一旦有了这个提议，你就可以通过将扰动应用到当前点上，得到一个新的候选点。所以，你可以把它想象成x1 = x0 + 扰动。

在传统的梯度下降法中，扰动是一个确定性计算出来的量，比如朝着梯度的方向迈出一步。但在Metropolis-Hastings中，扰动是随机生成的（有时会用梯度作为线索来决定随机移动的方向……但有时也完全是随机的，没有任何线索）。

当你得到了x1之后，你需要问自己：“我通过随机扰动x0做得对吗，还是把事情搞得一团糟？”这其中有一部分是关于是否保持在某些范围内，比如你的mybounds函数。另一部分则是看在新点上的目标函数值变得好多少或坏多少。

所以，有两种情况会拒绝x1的提议：第一，它可能违反了你设定的范围，成为一个不合适的点；第二，从Metropolis-Hastings的接受/拒绝评估步骤来看，它可能只是一个很糟糕的点，应该被拒绝。在这两种情况下，你都会拒绝x1，而是将x1 = x0，假装你仍然停留在同一个地方，重新尝试。

这和梯度类方法形成了对比，在梯度方法中，无论如何，你总会至少有一些移动（朝着梯度的方向迈出一步）。

好了，抛开这些不谈，让我们想想basinhopping函数是如何工作的。从文档中我们可以看到，典型的接受条件是通过take_step参数来访问的，文档中说：“默认的步进例程是坐标的随机位移，但其他步进算法可能对某些系统更好。”所以，即使不考虑你的mybounds范围检查器，这个函数也会随机位移坐标来生成新的尝试点。而且，由于这个函数的梯度只是常数1，它总是会朝着负梯度的方向迈出相同的大步（用于最小化）。

在实际操作中，这意味着提议的x1点总是会在区间[0,1]之外，而你的范围检查器总是会否决它们。

当我运行你的代码时，我看到这种情况一直在发生：

In [5]: spo.basinhopping(f1,x0,accept_test=mybounds,callback=print_fun,niter=200,minimizer_kwargs=minimizer_kwargs)
at minima -180750994.1924 accepted 0
[ -1.80746874e+08]
False
at minima -180746877.5530 accepted 0
[ -1.80746873e+08]
False
at minima -180746877.3896 accepted 0
[ -1.80750991e+08]
False
at minima -180750994.7281 accepted 0
[ -1.80746874e+08]
False
at minima -180746878.2433 accepted 0
[ -1.80746874e+08]
False
at minima -180746877.5774 accepted 0
[ -1.80746874e+08]
False
at minima -180746878.3173 accepted 0
[ -1.80750990e+08]
False
at minima -180750994.3509 accepted 0
[ -1.80750991e+08]
False
at minima -180750994.6605 accepted 0
[ -1.80746874e+08]
False
at minima -180746877.6966 accepted 0
[ -1.80746874e+08]
False
at minima -180746877.6900 accepted 0
[ -1.80750990e+08]
False
at minima -180750993.9707 accepted 0
[ -1.80750990e+08]
False
at minima -180750994.0494 accepted 0
[ -1.80750991e+08]
False
at minima -180750994.5824 accepted 0
[ -1.80746874e+08]
False
at minima -180746877.5459 accepted 0
[ -1.80750991e+08]
False
at minima -180750994.6679 accepted 0
[ -1.80750991e+08]
False
at minima -180750994.5823 accepted 0
[ -1.80750990e+08]
False
at minima -180750993.9308 accepted 0
[ -1.80746874e+08]
False
at minima -180746878.0395 accepted 0
[ -1.80750991e+08]
False

# ... etc.

所以它从来没有接受提议的点。输出并没有告诉你它找到了一个解决方案。它告诉你，随机扰动以探索可能的解决方案，结果总是得到一些看起来对优化器越来越好的点，但这些点却一直无法满足你的标准。它无法找到回到[0,1]的路径，以获取那些确实满足mybounds的点。