在使用
scipy.optimize.minimize
评估简单 LLS 问题的协方差矩阵时,使用方法 L-BFGS-B
时出现显着偏差。经过大量研究,我创建了详细的 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 在此配置中如此不同,以及如何调整此代码以获得协方差矩阵,就像我在前面的三个片段中所做的那样。
以下是我在研究过程中发现的一些相关问题,但没有回答当前问题:
为什么带有选项
minimize
的方法 method='L-BFGS-B'
返回如此不同的逆粗麻布矩阵,我应该如何修复它以获得正确的协方差矩阵。
正如 @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]])
这与其他结果相符。