高斯牛顿法

之前讨论了牛顿法和梯度下降法，前者基于海森矩阵对

设状态 $x\in\mathbb{R}^n$，残差 $r(x)\in\mathbb{R}^m$ 为向量函数，

代价函数为标量函数

$$F(x)=\frac{1}{2}\lVert r(x)\rVert^2 =\frac{1}{2}\sum_{i=1}^{m}r_i(x)^2 .$$

雅可比矩阵定义为

$$ J=\frac{\partial r}{\partial x}, \qquad J_{ij}=\frac{\partial r_i}{\partial x_j}, \qquad J\in\mathbb{R}^{m\times n}.$$

对代价函数展开求和

$$ F=\frac{1}{2}\left(r_1^2+r_2^2+\cdots+r_m^2\right).$$

对单个状态分量求偏导

$x_j$ 通过观测函数影响所有残差，故需对每一项求和：

$$ \frac{\partial F}{\partial x_j} =\frac{1}{2}\sum_{i=1}^{m}\frac{\partial\!\left(r_i^2\right)}{\partial x_j}.$$

对每一项用链式法则

$$ \frac{\partial\!\left(r_i^2\right)}{\partial x_j} =\underbrace{2r_i}_{\text{对外层平方求导}} \cdot \underbrace{\frac{\partial r_i}{\partial x_j}}_{\text{对内层用链式法则}}.$$

代回，系数 $\tfrac12$ 与 $2$ 相消

$$ \frac{\partial F}{\partial x_j} =\frac{1}{2}\sum_{i=1}^{m}2\,r_i\frac{\partial r_i}{\partial x_j} =\sum_{i=1}^{m}r_i\frac{\partial r_i}{\partial x_j}.$$

代入雅可比

$$ \boxed{\; \frac{\partial F}{\partial x_j} =\sum_{i=1}^{m}J_{ij}\,r_i \;}$$

即：第 $j$ 个梯度分量 $=$ 所有残差按“它们对 $x_j$ 的敏感度”加权求和。

写成矩阵形式

$$(J^{\top} r)_j=\sum_{i=1}^{m}(J^{\top})_{ji}r_i=\sum_{i=1}^{m}J_{ij}r_i$$

与第 4 步结果一致，故

$$\nabla F=J^{\top} r$$

也就是说，梯度可以表示为残差对状态的雅可比转置后乘以残差本身，这是由传递法则得到的结论，能够帮助我们得到后面下降的方向以及距离。

对刚才已经得到的F导数再求一次导

$$\frac{\partial^2F}{\partial x_j\partial x_k}=\sum_{i=1}^{m}\frac{\partial}{\partial x_k}\left(r_iJ_{ij}\right)$$

$$\frac{\partial}{\partial x_k}\left(r_iJ_{ij}\right)=\frac{\partial r_i}{\partial x_k}J_{ij}+r_i\frac{\partial J_{ij}}{\partial x_k}$$

代入并整理，得到两项

$$\frac{\partial^2F}{\partial x_j\partial x_k}=\sum_{i=1}^{m}J_{ij}J_{ik}+\sum_{i=1}^{m}r_i\frac{\partial^2r_i}{\partial x_j\partial x_k}$$

认出第一项就是 JᵀJ

$$(J^{\top}J)_{jk}=\sum_{i=1}^{m}J_{ij}J_{ik}$$

所以第二项被丢弃后（高斯牛顿近似：认为残差足够小，所以后面一项对曲率影响可以忽略）

$$\nabla^2F\approx J^{\top}J=H$$

之后直接带入牛顿法就能替代之前的二阶导数，用一阶导数去逼近了。

$$J^\top J\,\Delta x = -J^\top r$$

其中 $H$ 等于 $H=J^\top J$，注意到此时 $J$ 是残差 $r$ 关于自变量 $x$ 的雅可比矩阵（$m\times n$），而代价函数的梯度是 $g=J^\top r$。所以本质上该表达式仍然和之前牛顿法推导的 $H \Delta x = -g$ 是等价的，只不过这里用到了函数在某点处的二阶展开，本质上也就是导数在该点处的线性近似，就能顺理成章地使用“两个J的乘积”实现曲率的近似了。

这个方法的优点很明显，就是通过一阶导数去实现二阶导数的近似，从而避免了去计算曲率的复杂过程。当然同时，这个方法需要求解H的逆，在这个近似中 $H=J^TJ$显然半正定，无法保证总是有解，在鞍面奇异值点处不能保证是真正的最低点。再就是求出来的 $\Delta x$容易过大，导致无法收敛，也同样有可能出现类似梯度下降一样的左右震荡问题。

在生产实践中，我们一般采用LM法进行进行更加精准的下降。