为何scipy.optimize.minimize(newton-cg)在局部最大值处停滞？

我需要寻找一个二元函数的局部最小值，计划使用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参数改为hess=None、hess='cs'、hess='2-point'或hess='3-point'，结果依然相同。
• 使用'dogleg'、'trust-ncg'、'trust-krylov'、'trust-exact'或'trust-constr'等方法时，除nhev=1外，结果仍停留在x=[0,0]。

源码分析（newton-cg方法） #

调用minimize并指定method="newton-cg"时，会进入_minimize_newtoncg函数，关键流程如下：

1. 第2168行A = sf.hess(xk)：基于初始猜测值x0（即xk）计算海森矩阵，测试案例中该矩阵为：

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

   其中f_ij表示f对i和j的偏导数，满足f_xy = f_yx。
2. 第2183行Ap = A.dot(psupi)：计算海森矩阵与psupi（即f在xk处的负梯度）的乘积，即：

   Ap = [f_xx f_x + f_xy f_y, f_xy f_x + f_yy f_y]
   

疑似问题点 #

第2186行通过np.dot(psupi, Ap)计算曲率curv，推导可得：

curv = f_xx f_x² + 2 f_xy f_x f_y + f_yy f_y²

当初始值处于局部最大值时，f_x和f_y均为0，导致curv=0，触发循环中断，跳过xsupi、psupi等参数的更新，最终搜索方向pk=[0,0]，线搜索仍使用初始参数，从而终止在局部最大值。

需求背景 #

实际应用场景为：用宏自旋模型模拟磁滞回线，需从饱和状态开始，为每个外磁场步长寻找能量 landscape的局部最小值。饱和状态下，宏自旋与外磁场的夹角为0°，处于能量最小值；当外磁场减弱时，该饱和状态的局部最小值会转变为局部最大值。

• 最初尝试以上一步的角度及±微小增量作为初始值，预期只要不在局部最大值，结果应一致，再选择找到的最小最小值。但增量值（0.0001-0.01，角度范围为-π到π）对结果影响极大，因此放弃该方案。
• 后续计划仅在确认处于局部最大值时，根据梯度而非最终能量值选择方向。

更新 #

目前的解决方案：处于局部最大值或鞍点时，检查±增量位置的梯度，选择梯度最大的位置作为新初始值。该方法虽能运行，但增量值对结果的影响仍超出预期，需寻找适配多数场景的理想增量值。（无导数求解器效果也不理想）

内容的提问来源于stack exchange，提问作者legandro