Scipy curve_fit参数正确但协方差矩阵为无穷大问题排查

问题描述 #

测试scipy.optimize.curve_fit()拟合二次函数y=3x²+2的构造数据时，无论初始猜测如何，拟合参数都能得到预期的[3, 0, 2]，但当初始猜测的二次项系数a不接近0（尤其是大于1时），输出的协方差矩阵全为无穷大。

复现代码：

# python version: 3.9.7
# using a venv
# numpy version: 1.23.2
# scipy version: 1.9.0
import numpy as np
from scipy.optimize import curve_fit

# 构造二次函数y = 3*x**2 + 2的数据
x = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10], dtype=np.float64)
y = np.array([2, 5, 14, 29, 50, 77, 110, 149, 194, 245, 302], dtype=np.float64)

# 定义二次函数
def func(x, a, b, c):
    return a * x**2 + b * x + c

# 成功案例初始猜测
success = [0, 0, 0]
# 失败案例初始猜测
failure = [4, 0, 0]

popt, pcov = curve_fit(func, x, y, p0=failure)  # 切换p0可复现两种情况

print(popt)  # 预期输出 [3. 0. 2.]
print(pcov)  # 失败案例中输出全无穷大的协方差矩阵

原因分析 #

curve_fit默认使用Levenberg-Marquardt（lm）算法，计算协方差矩阵时需要对雅可比矩阵求逆。当前场景下的核心问题是：

• x范围为0-10，x²范围达到0-100，导致二次项a*x²的量级远大于一次项b*x（0-10）和常数项c（2）
• 当初始猜测的a值较大时，参数对应项的量级差异悬殊，雅可比矩阵出现病态（奇异），行列式接近0，无法正常求逆，最终协方差矩阵被填充为无穷大

虽然拟合参数能收敛到正确值，但数值稳定性问题导致协方差矩阵计算失败。

解决方法 #

1. 切换到更鲁棒的算法 #

trf或dogbox算法对参数缩放的敏感性更低，直接替换方法即可解决协方差矩阵异常问题：

popt, pcov = curve_fit(func, x, y, p0=failure, method='trf')

此时即使初始猜测为[4,0,0]，也能输出正常的协方差矩阵。

2. 对输入数据做缩放处理 #

将x缩放到相近量级的区间（比如0-1），让三个项的量级匹配，消除数值差异带来的稳定性问题：

x_scaled = x / 10  # 缩放x到0-1区间
popt_scaled, pcov_scaled = curve_fit(func, x_scaled, y, p0=failure)
# 还原参数：原二次项系数a = popt_scaled[0] / (10**2)
# 其他参数b = popt_scaled[1] / 10，c = popt_scaled[2]

3. 指定绝对误差参数 #

如果已知数据的误差范围，设置absolute_sigma=True并传入误差数组，让协方差计算基于绝对误差而非残差估计：

# 假设每个y值的绝对误差为1
popt, pcov = curve_fit(func, x, y, p0=failure, absolute_sigma=True, sigma=np.ones_like(y))

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