为什么我的scipy.optimize.minimize(method="newton-cg")函数卡在局部最大值上？

Question

我想找到一个依赖于两个变量的函数的局部最小值。为此，我打算使用scipy.optimize.minimize函数，并选择"newton-cg"方法，因为我可以通过数学公式计算雅可比矩阵和海森矩阵。

但是，当我从一个局部最大值开始猜测时，函数在第一次迭代时就成功结束了，停在了局部最大值上，尽管海森矩阵是负的。我写了一段简单的测试代码，能够重现这个问题：

import numpy as np
import scipy.optimize as o

def get_f_df(var):
    x = var[0]
    y = var[1]

    f = np.cos(x) + np.cos(y)
    df_dx = -np.sin(x)
    df_dy = -np.sin(y)

    return f, (df_dx, df_dy)

def hess(var):
    x = var[0]
    y = var[1]

    f_hess = np.zeros((2,2))
    f_hess[0,0] = -np.cos(x)
    f_hess[1,1] = -np.cos(y)
    return f_hess

min = o.minimize(get_f_df, (0, 0), jac=True, hess=hess, method="newton-cg")
print(min)

结果是：

message: Optimization terminated successfully.
 success: True
  status: 0
     fun: 2.0
       x: [ 0.000e+00  0.000e+00]
     nit: 1
     jac: [-0.000e+00 -0.000e+00]
    nfev: 1
    njev: 1
    nhev: 0

如果我使用 hess=None、hess='cs'、hess='2-point' 或 hess='3-point' 代替我自定义的海森矩阵函数，结果也是一样的。此外，如果我使用其他方法，比如 'dogleg'、'trust-ncg'、'trust-krylov'、'trust-exact' 或 'trust-constr'，基本上也会得到相同的结果，除了 nhev = 1，但结果在 x = [0,0] 时仍然是错误的。

要么我在这里做错了什么（这很可能），要么 minimize 函数，特别是 "newton-cg" 方法存在重大问题（这不太可能）。

关于后者，我也查看了源代码，想看看是否有什么问题，结果发现了一些奇怪的地方。不过，我并不完全理解整个代码，所以对我的担忧是否合理有些不确定：

让我们看看源代码

当调用 minimize 函数并指定方法为 "newton-cg" 时，它会跳转到 _minimize_newtoncg 函数（可以在这里查看源代码）。我想详细说明我认为这里发生了什么：

在第2168行，A = sf.hess(xk) 计算了海森矩阵，依赖于 xk，而 xk 最初是起始猜测 x0。在我的测试案例中，海森矩阵当然是

A = [[f_xx, f_xy], [f_xy, f_yy]]

其中 f_ij 是 f 关于 i 和 j 的导数。在我的情况下，f_xy = f_yx 也是成立的。

接下来，在第2183行，Ap = A.dot(psupi) 计算了海森矩阵 A 和 psupi 的乘积。psupi 基本上等于 b，也就是在 xk 处的 f 的负梯度。所以 Ap = A.dot(psupi) 的结果是

Ap = [f_xxf_x + f_xyf_y, f_xyf_x + f_yyf_y].

现在谈谈（可能的）问题

接下来，在第2186行，通过 np.dot(psupi, Ap) 计算曲率 curv。如上所述，psupi 是 f 的负梯度，因此结果是

curv = f_xxf_x² + 2 f_xyf_xf_y + f_yyf_y².

然而，所有这些导数都是在 xk 处计算的，而 xk 最初等于起始参数 x0。如果起始参数正好在一个局部最大值上，导数 f_x 和 f_y 就会等于 0。因此，curv = 0。这导致下一行的 for 循环中断，从而跳过更新 xsupi、psupi 和所有其他参数。因此，pk 变成了 [0,0]，并且调用了 _line_search_wolfe12，基本上使用所有起始参数。在这里，我对源代码的理解停止了，不过，我觉得在 curv = 0 后，事情就已经出错了，导致了 for 循环的中断。

编辑 - 我为什么需要这个：

由于我收到了些反馈，有人问我为什么不直接使用其他起始猜测，我想简单解释一下我的实际目标，也许这能帮助你帮助我。

我想使用宏自旋模型模拟磁滞回线。为此，我需要在每个外部磁场步长中找到能量景观的局部最小值，从饱和状态开始。在那里，磁宏自旋与外部磁场之间的角度为 0°。在饱和状态下，它们处于能量最小值。如果我减小外部磁场，我必须将之前的场步的角度作为新的起始猜测。随着外部磁场的减小，饱和状态下的能量景观中的局部最小值会转变为局部最大值。

起初，我对上一个场值的角度进行了局部最小化，并加上和减去一个小增量。我的想法是，只要不在局部最大值上，结果应该是相同的。然后我会选择找到最低最小值的那个值。不知道为什么，我选择的增量值对结果有巨大的影响。我的增量值通常在 0.0001 到 0.01 之间，而我的角度在 pi 的范围内（-3.141 到 3.141）。因此，我有点放弃了这个想法。

接下来，我想我会尝试检查一下我是否确实在局部最大值上，也许更考虑梯度，而不是最终的能量值来决定方向。如果那样有效，我会在这里更新。

更新：如果我处于局部最大值或鞍点，我现在会检查在加减增量值下的梯度，并选择梯度最高的位置作为新的起始猜测。这种方法有点有效，但确切的增量值对结果的影响仍然超出了我的预期。我想我需要尝试找到一个适合大多数情况的理想值。

（Matt Haberland 提出的无导数求解器似乎也没有真正帮助我。它们有点有效，但也有点无效。）

优化算法梯度下降数值优化雅可比矩阵局部最小值海森矩阵函数迭代磁滞回线

为什么我的scipy.optimize.minimize(method="newton-cg")函数卡在局部最大值上？

让我们看看源代码

现在谈谈（可能的）问题

编辑 - 我为什么需要这个：

1 个回答

撰写回答