Ceres Solvel非线性最小二乘问题的解算

发表于 2020-07-20 更新于 2024-04-26 分类于 Ceres Solvel ，最小二乘优化

简介

如果想要高效地使用Ceres Solver，需要掌握一定的非线性最小二乘解算基础知识。所以在这一部分将将要介绍Ceres中核心优化算法的工作原理。

设$x \in \mathbb{R}^n$是一个$n$维向量，并且$F(x) = [f_1(x),\cdots,f_m(x)] ^ T$是关于$x$的$n$维方程。那么我们关注下列优化问题

$$
argmin_x \frac{1}{2} \left | F(x) \right | ^ 2 \
L \le x \le U
$$
其中$L$和$U$分别是参数向量$x$的下限和上限。

因为对于一个一般的函数，求解全局最小值常常非常棘手，我们不得不关注局部最小值。$F(x)$的雅可比矩阵$J(x)$是一个$m \times n$的矩阵，其中$J_{ij}(x) = \partial _j f_i(x)$，函数的梯度向量$g(x) = \nabla \frac{1}{2} \left | F(x) \right | ^ 2 = J(x)^T F(x)$。

在计算非线性优化问题的一个通用的策略是，求解原问题的近似简化问题。在每一次循环中，根据近似问题的解可以确定向量$x$的修正值$\triangle x$。对于非线性最小二问题，我们可以通过线性化来建立近似问题，即$F(x + \triangle x) \approx F(x) + J(x) \triangle x$。那么上述问题就变成下列问题：

$$
min_{\triangle x} \frac{1}{2} \left |F(x) + J(x) \triangle x \right |^2
$$

这里官方教程跳了几步。要求全局最小值非常棘手，所以转而求局部最小值。而求局部最小值就是，从一个任意的起始点，观测四周的“更小值”。如果观测四周都比当前点大，那么当前点就是局部最小值点，算法达到收敛。否则，设这个新找到的最小值点为“当前点”，重复这一步骤。这也就是下文中”用$x \longleftarrow x + \triangle x$来更新“的含义。所以现在问题变成了，如何求解四周的点的值。即，给$x$赋予一个步长$\triangle x$，观察周围的$F(x + \triangle x)$，并且寻找最小值。

根据步长的控制方法，非线性优化可以分成两大类。

置信域法Trust Region （也有文献成为信任域法）置信域方法在搜索空间的子集内应用模型函数（通常是二次方程）来近似目标函数，这个空间被称为置信域。如果模型函数成功地使真正的目标函数最小化，则扩大置信域。否则收缩置信域并且再次尝试求解模型函数的优化问题。

线搜索法Line Search 线搜索方法首先找到一个下降的方向，目标函数将沿其下降。然后再确定步长，从而决定沿该方向到底移动多远。下降方向可以通过各种方法计算，如梯度下降法、牛顿法和拟牛顿法。步长可以精确或不精确地确定。

这里我们只对共轭梯度法进行介绍（因为我们用的库函数所采用的方法就是共轭梯度法）

对以下形式的线性最小二乘问题：

$$
min_{\triangle x} \frac{1}{2} \left |F(x) + J(x) \triangle x \right |^2
$$
令$H(x) = J(x) ^ TJ(x)$且$g(x) = - J(x) ^ {T}f(x)$。为了标记方便，我们舍弃对$x$的依赖，那么很容易看出求解上式等价于求解Normal Equation方程

$$
H\triangle x = g
$$
我们在程序中使用的计算上式的线性解算器是CGNR，该解算器在法向方程上使用共轭梯度求解器，但没有明确地构建Normal Equation方程。它利用下列关系：

$$
Hx = J ^ TJx = J^T(Jx)
$$
共轭梯度的收敛取决于调节器数量$\kappa(H)$，通常$H$的条件比较弱，并且必须使用预调节器才能获得合理的性能。目前只有JACOBI预处理器可以用于CGNR。它使用$H$的块对角线来预处理Normal Equation方程。

当用户选择CGNR作为线性求解器时，Ceres会自动从精确的步长算法切换到不精确的步长算法。

JACOBI预处理器

求解Normal Equation方程的收敛速度取决于$H$的特征值的分布。一个有用的上界是$\sqrt {\kappa H}$，$\sqrt {\kappa H}$是矩阵$H$的条件值。对于大多数$BA$问题来说，$\sqrt {\kappa H}$很大，直接应用共轭梯度法对该方程进行求解会导致糟糕的性能。

这个问题的解决方法是用一个预条件系统来替换$H\triangle x = g$方程。给定一个线性系统，$Ax = b$，预条件器$M$，预条件系统$M^{-1} Ax = M^{-1}b$。这个算法被称作预处理共轭梯度法（PCG），它最坏的情况现在取决于预条件矩阵$\sqrt{M^{-1}A}$的条件数。

使用预条件器$M$的计算成本是计算$M$以及任意向量$y$的乘积$M^{-1}y$，因此，要考虑两个相互竞争的因素：$H$的结构中有多少是由$M$占据的，以至于条件值$\kappa(HM^{-1})$，以及构造和使用$M$的计算成本。理想的预调节器是条件值$\kappa(M^{-1}A) = 1$。使用$M = A$可以实现这一效果，但是不是一个实际的选择，因为应用这个预调函数需要一个线性方程组，等价于无预先条件的问题。通常情况下，$M$的信息越多，使用的成本就越高。例如，基于不完全$Cholesky$分解的预条件器比$Jacobi$预条件器具有更好的收敛性能，但代价也更高。

所有的预调节器中最简单的是对角或雅可比的预调机。例如，$M = diag(A)$。对于像$H$这样的块状结构矩阵，可以将其推广到$Jacobi$的预调节器。Ceres实现的块$Jacobi$预调节器叫做JACOBI，当使用CGNR时，它指的是$H$的块对角线，即

$$
C_{i j}=\left{\begin{array}{cc}
A_{i i} & \text { if } i=j \
0 & \text { 其他 }
\end{array}\right.
$$

PCG算法的过程如下：假设$x_0 \in \mathbb{R}^n$是初始向量，

$$
\begin{aligned}
&r_{0}=b-A x_{0}\
&z_{0}=C^{-1} r_{0}\
&d_{0}=z_{0}\
&\text { For } k=0,1,2, \ldots\
&\alpha_{k}=\frac{z_{k}^{\top} r_{k}}{d_{k}^{\top} A d_{k}}\
&x_{k+1}=x_{k}+\alpha_{k} d_{k}\
&r_{k+1}=r_{k}-\alpha_{k} A d_{k}\
&z_{k+1}=C^{-1} r_{k+1}\
&\beta_{k+1}=\frac{z_{k+1}^{\top} r_{k+1}}{z_{k}^{\top} r_{k}}\
&d_{k+1}=x_{k+1}+\beta_{k+1} d_{k}
\end{aligned}
$$