非线性最小二乘
1. 前言
先考虑一个最简单的最小二乘问题： \[ \underset{x}{min} \frac{1}{2} \|f(\boldsymbol{x}) \|_{2}^{2} \] 这里自变量 \(\boldsymbol{x} \in \mathbb{R}^n\) ， \(f\) 是任意一个非线性函数，假设他有 \(m\) 维： \(f(\boldsymbol{x}) \in \mathbb{R}^m\) 。
考虑 \(f\) 的形式基本有两种解决方法：解析法和迭代法。
1.1 解析法
当 \(f\)
的数学形式比较简单的时候，或许问题可以直接用解析形式来求解。也就是让目标函数的导数为
0 ，然后求解 \(\boldsymbol{x}\)
的最优值： \[
\frac{\partial f}{\partial \boldsymbol{x}}=0 \nonumber
\] 解上述方程即可得到目标函数导数为 0
的极值，它可能是极大值、极小值或者鞍点处的值，依次比较即可得到最小值。
显然这种方法较为局限。
1.2 迭代法
对于不易直接求导的目标函数，我们可以用迭代的方法来求极小值。根据给定的初值，不断地更新当前的优化变量，使得目标函数朝着下降的方向变化。具体步骤列写如下：
- 给定初值 \(\boldsymbol{x}_0\)
；
- 对于第 \(k\) 次迭代，寻找一个增量
\(\Delta \boldsymbol{x}_k\)，使得 \(\|f(\boldsymbol{x}_k + \Delta
\boldsymbol{x}_k)\|_{2}^{2}\) 达到极小值；
- 如果 \(\Delta \boldsymbol{x}_k\)
足够小，则停止；
- 否则，令 \(\boldsymbol{x}_{k+1}=\boldsymbol{x}_k + \Delta \boldsymbol{x}_k\)，返回2。
这样问题就转化为如何确定增量 \(\Delta \boldsymbol{x}_k\)，使上述过程收敛，通常有下面的几种方法。
2. 一阶和二阶梯度法（最速下降法和牛顿法）
求解增量最直观的方法就是将 目标函数 在 \(\boldsymbol{x}\) 进行泰勒展开： \[ \|f(\boldsymbol{x}+\Delta\boldsymbol{x}) \|_{2}^{2}=\|f(\boldsymbol{x})\|_{2}^{2}+\boldsymbol{J}(\boldsymbol{x})\Delta \boldsymbol{x}+\frac{1}{2}\Delta \boldsymbol{x }^T \boldsymbol{H} \Delta \boldsymbol{x}. \]
这里 \(\boldsymbol{J}\) 是 \(\|f(\boldsymbol{x}) \|_{2}^{2}\) 关于 \(x\) 的导数（雅可比矩阵），\(\boldsymbol{H}\) 则是二阶导数（黑塞矩阵）。我们可以选择保留泰勒展开的一阶或二阶项，对应的求解方法则为一阶梯度或二阶梯度法。
如果保留一阶梯度，那么增量的方向为： \[
\Delta \boldsymbol{x}^{*}=-\boldsymbol{J}^T(\boldsymbol{x})
\] 这里的转置是为了满足矩阵乘法的维度匹配。
它的直观意义非常简单，只要我们沿着反向梯度方向前进即可。当然，我们还需要该方向上取一个步长
λ，求得最快的下降方式。这种方法被称为最速下降法。
如果保留二阶梯度，那么增量的解为： \[\label{二阶梯度法} \boldsymbol{H} \Delta \boldsymbol{x}^{*}= - \boldsymbol{J}^T \] 该方法称又为牛顿法。
2.1 小结
我们看到，一阶和二阶梯度法都十分直观，只要把函数在迭代点附近进行泰勒展开，并针对更新量作最小化即可。由于泰勒展开之后函数变成了多项式，所以求解增量时只需解线性方程即可，避免了直接求导函数为零这样的非线性方程的困难。
不过，这两种方法也存在它们自身的问题。最速下降法过于贪心，容易走出锯齿路线，反而增加了迭代次数。而牛顿法则需要计算目标函数的 \(\boldsymbol{H}\) 矩阵，这在问题规模较大时非常困难，我们通常倾向于避免 \(\boldsymbol{H}\) 的计算。所以，接下来我们详细地介绍两类更加实用的方法： Gauss-Newton 和 Levenberg-Marquadt 。
3. Gauss-Newton
Gauss Newton 是最优化算法里面最简单的方法之一。它的思想是将 \(f(\boldsymbol{x})\)
进行一阶的泰勒展开：
\[\label{局部近似}
f(\boldsymbol{x}+\Delta \boldsymbol{x}) \approx
f(\boldsymbol{x})+\boldsymbol{J}(\boldsymbol{x})\Delta\boldsymbol{x}
\] 这里 \(\boldsymbol{J}(\boldsymbol{x})\) 是 \(f(\boldsymbol{x})\) 关于 \(\boldsymbol{x}\) 的导数，实际上是一个 \(m\times m\)
的矩阵，也是一个雅可比矩阵。
这里是对 \(f(\boldsymbol{x})\) 在 \(\boldsymbol{x}\) 处进行泰勒展开，而不是对目标函数进行泰勒展开，需要和 梯度法 中进行区分对比！
根据前面的思想，我们当前的目标是为了寻找下降矢量 \(\Delta\boldsymbol{x}\)，使 \(\|f(\boldsymbol{x}+\Delta\boldsymbol{x})\|_{2}^{2}\)
达到最小。即
\[
\Delta \boldsymbol{x}^{*}=arg\ \underset{\Delta \boldsymbol{x}}{min}
\frac{1}{2}\|f(\boldsymbol{x})+\boldsymbol{J}(\boldsymbol{x})\Delta\boldsymbol{x}\|_{2}^{2}
\] 将目标函数展开化简得： \[
\begin{equation*}
\begin{aligned}
\frac{1}{2}\|f(\boldsymbol{x})+\boldsymbol{J}(\boldsymbol{x})\Delta\boldsymbol{x}\|_{2}^{2}&=\frac{1}{2}(f(\boldsymbol{x})+\boldsymbol{J}(\boldsymbol{x})\Delta\boldsymbol{x})^T(f(\boldsymbol{x})+\boldsymbol{J}(\boldsymbol{x})\Delta\boldsymbol{x})
\\
&=\frac{1}{2}(\|f(\boldsymbol{x})\|_{2}^{2} +
2f(\boldsymbol{x})^T\boldsymbol{J}(\boldsymbol{x}+\Delta\boldsymbol{x}^T
\boldsymbol{J}(\boldsymbol{x})^T\boldsymbol{J}(x)^T\Delta
\boldsymbol{x}))
\end{aligned}
\end{equation*}
\]
求上述目标函数对 \(\Delta \boldsymbol{x}\) 的导数，并令其为 0： \[\nonumber 2\boldsymbol{J}(\boldsymbol{x})^T \boldsymbol{J}(\boldsymbol{x})\Delta \boldsymbol{x}+2\boldsymbol{J}(\boldsymbol{x})^T f(\boldsymbol{x})=0 \]
可以得到： \[ \boldsymbol{J}(\boldsymbol{x})^T \boldsymbol{J}(\boldsymbol{x})\Delta \boldsymbol{x}=-\boldsymbol{J}(\boldsymbol{x})^T f(\boldsymbol{x}) \]
注意，我们要求解的变量是 \(\Delta
\boldsymbol{x}\)，因此这是一个线性方程组，我们称它为增量方程，也可以称为高斯牛顿方程
(Gauss Newton equations) 或者正规方程 (Normal
equations)。我们把左边的系数定义为 \(\boldsymbol{H}\)，右边定义为 \(\boldsymbol{g}\)，那么上式变为：
\[
\boldsymbol{H}\Delta \boldsymbol{x} = \boldsymbol{g}
\]
Gauss-Newton 的步骤可以写为：
- 给定初值 \(\boldsymbol{x}_0\)；
- 对于第 \(k\)
次迭代，求出当前的雅可比矩阵 \(\boldsymbol{J}(\boldsymbol{x}_k)\)
和误差\(f(\boldsymbol{x}_k)\)；
- 求解增量方程： \(\boldsymbol{H}\Delta
\boldsymbol{x}_k = \boldsymbol{g}\) ；
- 如果 \(\Delta \boldsymbol{x}_k\) 足够小，则停止；否则，令 \(\boldsymbol{x}_{k+1}=\boldsymbol{x}_k+\Delta \boldsymbol{x}_k\) ，返回 2。
求解增量方程是整个优化问题的核心所在。
3.1 小结
和 二阶梯度法（牛顿法）\(\eqref{二阶梯度法}\) 的增量解进行对比，可以发现 Gauss-Newton 中用 \(\boldsymbol{J}^{T}\boldsymbol{J}\) 作为牛顿法中二阶 Hessian 矩阵的近似，从而省略了计算 \(\boldsymbol{H}\) 的过程。
Gauss-Newton 的核心在于求解增量方程，但它要求 \(\boldsymbol{H}\) 是可逆且正定的，但实际中计算的 \(\boldsymbol{H}=\boldsymbol{J}^{T}\boldsymbol{J}\) 通常只是半正定的。因此，在使用 Gauss-Newton 的时候可能会出现 \(\boldsymbol{H}\) 是奇异或者病态的情况，导致算法不收敛。更严重的是，就算我们假设 \(\boldsymbol{H}\) 非奇异也非病态，如果我们求出来的步长 \(\Delta \boldsymbol{x}\) 太大，也会导致我们采用的局部近似 \(\eqref{局部近似}\) 不够准确，这样一来我们甚至都无法保证它的迭代收敛，哪怕是让目标函数变得更大都是有可能的。
尽管 Gauss-Newton 有这些缺点，但是它依然值得我们去学习，因为在非线性优化里，相当多的算法都可以归结为 Gauss-Newton 的变种。这些算法都借助了 Gauss-Newton 法的思想并且通过自己的改进修正 Gauss-Newton 法的缺点。例如一些线搜索方法 (line search method)，这类改进就是加入了一个标量 \(\alpha\)，在确定了 \(\Delta \boldsymbol{x}\) 进一步找到 \(\alpha\) 使得 \(\|f(\boldsymbol{x}+\alpha \Delta \boldsymbol{x})\|_{2}^{2}\) 达到最小，而不是像 Gauss-Newton 法那样简单地令 \(α = 1\)。
Levenberg-Marquadt 方法在一定程度上修正了这些问题，被称之为阻尼牛顿法(Damped Newton Method)，一般认为它比高斯牛顿更为鲁棒。尽管它的收敛速度可能会比高斯牛顿更慢，但是在 SLAM 里面却被大量应用。
4. Levenberg-Marquadt
前文提到当 \(\Delta \boldsymbol{x}\)
太大时，Gauss-Newton 方法中采用的近似二阶泰勒展开 \(\eqref{局部近似}\)
存在不够准确的问题，很自然的想到对 \(\Delta
\boldsymbol{x}\) 增加一个信赖区域 (Trust Region)，不能让 \(\Delta \boldsymbol{x}\)
因太大而使得近似不够准却。这种添加信赖区域的方法在非线性优化中被称为
信赖区域方法 (Trust Region Method)
。只有在规定的信赖区域中的 \(\Delta
\boldsymbol{x}\) 才认为近似是有效的。
那么如何确定这个信赖区域的范围呢？一个比较好的方法是根据我们的近似模型跟实际函数之间的差异来确定这个范围：如果差异小，我们就让范围尽可能大；如果差异大，我们就缩小这个近似范围。因此，考虑使用下式来判断泰勒近似是否够好。
\[
\rho=\frac{f(\boldsymbol{x+\Delta
\boldsymbol{x}})-f(\boldsymbol{x})}{\boldsymbol{J}(\boldsymbol{x})\Delta
\boldsymbol{x}}=\frac{实际函数下降值}{近似模型下降值}
\] 显然，当 \(\rho \rightarrow
1\) 时，近似效果最好；
当 \(\rho \ll 1\)
时，说明实际下降值远小于近似下降值，我们需要缩小近似范围，即减小 \(\Delta \boldsymbol{x}\) ；
当 \(\rho \gg 1\)
时，说明实际下降值远大于近似下降值，我们需要增大近似范围，即增大 \(\Delta \boldsymbol{x}\) 。
依此，我们可以改进 Gauss-Newton 非线性优化的框架：
- 给定初始值 \(\boldsymbol{x}_0\)，以及初始优化半径 \(\mu\)；
- 对于第 \(k\) 次迭代，求解：\[\label{L-M} \underset{\Delta \boldsymbol{x}_k}{min} \frac{1}{2}\|f(\boldsymbol{x}_k)+\boldsymbol{J}(\boldsymbol{x}_k)\Delta \boldsymbol{x}_k\|_{2}^{2},\ s.t.\|\boldsymbol{D}\Delta \boldsymbol{x}_k\|^{2}\le \mu\] 这里 \(\mu\) 是信赖区域的半径，\(\boldsymbol{D}\) 是一个调整矩阵
- 计算 \(\rho\)；
- 如果 \(\rho \gt \frac{3}{4}\)，则 \(\mu = 2\mu\)；如果 \(\rho \lt \frac{1}{4}\)，则 \(\mu = 0.5\mu\)；否则，\(\mu\) 保持不变。这里的阈值和扩大倍数都是经验值，可视情况调整；
- 如果 \(\rho\) 大于某阈值，认为近似可行。令 \(\boldsymbol{x}_{k+1}=\boldsymbol{x}_k+\Delta \boldsymbol{x}_k\)；
- 判断算法是否收敛。如果不收敛则返回2，否则结束。
关于 \(\eqref{L-M}\) 中 \(\boldsymbol{D}\) 的说明：
当 \(\boldsymbol{D}=\boldsymbol{I}\)
时，这里的约束条件可看作是将 \(\Delta
\boldsymbol{x}\) 约束在一个半径为 \(\mu\) 的球内；当 \(\boldsymbol{D} \ne \boldsymbol{I}\)
时，可以认为这是一个椭球。 在 L-M 算法中，取 \(\boldsymbol{D}\) 为非负对角矩阵（由 \(\boldsymbol{J}^{T}\boldsymbol{J}\)
的对角元素平方根构成），使得在梯度小的维度上约束范围更大一些。
根据\(\eqref{L-M}\) 求解增量， 拉格朗日法： \[ \nonumber L=\frac{1}{2}\|f(\boldsymbol{x}_k)+\boldsymbol{J}(\boldsymbol{x}_k)\Delta \boldsymbol{x}_k\|_{2}^{2}+\frac{\lambda}{2}\|\boldsymbol{D}\Delta\boldsymbol{x}_k\|_{2}^{2},\lambda \gt 0 \] 令偏导为0，可得增量方程： \[ \frac{\partial L}{\partial \Delta\boldsymbol{x}_k}=\boldsymbol{0} \Longrightarrow (\boldsymbol{H}+\lambda \boldsymbol{D}^{T}\boldsymbol{D})\Delta \boldsymbol{x}_k=\boldsymbol{g} \]
从结果来看，当 \(\lambda\) 较小时，\(\boldsymbol{H}\) 占主导地位，说明二次近似模型效果在该范围内是比较好的，此时 L-M 法接近于 G-N 法；当 \(\lambda\) 较大时，考虑 \(\boldsymbol{D}=\boldsymbol{I}\) 的简化形式，此时 \(\lambda \boldsymbol{D}^T \boldsymbol{D}=\lambda \boldsymbol{I}\) 占主导地位，说明二次近似的不够好，此时 L-M 法更接近于 一阶梯度下降法（最速下降法）。
L-M 的求解方式，可在一定程度上避免线性方程组的系数矩阵的非奇异和病态问题，提供更稳定更准确的增量 \(\Delta \boldsymbol{x}\)。
总而言之，非线性优化问题的框架，分为 Line Search 和 Trust Region
两类。
Line Search 先固定搜索方向，然后在该方向寻找步长，以最速下降法和
Gauss-Newton 法为代表。而 TrustRegion
则先固定搜索区域，再考虑找该区域内的最优点。此类方法以 L-M
为代表。实际问题中，我们通常选择 G-N 或 L-M 之一作为梯度下降策略。
5. 参考资料
[1] 视觉 SLAM 十四讲 —— 高翔