最小二乘问题

最小二乘问题

$$\min f\left(x\right)=\frac{1}{2}\sum _{i=1}^m{r}_i^2\left(x\right)=\frac{1}{2}r{\left(x\right)}^{T}r\left(x\right)x\in {\mathbb{R}}^n,m⩾n\text{\hspace{1em}(1)}$$
这里$r\left(x\right)={\left(r_1\left(x\right),r_2\left(x\right),\cdots ,r_m\left(x\right)\right)}^T$称为剩余函数，点$x$ 处剩余函数的值称为剩余量。若$r_i(x)$均为线性函数，则问题(1)为线性最小二乘问题，若至少有一个$r_i(x)$为非线性函数，则问题(1) 为非线性最小二乘问题。

f(x)的导数

设$J(x)$为$r(x)$的Jacobian矩阵
$$J\left(x\right)=\frac{\partial r}{\partial x}={\left[\nabla r_1\left(x\right),\dots ,\nabla r_m\left(x\right)\right]}^T\in {\mathbb{R}}^{m\times n}\text{\hspace{1em}(2)}$$
则$f(x)$的梯度为
$$g\left(x\right)=\nabla f\left(x\right)=\sum _{i=1}^m{r}_i\left(x\right)\nabla r_i\left(x\right)=J^T\left(x\right)r\left(x\right)\text{\hspace{1em}(3)}$$
$f(x)$的Hesse矩阵为
$$\begin{array}{rl}G\left(x\right)& ={\nabla }^{2}f\left(x\right)=\sum _{i=1}^m\nabla r_i\left(x\right)\nabla r_i{\left(x\right)}^T+\sum _{i=1}^m{r}_i\left(x\right){\nabla }^2{r}_i\left(x\right)\\ & =J^T\left(x\right)J\left(x\right)+S\left(x\right)\end{array}\text{\hspace{1em}(4)}$$
其中
$$S(x)=\sum _{i=1}^m{r}_i\left(x\right){\nabla }^2{r}_i\left(x\right)\text{\hspace{1em}(5)}$$
为便于讨论，我们采用以下记号：
$$\begin{gathered} J^{\ast }=J(x^{\ast }),J_k=J(x_k) \\ {S}^{\ast }=S(x^{\ast }),S_k=S(x_k) \end{gathered}$$

最小二乘问题的分类

在点$x^{\ast }$处，$\Vert S^{\ast }\Vert$的大小取决于剩余量与问题的非线性程度，对零剩余或线性最小二乘问题，$\Vert S^{\ast }\Vert =0$，随着剩余量的增大或$e_i(x)(i=1,\cdots ,m)$的非线性程度的增强，$\Vert S^{\ast }\Vert$的值变大。根据问题的这种特点，将算法分为小剩余算法和大剩余算法。小剩余算法处理$\Vert S^{\ast }\Vert$为零或不太大的问题，大剩余算法处理$\Vert S^{\ast }\Vert$较大的问题。

Newton法解最小二乘问题

$$f\left(x\right)=f\left(x_k\right)+\nabla f{\left(x_k\right)}^T\left(x-x_k\right)+\frac{1}{2}{\left(x-x_k\right)}^T{\nabla }^{2}f\left(x_k\right)\left(x-x_k\right)+O\left({\Vert x-x_k\Vert }^2\right)$$
使用二阶泰勒展开进行局部近似，这是一个二次型
$$q\left(x\right)=f\left(x_k\right)+\nabla f{\left(x_k\right)}^T\left(x-x_k\right)+\frac{1}{2}{\left(x-x_k\right)}^T{\nabla }^{2}f\left(x_k\right)\left(x-x_k\right)$$
二次型的极值可以通过令导数为0求得
$$q^{\prime }\left(x\right)=\nabla f\left(x_k\right)+{\nabla }^{2}f\left(x_k\right)\left(x-x_k\right)=0$$
令$d=x-x_k$为增量，代入$\nabla f(x),{\nabla }^{2}f(x)$得
$$\left(J_k^T{J}_k+S_k\right)d=-J_k^T{r}_k\text{\hspace{1em}(6)}$$
对最小二乘问题， Newton 方法的缺点是每次迭代都要求 $S_k$ ，即计算m个$n\times n$对称矩阵．显然，对一个算法而言， $S_k$ 的计算是一个沉重的负担．解决这个问题的方法是或者在 Newton 方程中忽略 $S_k$ ，或者用一阶导数信息近似$S_k$ 。而要忽略 $S_k$ ，则应在 $r_i(x)$接近于0或接近于线性时进行。这就是下面我们要讲的小剩余算法。

Gauss-Newton法

在Newton方程(6)中忽略$S_k$就得到Gauss-Newton(GN)方法。该方法也可以这样理解，在点$x_k$处线性化剩余函数$r_i(x_k+d)$，我们得到关于$d$的线性最小二乘问题
$$\underset{d\in {\mathbb{R}}^n}{\min }{q}_k(d)=\frac{1}{2}\Vert r_k+J_{k}d{\Vert }_2^2\text{\hspace{1em}(7)}$$
其中
$$\begin{array}{rl}{q}_k(d)& =\frac{1}{2}(J_{k}d+r_k{)}^T(J_{k}d+r_k)\\ & =\frac{1}{2}{d}^T{J}_k^T{J}_{k}d+d^T\left(J_k^T{r}_k\right)+\frac{1}{2}{r}_k^T{r}_k\end{array}\text{\hspace{1em}(8)}$$
这里$q_k(d)$是对$f(x_k+d)$的一种二次近似，它与$f(x_k+d)$的二次Taylor近似的差别在于二次项中少了$S_k$。
问题(7)的极小点$d_k$满足
$$J_k^T{J}_k{d}_k=-J_k^T{r}_k\text{\hspace{1em}(9)}$$
式(9)称为Gauss-Newton方程，由(9)式得到的方向$d_k$称为Gauss-Newton方向。

用 Gauss-Newton 方法求解最小二乘问题的算法如下

算法1 (Gauss-Newton 方法求解最小二乘问题）

1. 给定$x_0,\epsilon >0,k:=0$；
2. 若终止条件满足，则停止迭代；
3. 解$J_k^T{J}_{k}d=-J_k^T{r}_k$得$d_k$ ;
4. $x_{k+1}:=x_k+{\alpha }_k{d}_k$ ，其中${\alpha }_k$是一维搜索结果， $k:=k+1$，转2.

基本Gauss-Newton方法是指${\alpha }_k=1$的Gauss-Newton方法．带线搜索的Gauss-Newton方法称为阻尼Gauss-Newton 方法．

Gauss-Newton方法的优点在于它无须计算$r(x)$的二阶导数．另外，由(3)式和(9)式知
$$d_k^T{g}_k=d_k^T{J}_k^T{r}_k=-d_k^T{J}_k^T{J}_k{d}_k=-{\Vert J_k{d}_k\Vert }^2$$

这说明．当$J_k$满秩，$g_k$非零时，$d_k$是下降方向。

定理2（基本Gauss-Newton 方法的局部收敛性）
设 $r_i(x)\in C^2(i=1,\cdots ,m),x^{\ast }$是最小二乘问题(1)的最优解，且$J^{\ast T}{J}^{\ast }$正定。假设由基本Gauss-Newton法迭代产生的点列 { x k } \{x_k\} {xk​}收敛于$x^{\ast }$，则当$G(x)$与$J(x{)}^{T}J(x)$在$x^{\ast }$的邻域内Lipschitz连续时，有
$$\Vert h_{k+1}\Vert ⩽\Vert {\left(J^{\ast T}{J}^{\ast }\right)}^{-1}\Vert \Vert S^{\ast }\Vert \Vert h_k\Vert +O\left({\Vert h_k\Vert }^2\right)$$
其中$h_k=x_k-x^{\ast }$。
证明
因为$f\in C^2$，且$G(x)$ 在$x^{\ast }$的邻域内Lipschitz连续，当$x_k$充分接近$x^{\ast }$时，由Newton法收敛性的定理证明知
$$g\left(x_k+d\right)=g_k+G_{k}d+O\left({\Vert d\Vert }^2\right)$$
令$d=-h_k$，得
$$0=g^{\ast }=g_k-G_k{h}_k+O\left({\Vert h_k\Vert }^2\right)$$
将(3)(4)式代入上式得
$$J_k^T{r}_k-\left(J_k^T{J}_k+S_k\right)h_k+O\left({\Vert h_k\Vert }^2\right)=0\text{\hspace{1em}(10)}$$
因为$J^{\ast T}{J}^{\ast }$正定，当$x_k$充分接近$x^{\ast }$时，$J_k^T{J}_k$亦正定，我们用$(J_k^T{J}_k{)}^{-1}$左乘(10)式，由(8)式得
$$-d_k-h_k-{\left(J_k^T{J}_k\right)}^{-1}{S}_k{h}_k+O\left({\Vert h_k\Vert }^2\right)=0$$
因为
$$d_k+h_k=x_{k+1}-x_k+x_k-x^{\ast }=h_{k+1}$$
所以
$$\begin{array}{rl}{h}_{k+1}& =-{\left(J_k^T{J}_k\right)}^{-1}{S}_k{h}_k+O\left({\Vert h_k\Vert }^2\right)\\ \Vert h_{k+1}\Vert & ⩽\Vert {\left(J_k^T{J}_k\right)}^{-1}{S}_k\Vert \Vert h_k\Vert +O\left({\Vert h_k\Vert }^2\right)\\ & ⩽\Vert {\left(J_k^T{J}_k\right)}^{-1}{S}_k-{\left(J^{\ast T}{J}^{\ast }\right)}^{-1}{S}^{\ast }\Vert \Vert h_k\Vert +\Vert {\left(J^{\ast T}{J}^{\ast }\right)}^{-1}\Vert \Vert S^{\ast }\Vert \Vert h_k\Vert +O\left({\Vert h_k\Vert }^2\right)\end{array}\text{\hspace{1em}(11)}$$
在下面关于$S(x)$和$(J(z{)}^{T}J(z){)}^{-1}$在 $x^{\ast }$的邻域内Lipschitz连续的证明中，对于任意矩阵 $A(x)$，我们采用记号 $A_x=A(x)$．因为$G_x$和$J_x^T{J}_x$ 在$x^{\ast }$的邻域中Lipschitz连续，所以存在$\beta ,\gamma >0$，使得对$x^{\ast }$邻域内的任意两点 $x,y$ ，有
$$\begin{array}{rl}\Vert G\left(x\right)-G\left(y\right)\Vert & ⩽\beta \Vert x-y\Vert \\ \Vert J{\left(x\right)}^{T}J\left(x\right)-J{\left(y\right)}^{T}J\left(y\right)\Vert & ⩽\gamma \Vert x-y\Vert \end{array}$$
从而
$$\begin{array}{rl}\Vert S\left(x\right)-S\left(y\right)\Vert & =\Vert G\left(x\right)-a\left(y\right)-J{\left(x\right)}^{T}J\left(x\right)+J{\left(y\right)}^{T}J\left(Y\right)\Vert \\ & ⩽\Vert G\left(x\right)-G\left(y\right)\Vert +\Vert J{\left(x\right)}^{T}J\left(x\right)-J{\left(y\right)}^{T}J\left(y\right)\Vert \\ & ⩽\left(\beta +\gamma \right)\Vert x-y\Vert \end{array}$$
对$x^{\ast }$邻域内的任意点$x$，由$J^{\ast T}{J}^{\ast }$的正定性知，存在$\xi >0$，使得$\Vert (J_x^T{J}_x{)}^{-1}\Vert ⩽\xi$，从而
$$\begin{array}{rl}\Vert {\left(J_x^T{J}_x\right)}^{-1}-{\left(J_y^T{J}_y\right)}^{-1}\Vert & =\Vert {\left(J_x^T{J}_x\right)}^{-1}\left(J_y^T{J}_y-J{x}^T{J}_x\right){\left(J_y^{T}Jy\right)}^{-1}\Vert \\ & ⩽\Vert {\left(J_x^T{J}_x\right)}^{-1}\Vert \Vert {\left(J_y^T{J}_y\right)}^{-1}\Vert \Vert J_y^T{J}_y-J_x^T{J}_x\Vert \\ & ⩽\gamma {\xi }^2\Vert x-y\Vert \end{array}$$
所以$S_x$ 与$(J_x^T{J}_x{)}^{-1}$也在$x^{\ast }$的邻域内Lipschitz连续。
当$x_k$充分接近$x^{\ast }$时，有
$$\begin{array}{rl}& \Vert {\left(J_k^T{J}_k\right)}^{-1}{S}_k-{\left(J^{\ast T}{J}^{\ast }\right)}^{-1}{S}^{\ast }\Vert \\ & ⩽\Vert {\left(J_k^T{J}_k\right)}^{-1}{S}_k-{\left(J_k^T{J}_k\right)}^{-1}{S}^{\ast }\Vert +\Vert {\left(J_k^T{J}_k\right)}^{-1}{S}^{\ast }-{\left(J^{\ast T}{J}^{\ast }\right)}^{-1}{S}^{\ast }\Vert \\ & ⩽\left(\beta +\gamma \right)\Vert {\left(J_k^T{J}_k\right)}^{-1}\Vert \Vert h_k\Vert +\gamma {\xi }^2\Vert S^{\ast }\Vert \Vert h_k\Vert \\ & ⩽\left(\left(\beta +\gamma \right)\xi +\gamma {\xi }^2\Vert S^{\ast }\Vert \right)\Vert h_k\Vert \end{array}$$
所以
$$\Vert {\left(J_k^T{J}_k\right)}^{-1}{S}_k-{\left(J^{\ast T}{J}^{\ast }\right)}^{-1}{S}^{\ast }\Vert \Vert h_k\Vert ⩽\left(\left(\beta +\gamma \right)\xi +\gamma {\xi }^2\Vert S^{\ast }\Vert \right){\Vert h_k\Vert }^2$$
将上式代入(11)式可得
$$\Vert h_{k+1}\Vert ⩽\Vert {\left(J^{\ast T}{J}^{\ast }\right)}^{-1}\Vert \Vert S^{\ast }\Vert \Vert h_k\Vert +O\left({\Vert h_k\Vert }^2\right)$$
故定理结论成立。

该定理说明，若 $x_k\to x^{\ast }$，基本Gauss-Newton方法有如下两种情形的收敛速度：

• 二阶收敛速度．若$\Vert S(x^{\ast })\Vert ＝0$，即在零剩余问题或是线性最小二乘问题的情形，则方法在$x^{\ast }$附近具有Newton方法的收敛速度．
• 线性收敛速度．若$\Vert S(x^{\ast })\Vert \ne 0$，则方法的收敛速度是线性的，收敛速度随$S(x^{\ast })$的增大而变慢．

由此可见，基本Gauss-Newton方法的收敛速度是与$x^{\ast }$处剩余量的大小及剩余函数的线性程度有关的，即剩余量越小或剩余函数越接近线性，它的收敛速度就越快；反之就越慢，甚至对剩余量很大或剩余函数的非线性程度很强的问题不收敛.

LM方法

Gauss-Newton方法在迭代中会出现$J_k^T{J}_k$为奇异的情况，为了克服这个困难，提出LM (Levenberg-Marquardt)方法。修改Gauss-Newton方程为LM方程：
$$(J_k^T{J}_k+{\gamma }_{k}I)d=-J_k^T{r}_k\text{\hspace{1em}(12)}$$
其中${\gamma }_k⩾0$，使得$J_k^T{J}_k+{\gamma }_{k}I$正定，从计算角度出发，为保证该矩阵充分正定，${\gamma }_k$可能需要取得适当的大，$J_k^T{J}_k+{\gamma }_k$的正定性保证了得到的方向是下降方向。

（LM方程与信赖域问题的关系）$d_k$为信赖域子问题
$$\begin{array}{rl}& \underset{d}{\min }\frac{1}{2}{\Vert J_{k}d+r_k\Vert }^2\\ & s.t.{\Vert d\Vert }^2⩽{\Delta }_k^2,{\Delta }_k>0\end{array}\text{\hspace{1em}(13)}$$
的全局极小解的充分必要条件是，对满足(13)式的$d_k$，存在${\gamma }_k⩾0$，使得
$$\begin{array}{rl}& \left(J_k^T{J}_k+{\gamma }_{k}I\right)d_k=-J_k^T{r}_k\\ & {\gamma }_k\left({\Delta }_k^2-{\Vert d_k\Vert }^2\right)=0\end{array}\text{\hspace{1em}(14)}$$
证明：

必要性：

对于优化问题(13)，由有约束优化问题的最优性条件知，存在${\gamma }_k⩾0$，使得$d_k,{\gamma }_k$满足KKT条件。
Lagrange函数：
$$L(d,\gamma )=\frac{1}{2}{\Vert J_{k}d+r_k\Vert }^2-\frac{1}{2}\gamma \left({\Delta }_k^2-{\Vert d\Vert }^2\right)$$
KKT条件：
$$\begin{gathered} {\nabla }_{d}L\left(d_k,{\gamma }_k\right)=0\Rightarrow J_k^T{r}_k+\left(J_k^T{J}_k+{\gamma }_{k}I\right)d_k=0 \\ {\gamma }_k\left({\Delta }_k^2-{\Vert d_k\Vert }^2\right)=0\text{\hspace{1em}(互补性条件)} \end{gathered}$$
充分性：

因为$J_k^T{J}_k+{\gamma }_{k}I$半正定，所以方程(14)上式的解 $d_k$是
$${\tilde{q}}_k\left(d\right)=\frac{1}{2}{d}^T(J_k^T{J}_k+{\gamma }_{k}I)d+d^T\left(J_k^T{r}_k\right)+\frac{1}{2}{r}_k^T{r}_k$$
的全局极小点。由(8)式有
$${\tilde{q}}_k\left(d\right)=q_k\left(d\right)+\frac{1}{2}{\gamma }_k{\Vert d\Vert }^2$$
因为任给$d\in {\mathbb{R}}^n$，有${\tilde{q}}_k(d)⩾{\tilde{q}}_k(d_k)$
$$q_k\left(d\right)⩾q_k\left(d_k\right)+\frac{1}{2}{\gamma }_k\left({\Vert d_k\Vert }^2-{\Vert d\Vert }^2\right)$$
由(14)下式知，若${\gamma }_k=0$，有$q_k(d)⩾q_k(d_k)$；若${\gamma }_k\ne 0$，有$\Vert d_k{\Vert }^2={\Delta }_k^2$，所以
$$q_k(d)⩾q_k(d_k)+\frac{1}{2}{\gamma }_k({\Delta }_k^2-\Vert d{\Vert }^2)$$
这说明，对任意${\gamma }_k⩾0$和任意满足$\Vert d{\Vert }^2⩽{\Delta }_k^2$的$d$，$d_k$是问题(13)的全局最优解。

下面来考虑${\gamma }_k$的修正方法，它与信赖域半径${\Delta }_k$的修正是相关的。在信赖域方法中，从$x_k$到$x_k+d_k$，$f(x)$的实际减少量为
$$\Delta f_k=f(x_k)-f(x_k+d_k)$$
由(8)式给出的$f(x_k+d)$的二次近似函数$q_k(d)$的减少量为
$$\Delta q_k=q_k(0)-q_k(d_k)$$
这里$q_k(0)=f_k$，另外，由LM方程与$d_k^T{g}_k<0$知
$$\begin{array}{rl}\Delta q_k& =q_k\left(0\right)-q_k\left(d_k\right)\\ & =-\frac{1}{2}{d}_k^T{J}_k^T{J}_k{d}_k-d_k^T\left(J_k^T{r}_k\right)\\ & =\frac{1}{2}{d}_k^T\left(-J_k^T{J}_k{d}_k-{\gamma }_k{d}_k+{\gamma }_k{d}_k-2J_k^T{r}_k\right)\\ & =\frac{1}{2}{d}_k^T\left(-\left(J_k^T{J}_k+{\gamma }_{k}I\right)d_k+{\gamma }_k{d}_k-2J_k^T{r}_k\right)\\ & =\frac{1}{2}{d}_k^T\left({\gamma }_k{d}_k-g_k\right)>0\end{array}$$
其中$g_k=J_k^T{r}_k$
定义
$${\rho }_k=\frac{\Delta f_k}{\Delta q_k}$$
在第k步迭代，${\rho }_k$的值可以反映出$q_k(d_k)$近似$f(x_k+d_k)$的好坏。由LM方程知，${\gamma }_k$可以控制$\Vert d_k\Vert$的大小，从而可以控制信赖域的大小，${\gamma }_k$的修正应与信赖域方法中对${\Delta }_k$大小的修正相反。${\rho }_k$小于阈值，则说明近似效果差，应收缩信赖域，扩大$\gamma$，否则则相反。