使用 scipy.optimize.minimize 的 L-BFGS-B 方法时的协方差偏差

Question

在使用

scipy.optimize.minimize

评估简单 LLS 问题的协方差矩阵时，使用方法

L-BFGS-B

时出现显着偏差。经过大量研究，我创建了详细的 MCVE 来强调这个问题。

MCVE

那么好用

curve_fit

方法

假设我们想要将以下模型拟合到某些试验数据集：

import numpy as np
from scipy import optimize

def model(x, a, b, c):
    return a * x[:, 0]**2 + b * x[:, 0] + c

np.random.seed(12345)
X = np.linspace(-1, 1, 30).reshape(-1, 1)
y = model(X, 3, 2, 1)
s = 0.1 * np.ones(y.size)
n = s * np.random.normal(size=y.size)
yn = y + n

popt, pcov = optimize.curve_fit(model, X, yn, sigma=s, absolute_sigma=True)
popt, pcov

# (array([3.00876769, 2.00038216, 1.03068484]),
#  array([[ 3.29272571e-03, -1.11752553e-11, -1.17327007e-03],
#         [-1.11752553e-11,  9.35483883e-04,  2.98424026e-12],
#         [-1.17327007e-03,  2.98424026e-12,  7.51395072e-04]]))

curve_fit

方法直接提供协方差矩阵，将其视为预期结果。现在我们想使用不同的

optimize

方法重现这个协方差矩阵。

用基线方法检查

least_squares

使用

least_squares

函数（由

curve_fit

在后台使用）和 recipe

scipy

用于稳健地评估协方差，从而得到合规的结果：

def residuals_factory(x, y, sigma=1.):
    def wrapped(beta):
        return (y - model(x, *beta)) / sigma
    return wrapped

residuals = residuals_factory(X, yn, sigma=s)

sol0 = optimize.least_squares(residuals, x0=[1, 1, 1])
#     message: `xtol` termination condition is satisfied.
#     success: True
#      status: 3
#         fun: [-5.954e-01  9.965e-02 ... -3.855e-01  9.455e-01]
#           x: [ 3.009e+00  2.000e+00  1.031e+00]
#        cost: 16.317663438370303
#         jac: [[-1.000e+01  1.000e+01 -1.000e+01]
#               [-8.668e+00  9.310e+00 -1.000e+01]
#               ...
#               [-8.668e+00 -9.310e+00 -1.000e+01]
#               [-1.000e+01 -1.000e+01 -1.000e+01]]
#        grad: [ 1.070e-07 -1.632e-06  1.722e-07]
#  optimality: 1.632403261453419e-06
# active_mask: [ 0.000e+00  0.000e+00  0.000e+00]
#        nfev: 4
#        njev: 3

U, s, Vh = np.linalg.svd(sol0.jac, full_matrices=False)
tol = np.finfo(float).eps*s[0]*max(sol0.jac.shape)
w = s > tol
cov = (Vh[w].T/s[w]**2) @ Vh[w]

# array([[ 3.29272571e-03,  5.77282993e-12, -1.17327008e-03],
#        [ 5.77282993e-12,  9.35483875e-04,  1.14915683e-12],
#        [-1.17327008e-03,  1.14915683e-12,  7.51395089e-04]])

在这种情况下，我们可以重现协方差。

用

minimize

方法检查

使用

minimize

，我们可以通过声明卡方损失函数而不是残差来执行相同的操作。

def loss_factory(x, y, sigma=1.):
    def wrapped(beta):
        return 0.5 * np.sum(np.power((y - model(x, *beta)) / sigma, 2))
    return wrapped

loss = loss_factory(X, yn, sigma=s)

sol1 = optimize.minimize(loss, x0=[1, 1, 1])
#  message: Optimization terminated successfully.
#  success: True
#   status: 0
#      fun: 16.31766343837042
#        x: [ 3.009e+00  2.000e+00  1.031e+00]
#      nit: 5
#      jac: [ 1.907e-06  1.192e-06 -4.768e-07]
# hess_inv: [[ 3.293e-03  1.104e-11 -1.173e-03]
#            [ 1.104e-11  9.355e-04 -1.734e-12]
#            [-1.173e-03 -1.734e-12  7.514e-04]]
#     nfev: 32
#     njev: 8

我们再次能够得到协方差矩阵（参见

hess_inv

）。使用方法

BFGS

也可以按预期工作。

使用

L-BFGS-B

方法的问题

但是当我们使用

L-BFGS-B

方法时，返回的逆Hessian矩阵不符合之前计算的协方差矩阵：

sol2 = optimize.minimize(loss, x0=[1, 1, 1], method="L-BFGS-B")
#  message: CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH
#  success: True
#   status: 0
#      fun: 16.317663438454282
#        x: [ 3.009e+00  2.000e+00  1.031e+00]
#      nit: 8
#      jac: [-5.720e-05  2.917e-04 -4.182e-04]
#     nfev: 36
#     njev: 9
# hess_inv: <3x3 LbfgsInvHessProduct with dtype=float64>

sol2.hess_inv.todense()
# array([[ 0.03719202,  0.02371079, -0.02612526],
#        [ 0.02371079,  0.02215422, -0.01929008],
#        [-0.02612526, -0.01929008,  0.01984616]])

我想知道为什么逆 hessian 在此配置中如此不同，以及如何调整此代码以获得协方差矩阵，就像我在前面的三个片段中所做的那样。

相关问题

以下是我在研究过程中发现的一些相关问题，但没有回答当前问题：

• 如何使用 scipy.optimize.least_squares 计算标准差误差
• 在Scipy中，curve_fit如何以及为什么计算参数估计的协方差
• https://stats.stackexchange.com/questions/38115/how-do-i-interpret-the-covariance-matrix-from-a-curve-fit
• 如何从Python最小化函数获取Hessian矩阵？

问题

为什么带有选项

minimize

的方法

method='L-BFGS-B'

返回如此不同的逆粗麻布矩阵，我应该如何修复它以获得正确的协方差矩阵。

Answer 1

正如 @NickODell 指出的，方法

L-BFGS-B

计算 Hessian 的有限内存近似值。

底层函数

m

有一个

fmin_l_bfgs_b开关，但它在

minimize

签名中不可用，因此我们无法控制这个近似值的精确度。

作为解决方法，可以在最小化后计算 Hessian 矩阵：

from autograd import hessian

H = hessian(loss)
C = np.linalg.inv(H(sol2.x))

# array([[ 3.29272573e-03,  1.28301871e-19, -1.17327009e-03],
#        [ 2.44025127e-19,  9.35483871e-04, -6.92261149e-20],
#        [-1.17327009e-03, -3.24227332e-20,  7.51395089e-04]])

这与其他结果相符。