理解scipy盆地跳跃优化功能的例子

Posted 2023-03-12

【中文标题】理解scipy盆地跳跃优化功能的例子

【英文标题】：Example to understand scipy basin hopping optimization function

【发布时间】：2014-04-08 01:53:25

【问题描述】：

我在 scipy 中遇到了 basin hopping algorithm 并创建了一个简单的问题来了解如何使用它，但它似乎无法正常解决该问题。可能是我做错了什么。

代码如下：

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])

【问题讨论】：

尝试对边界外的 x 进行二次惩罚：scipy-optimize-leastsq-with-bound-constraints

【参考方案1】：

您正在测试的函数使用了一种称为 Metropolis-Hastings 的方法，该方法可以修改为一种称为模拟退火的过程，可以以随机方式优化函数。

其工作方式如下。首先你选择一个点，比如你的点x0。从那时起，您会生成一个随机扰动（这称为“提案”）。一旦有提议的扰动，您就可以通过将扰动应用于当前输出来获得新点的候选人。所以，你可以把它想象成x1 = x0 + perturbation。

在常规的旧梯度下降中，perturbation 项只是一个确定性计算的量，就像梯度方向上的一步。但在 Metropolis-Hastings 中，perturbation 是随机生成的（有时使用梯度作为随机去向的线索……但有时只是随机生成，没有任何线索）。

当你得到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.

所以它从不接受 posposal 点。输出并没有告诉您它找到了解决方案。它告诉您，探索可能的解决方案的随机干扰不断导致优化器看起来越来越好，但始终无法满足您的标准。它无法一路回到[0,1] 以获得确实满足mybounds 的积分。

【讨论】：

如果我们设置stepsize=0，也就是拒绝候选人，对吧？

那么，我们可以只使用 take_step 而不使用 mybounds 吗？

【参考方案2】：

按照您的编码，盆地跳跃的行为是将扰动与局部最小化结合起来。

由于局部优化部分，您的例程不断产生不可接受的结果。本质上，您使用的 BFGS 例程是完全不受约束的，因此它遵循梯度到负无穷大。然后，此结果会反馈到您的检查器中。

因此，无论您的流域跳跃点 x1 在哪里，BFGS 部分总是会变成一个巨大的负值。

您正在使用的基准函数x - 4 不是这里的理想目标。检查例如Rastrigin function。如果您确实需要优化线性函数，则有一整类算法可以做到这一点（参见***上的Linear Programming）。

【讨论】：

您能否澄清事实上不准确的部分，关于在 SciPy 中实现的盆地跳跃是做什么的？我不介意投反对票或其他任何事情，但很想得到纠正。据我了解，局部最小化部分是跟随负梯度的，而不是全局扰动部分。我最初发布的原因是我觉得区分不清楚，尽管我认为在说您的帖子不准确之前我应该​​更加小心。

OP 有一个自定义的接受/拒绝测试。如果他只是随机选择没有局部优化的点，他的优化器将保持在界限内。此外，如果您在 MH 测试和 OP 的边界测试之间进行逻辑与，则会发生同样的事情并且会收敛，因为优化器的自定义部分会拒绝该点。

这是我看到的机制...优化器选择x0 + delta。如果它不在界限内或它没有改善目标，它就会拒绝它。如果它改进并保持在界限内，那就成为新的x0。因此，即使分布均匀，在大量的抽奖中，每一个delta 都不太可能导致越界。我们可以同意吗？很难将这种差异解析到评论线程中。

不，这是不对的。该算法选择x0 + delta，然后，无论是否在边界内，从该点执行局部确定性优化器（在这种情况下运行到负无穷大）。在达到最大迭代限制后（或者在其他情况下，当它找到一个最优值时），然后它检查它是否在边界内以及它是否改进了适应度函数，并且 Metropolis 是否接受/拒绝（你可以看到它here）。

如果整个事情都不接受，那么自适应步进器向下调整 0.9 倍，您将从 [0.55, 1.45] 均匀抽取样本（再次从原点开始，仍然有 0.5 的机会不在区间内）。我确实同意，仅此一点，再加上足够的重复抽奖，就使得在 [0,1] 内没有此类抽奖的概率非常小。阅读代码时感到好奇，它甚至可能输出 OP 给出的答案，而不是总是输出初始点。

【参考方案3】：

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。