LM算法原理

非线性最小二乘法优化

高斯-牛顿法

参考文章：[优化] Gauss-Newton非线性最小二乘算法
算法流程如下图（来自参考文章）所示：

接下来本文使用的数学符号意义与上图一样。其中$x$是需要求解的参数，$f(x)$是一个残差向量。比如有一个优化问题，$y=asin(wt+b)+c$，给出m个数据$(t_i,y_i)(i=0,1,\cdots ,m-1)$，则
$$\begin{gathered} x=[a,w,b,c{]}^T \\ f(x)=[y_0-(asin(w{t}_0+b)+c),y_1-(asin(w{t}_1+b)+c),\cdots ,y_{m-1}-(asin(w{t}_{m-1}+b)+c){]}^T \end{gathered}$$
则$||f(x)|{|}^2$(向量二范数)就是最小二乘法的损失值。
设损失函数$l(x)=\frac{1}{2}||f(x)|{|}^2$。
另外，$J(x)$为$f(x)$的雅可比矩阵，假设$x$的长度为n，$f(x)$长度为m，则$J(X)$矩阵大小为(m,n)。
$H=J^{T}J$为$f(x)$的黑塞矩阵的近似矩阵。$B=-J^{T}f(x)$为损失函数$l(x)$（$\frac{1}{2}$只是为了求导后约掉$||f(x)|{|}^2$的指数2）的负梯度$-\frac{\partial l(x)}{\partial x}$。

最速下降法

参考文章：

• 【最优化】一文搞懂最速下降法
• 【最优化】为什么最速下降法中迭代方向是锯齿形的？

LM算法

在高斯-牛顿法中引入$\mu$得到LM算法

引入$\mu$的意义

• 高斯牛顿法的缺点
  • H有可能不可逆
    首先，$H=J^{T}J$为半正定对称矩阵（注：形如$A^{T}A$（A为任意矩阵）都是半正定对称矩阵，这个定理是奇异值分解的基础），可以分解为$H=Q\Lambda Q^T$，其中矩阵$Q$的每个列向量为$H$的特征向量，$\Lambda$为对角矩阵，对角元素为对应特征向量的特征值。
    因为$H$为半正定对称矩阵，因此特征值有可能为0，因此不可逆。因为若H可逆，则$H^{-1}=Q{\Lambda }^{-1}{Q}^T$，其中${\Lambda }^{-1}$对角元素为对应特征值$\lambda$的倒数$\frac{1}{\lambda }$，因此若特征值为0，则$H$不可逆。
  • 步长$\Delta x$可能过大，导致发散
    由高斯牛顿法的算法流程可知，其核心是在点$x_k$处利用$l(x)$的泰勒展开，用二次多项式$p_k(x)$（注：实际上$p_k(x)$不是真正泰勒展开的二次多项式，因为矩阵$H$只是黑塞矩阵的近似矩阵）近似$f(x)$。
    $$l(x_k+\Delta x)\approx p_k(x_k+\Delta x)=l(x_k)+(-B^T)\Delta x+\frac{1}{2}{\Delta x}^{T}H\Delta x$$
    然后求二次多项式$p_k(x)$的最小值点$x_{k+1}=x_k+\underset{\Delta x}{\mathrm{argmin}}{p}_k(x_k+\Delta x)$，然后$x_{k+1}$则是这一次迭代的结果。
    因此当$x_k$与$p_k(x)$的最小值点相距很远时，步长$\Delta x$会很大。但泰勒展开一般只在$x_k$的局部区域内能很好的近似原始函数$l(x)$，因此步长太大算法可能会发散（损失值不降反升）。
• 引入(非负数)$\mu$解决高斯牛顿法的缺点
  • 步长$\Delta x$太大的问题
    步长可能太大，那么一个自然的想法就是正则化。因此，修改损失函数为：
    $$p_k(x_k+\Delta x)=l(x_k)+(-B^T)\Delta x+\frac{1}{2}{\Delta x}^{T}H\Delta x+\frac{1}{2}\mu {\Delta x}^T\Delta x$$
    正则化系数$\mu$越大，则越能限制步长$\Delta x$的大小。
    求解$\underset{\Delta x}{\mathrm{argmin}}{p}_k(x_k+\Delta x)$的过程如下：
    (1) 求导：$\omega (\Delta x)=\frac{\partial p_k(x_k+\Delta x)}{\partial \Delta x}=(-B)+H\Delta x+\mu \Delta x=(-B)+(H+\mu I)\Delta x$
    (2) 令$\omega (\Delta x)=0$得：
    $$\underset{\Delta x}{\mathrm{argmin}}{p}_k(x_k+\Delta x)=(H+\mu I{)}^{-1}B$$
  • H不可逆的问题
    由上面可知现在$H$变成了$(H+\mu I)$，只要$\mu >0$，则$(H+\mu I)$一定可逆。因为：
    (1) 首先$(H+\mu I)$是对称矩阵（保证了$(H+\mu I)$有n个正交特征向量，n为$x$的长度,$(H+\mu I)$大小为(n,n)）。
    (2) 其次$(H+\mu I)$与$H$特征向量相同，并且：假设$Hx=\lambda x$，则$(H+\mu I)x=Hx+\mu x=(\lambda +\mu )x$。所以$(H+\mu I)$的特征值为$H$对应特征值加$\mu$。又因为$\lambda \ge 0$，所以若$\mu$大于0，则$(H+\mu I)$的特征值大于0。
    (3)结合(1)(2)得若$\mu >0$，则$(H+\mu I)$为对称正定矩阵，所以$(H+\mu I)$可逆。

如何自动调整$\mu$，LM与高斯牛顿法和最速下降法的关系，算法实现流程

• 如何自动调整$\mu$，LM与高斯牛顿法和最速下降法的关系
  参考文章：Levenberg–Marquardt算法学习
  • 其实信赖域法的本质就是看近似函数（比如这里就是泰勒展开的二阶形式）的损失值下降量$\Delta L_k$和实际损失函数的损失值下降量$\Delta F_k$的比值，如果$\frac{\Delta F_k}{\Delta L_k}$ 约等于1说明近似函数在步长${\Delta }_k$内与实际损失函数很近似，可以保持这个步长或者扩大步长，否则若$\frac{\Delta F_k}{\Delta L_k}$约等于0甚至是负数，就缩小步长。（需要保证$\Delta L_k>0$）
• 算法实现流程
  参考文章：A Brief Description of the Levenberg-Marquardt Algorithm Implemened by levmar

  • 注意： 里面的伪代码中有点错误，g应该是负梯度，也就是 $g:=-J^T{ϵ}_p$。

  • ${ϵ}^T{\Sigma }_y^{-1}ϵ$的作用
    参考文章 A Brief Description of the Levenberg-Marquardt Algorithm Implemened by levmar中提到了这样一段话：
    
    注意，这篇文章里的向量$x$（是本文中的真实值$y=[y_0,y_1,\cdots ,y_{m-1}{]}^T$）与本文的$x$意义不一样。因此下面本文用$y$代替这篇文章的$x$。${ϵ}^T{\Sigma }_y^{-1}ϵ$的作用是消除不同$y_i$有可能有不同量级的影响。
    我们假设${\Sigma }_y$为对角矩阵，也就是$y_i$之间相互独立，则对角值${\sigma }_i$为$y_i$的方差，${\sigma }_i$表示了$y_i$的变化范围（可以理解为量级）。量级越大，那么对应误差${ϵ}_i$的值变化范围也会大，因此在优化过程中会重点优化${ϵ}_i$。因此我们要避免这种由量级导致的误差过大或过小。因此算法以${ϵ}^T{\Sigma }_y^{-1}ϵ$作为损失值，代替${ϵ}^Tϵ$。
    若${\Sigma }_y$不是对角矩阵，但因为协方差矩阵和协方差矩阵的逆都是正定对称矩阵（只要没有互相关变量）。因此${\Sigma }_y^{-1}$可分解为$Q{\Lambda }^{-1}{Q}^T$。而${ϵ}^T{\Sigma }_y^{-1}ϵ$=$({ϵ}^{T}Q){\Lambda }^{-1}(Q^Tϵ)$，把$(Q^Tϵ)$当成新的随机变量。而$(Q^Tϵ)$的协方差矩阵为$\Lambda$，因此也实现了消除量级影响。

  • $\mu$初始值
    在参考文章A Brief Description of the Levenberg-Marquardt Algorithm Implemened by levmar的伪代码里，$\mu$的初始值如下图所示。其中$ma{x}_{i=1,\cdots ,m}(H_{ii})$(参考文章的$A$等于本文的$H$)。这其实是为了让$\mu$和$H$对角线上的值的数量级一致。因为我们有$H+\mu I$，因此$\mu$是加到$H$的对角线上的。
    

  • 参考文章建议的初始值 ：
    

非线性最小二乘法资料

• 《Methods for non-linear least squares problems》