非線性最小二乘法優(yōu)化

高斯-牛頓法

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

接下來(lái)本文使用的數(shù)學(xué)符號(hào)意義與上圖一樣。其中$x$是需要求解的參數(shù)，$f(x)$是一個(gè)殘差向量。比如有一個(gè)優(yōu)化問(wèn)題，$y=asin(wt+b)+c$，給出m個(gè)數(shù)據(jù)$(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$(向量二范數(shù))就是最小二乘法的損失值。
設(shè)損失函數(shù)$l(x)=\frac{1}{2}||f(x)|{|}^2$。
另外，$J(x)$為$f(x)$的雅可比矩陣，假設(shè)$x$的長(zhǎng)度為n，$f(x)$長(zhǎng)度為m，則$J(X)$矩陣大小為(m,n)。
$H=J^{T}J$為$f(x)$的黑塞矩陣的近似矩陣。$B=?J^{T}f(x)$為損失函數(shù)$l(x)$（$\frac{1}{2}$只是為了求導(dǎo)后約掉$||f(x)|{|}^2$的指數(shù)2）的負(fù)梯度$?\frac{?l(x)}{?x}$。

最速下降法

參考文章：

• 【最優(yōu)化】一文搞懂最速下降法
• 【最優(yōu)化】為什么最速下降法中迭代方向是鋸齒形的？

LM算法

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

引入$\mu$的意義

• 高斯牛頓法的缺點(diǎn)
  • H有可能不可逆
    首先，$H=J^{T}J$為半正定對(duì)稱矩陣（注：形如$A^{T}A$（A為任意矩陣）都是半正定對(duì)稱矩陣，這個(gè)定理是奇異值分解的基礎(chǔ)），可以分解為$H=Q\Lambda Q^T$，其中矩陣$Q$的每個(gè)列向量為$H$的特征向量，$\Lambda$為對(duì)角矩陣，對(duì)角元素為對(duì)應(yīng)特征向量的特征值。
    因?yàn)?span id="ozvdkddzhkzd" class="katex--inline">$H$為半正定對(duì)稱矩陣，因此特征值有可能為0，因此不可逆。因?yàn)槿鬑可逆，則$H^{?1}=Q{\Lambda }^{?1}{Q}^T$，其中${\Lambda }^{?1}$對(duì)角元素為對(duì)應(yīng)特征值$\lambda$的倒數(shù)$\frac{1}{\lambda }$，因此若特征值為0，則$H$不可逆。
  • 步長(zhǎng)$\Delta x$可能過(guò)大，導(dǎo)致發(fā)散
    由高斯牛頓法的算法流程可知，其核心是在點(diǎn)$x_k$處利用$l(x)$的泰勒展開(kāi)，用二次多項(xiàng)式$p_k(x)$（注：實(shí)際上$p_k(x)$不是真正泰勒展開(kāi)的二次多項(xiàng)式，因?yàn)榫仃?span id="ozvdkddzhkzd" class="katex--inline">$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$$
    然后求二次多項(xiàng)式$p_k(x)$的最小值點(diǎn)$x_{k+1}=x_k+\underset{\Delta x}{\mathrm{argmin}}{p}_k(x_k+\Delta x)$，然后$x_{k+1}$則是這一次迭代的結(jié)果。
    因此當(dāng)$x_k$與$p_k(x)$的最小值點(diǎn)相距很遠(yuǎn)時(shí)，步長(zhǎng)$\Delta x$會(huì)很大。但泰勒展開(kāi)一般只在$x_k$的局部區(qū)域內(nèi)能很好的近似原始函數(shù)$l(x)$，因此步長(zhǎng)太大算法可能會(huì)發(fā)散（損失值不降反升）。
• 引入(非負(fù)數(shù))$\mu$解決高斯牛頓法的缺點(diǎn)
  • 步長(zhǎng)$\Delta x$太大的問(wèn)題
    步長(zhǎng)可能太大，那么一個(gè)自然的想法就是正則化。因此，修改損失函數(shù)為：
    $$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$$
    正則化系數(shù)$\mu$越大，則越能限制步長(zhǎng)$\Delta x$的大小。
    求解$\underset{\Delta x}{\mathrm{argmin}}{p}_k(x_k+\Delta x)$的過(guò)程如下：
    (1) 求導(dǎo)：$\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不可逆的問(wèn)題
    由上面可知現(xiàn)在$H$變成了$(H+\mu I)$，只要$\mu >0$，則$(H+\mu I)$一定可逆。因?yàn)?#xff1a;
    (1) 首先$(H+\mu I)$是對(duì)稱矩陣（保證了$(H+\mu I)$有n個(gè)正交特征向量，n為$x$的長(zhǎng)度,$(H+\mu I)$大小為(n,n)）。
    (2) 其次$(H+\mu I)$與$H$特征向量相同，并且：假設(shè)$Hx=\lambda x$，則$(H+\mu I)x=Hx+\mu x=(\lambda +\mu )x$。所以$(H+\mu I)$的特征值為$H$對(duì)應(yīng)特征值加$\mu$。又因?yàn)?span id="ozvdkddzhkzd" class="katex--inline">$\lambda \ge 0$，所以若$\mu$大于0，則$(H+\mu I)$的特征值大于0。
    (3)結(jié)合(1)(2)得若$\mu >0$，則$(H+\mu I)$為對(duì)稱正定矩陣，所以$(H+\mu I)$可逆。

如何自動(dòng)調(diào)整$\mu$，LM與高斯牛頓法和最速下降法的關(guān)系，算法實(shí)現(xiàn)流程

• 如何自動(dòng)調(diào)整$\mu$，LM與高斯牛頓法和最速下降法的關(guān)系
  參考文章：Levenberg–Marquardt算法學(xué)習(xí)
  • 其實(shí)信賴域法的本質(zhì)就是看近似函數(shù)（比如這里就是泰勒展開(kāi)的二階形式）的損失值下降量$\Delta L_k$和實(shí)際損失函數(shù)的損失值下降量$\Delta F_k$的比值，如果$\frac{\Delta F_k}{\Delta L_k}$ 約等于1說(shuō)明近似函數(shù)在步長(zhǎng)${\Delta }_k$內(nèi)與實(shí)際損失函數(shù)很近似，可以保持這個(gè)步長(zhǎng)或者擴(kuò)大步長(zhǎng)，否則若$\frac{\Delta F_k}{\Delta L_k}$約等于0甚至是負(fù)數(shù)，就縮小步長(zhǎng)。（需要保證$\Delta L_k>0$）
• 算法實(shí)現(xiàn)流程
  參考文章：A Brief Description of the Levenberg-Marquardt Algorithm Implemened by levmar

  • 注意： 里面的偽代碼中有點(diǎn)錯(cuò)誤，g應(yīng)該是負(fù)梯度，也就是 $g:=?J^T{?}_p$。

  • ${?}^T{\Sigma }_y^{?1}?$的作用
    參考文章 A Brief Description of the Levenberg-Marquardt Algorithm Implemened by levmar中提到了這樣一段話：
    
    注意，這篇文章里的向量$x$（是本文中的真實(shí)值$y=[y_0,y_1,?,y_{m?1}{]}^T$）與本文的$x$意義不一樣。因此下面本文用$y$代替這篇文章的$x$。${?}^T{\Sigma }_y^{?1}?$的作用是消除不同$y_i$有可能有不同量級(jí)的影響。
    我們假設(shè)${\Sigma }_y$為對(duì)角矩陣，也就是$y_i$之間相互獨(dú)立，則對(duì)角值${\sigma }_i$為$y_i$的方差，${\sigma }_i$表示了$y_i$的變化范圍（可以理解為量級(jí)）。量級(jí)越大，那么對(duì)應(yīng)誤差${?}_i$的值變化范圍也會(huì)大，因此在優(yōu)化過(guò)程中會(huì)重點(diǎn)優(yōu)化${?}_i$。因此我們要避免這種由量級(jí)導(dǎo)致的誤差過(guò)大或過(guò)小。因此算法以${?}^T{\Sigma }_y^{?1}?$作為損失值，代替${?}^T?$。
    若${\Sigma }_y$不是對(duì)角矩陣，但因?yàn)閰f(xié)方差矩陣和協(xié)方差矩陣的逆都是正定對(duì)稱矩陣（只要沒(méi)有互相關(guān)變量）。因此${\Sigma }_y^{?1}$可分解為$Q{\Lambda }^{?1}{Q}^T$。而${?}^T{\Sigma }_y^{?1}?$=$({?}^{T}Q){\Lambda }^{?1}(Q^T?)$，把$(Q^T?)$當(dāng)成新的隨機(jī)變量。而$(Q^T?)$的協(xié)方差矩陣為$\Lambda$，因此也實(shí)現(xiàn)了消除量級(jí)影響。

  • $\mu$初始值
    在參考文章A Brief Description of the Levenberg-Marquardt Algorithm Implemened by levmar的偽代碼里，$\mu$的初始值如下圖所示。其中$ma{x}_{i=1,?,m}(H_{ii})$(參考文章的$A$等于本文的$H$)。這其實(shí)是為了讓$\mu$和$H$對(duì)角線上的值的數(shù)量級(jí)一致。因?yàn)槲覀冇?span id="ozvdkddzhkzd" class="katex--inline">$H+\mu I$，因此$\mu$是加到$H$的對(duì)角線上的。
    

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

非線性最小二乘法資料

• 《Methods for non-linear least squares problems》

總結(jié)

以上是生活随笔為你收集整理的LM算法原理的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。