非线性最小二乘法原理

一、一些基础概念[1]

1.最小二乘法定义
找到一个$x^{\ast }$，使得下面的式子到达局部最小值，
$$F(x)=\frac{1}{2}\sum _{i=1}^m{(f_i(x))}^2\text{\hspace{1em}(1-1)}$$这里的$f_i(x)$表示变量为$x$的函数，其中自变量$x$维度为$n$，函数个数$i=1,2,\cdots ,m$，最小二乘问题要求$m\ge n$，即方程的个数$\ge$未知数的个数。

一般而言，$f_i(x)=y_i(x)-y_i$，即$f_i(x)$由两个部分组成：含有自变量$x$的函数$y(x)$(估计值)+不含$x$的常数$y$(测量值)。而$f_i(x)$叫做残差项，式(1-1)的核心思想就是使残差的平方和最小。

2.全局最小值
给定目标函数$F$，找到一个$x$使得下面的式子最小，
$$x^{+}=\arg \underset{x}{\min }F(x)\text{\hspace{1em}(1-2)}$$
3.局部最小值
给定目标函数$F$，找到一个$x^{\ast }$使得,
$$F(x^{\ast })\le F(x)for\Vert x-x^{\ast }\Vert <\delta \text{\hspace{1em}(1-3)}$$这里的$\delta$是一个非常小的无限趋近于零的值。

4.局部最小值的必要条件
若$x^{\ast }$是一个局部最小值点，则有，
$$g^{\ast }=F^{{}^{\prime }}(x^{\ast })=0\text{\hspace{1em}(1-4)}$$这里的$g={\left[\begin{array}{lcc}\frac{\partial F(x)}{\partial x_1}& \cdots & \frac{\partial F(x)}{\partial x_n}\end{array}\right]}^T$，表示$F(x)$的梯度。

5.驻点
若有$x_s$使得，
$$g_s=F^{{}^{\prime }}(x_s)=0\text{\hspace{1em}(1-5)}$$则$x_s$为函数$F(x)$的驻点。
驻点的含义是在$x_s$处，$F^{{}^{\prime }}(x_s)=0$，$F(x)$的函数值在该点处既没有增大的趋势，也没有减小的趋势，那么驻点$x_s$有会有以下三种情况：

6.局部最小值点的充要条件
若$x_s$为驻点且$F^{{}^{\prime \prime }}(x_s)$正定，则$x_s$为局部最小值点。
这里的$F^{{}^{\prime \prime }}(x_s)=H=\left[\frac{{\partial }^{2}F(x)}{\partial x_i\partial x_j}\right]$,$H$为Hessian矩阵，表示函数$F(x)$的二阶导。

7.最小二乘问题的两种解的形式

<1>闭式解(closed-form solution)
能够直接给出优化问题解的解析形式，如下面要说的线性最小二乘。

<2>迭代解(iterative solution)
没有解析形式，只能从一个初始值$x_0$开始，一步步迭代求解得到$x_1,x_2,\cdots$。对于最优化问题，需要一步步迭代，并最终收敛到局部最小值点$x^{\ast }$，且迭代要满足以下条件：
$$F(x_{k+1})<F(x_k)\text{\hspace{1em}(1-6)}$$

二、下降方法

1. 下降方法是迭代解法的一种，其含义是：沿着一个下降的方向以一定的步长下降，一般这种下降方向就是朝着局部极小值点收敛的方向。所以下降法就要求解两个参数：

• 找到一个下降方向$h_d$;
• 确定步长$\alpha$使$F(x)$合适的减小。

2.下降法的步骤如下：
初值设置： 迭代次数$k=0$，迭代初值$x=x_0$，下降方向的找到的标志位$found=false$;
寻找步骤： 寻找下降方向$h_d=search\_direction(x)$和下降步长$\alpha =step\_length(x,h_d)$。
更新步骤： 更新解$x=x+\alpha h_d$和迭代次数$k=k+1$。
结束循环标志： 迭代次数达到设定值($k<k_{max}$)+没找到下降方向($found=true$,此时的$x$为驻点)。

3.寻找下降方向的原理
考虑目标函数$F(x)$在$x$处的泰勒展开，
$$\begin{array}{l}F(x+\alpha h)=F(x)+\alpha h^T{F}^{{}^{\prime }}(x)+O({\alpha }^2)\\ \\ \approx F(x)+\alpha h^T{F}^{{}^{\prime }}(x)\end{array}\text{\hspace{1em}(2-1)}$$这里式(2-1)里的$\alpha$要足够小，式(2-1)才能成立。

从式(2-1)里可以总结出：如果$h^T{F}^{{}^{\prime }}(x)<0$，则$h$一定是下降方向。这里很好证明，若$h^T{F}^{{}^{\prime }}(x)<0$，则式(2-1)里$F(x+\alpha h)\approx F(x)+负数$，会有$F(x+\alpha h)<F(x)$，符合式(1-6)里的定义。

根据这一原理，有两种下降方向：最速下降方向($h_{sd}$，Steepest Descent)和牛顿方向($h_n$)。

4.最速下降方向——最速下降法
含义：在点$x$处下降速度最快的方向。可以证明，最速下降方向为函数$F(x)$在$x$处的负梯度。
证明如下：
取式(2-1)中的最后一项，
$$h^T{F}^{{}^{\prime }}(x)=\Vert h\Vert \Vert F^{{}^{\prime }}(x)\Vert \cos \theta =\Vert F^{{}^{\prime }}(x)\Vert \cos \theta \text{\hspace{1em}(2-2)}$$其中$\theta$表示的是向量$h$和向量$F^{{}^{\prime }}(x)$的夹角。
那么，当式(2-2)中的$\theta =\pi$时，式(2-2)达到最小值。当$\alpha$取定时，$\Vert F(x+\alpha h)-F(x)\Vert$最大，即此时下降最快，且方向为负梯度方向，

$$h_{sd}=-\frac{F^{{}^{\prime }}(x)}{\Vert F^{{}^{\prime }}(x)\Vert }\text{\hspace{1em}(2-3)}$$

5.牛顿下降方向——牛顿法
牛顿方向考虑$F^{{}^{\prime }}(x)$在$x$处的泰勒展式，
$$\begin{array}{l}{F}^{{}^{\prime }}(x+h)=F^{{}^{\prime }}(x)+F^{{}^{\prime \prime }}(x)h+O(\Vert h\Vert )\\ \\ \approx F^{{}^{\prime }}(x)+F^{{}^{\prime \prime }}(x)h\end{array}\text{\hspace{1em}(2-4)}$$
因为局部最小值也是驻点的一种，所以在$x^{\ast }$附近，有$F^{{}^{\prime }}(x^{\ast })=0$。因此，令式(2-4)等于零，就可以得到牛顿下降方向，
$$h_n=-(F^{{}^{\prime \prime }}(x){)}^{-1}{F}^{{}^{\prime }}(x)\text{\hspace{1em}(2-5)}$$可以证明，式(2-5)所示的方向是一个下降方向。因为二阶导矩阵$H=F^{{}^{\prime \prime }}(x)$为正定矩阵，所以有$u^{T}Hu>0$，同理，对于一个非零向量$h^T$，有，
$$0<h_n^{T}H{h}_n=-h_n^T{F}^{{}^{\prime }}(x)\text{\hspace{1em}(2-6)}$$符合2.3中下降方向的判定条件。

6.最速下降法vs牛顿法

最速下降法使用的是最速下降方向，始终朝着下降速度最快的方向的下降，所以在距离最小值点$x^{\ast }$较远的地方，表现较好。并且可以证明，最速下降法的收敛速度为线性收敛，比较慢。
牛顿法使用了牛顿方向，从2.5的推算过程可以看出，它是在局部最小值附近，通过$F^{{}^{\prime }}(x^{\ast })=0$推出来的，所以其在局部最小值附近表现较好。并且也可以证明，它的收敛速度为平方收敛，速度快。但是用牛顿法求解最小值时，要求矩阵$H$正定，若$H$为非正定矩阵，则无法使用牛顿法。

所以，有人综合最速下降法和牛顿法，提出了hybrid方法。该方法描述如下：如果$H$正定，则取牛顿方向$h=h_n$；如果$H$非正定，则取最速下降方向$h=h_{sd}$。但在复杂的问题里，hybrid方法也很难被应用，因为在很多问题中，矩阵$H$一般计算不出来，所以这就要用到Gaussian-Newton法或者Levenberg-Marquardt法。

7.确定步长的方法1——Line Search[2]
当已知当前迭代变量$x$和下降方向$h$时，此时，需要做的是确定一个步长$\alpha$，使目标函数得到一个合适的下降。Line Search的做法就是把$x$和$h$当做常数，把$\alpha$当做待求解的变量，因此，有，
$$\phi (\alpha )=F(x+\alpha h)\text{\hspace{1em}(2-7)}$$对式(2-7)求导，有，

$${\phi }^{{}^{\prime }}(\alpha )=F(x)+\alpha h^T{F}^{{}^{\prime }}(x)\text{\hspace{1em}(2-8)}$$令式(2-8)中的$\alpha =0$，因为$h$为一个下降方向，所以有${\phi }^{{}^{\prime }}(0)=h^T{F}^{{}^{\prime }}(x)<0$。所以$\phi (\alpha )$的函数曲线大致如下，

对于上图，对$\alpha$的选择会产生三种情况：

• cond1: $\alpha$很小使得目标函数$F(x)$下降的幅度特别小，此时应该增大$\alpha$；
• cond2: $\alpha$很大，这样会造成$\phi (\alpha )\ge \phi (0)$，这不符合下降法的定义，即没有朝着下降的方向前进；
• cond3:$\alpha$非常接近函数的$\phi (\alpha )$的局部最小值点，此时合适可取。

精确的Line Search 方法如下，
$${\alpha }_e=\arg \underset{\alpha >0}{\min }F(x+\alpha h)\text{\hspace{1em}(2-9)}$$具体细节参考文献[2]里的Section2.5-2.6。
由于精确的Line Search方法非常耗时，而且当搜索方向$h$与方向$x-x\ast$相差很大时，就没有必要计算(2-9)里函数$\phi (\alpha )$真正的最小值。所以，对于Line Search就没有必要很精确，从而提出了soft line search。这种方法就是接受不满足cond1和cond2的$\alpha$值，具体细节参考文献[2]Section2.5。

8.确定步长的方法2——信赖域法(Trust Region)和阻尼法(Damped)
2.7中的Line Research方法要分两步走：先确定下降方向$h$后，再计算$\alpha$的大小。而这里的信赖域法和阻尼法则不用分步，它可以同时得到下降方向和步长。

它们的原理是在$x$的邻域内，用一个模型函数$L$来近似目标函数$F$：
$$F(x+h)\approx L(h)=F(x)+h^{T}c+\frac{1}{2}{h}^{T}B{h}^T\text{\hspace{1em}(2-10)}$$这里要注意式(2-10)的$L(h)$里的三项不是$F(x+h)$泰勒展式的前三项，因为，
$$F(x+h)\approx F(x)+h^{T}g+\frac{1}{2}{h}^{T}H{h}^T\text{\hspace{1em}(2-11)}$$这里的$c\ne g,B\ne H$。所以$L(h)$是$F(x)$在模型意义上的近似，它只是近似等于$F(x)$取泰勒展式前三项的结果。当然，这种近似只有在$h$很小时才成立。

那么，信赖域法的策略是：假设在当前变量$x$的一定范围内，即以$x$为中心，以$\Delta$为半径的范围内，模型$L(h)$对$F$的近似足够精确，那么此时的下降方向为，

$$h=h_{tr}=\arg \underset{\Vert h\Vert \le \Delta }{\min }L(h)\text{\hspace{1em}(2-12)}$$而阻尼法的策略是：给$L(h)$一个惩罚项$\frac{1}{2}\mu h^{T}h$(也叫正则项，机器学习里的$l_2$正则器就是这个)，以防止$h$过大造成的近似不精确，
h = h d m = arg ⁡ min ⁡ h { L ( h ) + 1 2 μ h T h } (2-13) h = {h_{dm}} = \arg \mathop {\min }\limits_h \{ L(h) + \frac{1}{2}\mu {h^T}h\}\tag{2-13} h=hdm​=arghmin​{L(h)+21​μhTh}(2-13)此时的步长选取策略为：若$h$是一个下降方向，直接取$\alpha =1$；若不是，取$\alpha =0$，本次迭代无效。

式(2-12)和式(2-13)中的$h$只能保证模型函数$L(h)$是下降的，由于$L(h)$只是$F(x)$的近似，所以，它的下降并不能保证函数$F$是不是真的在下降，所以，这里就提出了gain ratio的概念，这个参数是用来衡量模型下降品质的参数，定义如下：
$$\rho =\frac{F(x)-F(x+h)}{L(0)-L(h)}\text{\hspace{1em}(2-14)}$$显然，该式的分母始终为正，因为$h$会使$L(h)$下降。所以，当$\rho$为正时，表示分子也为正，此时$F(x)>F(x+h)$，符合式(1-6)，表示目标函数$F$也在下降，本次迭代有效；但当$\rho$为负时，则表示目标函数$F$在增大，不符合要求。

根据$\rho$值的大小，可以确定信赖域法的$\Delta$和阻尼法的$\mu$。

• 在信赖域法中，根据$\rho$值的不同，可以适当的调整$\Delta$的大小，使$F$获得一个合适的下降。通常，下面的这个策略被广泛使用，
  i f      ρ < 0.25            Δ = Δ / 2 e l s e      i f    ρ > 0.75            Δ = max ⁡ { Δ , 3 ∗ ∥ h ∥ }          (2-15) \begin{array}{l} if\;\;\rho < 0.25\\ \;\;\;\;\;\Delta = \Delta /2\\ else\;\;if\;\rho > 0.75\\ \;\;\;\;\;\Delta = \max \{ \Delta ,3*\left\| h \right\|\} \\ \;\;\;\; \end{array}\tag{2-15} ifρ<0.25Δ=Δ/2elseifρ>0.75Δ=max{Δ,3∗∥h∥}​(2-15)式(2-15)的含义是：当$\rho <0.25$时，说明$F$的下降或上升大小比$L$下降的四分之一还要小，表明该情况下的$L(h)$对$F$而言不是一个很好的近似，应该减小$\Delta$;当$\rho >0.75$时，说明$L$下降的幅度比较小，此时应该增大$\Delta$。
  对于式(2-15)里的阈值设定(0.25&0.75)以及除数因子(2)和乘法因子(3)的设定，信赖域法允许有较小的变动，但是选取这些参数的前提是：选取的参数最好不要使$\Delta$有震荡。

• 阻尼法中的$\mu$的选取原则与信赖域法中的$\Delta$选取同理。当$\rho$值很小时，说明$L$不能很好的近似$F$，所以要增大惩罚项，即增大$\mu$的值；反之，当$\rho$值很大时，说明$F$下降幅度小，要减小惩罚项，即减小$\mu$的值。基于此，马夸尔特提出了一种类似于式(2-15)的策略：
  $$\begin{array}{l}if\rho <0.25\\ \mu =\mu \ast 2\\ elseif\rho >0.75\\ \mu =\mu /3\\ \end{array}\text{\hspace{1em}(2-16)}$$对于式(2-16)里的阈值设定(0.25&0.75)以及除数因子(3)和乘法因子(2)的设定，阻尼法允许有较小的变动，但是选取这些参数的前提是：选取的参数最好不要使$\Delta$有震荡。
  除了式(2-16)，Nielsen提出了另外一种经验公式：
  i f      ρ > 0            μ = μ ∗ max ⁡ { 1 3 , 1 − ( 2 ρ − 1 ) 3 } ;      v = 2 e l s e                 μ = μ ∗ v ;      v = 2 ∗ v          (2-17) \begin{array}{l} if\;\;\rho > 0\\ \;\;\;\;\;\mu = \mu *\max \{ \frac{1}{3},1 - {(2\rho - 1)^3}\} ;\;\;v = 2\\ else\;\;\\ \;\;\;\;\;\mu = \mu *v;\;\;v = 2*v\\ \;\;\;\; \end{array}\tag{2-17} ifρ>0μ=μ∗max{31​,1−(2ρ−1)3};v=2elseμ=μ∗v;v=2∗v​(2-17)

9.阻尼法的进一步探讨——阻尼牛顿法
对于式（2-13），要想求其最小值，首先考虑$h$必须是$\phi (h)=L(h)+\frac{1}{2}\mu h^{T}h$的驻点，所以有，
$${\phi }^{{}^{\prime }}(h)=L^{{}^{\prime }}(h)+\mu h=0\text{\hspace{1em}(2-18)}$$再根据式（2-10）中$L(h)$的定义，有，
$$(B+\mu I)h_{dm}=-c\text{\hspace{1em}(2-19)}$$这里的$I$是单位阵，那么当$\mu$足够大时，$B+\mu I$肯定为正定矩阵，根据概念1.6，$h_{dm}$必为$L$的一个局部最小值点。

那么，再特殊化一点，取$c=F^{{}^{\prime }}(x),B=F^{{}^{\prime \prime }}(x)$，则式（2-19）可改写如下，
$$(F^{{}^{\prime \prime }}(x)+\mu I)h_{dn}=-F^{{}^{\prime }}(x)\text{\hspace{1em}(2-20)}$$由于该式与牛顿法只差一个含有$\mu$的阻尼项，所以，该式表示的方法叫做阻尼牛顿法(damped Newton)。该式可以从两个方面来看：

• 当$\mu$很小时，$(F^{{}^{\prime \prime }}(x)+\mu I)\approx F^{{}^{\prime \prime }}(x)$，则$h_{dn}=-(F^{{}^{\prime \prime }}(x){)}^{-1}{F}^{{}^{\prime }}(x)$，此时退化为牛顿法；
• 当$\mu$很大时，$(F^{{}^{\prime \prime }}(x)+\mu I)\approx \mu$，则$h_{dn}=-\frac{1}{\mu }{F}^{{}^{\prime }}(x)$，此时退化为了最速下降法。

10.信赖域法与阻尼法的一致性探讨
式(2-12)表示的信赖域法实际上等价于一个带约束的优化问题。
$$\underset{h}{\min }L(h)s.t.h^{T}h\le {\Delta }^2\text{\hspace{1em}(2-21)}$$对于学过机器学习里正则化理论的童鞋来说，这个式子是不是很眼熟。如果把目标函数$F(x)$当前一个未知的机器学习里所要估计的真实函数(实际上是求不出来的)，那么机器学习的目的就是用已有的数据取求出一个模型函数$L(h)$来尽可能的贴近$F(x)$，为了防止过拟合(过拟合的含义是：模型函数的表现很差，不能很好的近似$F(x)$)，通常可以采用$l_2$正则器来建立约束，尽可能的防止过拟合。而式(2-21)就是一个对$L(h)$进行正则化的过程。它的目的是为了让$L(h)$可以最大程度的近似$F(x)$，那么这样的情况下，对$F(x)$求最小值，就可以转换为对$L(h)$求最小值。而这里的$\Delta$影响的就是模型函数近似的品质好坏，而品质的好坏就是根据$\rho$判断出来的，所以$\Delta$的值是根据$\rho$的值不断变化的，即(2-15)里的策略。

式(2-21)其实是一个带约束的最优化问题，该类问题可以通过拉格朗日乘子法变成无约束的优化问题：
$$\underset{h}{\min }L(h)+\lambda (h^{T}h-{\Delta }^2)\text{\hspace{1em}(2-22)}$$显然，对式(2-22)求导然后等于零，可以求出其最小值。那么，式(2-22)对$h$求导有，

$$L^{{}^{\prime }}(h)+2\lambda h\text{\hspace{1em}(2-23)}$$再将上式反推回去，有，
$$L(h)+\lambda h^{T}h\text{\hspace{1em}(2-24)}$$式(2-22)和式(2-24)的导数式都是式(2-23)，所以可以得到式(2-22)中的最小值点其实与$\Delta$无关，只与拉格朗日乘子$\lambda$有关，若令$\lambda =\frac{1}{2}\mu$，则式(2-22)就变成了式(2-13)所示的阻尼法的求解形式。所以，信赖域法和阻尼法其实是一致的，它其实是一个最优化问题的两种表现形式。

三、线性最小二乘及其解法

1. 若概念1.1中的函数$f_i(x)$为关于$x$的线性函数，即$f_i(x)=a_{i}x-b_i$，其中估计值$y_i(x)=a_{i}x$(线性函数)，则概念1所示的问题为线性最小二乘问题。
对于线性最小二乘问题，可以直接闭式求解（closed-form），过程如下：
根据式(1-1)，写成矩阵形式，有，
$$\begin{array}{l}F(x)=(Ax-b{)}^T(Ax-b)\\ \\ =x^T{A}^{T}Ax-2x^T{A}^{T}b+b^{T}b\end{array}\text{\hspace{1em}(3-1)}$$这里的$A={\left[\begin{array}{lcc}{a}_1& \cdots & a_m\end{array}\right]}^T$，$b={\left[\begin{array}{lcc}-b_1& \cdots & -b_m\end{array}\right]}^T$。
对式(3-1)求导，有，
$$\frac{\partial F(x)}{\partial x}=2A^{T}Ax-2A^{T}b=0\text{\hspace{1em}(3-2)}$$则线性最小二乘的闭式解形式如下，
$$x=(A^{T}A{)}^{-1}{A}^{T}b\text{\hspace{1em}(3-3)}$$对于线性最小二乘的解法，还可以用QR分解来做。Matlab里的$A\backslash b$，就是用QR分解的方式来解的，它的求解精度要比式(1-8)高。

2. 举个用线性最小二乘拟合直线例子，若有以下六个点，坐标为：$(0,0),(2,0),(1,-1),(0,2),(-1,1),(1,1)$。对以上的点进行直线拟合，直线方程取$y=k_{1}x+k_2$。

在这个问题里，自变量为${\left[\begin{array}{lc}{k}_1& k_2\end{array}\right]}^T$。按概念1，估计值$y(x)=k_{1}x+k_2$，测量值$y=y$，且$A$和$b$分别为：

$$A={\left[\begin{array}{l}0,2,1,0,-1,1\\ 1,1,1,1,1,1\end{array}\right]}^T$$$$b={\left[0,0,-1,2,1,1\right]}^T$$matlab里的结果为$k_1{k}_2=[-0.4545,0.7273{]}^T$。
程序如下：

x=[0 2 1 0 -1 1]';
y=[0 0 -1 2 1 1]'
A=[x ones(6,1)];
b=y;
k1k2=A\b;

四、非线性最小二乘及其解法

1. 若概念1.1中的函数$f_i(x)$为关于$x$的非线性函数，即$f_i(x)=y_i(x)-y_i$，其中估计值$y_i(x)$为非线性函数(例如，$y_i(x)=cos(x)$)，测量值$y$为常数矩阵。则概念1.1所示的问题为非线性最小二乘问题。式(1-1)写成矩阵形式为，
$$F(x)=\frac{1}{2}\sum _{i=1}^m{(f_i(x))}^2=\frac{1}{2}f(x{)}^{T}f(x)\text{\hspace{1em}(4-1)}$$此时，对$f(x)$进行泰勒展开有，
$$f(x+h)=f(x)+J(x)h+O({\Vert h\Vert }^2)\text{\hspace{1em}(4-2)}$$上式中的$J(x)$为雅克比矩阵，且$J\in R^{m\times n}$。它的第$i$行和第$j$列表示如下，
$$J(x{)}_{ij}=\frac{\partial f_i(x)}{\partial x_j}\text{\hspace{1em}(4-3)}$$这表示第$i$个函数$f_i(x)$对第$j$维变量$x_j$求导。
注意到式(4-1)，若目标函数$F(x)$对$x_j$求导，则有，
$$\frac{\partial F(x)}{\partial x_j}=\sum _{i=1}^m{f}_i(x)\frac{\partial f_i(x)}{\partial x_j}\text{\hspace{1em}(4-4)}$$那么，$F(x)$的梯度可表示如下，
$$g=F^{{}^{\prime }}(x)=J(x{)}^{T}f(x)\text{\hspace{1em}(4-5)}$$同理，考虑$F(x)$的二阶Hessian矩阵，有，
$$\frac{{\partial }^{2}F(x)}{\partial x_j\partial x_k}=\sum _{i=1}^m(\frac{\partial f_i(x)}{\partial x_j}\frac{\partial f_i(x)}{\partial x_k}+f_i(x)\frac{{\partial }^2{f}_i(x)}{\partial x_j\partial x_k})\text{\hspace{1em}(4-6)}$$用雅克比矩阵表示上式为，
$$H=F^{{}^{\prime \prime }}(x)=J(x{)}^{T}J(x)+\sum _{i=1}^m{f}_i(x)f_i^{{}^{\prime \prime }}(x)\text{\hspace{1em}(4-7)}$$

非线性最小二乘一般是没有解析解的，所以只能迭代求解。而一般的优化方法也能用于非线性最小二乘的求解，但是效率比较慢。而高效的方法有高斯-牛顿法(Gaussian-Newton Method,G-N Method)，列文伯格-马夸尔特(Levenberg-Marquardt Method,L-M Method),Dog Leg方法等等。

2.高斯牛顿法(G-N Method)
前面说到牛顿法的原理是将原目标函数$F(x)$在$x$处进行一阶泰勒展开后得出来的；而高斯牛顿法则是对$f$进行一阶泰勒展开，
$$f(x+h)\approx l(h)=f(x)+Jh\text{\hspace{1em}(4-8)}$$根据式(4-1)，则有下式，
$$\begin{array}{l}F(x+h)\approx L(h)=\frac{1}{2}l(h{)}^{T}l(h)\\ \\ =\frac{1}{2}{f}^{T}f+h^T{J}^{T}f+\frac{1}{2}{h}^T{J}^{T}Jh\\ \\ =F(x)+h^T{J}^{T}f+\frac{1}{2}{h}^T{J}^{T}Jh\end{array}\text{\hspace{1em}(4-9)}$$所以，$L(h)$是对原函数$F(x)$的近似。到这一步就可以看出来了，这里的$L(h)$实际上就是2.8里所描述的信赖域法的模型函数。其实，就是对2.8节里的信赖域法取：$c=J^f，B=J^{T}J$的结果。

对$L(h)$进行求导，有，
$$L^{{}^{\prime }}(h)=J^{T}f+J^{T}Jh,L^{{}^{\prime \prime }}(h)=J^{T}J\text{\hspace{1em}(4-10)}$$根据牛顿法的推导，令$L^{{}^{\prime }}(h)=0$，有，
$$J^{T}J{h}_{gn}=-J^{T}f\text{\hspace{1em}(4-11)}$$式(4-11)中的$h_{gn}$必为一个下降方向。因为若(4-11)有解，则$J^{T}J$必正定，所以有，
$${h_{gn}}^T{F}^{{}^{\prime }}(x)={h_{gn}}^T{J}^{T}f=-{h_{gn}}^T{J}^{T}J{h}_{gn}<0\text{\hspace{1em}(4-12)}$$符合下降方向的判定条件。

高斯牛顿法与牛顿法相似，在接近局部最小值点$x^{\ast }$时收敛性较好。但是，若$J^{T}J$非正定，即$J^{T}J$为奇异矩阵，则式(4-11)无解，此时高斯牛顿法就失效了。

高斯牛顿法的实质： 在上一节讲线性最小二乘法的解法过程中，由于残差函数$f_i(x)$为线性函数，所以它有解析解，解的形式如下：
$$A^{T}Ax=-A^{T}b(这里为了方便解释，定义f_i(x)=a_{i}x+b_i，所以这里右边一项多了个负号)\text{\hspace{1em}(4-13)}$$比较式(4-13)和式(4-11)，
$$\begin{array}{l}线性最小二乘：A^{T}Ax=-A^{T}b\\ \\ 高斯牛顿法：J^{T}J{h}_{gn}=-J^{T}f\end{array}$$不难看出，其实两种解的形式是一致的，即把$A$看成$J$，把$b$看成$f$。造成这一现象的原因是：高斯牛顿法对$f$里的每一个函数进行了一阶泰勒近似，将原本的非线性函数$f_i(x)$变成了线性函数。当其线性化完成后，原本的非线性最小二乘问题也就成为了线性最小二乘，所以它们解的形式是一致的。

非线性函数通过一阶泰勒展开舍弃高阶项可以将其变成线性函数，即非线性函数在$x$处的一阶泰勒近似使其线性化。因为线性化函数是非常容易操作的，而非线性函数就非常复杂。所以，高斯牛顿法的实质就是讲非线性函数通过泰勒展开近似乘线性函数，然后进行求解。

其实这一思想，很多理论都有用到。举个滤波理论的例子：对于线性系统，我们可以用卡尔曼滤波(KF)进行状态估计；若系统为非线性系统，则就要通过泰勒展开使其线性化后再用KF求解，这就是扩展卡尔曼(EKF)的思想。

3.列文伯格-马夸尔特方法(L-M Method)
L-M方法解的形式如下，
$$(J^{T}J+\mu I)h_{lm}=-J^{T}f\text{\hspace{1em}(4-13)}$$可以看出L-M方法就是在G-N方法的基础上加了一项阻尼项，所以它刚被提出来的时候，又叫做阻尼高斯牛顿法(damped Gauss-Newton Method)。
看式(4-13)，若$\mu$的取值足够大，则$(J^{T}J+\mu I)$ 必为正定矩阵，式(4-13)一定有解。从这个方面来看，L-M方法其实就是在G-N方法的基础上加了一个阻尼项，使得G-N方法一定有解。
但是从另外一个方面看，L-M方法也是2.8节的阻尼法中，对模型函数$L(h)$取$\frac{1}{2}l(h{)}^{T}l(h)$的结果。那么，与阻尼法的分析类似，L-M方法里得到的下降方向$h_{lm}$其实也是最速下降方向和高斯牛顿方向的一个折中或者选择：

• 当$\mu$很大时，$J^{T}J+\mu I\approx \mu$，所以有$h_{lm}\approx -\frac{J^{T}f}{\mu }=-\frac{F^{{}^{\prime }}(x)}{\mu }$，此时$h_{lm}$近似为最速下降方向；
• 当$\mu$很小时，$J^{T}J+\mu I\approx J^{T}J$，所以有$J^{T}J{h}_{lm}=-J^{T}f$，此时$h_{lm}$近似为高斯牛顿方向。

同样，与阻尼法的分析类似，L-M方法迭代的品质是可以通过gain ratio来很衡量的，即，
$$\rho =\frac{F(x)-F(x+h_{lm})}{L(0)-L(h_{lm})}\text{\hspace{1em}(4-14)}$$这里证明一下分母肯定为正，这点在2.8节里有说明，但由于前面的$L(h)$没有具体形式，所以没证明。
$$\begin{array}{l}L(0)-L(h_{lm})=-{h_{lm}}^T{J}^{T}f-\frac{1}{2}{h_{lm}}^T{J}^{T}J{h}_{lm}\\ \\ =-\frac{1}{2}{h_{lm}}^T(2J^{T}f+(J^{T}J+\mu I-\mu I)h_{lm})\\ \\ =\frac{1}{2}{h_{lm}}^T(\mu h_{lm}-J^{T}f)\end{array}\text{\hspace{1em}(4-15)}$$上式中${h_{lm}}^T{h}_{lm}$必为正，$-h_{lm}{J}^{T}f$也必为正(根据式(4-12))，所以分母为正。

所以，关于L-M方法中的参数$\mu$选取，可以参考阻尼法中的式(2-16)和式(2-17)的策略。关于$\mu$的初值选取，从以上的分析中可以看出，$\mu$和$J^{T}J$是有关系的，所以通常$\mu$的初值设定为： μ 0 = τ ∗ max ⁡ { ( J T J ) i i } {\mu _0} = \tau * \max \{ {({J^T}J)_{ii}}\} μ0​=τ∗max{(JTJ)ii​}，这里的$\tau$由用户设定。

L-M方法迭代优化的停止条件：

• 若$x^{\ast }$为局部最小值点，则必有$F^{{}^{\prime }}(x^{\ast })=g(x\ast )=0$。根据这一特点，可以采用如下的策略：
  $${\Vert g\Vert }_{\infty }\le {\epsilon }_1\text{\hspace{1em}(4-16)}$$这里的${\epsilon }_1$是一个很小的接近于零的正数，有用户自己设置。上式的含义是：经过若干次迭代后，存在一个$x$的使得$\Vert g\Vert$趋近于零。
• 若迭代前后两次的$x$变化很小，也可以认为$x$到达了局部最小值点：
  $$\Vert x_{new}-x\Vert \le {\epsilon }_2(\Vert x\Vert +{\epsilon }_2)\text{\hspace{1em}(4-17)}$$
• 设置最大迭代次数：$k<k_{\max }$

具体的L-M方法流程：
[ 初 始 化 ] : k = 0 , v = 2 , x = x 0 A = J T J , g = J T f f o u n d = ( ∥ g ∥ ∞ ≤ ε 1 ) , μ = τ ∗ max ⁡ { a i i } [ 循 环 ] : w h i l e      ( n o t      f o u n d )      a n d      ( k < k max ⁡ )                        k = k + 1 ; S o l v e ( A + μ I ) h l m = − g                        i f ( ∥ h l m ∥ ≤ ε 2 ( ∥ x ∥ + ε 2 ) )                                  f o u n d = t r u e                        e l s e                                  x n e w = x + h l m                                  r h o = F ( x ) − F ( x n e w ) L ( 0 ) − L ( h l m )                                  i f      ρ > 0                                            x = x n e w                                            A = J T J , g = J T f                                            f o u n d = ( ∥ g ∥ ∞ ≤ ε 1 )                                            μ = μ ∗ max ⁡ { 1 3 , 1 − ( 2 ρ − 1 ) 3 } , v = 2                                  e l s e                                            μ = μ ∗ v , v = 2 ∗ v \begin{array}{l} [初始化]:\\ k = 0,v = 2,x = {x_0}\\ A = {J^T}J,g = {J^T}f\\ found = ({\left\| g \right\|_\infty } \le {\varepsilon _1}),\mu = \tau * \max \{ {a_{ii}}\} \\\\ [循环]:\\ while\;\;(not\;\;found)\;\;and\;\;(k < {k_{\max }})\\ \;\;\;\;\;\;\;\;\;\;\;k = k + 1;Solve(A + \mu I){h_{lm}} = - g\\ \;\;\;\;\;\;\;\;\;\;\;if(\left\| {{h_{lm}}} \right\| \le {\varepsilon _2}(\left\| x \right\| + {\varepsilon _2}))\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;found = true\\ \;\;\;\;\;\;\;\;\;\;\;else\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{x_{new}} = x + {h_{lm}}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;rho = \frac{{F(x) - F({x_{new}})}}{{L(0) - L({h_{lm}})}}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;if\;\;\rho > 0\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;x = {x_{new}}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;A = {J^T}J,g = {J^T}f\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; found = ({\left\| g \right\|_\infty } \le {\varepsilon _1})\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\mu = \mu * \max \{ \frac{1}{3},1 - {(2\rho - 1)^3}\} ,v = 2\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;else\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\mu = \mu * v,v = 2 * v \end{array} [初始化]:k=0,v=2,x=x0​A=JTJ,g=JTffound=(∥g∥∞​≤ε1​),μ=τ∗max{aii​}[循环]:while(notfound)and(k<kmax​)k=k+1;Solve(A+μI)hlm​=−gif(∥hlm​∥≤ε2​(∥x∥+ε2​))found=trueelsexnew​=x+hlm​rho=L(0)−L(hlm​)F(x)−F(xnew​)​ifρ>0x=xnew​A=JTJ,g=JTffound=(∥g∥∞​≤ε1​)μ=μ∗max{31​,1−(2ρ−1)3},v=2elseμ=μ∗v,v=2∗v​

参考文献：
[1] METHODS FOR NON-LINEAR LEAST SQUARES PROBLEMS
[2] Unconstrained Optimization, 3rd Edition.