Ceres(5): Solver

转载地址：

http://blog.csdn.net/HUAJUN998/article/details/76271938

http://ceres-solver.org/nnls_solving.html

定义 x∈Rn是一个n-维的向量 （向量存储的都是变量），F(x)=[f1(x),...,fm(x)]⊤ 是一个m维的关于x的方程，我们的目的是解决优化问题 

(1)argminx12∥F(x)∥2 .L≤x≤U

由于对F(x)的全局最小化是一个棘手的问题，我们将不得不满足于寻找局部最小值。

F(x) 的Jacobian J(x) 是一个m×n的矩阵，Jij(x)=∂jfi(x) ,梯度向量 g(x)=∇12|F(x)|2=J(x)⊤F(x).

解决非线性优化问题的一般策略是，解决原始问题的一系列逼近问题。在每次迭代中，这个近似都得到了解决，以确定Δx修正向量x。对于非线性最小二乘，可以利用线性化F(x+Δx)≈F(x)+J(x)Δx来构造一个近似，这就引出了下面的线性最小二乘问题: 

(2)minΔx12∥J(x)Δx+F(x)∥2

Trust Region Methods

整体的大致框架如下所示：

μ是置信域半径，D(x)是一个定义在F(x)域上的矩阵， ρ用来衡量，线性模型在变量进行 delta_x的变化后，

预测的目标函数的变化量，和实际目标函数(非线性模型)的变化量之间的近似程度,

delta_x的解依赖于置信域半径 u 的大小,所以用 rho 来调整 u。

 在置信域算法中的关键步骤是解决下面的约束优化问题

Ceres实现了两个置信域算法L-M和Dogleg。如果有边界约束，那么每一个都被加了一个线性搜索。用户可以通过设置Solver::Options::trust_region_strategy_type.来选择它们。

Levenberg-Marquardt 

L-M算法是解决非线性最小二乘问题最流行的方法，它也是最先开发出来的置信域方法。Ceres实现了L-M算法的一个精确step和一个不精确step版本。

可以证明，解决(3)可以通过解决具有如下形式的非约束优化来获得 

argminΔx12∥J(x)Δx+F(x)∥2+λ∥D(x)Δx∥2

λ 是拉格朗日乘数与 μ 逆相关 

argminΔx12∥J(x)Δx+F(x)∥2+1μ∥D(x)Δx∥2

                                                                         （4）

矩阵D(x)是一个非负的对角矩阵，通常是矩阵J(x)⊤J(x)的对角线的平方根。

/*****************************/

/*****************************/

Ceres提供了两类方法来解(5)式：分解法，迭代法。

分解方法是精确的求解(4)，使用 Cholesky 或 QR分解的方法，得到精确的step。

迭代的方法：比如牛顿法。这类方法有两个要素：

一：一个简单的方法来大致给出线性化的系统模型的解；

二：迭代的终止条件，一个典型的比如：

（6）∥H(x)Δx+g(x)∥≤ηk∥g(x)∥.

 0<ηk<1is known as the forcing sequence.

当用户选择基于分解的linear solver时，使用精确step L-M算法。当选择迭代的linear solver时，将使用不精确的step L-M算法。

Dogleg

关键思想是计算两个向量 

ΔxGauss-NewtonΔxCauchy=argminΔx12∥J(x)Δx+f(x)∥2.=−∥g(x)∥2∥J(x)g(x)∥2g(x).

Dogleg methods finds a vector Δx defined by ΔxGauss-Newtonand ΔxCauchy that solves the trust region problem.

Ceres支持两种变体，可以通过设置Solver::Options::dogleg_type选择。

TRADITIONAL_DOGLEG 就像 Powell 所描述的，用高斯-牛顿和柯西矢量构造两条线段，发现沿着这条线最远处的点，形状就像一条狗腿(因此得名)。关于具体的推理和计算的更多细节，请参见 K. Madsen, H.B. Nielsen, and O. Tingleff, Methods for nonlinear least squares problems, 2004.

SUBSPACE_DOGLEG 是一种更复杂的方法 ，它考虑了由这两个向量形成的2维子空间，并在这个子空间中找到最小化置信域问题的点。R.H. Byrd, R.B. Schnabel, and G.A. Shultz, Approximate solution of the trust region problem by minimization over two dimensional subspaces, Mathematical programming, 40(1):247–263, 1988.

对L-M来说，Dogleg的关键优势在于，如果对μ的某一特定选择计算的step并没有导致目标函数值的明显减少，那么L-M得用更小μ值从头开始解决线性近似，而Dogleg只需要计算高斯-牛顿和柯西矢量之间的值，它们都不依赖于 μ的值。

Non-monotonic(非单调) Steps

基本的置信域算法是下降算法，即在迭代中选取  X=X+delta_x   的时候，只有使目标函数减小的点才会被接受。

放宽这一要求，可以使算法在长期内更有效，代价是在局部过程中，增加了目标函数的值。

将 Solver::Options::use_nonmonotonic_steps设为 true ，将使用 由Conn, Gould & Toint 提出的非单调置信域算法。

在优化过程中，尽管目标函数的值可能大于遇到的最小值，但是返回给用户的最终参数是与所有迭代的最小cost对应的参数。

Line Search Methods

目前，Ceres的线性搜索方法不能处理边界约束，因此只能用于解决不受约束的问题。

线性搜索方法

1. 给定一个初始点x
2. Δx=−H−1(x)g(x)
3. argminμ12∥F(x+μΔx)∥2
4. x=x+μΔx
5. 跳转到2

H(x)是目标函数Hessian的近似，g(x)是x处的梯度。依靠对H(x)的选择我们得到各种不同的搜索方向Δx。 
第4步，这是一个一维的优化或者是沿着delta_x 方向的线性搜索。

不同的线搜索算法在选择搜索方向Δx和沿着Δx的一维优化方法上有不同的选择。H(x)的选择是这些方法中计算复杂度的主要来源。

目前，Ceres Solver支持4种搜索方向的选择，都是针对大尺度的问题。

1. STEEPEST_DESCENT 这对应于H(x)是单位矩阵。除了最简单的问题，这并不是一个好的搜索方向。它只是为了完整性而被包含在这里。
2. NONLINEAR_CONJUGATE_GRADIENT 共轭梯度法在非线性函数中的generalization。generalization可以以多种不同的方式实现，从而产生不同的搜索方向。Ceres Solver目前支持FLETCHER_REEVES, POLAK_RIBIEREand HESTENES_STIEFEL方向。
3. BFGS Secant的方法在多维中的generalization，在这种情况下，得到了一种完全的、密集的逆Hessian的近似，并用来计算拟牛顿step。BFGS是目前已知的最著名的拟牛顿算法。
4. LBFGS 用有限的内存近似完整的BFGS方法，在此方法中，最后的M次迭代被用来估计Hessian的逆。

目前，Ceres Solver同时支持基于回溯和插值的 Armijo线性搜索算法，以及一种基于分割/缩放插值Wolfe condition线性搜索算法。然而，注意，为了保证BFGS和LBFGS方法的基础假设，应该使用Wolfe线性搜索算法。

LinearSolver

回想一下，在上面描述的两个置信域方法中，关键的计算是一个具有如下形式的线性最小二乘问题 

(7)minΔx12∥J(x)Δx+f(x)∥2.

让H(x)=J(x)⊤J(x),g(x)=−J(x)⊤f(x)。很容易看出解(7)等价于解 normal equations （线性最小二乘的解）

(8)HΔx=g

Ceres提供了许多不同的方案解决(8)。

DENSE_QR

对于较小的问题(几百个参数和几千个残差)，相对稠密的Jacobians，DENSE_QR 是一个选择。让J=QR是J的QR分解。Q是一个正交矩阵，R是一个上三角矩阵。然后我们能得到（8）的解 

Δx∗=−R−1Q⊤f

Ceres使用了Eigen的稠密QR分解步骤。

DENSE_NORMAL_CHOLESKY & SPARSE_NORMAL_CHOLESKY

大型非线性最小二乘问题通常是稀疏的。在这种情况下，使用稠密的QR分解是低效的。让H=R⊤R是normal equations 的Cholesky 分解，R是一个上三角矩阵，对（8）的解是 

Δx∗=R−1R−⊤g.

细心的读者会注意到，H的Cholesky分解中的R和J的QR分解中的R是相同的上三角矩阵。由于Q是正交矩阵，J=QR可以得到J⊤J=R⊤Q⊤QR=R⊤R，有两种不同的Cholesky分解——稀疏和稠密。

DENSE_NORMAL_CHOLESKY

Ceres usesEigen‘s dense LDLT factorization routines.

SPARSE_NORMAL_CHOLESKY

Ceres uses the sparse Cholesky factorization routines in Professor Tim Davis’ SuiteSparse or CXSparse packages

or the sparse Cholesky factorization algorithm in Eigen

CGNR

对于一般的稀疏问题来说，如果这个问题对CHOLMOD来说太大，或者一个稀疏的线性代数库没有链接到Ceres，那么另一个选择就是CGNR方案。

这个方案是对normal equations使用共轭梯度求解，但不显式构造normal equations。它利用 

Hx=J⊤Jx=J⊤(Jx)

共轭梯度的收敛依赖于调节数κ(H)。通常情况下，H的条件很差（poorly conditioned），必须使用Preconditioner才能获得合理的性能。

目前只有JACOBI preconditioner 可与CGNR一起使用。它使用了“H”的块对角特性来precondition the normal equations。

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

DENSE_SCHUR & SPARSE_SCHUR

虽然可以使用SPARSE_NORMAL_CHOLESKY来解决BA问题，但BA问题有一个特殊的结构，可以构造一个更有效的方案解决(8)。

假设SfM问题由p个摄像机和q个点组成，而变量向量x有块结构x=[y1,...,yp,z1,...,zq]。y和z对应的是相机和点参数。此外，让摄像机块大小为c，并且点块大小为s(对于大多数问题，c=6-9和s=3)。Ceres并没有对这些块大小施加任何恒定性要求，但是选择它们是常量，简化了说明。

BA问题的一个关键特性是，没有一个fi包含两个或多个点块。这就意味着矩阵的H的形式是 

(9)H=[BE⊤EC] ,

B∈Rpc×pc是一个有p个c×c大小的稀疏矩阵块，C∈Rqs×qs是一个有q个s×s大小的对角矩阵块，E∈Rpc×qs是一个对每个观察有c×s大小的稀疏矩阵块。现在我们来分割块Δx=[Δy,Δz],g=[v,w]重新分析(8)作为块结构化线性系统 

(10)[BE⊤EC][ΔyΔz]=[vw] ,

并应用高斯消元。正如我们上面所提到的，C是一个块对角矩阵，对角块大小s×s。因此，计算C的逆矩阵通过每个块的逆运算是很简单的。这让我们可以消除Δz通过观察Δz=C−1(w−E⊤Δy)。得到 

(11)[B−EC−1E⊤]Δy=v−EC−1w .

其中 

S=B−EC−1E⊤

is theSchur complement of C in H.它也被称为 reduced camera matrix，因为参与(11)的唯一变量是与摄像机相对应的变量。S∈Rpc×pc是一个块结构的对称正定矩阵，块的大小是c×c。块Sij对应于一对图像i和j的值是非0的，当且仅当两张图像观察到至少同一个点时。

现在，(10)可以通过先构造S，求解Δy，然后用Δy来求Δz的值。因此，n×n,n=pc+qs线性系统的解被简化为块对角矩阵C的逆矩阵，几个矩阵与矩阵，矩阵与向量相乘，以及块稀疏pc×pc线性系统(11)的解。对于大多数问题，相机的数量要远小于点的数量，p≪q,因此求解（11）比（10）简单的多，这就是 Schur complement的trick。

对于(11)式仍有其他的求解途径。对于对称正定系统的求解完全可以通过Cholesky分解。根据系统的对称正定矩阵的结构，有两种方式。

第一种，直接分解，我们存储和分解S当作一个稠密的矩阵。这个方法有O(p2)的空间复杂度和O(p3)的时间复杂度，只有在相机数目只有几百时是可行的。Ceres实现这个策略叫做DENSE_SCHUR。

但是，当S是典型的相当稀疏的矩阵时，这个时候因为大多数图片只能看到一小部分场景。这就让我们有了第二个选择：稀疏直接法。这些方法将S作为一个稀疏矩阵，使用行和列的重新排序算法来最大限度地增加Cholesky 分解的稀疏性，并将它们的计算工作重点放在分解的非零部分上。稀疏直接法，依赖于Schur complement的稀疏结构，使BA明显加速，相比基于稠密分解的算法。Ceres实现这个策略叫做SPARSE_SCHUR 。

ITERATIVE_SCHUR

另一种求解BA问题的途径是，应用 Preconditioned Conjugate Gradients 对 reduced camera matrix S，而不是对H。这样做的一个原因是S要比H小的多，但是更重要的是，你可以看到κ(S)≤κ(H)。Ceres对S实现Conjugate Gradients的方法叫做ITERATIVE_SCHUR 。当用户选择ITERATIVE_SCHUR作为线性求解方式时，Ceres自动将精确求解转换为不精确求解。

当应用Conjuagate Gradients的关键是，估计矩阵向量乘积Sx,x是任意向量。有两种方法可以对该乘积估计，并且可以通过Solver::Options::use_explicit_schur_complement来控制。根据要解决的问题，这两种方法的性能差异可能相当大。

    Implicit（隐式的）  这是默认的。隐式求值适用于大问题。

    Explicit （显式的） 矩阵向量乘积求值的复杂性与雅可比矩阵中的非零数有关。对于小到中型的问题，构建Schur complement的成本足够小，可以在内存中显式地构造它，并使用它来评估Sx乘积。

Preconditioner略

Ordering

在线性求解器中，消除变量的顺序对算法的效率和准确性有很大的影响。例如，当做稀疏的Cholesky 分解时，对于很多矩阵，一个好的消元排序将会给出一个带有O(n)的存储，一个坏的消元顺序将导致一个密集的分解。(如果上面先解 delta_z,再解 delta_y,你会发现求逆的矩阵会非常大，且很难求)

ΔzΔz

Ceres允许用户向解析器提供各种关于变量消除顺序的提示。这可以从没有任何提示（这时解析器自己来决定最佳的顺序），到一个精确的变量被消除的顺序，以及在其它们之间的各种可能性。

ParameterBlockOrdering 类的实例用于向Ceres传递这些信息。

排序是参数块的有序分割。每个参数块仅属于一个组，每个组都有一个与之关联的唯一整数，这决定了它们在组中的顺序。我们称这些组为消除组 （Elimination Groups）。

Ceres确保最低编号的消除组的参数块首先被消除，然后消除下一个最低编号的消除组，依次下去。在每个消除组中，Ceres自由地按照它所选择的参数块顺序来排序。

下面这些类的APIs 控制着整个解析器的行为，性能和各种详细的报告，是需要我们熟练掌握的。详见http://blog.csdn.net/HUAJUN998/article/details/76341292

http://ceres-solver.org/nnls_solving.html#solver-options

classSolver::Options

classParameterBlockOrdering

classIterationSummary

classCRSMatrix

classSolver::Summary