非线性最小二乘问题的求解方法（二）

在上一篇文章中已经介绍了非线性最小二乘问题的常规解法，对牛顿法、梯度法、线性搜索等进行了梳理，接下来将会介绍一些常用的非线性最小二乘问题的求解方法，它们在很多情形下能够有更好的表现。

3.非线性最小二乘问题的求解方法

这里仍再次给出非线性最小二乘问题：对于给定的函数$f:R^n\to R^m$，求解使得$||f(x)||$取得最小值时的$x$值，其中$m>=n$。这个问题等价于求解：

$x^{\ast }=argmi{n}_x${$F(x)$}
其中$F(x)=\frac{1}{2}{\sum }_{i=1}^m{f}_i(x{)}^2=\frac{1}{2}{\sum }_{i=1}^m||f_i(x)|{|}^2=\frac{1}{2}f(x{)}^{T}f(x)$

假设$f(x)$存在二阶连续偏导数，那么根据泰勒展开式，有：

$f(x+h)=f(x)+J(x)+O(||h||{)}^2$，其中矩阵$J\in R^{m\times n}$为雅可比矩阵，由$f(x)$的一阶偏导数组成，$(J(x){)}_{ij}=\frac{\partial f_i}{\partial x_j}(x)$

对于函数$F:R^n\to R$，根据$F(x)=\frac{1}{2}f(x{)}^{T}f(x)$，有$\frac{\partial F}{\partial x_j}(x)={\sum }_{i=1}^m{f}_i(x)\frac{\partial f_i}{\partial x_j}(x)$，因此

$F^{\prime }(x)=J(x{)}^{T}f(x)$

这里需要先计算函数$F$的海塞矩阵，该矩阵的$(j,k)$位置的元素为

$\frac{{\partial }^{2}F}{\partial x_j\partial x_k}(x)={\sum }_1^m(\frac{\partial f_i}{\partial x_j}(x)\frac{\partial f_i}{\partial x_k}(x)+\frac{{\partial }^2{f}_i(x)}{\partial x_j\partial x_k}(x))$

因此$F^{\prime \prime }(x)=J(x{)}^{T}J(x)+{\sum }_1^m{f}_i(x)f_i^{\prime \prime }(x)$
对于线性最小二乘问题$f(x)=b-Ax$，容易求得$F^{\prime }(x)=-A^T(b-Ax)$，求解一阶导数为零时的解$x^{\ast }$即为所求的解，这种方法称之为normal equations：

$(A^{T}A)x^{\ast }=A^{T}b$

上述问题可简化为$A{x}^{\ast }\approx b$，通过求解正交矩阵$Q$使得$Q^{T}A=\left[\begin{array}{c}R\\ 0\end{array}\right]$，然后回代法求解$R{x}^{\ast }=(Q^{T}b{)}_{1:n}$得到$x^{\ast }$，其中$R\in R^{n\times n}$为上三角矩阵。这种方法称为正交变换法。

3.1 Gauss-Newton方法

Gauss-Newton法在$x$的邻域内对函数$f$进行线性逼近：当$||h||$足够小时，根据泰勒展开式，有：

$f(x+h)\approx l(h)\equiv f(x)+J(x)h$

而对于$F(x)=\frac{1}{2}f(x{)}^{T}f(x)$则有：

$F(x+h)\approx L(h)\equiv \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$，其中$f=f(x)$，$J=J(x)$

Gauss-Newton法求解的迭代量$h_{gn}$可以使得$L(h)$取得极小值，即$h_{gn}=argmi{n}_h${$L(h)$}。很容易求得$L$的导数和海塞矩阵：

$L^{\prime }(h)=J^{T}f+J^{T}Jh$，$L^{\prime \prime }(h)=J^{T}J$

由以上可知，$L^{\prime }(0)=F^{\prime }(x)$；并且海塞矩阵$L^{\prime \prime }(h)$与$h$无关，且为对称矩阵。如果矩阵$J$满秩，则易证矩阵$L^{\prime \prime }(h)$为正定矩阵，此时表明$L(h)$具有唯一最小点，可以通过求解以下方程：

$(J^{T}J)h_{gn}=-J^{T}f$

该解为函数$F$的下降方向，因为：

$h_{gn}^T{F}^{\prime }(x)=h_{gn}^T(J^{T}f)=-h_{gn}^T(J^{T}J)h_{gn}<0$

因此可以通过上一章所述方法求解，通常迭代计算为， 求解$(J^{T}J)h_{gn}=-J^{T}f$，然后$x:=x+\alpha h_{gn}$，其中$\alpha$通过线性搜索法求解，经典的Gauss-Newton法在所有迭代步骤中取$\alpha =1$。当{$x|F(x)<=F(x_0)$}有界并且在所有迭代步骤中雅可比矩阵$J(x)$始终满秩时，结合线性搜索法的Gauss-Newton法可以被证明收敛。在上一篇文章提到，牛顿法求解优化问题时，最终阶段接近平方收敛，但是Gauss-Newton法却不一般不具备这样的性质。牛顿法和Gauss-Newton法两种方法求解下降方向的方程分别为：

牛顿法：$F^{\prime \prime }(x)h_n=-F^{\prime }(x)$ \quad \quad \quad \quad Gauss-Newton法：$L^{\prime \prime }(h)h_{gn}=-L^{\prime }(0)$

上述两个方程的等式右侧是相等的，但是等式左侧并不相等：

$F^{\prime \prime }(x)=L^{\prime \prime }(h)+{\sum }_{i=1}^m{f}_i(x)f_i^{\prime \prime }(x)$

因此，如果$f(x^{\ast })=0$，那么当$x$趋近于$x^{\ast }$时，有$L^{\prime \prime }(h)\approx F^{\prime \prime }(x)$，此时Gauss-Newton法最终收敛也接近平方收敛。如果函数{$f_i$}的曲率较小（即二次导数接近于0）或者函数{$|f_i(x^{\ast })|$}较小时，我们使用Gauss-Newton法能得到超线性收敛，但是总体来看Gauss-Newton法通常是线性收敛。需要注意的是，函数$F(x^{\ast })$的值对于收敛速度具有较大影响。
对于非线性最小二乘问题，无论是否采用线性搜索的方法，Gauss-Newton法都有可能求解失败。相对而言，牛顿法通过求解二阶导数往往能获得二次收敛，Gauss-Newton法虽然只能获得线性收敛，但效果在有些情况下表现不错。在3.2和3.3节将会介绍两种求解最优解的方法，它们对于求解全局最优解具有更优异的表现；而3.4节将会给出Gauss-Newton法的改进方法，从而能够获得超线性收敛。

3.2 Levenberg-Marquardt方法

Levenberg-Marquardt法其实是添加了阻尼项的Gauss-Newton法，迭代求解$h_{lm}$可通过求解以下方程得出：

$(J^{T}J+\mu I)h_{lm}=-g$，其中$g=J^T$，$\mu >=0$

其中$J=J(x)$，$f=f(x)$。阻尼参数$\mu$对迭代求解的影响如下：>

$1^o$对于所有$\mu >0$，系数矩阵为正定型，这就保证了$h_{lm}$为函数$L(x)$的下降方向；
$2^o$当$\mu$相当大时，可知$h_{lm}\approx -\frac{1}{\mu }g=-\frac{1}{\mu }{F}^{\prime }(x)$，此时表明$x$沿着梯度下降方向移动了较小的步长，在当前迭代值距离最终求解值较远时，这种迭代策略可以对整个迭代过程取得较好的效果；
$3^o$当$\mu$非常小时，可知$h_{lm}\approx h_{gn}$，在迭代的最终阶段，即$x$趋近于$x^{\ast }$时，是一种优良策略。如果$F(x^{\ast })=0$或$F(x^{\ast })$非常小时，那么这种方法可获得近似于平方收敛。

因此，阻尼系数同时影响了迭代的下降方向和步长，因此我们可以放弃精确的线性搜索求解，选择通过控制阻尼系数来调整迭代过程。在迭代开始时，$\mu$的初始值选取应该与矩阵$A_0=J(x_0{)}^{T}J(x_0)$中元素大小有关，取${\mu }_0=\tau \cdot ma{x}_i${$a_{ii}^{(0)}$}，其中$\tau$由算法工作人员自己决定。迭代过程可根据收益比进行调整：

$\rho =\frac{F(x)-F(x+h_{lm})}{L(0)-L(h_{lm})}$

上式分母是线性模型$L$的变化量：

$L(0)-L(h_{lm})=-h_{lm}{J}^{T}f-\frac{1}{2}{h}_{lm}^T{J}^{T}J{h}_{lm}=\frac{1}{2}{h}_{lm}^T(2g+(J^{T}J+\mu I-\mu I)h_{lm})=\frac{1}{2}{h}_{lm}^T(\mu h_{lm}-g)$

注意到，$h_{lm}^T{h}_{lm}$和$-h_{lm}^{T}g$项都是正数项，因此确保了该变化值$L(0)-L(h_{lm})$必定为正值。当$\rho$较大时，表明模型$L(h_{lm})$是函数$F(x+h_{lm})$的近似估计，所以我们可以通过减少阻尼系数$\mu$，从而在下一步迭代中使得迭代量$h_{lm}$更接近Gauss-Newton法的$h_{gn}$；当$\rho$较小甚至为负值时，表明模型$L(h_{lm})$对函数$F(x+h_{lm})$的近似估计效果较差，因此需要增加阻尼系数，从而使得迭代方向更接近梯度法下降方向并减小迭代步长。
算法终止迭代时要求在全局最小值点满足$F^{\prime }(x^{\ast })=g(x^{\ast })=0$，并且在当前迭代点处$x$的迭代变化非常小，并且为了避免无限迭代应限制迭代次数$k$，因此可以取以下迭代终止条件：

$||g|{|}_{infin}<={\epsilon }_1;||x_{new}-x||<={\epsilon }_2(||x||+{\epsilon }_2);k>=k_{max}$

其中${\epsilon }_1,{\epsilon }_2$为算法工作人员选定的较小的正数，最大迭代次数$k_{max}$应为某选定的正整数。不等式$||x_{new}-x||<={\epsilon }_2(||x||+{\epsilon }_2)$对于$x$取较大值到接近0时，变化从相对步长${\epsilon }_2$逐渐变化到绝对步长${\epsilon }_2^2$，因此对于不同的$x$值均可以有效平衡和限制。
由上述可知，Gauss-Newton法将最小二乘问题近似为线性问题：$f(x)+J(x)h\approx 0$，类似的，L-M法是这个线性问题的normal equations法的近似，即

$\left[\begin{array}{c}f(x)\\ 0\end{array}\right]+\left[\begin{array}{c}J(x)\\ \sqrt{\mu I}\end{array}\right]h\approx 0$

上述方程可通过正交变换进行精确求解，但是由于$h_{lm}$只是迭代的一步，所以并不需要精确求解，而通过normal equations法求解上述方程计算更简单，因此在实际问题中应用较为广泛。
L-M法的一般计算过程如下：

$begin$
\quad $k:=0;\nu :=2;x:=x_0$
\quad $A:=J(x{)}^{T}J(x);g:=J(x{)}^{T}f(x)$
\quad $found:=(||g|{|}_{\infty }<={\epsilon }_1)$
\quad $while(notfound)and(k<k_{max})$
\quad {
\quad \quad $k:=k+1;Solve(A+\mu I)h_{lm}=-g$
\quad \quad $if(||h_{lm}||<={\epsilon }_2(||x||+{\epsilon }_2)$
\quad \quad {
\quad \quad \quad $found:=true$
\quad \quad }
\quad \quad $else$
\quad \quad {
\quad \quad \quad $x_{new}:=x+h_{lm}$
\quad \quad \quad $\rho :=(F(x)-F(x_{new}))/(L(0)-L(h_{lm}))$
\quad \quad \quad $if(\rho >0)$
\quad \quad \quad {
\quad \quad \quad \quad $x:=x_{new}$
\quad \quad \quad \quad $A:=J(x{)}^{T}J(x);g:=J(x{)}^{T}f(x)$
\quad \quad \quad \quad $found:=(||g|{|}_{\infty }<={\epsilon }_1)$
\quad \quad \quad \quad $\mu :=\mu \ast max${$\frac{1}{3},1-(2\rho -1{)}^3$};$\nu :=2$
\quad \quad \quad }
\quad \quad \quad $else$
\quad \quad \quad {
\quad \quad \quad \quad $\mu :=\mu \ast \nu ;\nu :=2\ast \nu$
\quad \quad \quad }
\quad \quad }
\quad }
$end$

3.3 Powell’s Dog Leg方法

L-M方法通过添加阻尼项，结合了Gauss-Newton法和梯度下降法的优势。Powell’s Dog Leg法则是信赖域法的改进：
给定函数$f:R^n\to R^m$，Gauss-Newton法的迭代量$h_{gn}$是线性问题$J(x)h\approx -f(x)$的最小二乘解法，可通过求解normal equations：$(J(x{)}^{T}J(x))h_{gn}=-J(x{)}^{T}f(x)$获得。
而梯度法的下降方向由以下等式给出：$h_{sd}=-g=-J(x{)}^{T}f(x)$，迭代步长则需要根据线性模型求解：由$f(x+\alpha h_{sd})\approx f(x)+\alpha J(x)h_{sd}$可得$F(x+\alpha h_{sd})\approx \frac{1}{2}||f(x)+\alpha J(x)h_{sd}|{|}^2=F(x)+\alpha h_{sd}^{T}J(x{)}^{T}f(x)+\frac{1}{2}{\alpha }^2||J(x)h_{sd}|{|}^2$。随着$\alpha$变化使得$F(x+\alpha h_{sd})$取得最小值时，$F(\alpha )$的一阶导数为零，因此：$\alpha =-\frac{h_{sd}^{T}J(x{)}^{T}f(x)}{||J(x)h_{sd}||}=\frac{||g|{|}^2}{||J(x)g|{|}^2}$。
综上，我们下一步迭代有两种方案，一种是根据梯度法迭代当前点$x$移动$a=\alpha h_{sd}$，另一种是根据Gauss-Newton法从当前点$x$移动迭代量$b=h_{gn}$。而Powell’s Dog Leg方法根据信赖域半径$\Delta$采用以下的策略来选择迭代方案：

$if(||h_{gn}||<=\Delta )$
$h_{dl}:=h_{gn}$
$elseif(||\alpha h_{sd}||>=\Delta )$
$h_{dl}=(\Delta /||h_{sd}||)h_{sd}$
else
$h_{dl}=\alpha h_{sd}+\beta (h_{gn}-\alpha h_{sd})$，其中适当选择$\beta$使得$||h_{dl}||=\Delta$

根据上面对$a$和$b$的定义，并定义$c=a^T(b-a)$，可得：

$\psi (\beta )\equiv ||a+\beta (b-a)|{|}^2-{\Delta }^2=||b-a|{|}^2{\beta }^2+2c\beta +||\alpha |{|}^2-{\Delta }^2$

接下来讨论上述二次多项式的一个根（即多项式值为0时自变量取值），注意到当$\beta \to -\infty$时，$\psi \to +\infty$；$\psi (0)=||a|{|}^2-{\Delta }^2<0$；$\psi (1)=||h_{gn}|{|}^2-{\Delta }^2>0$，因此函数$\psi$必定存在一个负根并且在$[0,1]$上存在一个根。$[0,1]$上的根可由以下方法求解：

$if(c<=0)$
$\beta =({\Delta }^2-||a|{|}^2)/(c+\sqrt{c^2+||b-a|{|}^2({\Delta }^2-||a|{|}^2)})/||b-a|{|}^2$
else
$\beta =({\Delta }^2-||a|{|}^2)/(c+\sqrt{c^2+||b-a|{|}^2({\Delta }^2-||a|{|}^2)})$

收益比$\rho$的定义与L-M法中相似，为：

$\rho =\frac{F(x)-F(x+h_{dl}}{L(0)-L(h_{dl})}$

同样的，其中$L$为线性模型

$L(h)=\frac{1}{2}||f(x)+J(x)h|{|}^2$

在L-M法中采用收益比$\rho$来调整阻尼系数的大小，而在Powell’s Dog Leg方法中则采用收益比$\rho$来调整信赖域的半径。当$\rho$较大时表明当前线性模型近似效果较好，因此可以增加信赖域半径$\Delta$从而采取更大的步长，此时迭代方向更接近Gauss-Newton法的迭代方向；当$\rho$较小时（甚至负值）表明我们应当减小信赖域半径$\Delta$，采取更小的迭代步长，此时迭代方向更接近梯度下降法的迭代方向。总结如下：
$1^o$初始化。其中$x_0$和${\Delta }_0$应当自行给定；
$2^o$除3.2节迭代终止条件外，补充$||f(x)|{|}_{\infty }<={\epsilon }_3$，从而保证在$m=n$这种特殊的非线性系统模型下保证$f(x^{\ast })=0$；
$3^o$当$m=m$时，此时$J(x)h=-f(x)$而不是近似等于，而且不需使用normal equations进行求解；
$4^o$对于信赖域半径$\Delta$取值不同的三种情况，令$h1=h_{dl},h2=h_{gn},h3=\frac{-\Delta }{||g||}g$，有：
$L(0)-L(h_{dl})=\{\begin{array}{ll}F(x)& \text{h1=h2}\\ \frac{\Delta (2||\alpha g||-\Delta )}{2\alpha }& \text{h1=h3}\\ \frac{1}{2}\alpha (1-\beta {)}^2||g|{|}^2+\beta (2-\beta )F(x)& \text{otherwise}\end{array}$
$5^o$采用下列方法更新信赖域半径：

if($\rho <0.25$)
$\Delta :=\Delta /2$
else
$\Delta :$=max{$\Delta ,3\ast ||h||$}

$6^o$附加迭代终止条件：$\Delta <={\epsilon }_2(||x||+{\epsilon }_2)$，此时下一步迭代中条件$||x_{new}-x||<={\epsilon }_2(||x||+{\epsilon }_2)$必然满足。算法具体过程如下：

$begin$
\quad $k:=0;x:=x_0;\Delta :={\Delta }_0;g:=J(x{)}^{T}f(x)$
\quad $found:=(||f(x)|{|}_{\infty }<={\epsilon }_3)or(||g|{|}_{\infty }<={\epsilon }_1)$
\quad $while(notfound)and(k<k_{max})$
\quad {
\quad \quad $k:=k+1;\alpha =\frac{||g|{|}^2}{||J(x)g|{|}^2}$
\quad \quad $h_{sd}:=-\alpha g;SolveJ(x)h_{gn}\approx -f(x)$
\quad \quad 更新计算$h_{dl}$（方法见上文本节部分）
\quad \quad $if(||h_{dl}||<={\epsilon }_2(||x||+{\epsilon }_2)$
\quad \quad {
\quad \quad \quad $found:=true$
\quad \quad }
\quad \quad $else$
\quad \quad {
\quad \quad \quad $x_{new}:=x+h_{dl}$
\quad \quad \quad $\rho :=(F(x)-F(x_{new}))/(L(0)-L(h_{dl}))$
\quad \quad \quad $if(\rho >0)$
\quad \quad \quad {
\quad \quad \quad \quad $x:=x_{new};g:=J(x{)}^{T}f(x)$
\quad \quad \quad \quad $found:=(||f(x)|{|}_{\infty }<={\epsilon }_3)or(||g|{|}_{\infty }<={\epsilon }_1)$
\quad \quad \quad }
\quad \quad \quad $if(\rho >0.75)$
\quad \quad \quad {
\quad \quad \quad \quad $\Delta :$=max{$\Delta ,3\ast ||h||$}
\quad \quad \quad }
\quad \quad \quad $elseif(\rho <0.25)$
\quad \quad \quad {
\quad \quad \quad \quad $\Delta :=\Delta /2;found:=(\Delta <={\epsilon }_2(||x||+{\epsilon }_2))$
\quad \quad \quad }
\quad \quad }
\quad }
$end$

3.4 基于L-M和Quasi-Newton的混合方法

Madsen在1988年提出了一种基于L-M（如果$F(x^{\ast })=0$可获得平方收敛，否则接近线性收敛）和Quasi-Newton的混合方法（即使$F(x^{\ast })\ne 0$也能获得超线性收敛。这种迭代方法在初始阶段选择L-M方法，当$F(x^{\ast })$明显不等于零时，选择Quasi-Newton法，如果条件允许，可能又调整为L-M方法。选择Quasi-Newton法的条件为：在连续三次成功迭代过程中，满足$||F^{\prime }(x)|{|}_{\infty }<0.02\ast F(x)$。这个条件可以理解为当前迭代值$x$接近使得$F^{\prime }(x^{\ast })=0$的$x^{\ast }$，并且显然$F(x^{\ast })\ne 0$。如前所述，此时可能获得较慢的线性收敛。Quasi-Newton法从当前迭代点$x$出发，基于海塞矩阵$F^{\prime \prime }(x)$的近似矩阵$B$求解迭代量$h_{qn}$，通过以下方程求解：

$B{h}_{qn}=-F^{\prime }(x)$

以上方程是牛顿法的方程的近似。近似矩阵$B$在迭代中根据BFGS策略更新：在迭代过程中矩阵$B$始终是对称的（和$F^{\prime \prime }(x)$一样）和正定型，这就保证了$h_{qn}$为函数$F(x)$的下降方向。矩阵$B$的初始值取为对称正定型矩阵$B_0=I$，并且BFGS策略在每次迭代时将有秩为2的矩阵加到当前矩阵B上。Madsen在1988年提出的策略为：

$h:=x_{new}-x$；$y:J_{new}^T{J}_{new}h+(J_{new}-J{)}^{T}f(x_{new})$
$if(h^{T}y>0)$
$v:=Bh$；$B:=B+(\frac{1}{h^{T}y}y)y^T-(\frac{1}{h^{T}v}v)v^T$

其中：$J=J(x)$，$J_{new}=J(x_{new})$。由于当前$B$矩阵为正定型并且仅有$h_{T}y>0$时才会改变，经过上式变换可保证迭代后的矩阵$B$也是正定型（具体证明过程未给出，笔者将在后面的文章中进行推导补充）。Quasi-Newton法在迭代的全过程并不稳健，也不能确保所求迭代方向为下降方向。目标解$x^{\ast }$使得$F(x^{\ast })=0$，并且最终迭代阶段，好的求解模型迭代使得$||F^{\prime }(x)||$下降较快。如果$||F^{\prime }(x)||$下降较慢，则应选择L-M算法。总结如下：
$1^o$初始化。其中${\mu }_0$可由${\mu }_0=\tau \cdot ma{x}_i${$a_{ii}^{(0)}$}确定，迭代终止条件参考3.2节；
$2^o$在迭代过程中，由于对于同样的$x$，函数值$f$和雅可比矩阵$J$等是确定的，因此计算过程中这些值可以存储下来，避免同一步迭代中重复计算；
$3^o$迭代过程中，L-M和Q-N算法可为计算海塞矩阵的近似矩阵提供信息；
$4^o$迭代过程中，收益比$\rho$同样用来更新$\mu$值，参考L-M算法部分；
$5^o$注意更换迭代计算方法的时机。在迭代开始时，参数$count$需要被初始化为零；
$6^o$在3次连续迭代中，如果都满足$||F^{\prime }(x)|{|}_{\infty }<0.02\ast F(x)$，其中$\rho >0$，即每次$x$发生了改变，就应该改为Q-N迭代过程；
$7^o$结合信赖域算法与Quasi-Newton算法，当从L-M算法调整为Quasi-Newton算法时，信赖域半径$\Delta$初始化为$max${$1.5{\epsilon }_2(||x||+{\epsilon }_2),\frac{1}{5}||h_{lm}||$}；
$8^o$采用以下方法更新$\Delta$：

if($\rho <0.25$)
$\Delta :=\Delta /2$
else
$\Delta :$=max{$\Delta ,3\ast ||h||$}

$9^o$当$F^{\prime }$接近零时，允许函数$F$有极小的增大，即$\delta =\sqrt{{\epsilon }_M}$，其中${\epsilon }_M$为计算机的$unitroundoff$；
$10^o$对于Q-N迭代过程，下降的速度可能不够快

$begin$
\quad $k:=0;x:=x_0;\mu :={\mu }_0;B:=I$
\quad $found:=(||F^{\prime }(x)|{|}_{\infty }<={\epsilon }_1);method:=L-M$
\quad $while(notfound)and(k<k_{max})$
\quad {
\quad \quad $k:=k+1;$
\quad \quad $case(method)$
\quad \quad {
\quad \quad \quad $L-M:$
\quad \quad \quad \quad $[x_{new},found,better,method,\cdots ]:=LMstep(x,\cdots )$
\quad \quad \quad $Q-N:$
\quad \quad \quad \quad $[x_{new},found,better,method,\cdots ]:=QNstep(x,B,\cdots )$
\quad \quad }
\quad \quad 采用Madsen的策略（见上文）更新B
\quad \quad $if(better)$
\quad \quad {
\quad \quad \quad x:=x_{new};
\quad \quad }
\quad }
$end$

其中L-M迭代方法的步骤如下：

$[x_{new},found,better,method,\cdots ]:=LMstep(x,\cdots )$
$begin$
\quad $x_{new}:=x;method:=L-M$
\quad $Solve(J(x{)}^{T}J(x)+\mu I)h_{lm}=-F^{\prime }(x)$
\quad $if(||h_{lm}||<={\epsilon }_2(||x||+{\epsilon }_2))$
\quad {
\quad \quad $found:=true$
\quad }
\quad $else$
\quad {
\quad \quad $x_{new}:=x+h_{lm}$
\quad \quad $\rho =(F(x)-F(x_{new}))/((L(0)-L(h_{lm}))$
\quad \quad $if(\rho >0)$
\quad \quad {
\quad \quad \quad $better:=true;found:=(||F^{\prime }(x_{new})|{|}_{\infty }<={\epsilon }_1)$
\quad \quad \quad $if(||F^{\prime }(x_{new})|{|}_{\infty }<0.02\ast F(x_{new})$
\quad \quad \quad {
\quad \quad \quad \quad $count:=count+1$
\quad \quad \quad \quad $if(count=3)$
\quad \quad \quad \quad {
\quad \quad \quad \quad \quad method:=Q-N
\quad \quad \quad \quad }
\quad \quad \quad }
\quad \quad \quad $else$
\quad \quad \quad {
\quad \quad \quad \quad count:=0
\quad \quad \quad }
\quad \quad }
\quad \quad $else$
\quad \quad {
\quad \quad $count:=0;better:=false$
\quad \quad }
\quad }
$end$

Q-N迭代方法的步骤如下：

$[x_{new},found,better,method,\cdots ]:=QNstep(x,\cdots )$
$begin$
\quad $x_{new}:=x;method:=Q-N;better:=false$
\quad $SolveB{h}_{qn}=-F^{\prime }(x)$
\quad $if(||h_{qn}||<={\epsilon }_2(||x||+{\epsilon }_2))$
\quad {
\quad \quad $found:=true$
\quad }
\quad $else$
\quad {
\quad \quad $if(||h_{qn}||>\Delta )$
\quad \quad {
\quad \quad \quad $h_{qn}:=(\Delta /||h_{qn})\ast h_{qn}$
\quad \quad }
\quad \quad $x_{new}:=x+h_{lm}$
\quad \quad $if(||F^{\prime }(x_{new})|{|}_{\infty }<={\epsilon }_1)$
\quad \quad {
\quad \quad \quad $found:=true$
\quad \quad }
\quad \quad $else$
\quad \quad {
\quad \quad \quad $better:=(F(x_{new})<F(x))or((F^{\prime }(x_{new})<=(1+\delta )F(x))and(||F^{\prime }(x_{new})|{|}_{\infty }<||F^{\prime }(x)|{|}_{\infty }))$
\quad \quad \quad $if(F^{\prime }(x_{new})|{|}_{\infty }>=||F^{\prime }(x)|{|}_{\infty }$
\quad \quad \quad {
\quad \quad \quad \quad method:=L-M
\quad \quad \quad }
\quad \quad }
\quad }
$end$

3.5 L-M方法的割线型

本文所述的最小二乘问题中，假设函数$f$均为可微函数，即存在雅可比矩阵：$J(x)=\left[\frac{\partial f_i}{\partial x_j}\right]$。在很多实际情形中，由于函数$f$的具体形式未给出（即黑盒子），因此难以给出雅可比矩阵$J$中元素的计算公式，割线型L-M方法为这种问题提供了一种解决思路。对于这类问题，最简单的改进措施是采用差分法得到矩阵$B$来替换雅可比矩阵$J(x)$，其中矩阵$B_{ij}$（即矩阵$B$的第$(i,j)$项）可以通过有限差分近似进行求解：

$B_{ij}=\frac{\partial f_i}{\partial x_j}(x)\approx \frac{f_i(x+\delta e_j)-f_i(x)}{\delta }\equiv b_{ij}$，其中$e_j$为沿着第$j$坐标轴方向的单位向量，$\delta$为适当较小的实数

这种方法每次迭代一次$x$时需要计算$n+1$次函数$f$的值，并且由于$\delta$可能比距离$||x-x^{\ast }||$更小，所以在全局搜索的表现上这种方法并不能获得比计算函数$f(x)$更多的信息，因此需要对方法进行改进以提高效率。
对于函数$f:R^n\to R^m$的近似线性模型$f(x+h)\approx l(h)\equiv f(x)+J(x)h$，替换其中的雅可比矩阵$J(x)$，可得到：

$f(x+h)\approx \lambda (h)\equiv f(x)+Bh$，其中矩阵$B$是雅可比矩阵$J(x)$的近似

在下一步迭代中，对应$x_{new}$需要更新矩阵$B_{new}$，即：

$f(x_{new}+h)\approx \lambda (h)\equiv f(x_{new})+B_{new}h$

进一步地，如果模型能够保证迭代过程$h=x-x_{new}$，即：

$f(x)=f(x_{new})+B_{new}(x-x_{new})$

上式为我们提供了矩阵$B_{new}$的$m\cdot n$个元素中的$m$个等式，显然需要补充更多的条件。Broyden在1965年提出如下补充：

$B_{new}v=Bv$，其中$v$为满足$v\perp (x-x_{new})$的任意解

以上条件很容易被以下Broyden’s Rank One Update迭代过程满足：

$B_{new}=B+\mu h^T$，其中$h=x_{new}-x,u=\frac{1}{h^{T}h}(f(x_{new})-f(x)-Bh)$

注意到当$n=1$时条件 $f(x)=f(x_{new})+B_{new}(x-x_{new})$和条件$B_{new}v=Bv$是相一致的。这里加ing这种迭代方式成为广义割线型。L-M算法结合这种迭代方式进行改进后，大概形式如下：

$Solve(B^{T}B+\mu I)h_{slm}=-B^{T}f(x)$
$x_{new}:=x+h_{slm}$
采用Broyden’s Rank One Update迭代过程更新$B$
根据L-M算法更新$\mu$和$x$（见3.2节）

Powell已经证明了如果迭代过程中自变量$x_0,x_1,x_2,\cdots ,$收敛于$x^{\ast }$并且迭代步长{$h_k\equiv x_k-x_{k-1}$}满足列向量{$h_{k-n+1},\cdots ,h_k$}为线性无关(遍布整个空间$R^n$（其中k>n))，那么近似矩阵{$B_k$}收敛于雅可比矩阵$J(x^{\ast })$，而与初始矩阵$B_0$的选取无关。在实际情形中，前$n$步迭代往往并不满足上述假设（书中描述为未遍布整个空间$R^n$），因此可能在迭代若干次后，矩阵$B$是雅可比矩阵的较差近似估计，这就导致$-B^{T}f(x)$不一定是较好的下降方向甚至都不是下降方向。在这种情形下，变量$x$会一直保持不变而$\mu$则逐渐增大。此时虽然矩阵$B$发生了变化但是可能仍然是较差的近似估计，从而会导致$\mu$值的进一步增加。最终迭代过程可能因为$h_{slm}$过小（满足迭代终止条件之一）而终止迭代，尽管此时$x$距离$x^{\ast }$仍然较远。
针对这个问题，目前已经提出了许多方法，如偶尔采用差分法重新计算矩阵$B$。本文给出了一种基于循环坐标迭代过程的路线：由$h$控制当前迭代量，用$j$表示当前的坐标数。当$h$和向量$e_j$之间的夹角$\theta$过大，那么对于雅可比矩阵的第$j$列就需要用差分法重新计算近似值。准确地来说，条件为：

$cos\theta =\frac{|h^T{e}_j|}{||h||\cdot ||e_j||}<\gamma \rightleftharpoons |h_j|<\gamma ||h||$

实验表明选择$\gamma =0.8$具有较好的效果（悲观估计），在这种选择下每步迭代中大概需要计算$f$的两列向量。整个迭代方法总结如下：
$1^o$初始化。$x_0$为输入值，矩阵$B_0$为输入值或根据差分法计算，迭代终止条件及其参数参考3.2节，矩阵$B_0$差分法计算中$\delta$也为输入值；
$2^o$$mod(j,n)表示$j$被$n$除后的余数；
$3^o$迭代参数$\eta$的计算标准为：如果$x_j=0$，那么$\eta :={\delta }^2$；否则$\eta :=\delta |x_j|$；
$4^o$在迭代过程中，近似矩阵$B$每一步迭代都会发生改变，而$x$仅在下降方向满足时才会改变。因此尽管函数$f(x)$保持不变时，梯度近似估计$g$的值也会发生改变。
迭代过程如下：

$begin$
\quad $k:=0;x:=x_0;B:=B_0;j:=0$
\quad $g:=B^{T}f(x);found:=(||g|{|}_{\infty }<={\epsilon }_1)$
\quad $while(notfound)and(k<k_{max})$
\quad {
\quad \quad $k:=k+1;Solve(B^{T}B+\mu I)h=-g$
\quad \quad $if(||h||<={\epsilon }_2(||x||+{\epsilon }_2))$
\quad \quad {
\quad \quad \quad $found:=true$
\quad \quad }
\quad \quad $else$
\quad \quad {
\quad \quad \quad $j:=mod(j,n)+1$
\quad \quad \quad $if(|h_j|<0.8||h||)$
\quad \quad \quad {
\quad \quad \quad \quad $根据B_{new}=B+\mu h^T迭代计算B$
\quad \quad \quad }
\quad \quad \quad $if(F(x_{new})<F(x))$
\quad \quad \quad {
\quad \quad \quad \quad $x:=x_{new}$
\quad \quad \quad }
\quad \quad \quad $g:=B^{T}f(x);found:=(||g|{|}_{\infty }<={\epsilon }_1)$
\quad \quad }
\quad }
$end$

在许多应用中，$m$和$n$较大，但是每个函数$f_i(x)$仅仅依赖于$x$的若干元素。在这种情况下，一阶偏导数中$\frac{\partial f_i}{\partial x_j}(x)$很多都为零，即矩阵$J(x)$为稀疏矩阵。Nielen(1997)介绍了一系列在L-M算法中利用矩阵稀疏性的方法。在Broyden’s Rank One Update迭代公式中，正常情况下列向量$h$和$u$中元素都是非零项，所以矩阵$B_{new}$将会是稠密矩阵。本文不对这部分内容进一步介绍，具体内容可参考Gill(1984)和Toint(1987)的著作。

3.6 Dog Leg方法的割线型

采用割线型来对雅可比矩阵进行近似的方法也同样可以用到Dog Leg方法中。这里，我们考虑非线性系统问题中$m=n$的特殊情况。Broyden不仅给出了以下公式来更新雅可比矩阵的近似矩阵：

$B_{new}=B+(\frac{1}{h^{T}h}(y-Bh))h^T$，其中$h=x_{new}-x$，$y=f(x_{new})-f(x)$

并且还给出了雅可比矩阵的逆矩阵近似$D\approx J(x{)}^{-1}$的更新公式：

$D_{new}=D+(\frac{1}{h^{T}Dy}(h-Dy))(h^{T}D)$，其中$h=x_{new}-x$，$y=f(x_{new})-f(x)$

需要注意的是，在矩阵$D$的迭代公式中分母可能很小甚至等于0，因此如果$||h^{T}Dy||<\sqrt{{\epsilon }_M}||h||$，那么就不应用上述方法更新，而应用$D=B^{-1}$进行计算。
基于上述近似矩阵，梯度下降法的下降方向$h_{sd}$和Gauss-Newton法的迭代步长$h_{gn}$可以近似为：

$h_{ssd}=-B^{T}f(x)$，并且$h_{sgn}=-Df(x)$

Dog Leg算法简单改进后就可以用到上述的近似量，初始矩阵$B=B_0$可以通过不同的近似方法得到，并且$D_0$可以近似为$B_0^{-1}$。显然可知当前情况下，矩阵$B$和矩阵$D$满足$BD=I$。参数$\alpha$可以由下式求得：

$\alpha =-\frac{h^T{B}^{T}f(x)}{||Bh||}=\frac{||g|{|}^2}{||Bg|{|}^2}$

和L-M算法的割线型类似，这种方法也需要额外计算来进行迭代从而保证矩阵$B$和$D$是雅可比矩阵和其逆矩阵的近似估计效果。我们发现采用3.6节中$cos\theta =\frac{|h^T{e}_j|}{||h||\cdot ||e_j||}<\gamma \rightleftharpoons |h_j|<\gamma ||h||$公式在这里也有较好的表现。
每一次迭代计算$D_{new}$大概需要$10n^2$次浮点运算，并且计算$h_{ssd}$和$h_{sgn}$以及$\alpha$大概需要$6n^2$次浮点计算，因此结合自由梯度法的Dog Leg方法每次迭代除计算$f(x_{new})外需要$16n^2$次浮点运算，而Dog Leg算法每次迭代除计算$f(x_{new})$和$J(x_{new})$外还需$\frac{2}{3}{n}^3+6n^2$次浮点运算。因此当$n$很大时，结合自由梯度法的Dog Leg方法效率更高。但是考虑到迭代次数通常非常大，并且如果雅可比矩阵可以计算时，梯度法的Dog Leg方法通常收敛更快。

3.7 总结

文中介绍了一系列用于解决非线性最小二乘问题的算法，可在任何一个优秀的程序库中找到，其具体应用可以在GAMS(Guide to Available Mathematical Software)中找到，网址为：$http://gams.nist.gov$。本文中的示例可以在MATLAB中进行计算，相关的程序可以toolbox immoptibox中找到，详见：$http:/www.imm.dtu.dk/~hbn/immoptibox。最后，需要提及的是，针对具体问题，有时候重新推导公式可能会使得问题简化。（博客仅仅将文中的思路和文字进行了翻译，舍弃了大量的examples，但是这些example对于深刻理解和熟悉算法是特别重要的，感兴趣的作者可移步原文）。