起源
L-M是参数估算和数值拟合都经常会用到的一种数学方法，在研究全景图像时，就留意到Position Adjust的Optimize阶段，诸多Optimize工具都用到L-M来估计准确的Position;
BundleAdjust的思路
BundleAdjust是论文[1]提出的对camera parameters进行统一估算的算法，其核心思路就是选定sphere坐标系到二维平面坐标系的转换矩阵后，根据Invariant Features得到的匹配点对，构造代价函数：
公式14给定的是单个匹配点在经过初始化的转换矩阵映射后的点的距离，而17式对Inner Point和Outlier进行了区分（非常有价值的思路）；论文后续的讨论都是对此式的求解以及降低计算量的方法；当然，解法的落点也是L-M算法
上述就是研究L-M算法的实际问题背景
基础概念
从L-M算法的定义论述，就可以知道L-M绕不开牛顿法和梯度下降法，为了对L-M算法有完整的直觉上的理解，得益于这么多次搬家都没有扔掉的数学分析和数值分析课本，从[2]和[3]可以得到基础概念的理解，微分方程/泰勒展开/雅可比行列式/Hessian矩阵；
[4]的Chapter 10，讲述了非线性方程组的数值解法，里面涉及到的一些算法和思路从一元方程的求解（Chapter 2）出发会更好理解，而且总体的思路上一元方程的求解和非线性方程组的解法一脉相承；这里面我觉得最核心的思想是不动点迭代
有前面的基础后基本能理解L-M的算法思路，[5]里面大量讨论了迭代算法的原理和细节，linear search/雅可比矩阵的替代B矩阵/割线替代等
下面是个人一些理解的总结
一、面对的问题
1 非线性方程组求解
f1(x1,x2,...,xn) = 0
f2(x1,x2,...,xn) = 0
...
fn(x1,x2,...,xn) = 0
2 最小化问题
最小化问题最典型的就是最小二乘问题
3 问题的统一
非线性方程组的求解可以转换为最小化问题；如1所示的方程组的解可以通过求解下面方程G的最小值来得到；
G(x1,x2,...,xn) = f1(x1,x2,...,xn) ^2 + f1(x1,x2,...,xn) ^2 + ... + fn(x1,x2,...,xn) ^2
这个问题的统一涉及到对多元函数梯度和雅可比行列式的区别的理解，后续会讲到
数值迭代方法
1 数值求解
对一元以及多元非线性方程（组）直接求解是很困难的，印象里有个数学家证明过，对于某些高元非线性方程是无法直接求解的，只能通过数值求解；即通过计算的方法来迭代求解；
对于任意给定的函数求最值或者极值或者方程（组）求解，我们已知的信息就是函数或者方程本身；我们想要得到的是一组(x1,x2,...xn)值，使得函数（方程）在该值上具备某些特性（等于0，具备最值或者极值），那么我们就要尝试构造函数或者方程跟(x1,x2,...xn)之间的关系；
那么不动点迭代的出现，其基本思想就非常清晰了
2 不动点迭代
p=G(p)
如果给定点p，满足上式，即p为函数G(x)的不动点；这个概念很容易推广到多元函数上
观察这个式子，左边是我们想要求解的项，右边是函数本身，该式子直接构造了求解项和已知项的联系；
并且，该式能很容易地跟方程求解问题关联，G(x)的不动点就是函数 F(x)=x - G(x)的解，那么如果我们需要求解F(x)，只需要构造一个G(x)使得 F(x)=x - G(x)；就可以将F(x)的求解问题转换为G(x)求不动点的问题
x=G(x)如何通过数值方法求解
构造不动点方程的最大优势是，通过数值方法求解不动点几乎是天然的，观察上式它直接给出了不动点和函数直接的关系，假设给定一个x0，如果G(x0)=x0则结束，如果不等，则x1= G(x0)；通过这样的更替，我们天然地能够得到一个x的迭代序列
这就是不动点迭代
3 不动点迭代的理论基础
x= G(x)
上式是否一定有不动点以及对于任意给定的F(x)，我们可以构造的不动点函数G(x)有很多个（参考[4] Chapter 2 EXAMPLE 3）那是否任意构造的G(x)都能够收敛到真正的不动点？答案是否定的
那[4] Chapter 2 Theorem 2.2给出了不动点一定存在的理论前提
[4] Chapter 2 Theorem 2.3给出了如何考察构造出来的G(x)是否能真正收敛到不动点的理论基础；
这部分偏理论，不展开讲
4 泰勒展开
上式就是泰勒展开式，会发现一阶泰勒展开非常适合构造不动点迭代，并且很容易推导出来，一阶泰勒展开具备收敛到不动点的性质。而且其误差能进行定量的分析，即求解F(x) = 0的问题可以简化为求解如下的不动点问题：
x = F(x) / F'(x)
迭代公式：
p(n) = p(n-1) - F(p(n-1))/ F'(p(n-1))
这就是牛顿法
5 牛顿法
牛顿法需要知道原函数的一阶导数，在数值方法里面，求解导数很多时候是不方便的（也不通用，这意味着需要针对每一个函数分析其导函数）因此后续多采用的还是割线法以及 method of False Position等变种算法，其通过提供两个初始点来近似F'(x)；可参考[4] Figure 2.10，有直观的图像表示
上述是一元函数的牛顿法论述，在多元函数中，讨论高阶导数变的困难(阅读 [3] 第二十三章，对n维空间/向量函数/雅可比/黑塞矩阵等有具体的定义)，一般最多讨论到二阶导数，通过偏导数，我们得到雅可比矩阵和黑塞矩阵
而因为讨论对象从一元x变换成多元 (x,y,z,....)，其迭代公式中需要将(x1,y1,z1,...)迭代到下一个 (x2,y2,z2,...)；需要的要么是(dx,dy,dz,....)，要么就是变换矩阵A
这两个对应的多元函数中的梯度法和牛顿法
6 梯度法
记得在前面提到的“问题的统一”么，来谈多元函数的梯度和雅可比行列式
在一元函数中导数确定的方向就是梯度方向；
那么多元函数涉及的问题会比较复杂，一般研究的对象会有两种：
1 F(x)的最小化, 典型的像[5] 一直研究的就是Least Squares Problems
2 多元函数线性方程组
f1(x1,x2,...xn) = 0
f2(x1,x2,...xn) = 0
....
fn(x1,x2,...xn) = 0
仔细看这两种问题，一个是求平方和的最小，一个是方程组的解；会发现这两者其实是一个问题；假如方程组有解，那么当一组(x1,x2,...xn)为方程组的解时，F(x)取最小值；如果一组(x1,x2,...xn)使得F(x)取得最小值，即使该最小值大于0，那么这组(x1,x2,...xn)也是最接近方程组的解的一组值。
所以在讨论多元函数相关概念的时候需要区别是针对方程组还是针对F(x)，一般地在多元函数的讨论中，我们说梯度，针对的是F(x)，而当我们讨论雅可比行列式的时候，直接针对的是方程组；而当我们讨论黑塞矩阵，我们针对的是F(x)；
所以梯度下降直接利用F(x)函数对(x1,x2,...,xn)里面的各个分量的偏导数得到梯度方向；而雅可比从组成F(x)的函数分量对对(x1,x2,...,xn)里面的各个分量的偏导数得到
梯度下降是线性收敛，而牛顿法是Quadratic convergence （见 [5] 2.3式以及2.2节）
[4] Chapter 10(简写为 [4]-C10) 和[5] 都介绍牛顿法，[4]-C10直接引入的问题是非线性方程组的解，而[5]直接引入的问题是最小二乘的最小值；
在思路上，[4]-C10引入的还是不动点迭代， EXAMPLE2 是一个具体的构造 P = G(P)的例子；然后引入不动点存在以及收敛的理论基础(Theorem 10.4 以及 Theorem 10.6)，最后介绍牛顿法时水到渠成；
[5] 在思路上从下降法开始介绍，下降法思路总体分为两部，寻找下降方向，确定下降步长，然后迭代
7 在L-M之前
[5] 中对L-M的介绍是，“Levenberg (1944) and later Marquardt (1963) suggested to use a damped Gauss-Newton method”；所以理解L-M需要理解什么是 Gauss-Newton method 以及什么是damped method；
7.1 牛顿法和高斯-牛顿法的区分
牛顿法的思路很直接，利用函数F(x)如果在x处取得最小值，则F'(x)=0的性质来确定下降方向；
F'(X + h) = F'(X) + F''(X)h + 高阶极小项 = 0
当h足够小的时候，忽略高阶极小项；这样直接得到下降方向
h = -F''(X) / F'(X)
如果F''(X)是正定矩阵，很容易证明h一定是下降方向（下降方向的定义 [5] 中 Definition 2.6）；
高斯牛顿法的思路是展开F(X)的各个组成子方程的泰勒展开；研究对象是最小二乘法的最小值定义形式的F(X):
F(X) = 0.5 * ||f(x)||^2
f(x+h) = f(x) + J(x)h + 高阶极小项
当h足够小，忽略高阶极小项；然后将该式子带入F(X)后计算得到([5] 3.7式)；
F(X + h) = L(h) = F(X) + h(t) * J(t) * f + 0.5 * h(t) * J(t) * J * h
J(t)是J的转置，h(t)是 h的转置
研究L(h)函数， L'(h) = J(t) * f + J(t) * J * h L''(h) = J(t) * J
可以看到L''(h)的值跟h没有关系，而且如果J是满秩，L''(h)恒为正定矩阵，这说明L(h)有唯一的最小值，并且可以通过 L'(h) = 0 来确定下降方向；
h = - J(t) * f / ( J(t) * J)
高斯牛顿法和牛顿法研究的对象是不同的，但是思路是一致的，在高斯牛顿法中，可以认为高斯牛顿法将F(X)转换L(h)后，对L(h)的研究用的就是牛顿法；
高斯牛顿法和牛顿法的区别还在于下降速度上，牛顿法是二次收敛速度，高斯牛顿法在特定情况下可以达到二次收敛速度，以及一些特定情况下可以达到超线性收敛，但是一般地，只能得到线性收敛速度（详细见 [5] 3.12式及后面的段落讨论）
7.2 damped method 阻尼法
“Levenberg (1944) and later Marquardt (1963) suggested to use a damped Gauss-Newton method” 高斯牛顿法前面讨论了，那么阻尼法是什么？具体的公式见于 [5] 2.4章节
此处记录下个人的理解：
之前提到，[5] 在思路上从下降法开始介绍，下降法思路总体分为两部，寻找下降方向，确定下降步长，然后迭代
那么不管是梯度法，牛顿法，高斯牛顿法都是确定下降方向的方法，对于步长没有讨论，那么2.4提到 Trust Region 和 damped method两种方法，其核心思想是同时考虑下降方向和步长
如何做呢？通过研究F(X)的二阶泰勒展开式的近似函数的性质：
L(h) = F(X) + h(t) * c + 0.5 * h(t) * B * h
L(h)是一个关于h的函数，然后研究L(h)在一定步长情况下何时取得极小值；
那么damped method就是研究 M(h) = L(h) + 0.5 * u * h(t) * h , u >= 0 何时取得极小值
当其取得极小值时，如果 F(X + h) < F(X) 则说明h是下降方向，则迭代F(x) 否则迭代u
那么L(h)何时取得极小值？ M'(h) = L'(h) + uh = 0 时取得；从而得到下降方向
h = -c / (B + u * I )
对于u的考量可以通过 p = (F(X) - F(X + h) ) / ( L(0) - L(h) ) 来考量，详细讨论见[5] 中 2.18,2.19,2.20式子的讨论
而从该数学公式上看，p越接近1说明L函数越能近似F函数，即在h方向上，L函数能很好地近似F函数，当p越接近0，甚至为负值时，则说明L函数不能很好地近似F函数；在[5] 3.2节3.14式子后继续讨论了gain ratio
对u的考量就是讨论当h不是下降方向的时候迭代u的方法以及该方法的数学理论基础，此处不赘述
8 L-M
[5] 中对L-M的介绍是，“Levenberg (1944) and later Marquardt (1963) suggested to use a damped Gauss-Newton method”；
所以L-M其实是一种 damped Gauss-Newton method；所以如果理解了Gauss-Newton method以及 Damped Method，那么L-M其实很好理解，L-M定义的下降方向如下：
(J(t) * J + u * I) * h = -J(t) * f
讨论下L-M的优点，首先对所有 u > 0 ，系数矩阵 (J(t) * J + u * I) 是正定矩阵，这确保了L-M方法确定的h方向一定是下降方向；另外就是当u较大的时候，L-M接近梯度下降方向，u较小的时候接近高斯-牛顿方向
9 割线法
牛顿法和L-M都有割线法；所有的割线法解决的问题都是替代雅可比矩阵；雅可比矩阵在实际问题中不好计算，因为我们面对的方程是多样的（各种模型组成的函数差异），如果都需要其偏导函数，那么无疑增加来解决问题的难度以及丧失了通用性，割线法只用到了函数本身，使得问题的解决能够通用并能减少计算量
9.1 牛顿法的割线替代
[4] 2.3讲到Newton‘s Method的两种割线方法： Secant Method 和 Method of False Method， 思路类似，使用两个初始点，然后用F(xn) - F(xn-1) / (xn - (xn-1) ) 来近似导数；此处需要注意，随着xn序列的更新，实际上导数的近似 F(xn) - F(xn-1) / (xn - (xn-1) ) 的值也在跟着更新
9.2 L-M 的割线替代B矩阵
如下所示的L-M方法中最核心的计算就是雅可比矩阵，那么我们同理假设有一个矩阵B来近似J，那么核心就是如何得到初始的B0以及如何根据xn的序列来更新B
(J(t) * J + u * I) * h = -J(t) * f
10 OpenCV LMSolver 解析
L-M在OpenCV中的实现有两个，一个是CvLevMarq，一个是LMSolver； 建议使用LMSolver，根据注释，LMSolver 是Matlab的LMSolve包转换成C++语言，“translation to C++ from the Matlab's LMSolve package by Miroslav Balda”；使用LMSolver的时候只需要实现LMSolver::Callback；计算每次迭代的雅可比矩阵以及error；
11 Bundle Adjust的 LMSolver 实现
OpenCV的BundleAdjusterBase使用 CvLevMarq 来进行L-M求解，而且雅可比矩阵计算时，每个方程使用了7参数；我们尝试使用LMSolver来进行重写，并且简化方程参数；关于映射模型，更详细的介绍参见：
Opencv2.4.9源码分析——Stitching（二）
Opencv2.4.9源码分析——Stitching（三）
[参考资料]
[1] Automatic Panoramic Image Stitching using Invariant Features
[2] 数学分析（上）第三版 华东师范大学数学系 第五 导数和微分、第六章 微分中值定理及其应用
[3] 数学分析（下）第三版 华东师范大学数学系 第二十三章 流形上微积分学初阶
[4] NUMERICAL ANALYSIS (Seventh Edition) Chapter 2 Solution of Equations in one Variable, Chapter 10 Numerical Solutions of Nonlinear Systems of Equations
[5] METHODS FOR NON-LINEAR LEAST SQUARES PROBLEMS 2nd Edition, April 2004,Informatics and Mathematical Modeling Technical University of Denmark
[6] https://blog.csdn.net/zhaocj/article/details/78809143
[7] https://blog.csdn.net/zhaocj/article/details/78799194
