scipy.optimize.minimize参数方差求解：Hessian逆矩阵结果异常

嘿，这个问题我太熟了！咱们捋清楚为啥你的结果差了100倍，以及正确的计算方式：

首先得明确：你用Hessian逆对角线得到的结果差100倍，大概率是没考虑残差的方差估计和自由度缩放，也就是你看到的“测量次数缩放”的核心逻辑。咱们一步步拆解：

核心原理：协方差矩阵的推导 #

非线性最小二乘的目标是最小化残差平方和 ( S(\theta) = \sum_{i=1}^n (y_i - f(x_i, \theta))^2 )，其中 ( \theta ) 是待拟合参数，( n ) 是测量次数（数据点数量），( p ) 是参数个数。

参数的协方差矩阵正确公式是：
[ \text{Cov}(\hat{\theta}) = \sigma^2 \cdot (J^T J)^{-1} ]
这里的关键项：

• ( J ) 是残差的Jacobian矩阵（每行对应一个数据点对参数的偏导）
• ( \sigma^2 ) 是噪声的方差估计，计算方式为 ( \sigma^2 = \frac{S(\hat{\theta})}{n - p} )，其中 ( n-p ) 是自由度（测量次数减去参数个数，用来修正过拟合的偏差）

为啥你的结果大了100倍？ #

你直接用Hessian逆的对角线，大概率没做两件事：

1. 没除以自由度 ( n-p )：如果你的数据点数量比参数多100个（比如n-p=100），那直接用 ( S(\hat{\theta}) ) 而不是 ( S(\hat{\theta})/(n-p) ) 来估计 ( \sigma^2 )，结果就会大100倍，刚好对应你的情况。
2. 没考虑优化目标的缩放因子：很多优化器（包括scipy的minimize）会把目标函数定义为 ( 0.5 \times S(\theta) )（这样梯度和Hessian的形式更简洁），这时候目标函数的Hessian ( H \approx J^T J )，所以 ( H^{-1} = (J^T J)^{-1} )；但如果你自己定义的目标函数是原始的 ( S(\theta) )（没有乘0.5），那 ( H \approx 2J^T J )，这时候 ( (J^T J)^{-1} = H^{-1}/2 )，协方差矩阵需要乘以2才能修正。

用scipy.optimize.minimize的实操步骤 #

我给你举个具体的代码例子，比如拟合非线性函数 ( y = a e^{-b x} + c )：

1. 生成模拟数据 #

import numpy as np
from scipy.optimize import minimize

# 真实参数
true_params = [2.5, 0.3, 0.8]
a_true, b_true, c_true = true_params

# 生成带噪声的数据
x = np.linspace(0, 10, 100)
y = a_true * np.exp(-b_true * x) + c_true + np.random.normal(0, 0.1, size=len(x))

2. 定义模型和目标函数 #

这里我们用带0.5因子的目标函数（和scipy优化器的习惯一致）：

def model(params, x):
    a, b, c = params
    return a * np.exp(-b * x) + c

def objective(params, x, y):
    residuals = y - model(params, x)
    return 0.5 * np.sum(residuals ** 2)  # 带0.5因子的残差平方和

# 定义Jacobian（可选，但能让Hessian逆更准确）
def jacobian(params, x, y):
    a, b, c = params
    n = len(x)
    jac = np.zeros((n, 3))
    jac[:, 0] = np.exp(-b * x)
    jac[:, 1] = -a * x * np.exp(-b * x)
    jac[:, 2] = np.ones(n)
    residuals = y - model(params, x)
    return jac.T @ residuals  # 目标函数的梯度是J^T r

3. 拟合参数 #

initial_guess = [2, 0.2, 0.5]
res = minimize(objective, initial_guess, args=(x, y), jac=jacobian, method='L-BFGS-B')
theta_hat = res.x

4. 计算参数协方差和方差 #

# 计算残差平方和
residuals = y - model(theta_hat, x)
S = np.sum(residuals ** 2)

# 自由度
n = len(y)
p = len(theta_hat)
df = n - p

# 估计噪声方差σ²
sigma_sq = S / df

# 获取Hessian逆（这里因为目标函数带0.5因子，Hessian≈J^T J，所以hess_inv≈(J^T J)^{-1}）
hess_inv = res.hess_inv.todense()  # L-BFGS-B返回的是稀疏矩阵，转成稠密矩阵

# 协方差矩阵
cov_matrix = sigma_sq * hess_inv

# 参数方差（协方差矩阵的对角线）
param_variances = np.diag(cov_matrix)
param_stds = np.sqrt(param_variances)

print("拟合参数：", theta_hat)
print("参数方差：", param_variances)
print("参数标准差：", param_stds)

额外提示 #

如果你用scipy.optimize.curve_fit，它会直接返回pcov（已经处理好的协方差矩阵），本质上就是做了上述的自由度缩放和Hessian修正，你可以用它来验证你的计算结果是否正确。

最后再强调一遍：核心是用自由度(n-p)缩放残差平方和得到σ²，这就是你看到的“测量次数缩放”的本质——不是单纯除以测量次数，而是除以自由度，当参数个数远小于测量次数时，两者近似相等。

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