非線性最小二乘法優化

高斯-牛頓法

參考文章：[優化] Gauss-Newton非線性最小二乘算法
算法流程如下圖（來自參考文章）所示：

接下來本文使用的數學符號意義與上圖一樣。其中$x$是需要求解的參數，$f(x)$是一個殘差向量。比如有一個優化問題，$y=asin(wt+b)+c$，給出m個數據$(t_i,y_i)(i=0,1,?,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),?,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{?l(x)}{?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{?p_k(x_k+\Delta x)}{?\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,?,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,?,m}(H_{ii})$(參考文章的$A$等于本文的$H$)。這其實是為了讓$\mu$和$H$對角線上的值的數量級一致。因為我們有$H+\mu I$，因此$\mu$是加到$H$的對角線上的。
    

  • 參考文章建議的初始值 ：
    

非線性最小二乘法資料

• 《Methods for non-linear least squares problems》

總結

以上是生活随笔為你收集整理的LM算法原理的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。