高斯牛顿法 | Gauss-Newton Method

非线性最小二乘问题常用的优化算法有 Newton 法、Gauss-Newton 法 (GN)、Levenberg-Marquardt 算法 (LM) 和 Dog-leg 算法。本文介绍了其常用解法之一的高斯牛顿法（Gauss-Newton Method）的方法和数学推导过程。相比于牛顿法，高斯牛顿法避免了计算海森矩阵，在计算效率和实现难度上都具有优势。

1. 非线性最小二乘问题

非线性最小二乘问题是指最小化目标函数：

其中 \mathbf{x} \in \mathbb{R}^n 是待优化的参数向量，\mathbf{f}(\mathbf{x}) = [f_1(\mathbf{x}), f_2(\mathbf{x}), \cdots, f_m(\mathbf{x})]^{\top} 是残差向量，每个 f_i(\mathbf{x}) 都是非线性函数。这种形式在许多实际问题中都很常见，比如曲线拟合、相机标定等。

2. Gauss-Newton 法

高斯牛顿法的核心思想是通过对残差函数进行一阶泰勒展开，将非线性最小二乘问题转化为线性最小二乘问题。

2.1 一阶泰勒展开

对残差函数 \mathbf{f}(\mathbf{x}) 在当前点 \mathbf{x} 处进行一阶泰勒展开：

其中 \mathbf{J} 是残差函数 \mathbf{f}(\mathbf{x}) 的雅克比矩阵：

2.2 矩阵维度分析

在进行后续推导之前，我们先明确各个矩阵的维度。假设有 n 个参数要优化，m 个观测值（残差），则：

残差向量 \mathbf{f}(x) 的维度为 m \times 1：
\mathbf{f}(x) = \begin{bmatrix} f_1(x) \\ f_2(x) \\ \vdots \\ f_m(x) \end{bmatrix} \in \mathbb{R}^{m \times 1}

雅可比矩阵 \mathbf{J} 的维度为 m \times n：
\mathbf{J}(x) = \begin{bmatrix}\frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \cdots & \frac{\partial f_1}{\partial x_n} \\\frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} & \cdots & \frac{\partial f_2}{\partial x_n} \\\vdots & \vdots & \ddots & \vdots \\\frac{\partial f_m}{\partial x_1} & \frac{\partial f_m}{\partial x_2} & \cdots & \frac{\partial f_m}{\partial x_n}\end{bmatrix} \in \mathbb{R}^{m \times n}

参数增量 \Delta \mathbf{x} 的维度为 n \times 1：
\Delta \mathbf{x} = \begin{bmatrix} \Delta x_1 \\ \Delta x_2 \\ \vdots \\ \Delta x_n \end{bmatrix} \in \mathbb{R}^{n \times 1}

2.3 最小二乘损失函数推导

将一阶泰勒展开代入原目标函数，得到关于增量 \Delta \mathbf{x} 的近似二次型：

这里需要说明交叉项 \Delta \mathbf{x}^{\top} \mathbf{J}^{\top} \mathbf{f} 的推导。根据矩阵维度分析：

• \mathbf{f}^{\top}\mathbf{J}\Delta \mathbf{x} 的计算过程：
  • \mathbf{f}^{\top}: 1 \times m
  • \mathbf{J}: m \times n
  • \Delta \mathbf{x}: n \times 1
  • 最终结果为 1 \times 1 的标量
• \Delta \mathbf{x}^{\top}\mathbf{J}^{\top}\mathbf{f} 的计算过程：
  • \Delta \mathbf{x}^{\top}: 1 \times n
  • \mathbf{J}^{\top}: n \times m
  • \mathbf{f}: m \times 1
  • 最终结果为 1 \times 1 的标量

由于这两个表达式都得到标量，且标量的转置等于其本身，因此：
\mathbf{f}^{\top}\mathbf{J}\Delta \mathbf{x} = \Delta \mathbf{x}^{\top}\mathbf{J}^{\top}\mathbf{f}\tag{5}

2.4 求解增量方程

为了求解最优的 \Delta \mathbf{x}，我们需要：
\Delta \mathbf{x} \equiv \arg \min _{\Delta \mathbf{x}}\left\{L(\Delta \mathbf{x})\right\}\tag{6}

对 L(\Delta \mathbf{x}) 求关于 \Delta \mathbf{x} 的偏导数时，需要用到以下矩阵求导法则：

• 标量对向量的求导：
  • 若 \mathbf{a} 是 n \times 1 维向量，\mathbf{x} 是 n \times 1 维向量，则：
    \frac{\partial (\mathbf{a}^{\top}\mathbf{x})}{\partial \mathbf{x}} = \mathbf{a}
  • 对称的情况，若 \mathbf{x} 是 n \times 1 维向量，\mathbf{a} 是 n \times 1 维向量，则：
    \frac{\partial (\mathbf{x}^{\top}\mathbf{a})}{\partial \mathbf{x}} = \mathbf{a}

  这两个公式是等价的，因为对于标量来说 \mathbf{a}^{\top}\mathbf{x} = \mathbf{x}^{\top}\mathbf{a}。

• 二次型对向量的求导：\frac{\partial (\mathbf{x}^{\top}A\mathbf{x})}{\partial \mathbf{x}} = (A + A^{\top})\mathbf{x}

根据公式 (4)，L(\Delta \mathbf{x}) 包含三项：
L(\Delta \mathbf{x}) = \frac{1}{2} \mathbf{f}^{\top} \mathbf{f}+\Delta \mathbf{x}^{\top} \mathbf{J}^{\top} \mathbf{f}+\frac{1}{2} \Delta \mathbf{x}^{\top} \mathbf{J}^{\top} \mathbf{J} \Delta \mathbf{x}\tag{7}

分别对这三项求导：

1. 第一项：
\frac{1}{2} \mathbf{f}^{\top} \mathbf{f} 中 \mathbf{f} = \mathbf{f}(\mathbf{x}) 是 \mathbf{x} 的函数与 \Delta \mathbf{x} 无关，因此对 \Delta \mathbf{x} 求导为 0

2. 第二项：
\Delta \mathbf{x}^{\top} \mathbf{J}^{\top} \mathbf{f} 是线性项，其中 \mathbf{J}^{\top} \mathbf{f}\in \mathbb{R}^{n \times 1}，根据标量对向量的求导法则：
\frac{\partial (\Delta \mathbf{x}^{\top} \mathbf{J}^{\top} \mathbf{f})}{\partial \Delta \mathbf{x}} = \mathbf{J}^{\top} \mathbf{f}\tag{8}

3. 第三项：
\frac{1}{2} \Delta \mathbf{x}^{\top} \mathbf{J}^{\top} \mathbf{J} \Delta \mathbf{x} 是二次型，令 A = \frac{1}{2}\mathbf{J}^{\top} \mathbf{J}，根据二次型求导法则：
\begin{aligned}\frac{\partial (\frac{1}{2} \Delta \mathbf{x}^{\top} \mathbf{J}^{\top} \mathbf{J} \Delta \mathbf{x})}{\partial \Delta \mathbf{x}} &= \frac{1}{2}[(\mathbf{J}^{\top} \mathbf{J}) + (\mathbf{J}^{\top} \mathbf{J})^{\top}]\Delta \mathbf{x} \\&= \frac{1}{2}[2\mathbf{J}^{\top} \mathbf{J}]\Delta \mathbf{x} \\&= \mathbf{J}^{\top} \mathbf{J}\Delta \mathbf{x}\end{aligned}\tag{9}

这里用到了 \mathbf{J}^{\top} \mathbf{J} 是对称矩阵的性质，即 (\mathbf{J}^{\top} \mathbf{J})^{\top} = \mathbf{J}^{\top} \mathbf{J}。

将三项的导数相加并令其为零：
\begin{aligned}\frac{\partial L(\Delta \mathbf{x})}{\partial \Delta \mathbf{x}} &= 0 + \mathbf{J}^{\top} \mathbf{f} + \mathbf{J}^{\top} \mathbf{J}\Delta \mathbf{x} = 0\end{aligned}\tag{10}

最终得到高斯牛顿法的增量方程，也可以称为高斯牛顿方程（Gauss-Newton equation）或正规方程（Normal equation）：
\underbrace{\left(\mathbf{J}^{\top} \mathbf{J}\right)}_{\mathbf{H}} \Delta \mathbf{x}_{\mathbf{g n}}=\underbrace{-\mathbf{J}^{\top} \mathbf{f}}_{\mathbf{g}}\tag{11}

把左右系数分别定义为 \mathbf{H} 和 \mathbf{g}，则上述方程简写为：
\mathbf{H} \Delta \mathbf{x}_{\mathbf{g n}}=\mathbf{g}\tag{12}

高斯牛顿法的一个重要特点是不需要求解二阶的海森矩阵，而是用 \mathbf{J}^{\top} \mathbf{J} 作为海森矩阵的近似，这使得其计算量比牛顿法更小。此外，由于 \mathbf{J}^{\top} \mathbf{J} 是对称正定矩阵（当 \mathbf{J} 满秩时），这保证了增量方程有唯一解。

需要注意的是，高斯牛顿法也存在一些局限性：

• 当 \mathbf{J} 不满秩时，\mathbf{J}^{\top} \mathbf{J} 不可逆，算法可能失效
• 由于使用了一阶泰勒展开的近似，在远离最优解的区域收敛性可能不好
• 没有步长选择策略，可能会出现震荡或发散

这些问题在后续的 Levenberg-Marquardt 算法中得到了改进。

3. 算法实现与迭代过程

高斯牛顿法的迭代过程可以分为三个主要步骤：

3.1 初始化

• 选择初始参数估计值 \mathbf{x}_0
• 设置迭代终止条件：
  • 最大迭代次数 N_{\max}
  • 参数更新量阈值 \epsilon_{\Delta x}
  • 残差变化阈值 \epsilon_f

3.2 迭代计算

对于第 k 次迭代：

1. 计算雅可比矩阵
\mathbf{J}(\mathbf{x}_k) = \begin{bmatrix}\frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \cdots & \frac{\partial f_1}{\partial x_n} \\\frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} & \cdots & \frac{\partial f_2}{\partial x_n} \\\vdots & \vdots & \ddots & \vdots \\\frac{\partial f_m}{\partial x_1} & \frac{\partial f_m}{\partial x_2} & \cdots & \frac{\partial f_m}{\partial x_n}\end{bmatrix}\tag{13}

2. 求解增量方程
\mathbf{J}(\mathbf{x}_k)^{\top}\mathbf{J}(\mathbf{x}_k)\Delta\mathbf{x}_k = -\mathbf{J}(\mathbf{x}_k)^{\top}\mathbf{f}(\mathbf{x}_k)\tag{14}

可以通过各种方法求解这个方程，得到 \Delta\mathbf{x}_k

3. 更新参数
\mathbf{x}_{k+1} = \mathbf{x}_k + \Delta\mathbf{x}_k\tag{15}

3.3 终止判断

算法在满足以下任一条件时终止：

• 残差变化小于阈值：
  \|\mathbf{f}(\mathbf{x}_{k+1})\| - \|\mathbf{f}(\mathbf{x}_k)\| < \epsilon_f\tag{16}
• 参数更新量小于阈值：
  \|\Delta\mathbf{x}_k\| < \epsilon_{\Delta\mathbf{x}_k}\tag{17}
• 达到最大迭代次数：
  k \geq N_{\max}\tag{18}

参考文献

1. 贺一家, 高翔, 崔华坤.《从零开始手写 VIO》讲座
2. Gauss–Newton algorithm - Wikipedia
3. Nocedal, J., & Wright, S. (2006). Numerical optimization. Springer Science & Business Media.
4. Madsen, K., Nielsen, H. B., & Tingleff, O. (2004). Methods for non-linear least squares problems.