Scipy optimize.minimize使用L-BFGS-B时的协方差偏差问题

问题原因 #

• L-BFGS-B是带边界约束的优化算法，它维护的有限内存近似逆海森矩阵，是基于无约束区域的曲率信息构建的。即便你没给LLS问题加显式约束，算法内部的边界处理逻辑也可能干扰逆海森矩阵的近似精度——毕竟LLS本质是无约束问题，带约束的优化器的曲率近似逻辑不匹配。
• scipy中minimize返回的hess_inv对L-BFGS-B而言，只是迭代过程中维护的近似矩阵，并非严格的真实逆海森矩阵（或其缩放版本）。而curve_fit、least_squares、BFGS这些方法，要么是针对最小二乘问题推导了精确的协方差计算逻辑（比如curve_fit用雅可比伪逆），要么是无约束优化下的逆海森近似更贴合真实情况。

修正方案 #

方法1：手动计算协方差矩阵 #

LLS问题的协方差矩阵可以通过雅可比矩阵直接推导，不需要依赖优化器返回的近似结果，步骤如下：

1. 定义残差函数，返回每个样本的残差值（即观测值与模型预测值的差）。
2. 用L-BFGS-B得到最优参数。
3. 计算残差在最优参数处的雅可比矩阵（可通过数值微分工具实现）。
4. 按公式计算协方差矩阵：cov = (J.T @ J)^{-1} * sigma²，其中sigma²是残差的方差估计，公式为np.sum(残差²) / (样本数 - 参数个数)。

示例代码：

import numpy as np
from scipy.optimize import minimize
from numdifftools import Jacobian

# 构造简单线性最小二乘问题
x = np.linspace(0, 10, 100)
true_params = [2.5, 1.2]
y = true_params[0] * x + true_params[1] + np.random.normal(0, 0.5, size=len(x))

# 定义残差和目标函数
def residuals(params):
    return y - (params[0] * x + params[1])

def objective(params):
    return np.sum(residuals(params)**2)

# L-BFGS-B优化得到最优参数
opt_result = minimize(objective, x0=[1, 0], method='L-BFGS-B')
params_opt = opt_result.x

# 计算雅可比矩阵
jac = Jacobian(residuals)(params_opt)

# 计算残差方差估计
sigma_sq = np.sum(residuals(params_opt)**2) / (len(y) - len(params_opt))

# 计算协方差矩阵
cov_matrix = np.linalg.inv(jac.T @ jac) * sigma_sq

print("手动计算的协方差矩阵:\n", cov_matrix)

方法2：切换无约束优化方法 #

如果问题本身不需要参数边界限制，直接用BFGS方法替代L-BFGS-B，它返回的hess_inv更接近真实逆海森矩阵，只需乘以残差方差就能得到正确的协方差矩阵。

方法3：改用least_squares工具 #

scipy.optimize.least_squares是专门针对最小二乘问题的优化器，支持边界约束（如果有需求），会直接返回雅可比矩阵，可基于此精确计算协方差矩阵，比L-BFGS-B更适配LLS场景。

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