非线性最小二乘：牛顿法、高斯—牛顿与 LM

1. 问题形式

非线性最小二乘问题写作：

\[ \min_x F(x),\qquad F(x)=\frac12\|r(x)\|^2 =\frac12\sum_{i=1}^m r_i(x)^2, \]

其中 \(x\in\mathbb R^n\) 是待估参数，\(r(x)\in\mathbb R^m\) 是残差向量。与线性最小二乘不同，\(r(x)\) 对 \(x\) 是非线性的，不能一次整理成固定的 \(Ax=b\) 后直接求解。

例如，用

\[ y=ae^{bt} \]

拟合数据 \((t_i,y_i)\)。令 \(x=[a,b]^T\)，第 \(i\) 个残差为

\[ r_i(a,b)=ae^{bt_i}-y_i, \]

目标就是最小化这些残差的平方和。

2. 从牛顿法开始

对一般优化问题

\[ \min_x f(x), \]

在当前点 \(x_k\) 附近，用二阶泰勒展开建立局部模型：

\[ m(\Delta x) =f(x_k)+\nabla f(x_k)^T\Delta x +\frac12\Delta x^T H(x_k)\Delta x, \]

其中 \(H=\nabla^2f\) 是 Hessian。此时 \(x_k\) 已固定，未知量是“下一步往哪里走、走多远”的位移 \(\Delta x\)。

对 \(\Delta x\) 求梯度并令其为零：

\[ \nabla_{\Delta x}m =\nabla f+H\Delta x=0, \]

得到牛顿方程

\[ H\Delta x=-\nabla f, \]

再更新

\[ x_{k+1}=x_k+\Delta x. \]

实践中应解线性方程 \(H\Delta x=-\nabla f\)，而不是显式计算 \(H^{-1}\)。

为什么一阶展开不够

一阶模型

\[ f(x+\Delta x)\approx f(x)+\nabla f(x)^T\Delta x \]

是一个平面，只能给出下降方向，不能独立确定有限步长。例如 \(f(x)=x^2\) 在 \(x=10\) 处的一阶模型为

\[ 100+20\Delta x. \]

若直接最小化这个线性模型，会得到 \(\Delta x\to-\infty\)。因此梯度下降需要另选学习率，而牛顿法用 Hessian 描述局部曲率，从二次模型中同时得到方向和步长。

3. 高斯—牛顿推导

高斯—牛顿专门利用最小二乘的平方结构。它不对目标函数 \(F\) 做二阶展开，而只对残差做一阶展开：

\[ r(x_k+\Delta x) \approx r(x_k)+J(x_k)\Delta x, \]

其中

\[ J=\frac{\partial r}{\partial x} \]

是残差关于参数的 Jacobian。简记当前的 \(r(x_k)\) 和 \(J(x_k)\) 为 \(r\)、\(J\)，局部问题变为

\[ \min_{\Delta x}\frac12\|r+J\Delta x\|^2. \]

虽然残差只做了一阶展开，但代入平方损失后会自然产生二次项：

\[ \frac12(r+J\Delta x)^T(r+J\Delta x) =\frac12r^Tr +\Delta x^TJ^Tr +\frac12\Delta x^TJ^TJ\Delta x. \]

对 \(\Delta x\) 求梯度并令其为零：

\[ J^T(r+J\Delta x)=0, \]

因此高斯—牛顿方程为

\[ \boxed{J^TJ\Delta x=-J^Tr}. \]

求得增量后更新：

\[ \boxed{x_{k+1}=x_k+\Delta x}. \]

每轮重新计算 \(r\) 和 \(J\)，直到步长、梯度或目标函数变化足够小。

整个过程可以概括为：

\[ x_k \longrightarrow (r,J) \longrightarrow \min_{\Delta x}\|r+J\Delta x\|^2 \longrightarrow \Delta x \longrightarrow x_{k+1}. \]

本质上，每次迭代都把一个非线性问题局部变成线性最小二乘问题。

4. 高斯—牛顿与牛顿法的关系

对

\[ F(x)=\frac12r(x)^Tr(x), \]

其梯度为

\[ \boxed{\nabla F=J^Tr}, \]

精确 Hessian 为

\[ \boxed{ \nabla^2F =J^TJ+\sum_{i=1}^m r_i\nabla^2r_i }. \]

将其代入牛顿方程，本应得到

\[ \left(J^TJ+\sum_i r_i\nabla^2r_i\right)\Delta x=-J^Tr. \]

高斯—牛顿忽略残差的二阶导数项，使用近似

\[ \boxed{\nabla^2F\approx J^TJ}, \]

于是得到

\[ J^TJ\Delta x=-J^Tr. \]

因此二者不是同一种方法，但关系很近：

\[ \boxed{ \text{高斯—牛顿} =\text{针对非线性最小二乘的近似牛顿法} }. \]

当解附近的残差较小，即 \(r_i\approx0\) 时，被忽略的

\[ \sum_i r_i\nabla^2r_i \]

通常也较小，因此 \(J^TJ\) 对真实 Hessian 的近似较好。即使残差不小，只要各 \(r_i\) 的非线性较弱，这个近似也可能有效。

5. 示例：指数曲线拟合

模型与残差为

\[ y=ae^{bt},\qquad r_i(a,b)=ae^{bt_i}-y_i. \]

两个偏导数是

\[ \frac{\partial r_i}{\partial a}=e^{bt_i}, \qquad \frac{\partial r_i}{\partial b}=at_i e^{bt_i}, \]

所以

\[ J= \begin{bmatrix} e^{bt_1} & at_1e^{bt_1}\\ e^{bt_2} & at_2e^{bt_2}\\ \vdots & \vdots\\ e^{bt_m} & at_me^{bt_m} \end{bmatrix}. \]

在当前参数 \((a,b)\) 处计算 \(r\) 和 \(J\)，求解

\[ J^TJ \begin{bmatrix} \Delta a\\ \Delta b \end{bmatrix} =-J^Tr, \]

然后更新

\[ a\leftarrow a+\Delta a, \qquad b\leftarrow b+\Delta b. \]

一个用于理解算法的 NumPy 实现如下：

import numpy as np

t = np.array([0.0, 1.0, 2.0, 3.0])
y = np.array([2.0, 3.1, 5.5, 8.8])
x = np.array([1.0, 0.5])  # [a, b]

for _ in range(20):
    a, b = x
    r = a * np.exp(b * t) - y
    J = np.column_stack([
        np.exp(b * t),
        a * t * np.exp(b * t),
    ])

    # 直接解 J @ dx ≈ -r，避免构造逆矩阵。
    dx = np.linalg.lstsq(J, -r, rcond=None)[0]
    x += dx

    if np.linalg.norm(dx) < 1e-8:
        break

数值计算时，通常不应显式使用

np.linalg.inv(J.T @ J) @ (-J.T @ r)

因为构造 \(J^TJ\) 会放大条件数，而显式求逆也没有必要。直接用 QR 或 SVD 解 \(J\Delta x\approx-r\) 更稳定。

6. LM：给高斯—牛顿加入阻尼

Levenberg—Marquardt（LM）的典型形式为

\[ \boxed{(J^TJ+\lambda I)\Delta x=-J^Tr}, \]

其中 \(\lambda I\) 是阻尼项。

当 \(\lambda\) 很小时，LM 接近高斯—牛顿；当 \(\lambda\) 很大时，

\[ \Delta x\approx-\frac1\lambda J^Tr, \]

而 \(J^Tr\) 正是目标函数的梯度，所以此时更接近小步长梯度下降。

LM 的核心作用是：局部模型不可靠或 \(J^TJ\) 病态时走得更保守；局部模型可靠时逐渐采用快速的高斯—牛顿步。实际算法会根据试探步带来的真实下降与模型预测下降是否一致，动态调整阻尼，而不是固定使用一个 \(\lambda\)。

7. 初始值与收敛

高斯—牛顿及 LM 都是局部优化方法。它们依赖

\[ r(x+\Delta x)\approx r(x)+J(x)\Delta x \]

在当前点附近足够准确。因此：

• 初值较好时，通常能快速精修参数；
• 初值较差时，可能走向错误的局部极小值、震荡或发散；
• 非线性越强，对初值往往越敏感；
• LM 和信赖域方法能扩大可收敛区域，但不能消除初值问题。

常见工程流程是：

\[ \boxed{ \text{粗初始化或全局搜索} \longrightarrow \text{高斯—牛顿、LM 或信赖域方法精修} }. \]

8. 值得保留的问题与答案

Q1：只对残差做一阶展开，为什么能确定步长？

因为残差被放进平方损失后，一阶模型 \(r+J\Delta x\) 的平方会产生 \(\Delta x^TJ^TJ\Delta x\)，自然形成一个关于 \(\Delta x\) 的二次模型。高斯—牛顿正是利用了这个特殊结构。

Q2：“对 \(\Delta x\) 求导”是什么意思？

当前点 \(x_k\) 已经固定，\(f(x_k)\)、梯度、Hessian、\(r(x_k)\) 和 \(J(x_k)\) 都视为常量。把下一步位移 \(\Delta x\) 当成局部模型中的未知变量，对它求梯度并令梯度为零，就是在寻找令局部模型最小的那一步。

Q3：LM 为什么通常比裸高斯—牛顿稳？

阻尼项 \(\lambda I\) 能改善线性系统的可解性，并限制不可靠的大步。LM 会根据局部模型的预测效果调节 \(\lambda\)，在较稳的梯度型步伐和较快的高斯—牛顿步之间过渡。

Q4：实际计算该使用什么 API 和算法？

通用非线性最小二乘可使用 scipy.optimize.least_squares。只需提供残差函数，也可以额外提供 Jacobian：

import numpy as np
from scipy.optimize import least_squares

t = np.array([0.0, 1.0, 2.0, 3.0])
y = np.array([2.0, 3.1, 5.5, 8.8])

def residual(x):
    a, b = x
    return a * np.exp(b * t) - y

def jac(x):
    a, b = x
    return np.column_stack([
        np.exp(b * t),
        a * t * np.exp(b * t),
    ])

result = least_squares(
    residual,
    x0=[1.0, 0.5],
    jac=jac,
    method="trf",
)

普通曲线拟合还可使用更高层的 scipy.optimize.curve_fit。NumPy 提供线性代数工具，但不提供完整的通用非线性最小二乘优化器；scikit-learn 的重点则是现成的机器学习估计器，而不是任意残差函数的非线性最小二乘。

least_squares 提供三种主要算法：

信赖域方法的共同思想是：只在当前点附近的有限区域内信任局部模型。若一步的实际效果符合预测，就扩大信赖域；若效果较差，就缩小信赖域。

一个够用的选择原则是：有边界或没有明确偏好时用 trf；无边界、小中型且初值较好时可试 lm；dogbox 可作为小型带边界问题的备选。

Q5：好初值是不是必要条件？

严格说不是每次都“必须”，但它显著影响能否收敛以及收敛到哪个解。高斯—牛顿尤其适合“已经大致知道答案在哪里，再把参数估准”，不适合单独承担全局搜索。

9. 核心记忆

\[ \boxed{ \text{非线性残差} \xrightarrow{\text{一阶线性化}} r+J\Delta x \xrightarrow{\text{平方损失}} \text{二次局部模型} \xrightarrow{\text{最小化}} J^TJ\Delta x=-J^Tr } \]

牛顿法直接用二阶模型；高斯—牛顿借助最小二乘结构，用 \(J^TJ\) 近似 Hessian；LM 再通过阻尼控制高斯—牛顿步的激进程度。