4. 参数估计的最小二乘原理

在非线性曲线拟合中，我们通常已经选定了模型形式，例如指数衰减 y = a·exp(-b·x) + c、Logistic 生长曲线或双指数模型，剩下的核心任务就是确定模型中的未知参数。最小二乘法给出了一个清晰而通用的参数估计准则：让模型预测值与观测值之间的“偏差总量”尽可能小。理解这一准则的目标函数、几何含义与最优性条件，是后续使用高斯-牛顿法、Levenberg-Marquardt 算法进行迭代求解的理论基础。

一、目标函数：残差平方和

设有 n 组观测数据 (x_i, y_i)，i = 1, 2, …, n，选定的非线性模型为 f(x; θ)，其中参数向量 θ = (θ₁, θ₂, …, θ_p)。对第 i 个观测点，残差定义为观测值减去模型预测值：

r_i(θ) = y_i − f(x_i; θ)

注意残差带符号：模型预测偏高时残差为负，预测偏低时残差为正。若直接把残差相加，正负会相互抵消，无法反映整体偏差大小；若取绝对值，函数在零点不可导，不利于优化。因此最小二乘法采用残差平方和作为目标函数：

S(θ) = Σ_{i=1}^{n} r_i(θ)² = Σ_{i=1}^{n} [y_i − f(x_i; θ)]²

参数估计问题于是转化为一个无约束优化问题：求 θ̂ = argmin S(θ)。平方运算带来两个好处：一是消除符号，使正负偏差同等受罚；二是对大残差施加更重的惩罚，因为偏差翻倍，平方项变为四倍。这也意味着最小二乘对离群点较为敏感——一个远离主体的异常点可能显著拉动拟合曲线，实际分析中常需先检查数据质量或改用稳健损失函数。

若各观测点的测量精度不同，可引入权重 w_i > 0，得到加权最小二乘目标：

S(θ) = Σ_{i=1}^{n} w_i · [y_i − f(x_i; θ)]²

权重通常取测量误差方差的倒数，精度高的点权重大。当 w_i = 1 时退化为普通最小二乘。下文为叙述简洁，以普通最小二乘为主，加权情形只需在相应公式中把 r_i² 换成 w_i·r_i²。

二、几何解释：在观测空间中的距离最小化

把 n 个观测值排成一个向量 y = (y₁, y₂, …, y_n)ᵀ，它对应 n 维观测空间中的一个点。当参数 θ 在 p 维参数空间中变化时，模型在所有观测点上的预测值构成向量 f(θ) = (f(x₁; θ), …, f(x_n; θ))ᵀ。随着 θ 变化，f(θ) 在观测空间中扫出一张 p 维的曲面（称为模型流形）。

残差向量 r(θ) = y − f(θ) 正是从模型流形上的点指向观测点 y 的向量，其欧氏长度的平方就是 S(θ)。因此，最小二乘估计的几何意义是：在模型流形上寻找距离观测点最近的那个点，对应的参数就是最优参数。若模型是线性的，流形退化为一个平面（或超平面），最近点就是 y 在该平面上的正交投影，残差向量与平面垂直——这正是线性最小二乘“正规方程”的几何来源。非线性情形下流形弯曲，最近点处残差向量应与流形在该点的切空间正交，这一直觉将直接导出下面的最优性条件。

三、最优性条件：梯度为零与正交性

目标函数 S(θ) 是 p 个参数的标量函数。若 θ̂ 是局部极小点且 S 在该点可微，则一阶必要条件为梯度为零：

∂S/∂θ_j = −2 · Σ_{i=1}^{n} r_i(θ) · ∂f(x_i; θ)/∂θ_j = 0,   j = 1, …, p

写成矩阵形式，记残差向量 r(θ)、雅可比矩阵 J(θ)，其第 i 行第 j 列为 ∂f(x_i; θ)/∂θ_j，则梯度条件等价于：

J(θ̂)ᵀ · r(θ̂) = 0

这正是几何解释的代数表达：残差向量与雅可比矩阵的每一列都正交，而雅可比列张成的正是模型流形在 f(θ̂) 处的切空间。因此“残差垂直于切空间”与“Jᵀr = 0”是同一件事的两种说法。对线性模型，J 为常数矩阵，该式直接化为正规方程 JᵀJ·θ = Jᵀy，可一步求解；对非线性模型，J 依赖 θ，方程非线性，只能迭代逼近。

需要强调，Jᵀr = 0 只是局部极小的必要条件，满足它的点可能是局部极大或鞍点。要确认极小，还需二阶条件：海森矩阵近似正定。对最小二乘问题，海森矩阵为

H ≈ 2 · JᵀJ

这里略去了含二阶导数与残差乘积的项（高斯-牛顿近似）。只要 J 列满秩，JᵀJ 即正定，局部极小成立。这也解释了为什么参数可辨识性如此重要：若 J 列相关（例如模型存在参数冗余或数据覆盖的 x 范围太窄），JᵀJ 接近奇异，优化会沿某些方向“平坦”，参数估计不稳定、方差巨大。

四、一个可手算的示例

考虑模型 f(x; a) = a·x²，只有单个参数 a，数据为 (1, 2.1)、(2, 7.8)、(3, 18.5)。残差 r_i = y_i − a·x_i²，目标函数：

S(a) = (2.1 − a)² + (7.8 − 4a)² + (18.5 − 9a)²

对 a 求导并令其为零：

dS/da = −2(2.1 − a) − 8(7.8 − 4a) − 18(18.5 − 9a) = 0
整理得：2·(1²+2²+3²)·a = 2·(1·2.1 + 2·7.8 + 3·18.5)
即 28a = 155.4，a ≈ 5.55

该模型对参数是线性的，所以一步得到闭式解。若把模型改为 f(x; a, b) = a·exp(b·x)，导数中会出现 a·x·exp(bx)，方程无法直接解出，必须借助迭代。此时可以这样验证迭代结果是否收敛到正确的最优解：检查 Jᵀr 的各分量是否接近零（例如绝对值小于 1e-6），并观察连续两次迭代的 S 值下降是否已小于设定容差。若 S 不再下降而 Jᵀr 仍明显非零，往往说明迭代停滞或陷入局部极小，需要更换初值或检查模型设定。

五、实践中的几个要点

• 残差定义要统一。本文采用 r = y − f，则梯度为 −2Jᵀr；若采用 f − y，符号相反但平方和不变，最优解一致。代码实现时务必与推导保持一致，避免符号错误导致迭代发散。
• 参数尺度影响收敛。若 a 的量级为 10³ 而 b 为 10⁻³，雅可比列范数差异巨大，JᵀJ 条件数很差。可先对参数做归一化，或直接使用 Levenberg-Marquardt 方法中的阻尼项改善。
• 初值决定能否找到全局最优。非线性最小二乘通常只保证局部最优。对指数类模型，可先对数据取对数做线性回归获得初值；对多峰模型，可网格扫描初值后取 S 最小者。
• 离群点会放大影响。由于平方惩罚，单个大残差点的贡献可能主导 S。拟合后应绘制残差图，检查是否存在系统性弯曲（模型形式不当）或孤立大残差（异常点）。

综上，最小二乘原理可以概括为三层：目标函数把“拟合好坏”量化为残差平方和；几何上它是在模型流形上寻找距观测点最近的点；最优性条件则要求残差与切空间正交，即 Jᵀr = 0，并辅以 JᵀJ 正定的二阶确认。掌握这三层含义，后续学习各类非线性优化算法时，就能清楚每一步迭代究竟在朝哪个方向、依据什么准则推进。