2. 非线性最小二乘基础

在参数拟合问题中，我们通常有一组观测数据 \((t_i, y_i)\)，以及一个由待定参数向量 \(\boldsymbol{\theta}\) 决定的模型函数 \(f(t;\boldsymbol{\theta})\)。当模型对参数是非线性的（例如指数衰减 \(y=Ae^{-kt}\)、双指数、Logistic 曲线、多峰高斯叠加等），就不能像线性最小二乘那样用一次矩阵求解直接得到解析解，而必须借助迭代优化。本节要建立的，正是这条迭代求解链路的数学骨架：残差如何定义、目标函数如何构建、梯度与雅可比矩阵如何推导、正规方程如何形成，以及高斯-牛顿法如何利用这些量完成一步参数更新。

残差与目标函数

先明确符号约定：设共有 \(m\) 个观测点，参数向量 \(\boldsymbol{\theta}\in\mathbb{R}^{n}\)，第 \(i\) 个点的模型预测值为 \(f_i(\boldsymbol{\theta})=f(t_i;\boldsymbol{\theta})\)。定义残差为“观测减预测”：

\[ r_i(\boldsymbol{\theta}) = y_i - f_i(\boldsymbol{\theta}),\qquad i=1,\dots,m \]

把全部残差堆成列向量 \(\mathbf{r}(\boldsymbol{\theta})\in\mathbb{R}^{m}\)。最小二乘的目标是让残差平方和最小，即构建目标函数：

\[ S(\boldsymbol{\theta}) = \sum_{i=1}^{m} r_i^2(\boldsymbol{\theta}) = \mathbf{r}^{\mathsf{T}}\mathbf{r} \]

这里采用“观测减预测”的约定，因此残差符号与后续梯度方向必须保持一致；若改用“预测减观测”，则梯度、雅可比和更新式中的符号会整体翻转，但最终迭代结果不变。本文与图示统一采用“观测减预测”。

需要注意，非线性最小二乘的 \(S(\boldsymbol{\theta})\) 一般不是参数的二次型，因此可能存在多个局部极小值，这正是本课程后续要专门处理局部极值的原因。此外，各观测点的噪声水平可能不同，工程上常引入权重 \(w_i\)（如 \(w_i=1/\sigma_i^2\)），得到加权目标 \(S_w=\sum w_i r_i^2\)，其推导结构与无权情形完全一致，只是雅可比矩阵的每一行要乘以 \(\sqrt{w_i}\)。

梯度与雅可比矩阵

要对 \(S\) 做迭代下降，首先需要它的一阶导数。对 \(S=\sum r_i^2\) 关于第 \(j\) 个参数求偏导，用链式法则：

\[ \frac{\partial S}{\partial \theta_j} = 2\sum_{i=1}^{m} r_i\frac{\partial r_i}{\partial \theta_j} = -2\sum_{i=1}^{m} r_i\frac{\partial f_i}{\partial \theta_j} \]

其中负号来自 \(r_i=y_i-f_i\)。定义雅可比矩阵 \(\mathbf{J}\in\mathbb{R}^{m\times n}\)，其元素为模型对参数的偏导：

\[ J_{ij} = \frac{\partial f_i}{\partial \theta_j} \]

于是梯度可紧凑写成 \(\nabla S = -2\mathbf{J}^{\mathsf{T}}\mathbf{r}\)。雅可比矩阵是整条链路的核心：它把“参数空间的一小步”映射为“预测值的变化”，从而决定残差如何变化。对常见模型，雅可比可以解析写出。以单指数衰减 \(f(t;A,k)=Ae^{-kt}\) 为例，参数 \(\boldsymbol{\theta}=(A,k)\)，则

∂f/∂A = e^(-k t)
∂f/∂k = -A t e^(-k t)

若模型过于复杂难以手写偏导，也可用有限差分近似 \(J_{ij}\approx[f_i(\theta_j+h)-f_i(\theta_j)]/h\)，但会引入截断误差，且 \(h\) 的选取需与参数尺度匹配。

正规方程与高斯-牛顿法原理

有了梯度，最朴素的做法是梯度下降：\(\boldsymbol{\theta}_{k+1}=\boldsymbol{\theta}_k-\alpha\nabla S\)。但梯度下降只利用一阶信息，步长 \(\alpha\) 难选，在参数尺度差异大或谷底狭长时收敛很慢。高斯-牛顿法改用二阶近似，思路是把非线性模型在当前参数处做一阶泰勒展开：

\[ f_i(\boldsymbol{\theta}+\Delta\boldsymbol{\theta}) \approx f_i(\boldsymbol{\theta}) + \sum_j J_{ij}\Delta\theta_j \]

代入残差，得到线性化残差 \(\mathbf{r}(\boldsymbol{\theta}+\Delta\boldsymbol{\theta})\approx \mathbf{r}-\mathbf{J}\Delta\boldsymbol{\theta}\)。于是目标函数近似为

\[ S(\boldsymbol{\theta}+\Delta\boldsymbol{\theta}) \approx \|\mathbf{r}-\mathbf{J}\Delta\boldsymbol{\theta}\|^2 \]

这是一个关于 \(\Delta\boldsymbol{\theta}\) 的二次型，对其求导并令导数为零，即得正规方程：

\[ (\mathbf{J}^{\mathsf{T}}\mathbf{J})\,\Delta\boldsymbol{\theta} = \mathbf{J}^{\mathsf{T}}\mathbf{r} \]

解出 \(\Delta\boldsymbol{\theta}\) 后更新 \(\boldsymbol{\theta}\leftarrow\boldsymbol{\theta}+\Delta\boldsymbol{\theta}\)，重复“计算残差与雅可比—解正规方程—更新参数”的循环，直到满足收敛判据。矩阵 \(\mathbf{J}^{\mathsf{T}}\mathbf{J}\) 是 \(S\) 的 Gauss-Newton 近似 Hessian（忽略了残差二阶导项 \(\sum r_i\partial^2 f_i/\partial\theta_j\partial\theta_k\)），当残差较小时该近似很好，收敛接近二次速度；当残差较大或模型强非线性时，近似变差，可能步长过大甚至发散。

工程实现中通常不显式求逆，而是用 Cholesky 分解或 QR 分解求解正规方程，数值更稳定。若 \(\mathbf{J}^{\mathsf{T}}\mathbf{J}\) 接近奇异（参数不可辨识或高度相关），可加阻尼项 \((\mathbf{J}^{\mathsf{T}}\mathbf{J}+\lambda\mathbf{I})\Delta\boldsymbol{\theta}=\mathbf{J}^{\mathsf{T}}\mathbf{r}\)，这正是 Levenberg-Marquardt 方法的雏形：\(\lambda\) 大时退化为小步梯度下降，\(\lambda\) 小时接近高斯-牛顿。

一个可手算的示例

设模型 \(f(t;\theta)=\theta t\)（单参数线性模型，用于验证链路），观测为 \((t,y)=(1,2.1),(2,3.9),(3,6.2)\)，初值 \(\theta_0=1.5\)。此时 \(J_{i1}=t_i\)，\(\mathbf{J}=[1,2,3]^{\mathsf{T}}\)。预测为 \(1.5,3.0,4.5\)，残差 \(\mathbf{r}=[0.6,0.9,1.7]^{\mathsf{T}}\)。计算 \(\mathbf{J}^{\mathsf{T}}\mathbf{J}=1+4+9=14\)，\(\mathbf{J}^{\mathsf{T}}\mathbf{r}=0.6+1.8+5.1=7.5\)，故 \(\Delta\theta=7.5/14\approx0.536\)，更新后 \(\theta_1\approx2.036\)。再迭代一次残差显著减小，逐步逼近最小二乘解 \(\hat\theta=\sum t_iy_i/\sum t_i^2=35.1/14\approx2.507\)。该例虽为线性，但完整走通了“残差—雅可比—正规方程—更新”的每一步，便于对照检查符号与矩阵维度是否正确。