数值最优化—非线性方程组与最小二乘问题

一、参考

《数值最优化算法与理论》

二、非线性方程组的局部算法

1. 局部Newton法

设$F:R^n\to R^n$连续可微。考察如下非线性方程组：
$$F(x)=0$$
若$F$是某个函数$f:R^n\to R^n$的梯度，由 数值最优化—无约束问题的下降算法与线性搜索：无约束问题解的一阶必要条件 可知，上面方程表示无约束问题$minf(x)$的最优解的一阶必要条件。

记$F^{\prime }(x)$为函数$F$在$x$处的$Jacobi$矩阵。求解上面方程的局部Newton法的迭代格式如下：$x^{(k+1)}=x^{(k)}+d^{(k)},d^{(k)}$是线性方程组
$$F^{\prime }(x^{(k)})d+F(x^{(k)})=0$$
的解。

经过与 数值最优化—无约束问题最速下降法和Newton法：Newton法 比较，不难发现，若$F$是函数$f:R^n\to R$的梯度，则求解非线性方程组的局部Newton法与求解无约束问题$minf(x)$的古典Newton法一致。因此，求解非线性方程组的Newton法是求解无约束最优化问题Newton法的一种推广。

2. 局部拟Newton法

局部拟Newton法的迭代格式为：$x^{(k+1)}=x^{(k)}+d^{(k)},d^{(k)}$是线性方程组：
$$B_k{d}^{(k)}+F(x^{(k)})=0$$
的解，其中，矩阵$B_k$是$F^{\prime }(x^{(k)})$的近似，它满足下面的拟Newton（割线）方程：
$$B_{k+1}{s}^{(k)}=y^{(k)}$$
其中$y^{(k)}=F(x^{(k+1)})-F(x^{(k)}),s^{(k)}=x^{(k+1)}-x^{(k)}$。注意到$y^{(k)}\approx F^{\prime }(x^{(k1)})s^{(k)}$，因此，$B_{k+1}$与$F^{\prime }(x^{(k+1)})$沿方向$s^{(k)}$很接近。

三、非线性方程组的全局化算法

1. 阻尼Newton法

设$F(x)=(F_1(x),F_2(x),\cdot \cdot \cdot ,F_n(x){)}^T$。引入函数：
$$\theta (x)=\frac{1}{2}||F(x)|{|}^2=\frac{1}{2}\sum _{i=1}^n{F}_i(x{)}^2$$
我们称$\theta$是方程组的残量函数活模函数。经计算，容易得到$\theta$的梯度具有如下形式：
$$\nabla \theta (x)=\sum _{i=1}^n{F}_i(x)\nabla F_i(x)=F^{\prime }(x{)}^{T}F(x)$$
设$\bar{d}$是下面的线性方程组的解：
$$F^{\prime }(x)d+F(x)=0$$
则有：
$$\nabla \theta (x{)}^T\bar{d}=-||F(x)|{|}^2<0$$
即$\bar{d}$是函数$\theta$在$x$处的一个下降方向。利用此性质，我们可以构造求解非线性方程的线性搜索型Newton法。称之为阻尼Newton法。如下：

1. 给定初始点$x^{(0)}\in R^n$，常数$\rho \in (0,1),{\sigma }_1\in (1,1/2)$，精度$ϵ>0$，令$k=0$。
2. 若$||F(x^{(k)})||\le ϵ$，则终止算法，得解$x^{(k)}$。否则，转3。
3. 解线性方程序组 $F^{\prime }(x^{(k)})d+F(x^{(k)})=0$，得方向$d^{(k)}$。
4. 确定步长${\alpha }_k$ 为集合 { ρ i ∣ i = 0 , 1 , 2 ⋅ ⋅ ⋅ } \{\rho ^i | i=0,1,2···\} {ρi∣i=0,1,2⋅⋅⋅}中使得下面的不等式成立的最大者：
   $$\theta (x^{(k)}+{\alpha }_k{d}^{(k)})\le (1-2{\sigma }_1{\alpha }_k)\theta (x^{(k)})$$
5. 令$x^{(k+1)}=x^{(k)}+{\alpha }_k{d}^{(k)}，k=k+1$。转2。

2. 信赖域算法

类似于求解无约束问题的信赖域算法，我们构建求解非线性方程组的信赖域算法，其子问题为：
$$\{\begin{array}{l}min{m}_k(d)\stackrel{△}{=}\frac{1}{2}||F(x^{(k)})+F^{\prime }(x^{(k)})d|{|}^2\\ s.t.||d||\le {\Delta }_k\end{array}$$
其中${\Delta }_k>0$是信赖域半径。记该问题的解为$d^{(k)}$。函数$m_k$是$\theta$的近似函数，它由在$\theta$中对$F$作线性近似产生。Newton-信赖域算法如下：

1. 取初始点$x^{(0)}\in R^n,\bar{\Delta }>0,{\Delta }_0\in (0,\bar{\Delta }),\eta \in [0,\frac{1}{4})$，精度$ϵ>0$，令$k=0$。
2. 若$||\theta (x^{(k)})||\le ϵ$，则终止算法，得解$x^{(k)}$。否则，转3。
3. 求解上面的信赖域子问题，得解$d^{(k)}$。
4. 计算$r_k$。若$r_k>\frac{3}{4}$，则令 Δ k + 1 = m i n { 2 Δ k , Δ ˉ } \Delta _{k+1} = min\{2\Delta _k,\bar \Delta \} Δk+1​=min{2Δk​,Δˉ}；若$\eta <r_k<\frac{1}{4}$，则令${\Delta }_{k+1}=\frac{1}{2}{\Delta }_k$；若$\frac{1}{4}\le r_k\le \frac{3}{4}$，则令${\Delta }_{k+1}={\Delta }_k$。
5. 若$r_k\le \eta$，令$x^{(x+1)}=x^{(k)},k=k+1$，转3。否则令$x^{(k+1)}=x^{(k)}+d^{(k)}，k=k+1$，转2。

注：4.中的常数1/4，3/4, 1/2是根据经验选取的。实际计算时，可根据问题对它们进行调整。

四、最小二乘问题

如下特殊形式的无约束问题：
$$minf(x)=\frac{1}{2}\sum _{i=1}^m{F}_i(x{)}^2$$
其中，$F_i:R^n\to R(i=1,2,\cdot \cdot \cdot ,m)$连续可微。称该问题为最小二乘问题。非线性最小二乘问题包含非线性方程组作为其特殊情形，即$m=n$。且该问题的最优解处的目标函数值为0。当$F_i(i=1,2,\cdot \cdot \cdot ,m)$都是线性函数时，该问题称为线性最小二乘问题。

1. 线性最小二乘问题

设$F_i(x)(i=1,2,\cdot \cdot \cdot ,m)$为线性函数：
$$F_i(x)=a_i^{T}x-b_i(i=1,2,\cdot \cdot \cdot ,m)$$
其中，$a_i\in R^n,b_i\in R(i=1,2,\cdot \cdot \cdot ,m)$。考虑如下线性最小二乘问题：
$$minf(x)=\frac{1}{2}\sum _{i=1}^m(a_i^{T}x-b_i{)}^2$$
容易证明，这是一个凸二次规划问题。（参考：数值最优化—概述）。

由无约束问题解的最优性条件，该问题等价于下面的线性方程组：
$$\nabla f(x)=\sum _{i=1}^m(a_i^{T}x-b_i)a_i=0$$
若记$A=(a_1,a_2,\cdot \cdot \cdot ,a_m{)}^T,b=(b_1,b_2,\cdot \cdot \cdot ,b_m{)}^T$，则上面的线性方程组可写为：
$$A^{T}Ax-A^{T}b=0$$
此线性方程组称为该问题的正规方程组。由于正规方程组的系数矩阵的秩与其增广矩阵的秩相等，因此方程组有解。

2. 非线性最小二乘问题

当问题中的至少有一个$F_i$是非线性函数，问题称为非线性最小二乘问题。设$F$二次连续可微。经直接计算可得：
$$\begin{gathered} \nabla f(x)=\sum _{i=1}^m{F}_i(x)\nabla F_i(x)=F^{\prime }(x{)}^{T}F(x) \\ {\nabla }^{2}f(x)=F^{\prime }(x{)}^T{F}^{\prime }(x)+\sum _{i=1}^m{F}_i(x){\nabla }^2{F}_i(x) \end{gathered}$$
其中，
$$F^{\prime }(x)=(\nabla F_1(x),\nabla F_2(x),\cdot \cdot \cdot ,\nabla F_m(x){)}^T$$
表示函数$F$在$x$处的Jacobi矩阵。

非线性最小二乘问题是一个无约束最优化问题，因此，我们可以利用求解无约束最优化问题的算法求解。注意到该问题的特殊性，我们可以针对该问题设计独特的算法。下面介绍求解非线性最小二乘问题的两种常用算法：Gauss-Newton法和Levenberg-Marquardt算法。

① Gauss-Newton法

求解非线性最小二乘问题的Gauss-Newton法中的$d^{(k)}$是下面的线性方程组的解：
$$F^{\prime }(x^{(k)}{)}^T{F}^{\prime }(x^{(k)})d+F^{\prime }(x^{(k)}{)}^T{F}^{\prime }(x^{(k)})=0$$
或等价地，$d^{(k)}$是线性最小二乘问题：
$$min\frac{1}{2}||F(x^{(k)})+F^{\prime }(x^{(k)})|{|}^2$$
的解。

由于$d^{(k)}$是$f$在$x$处的一个下降方向，因此，可以用下降算法构造求解非线性最小二乘问题的算法如下：

1. 给定初始点$x^{(0)}\in R^n$，常数$\rho \in (0,1),{\sigma }_1\in (0,1/2)$，精度$ϵ>0$，令$k=0$。
2. 若$||\nabla f(x^{(k)})|{|}^2\le ϵ$，则终止算法，得解$x^{(x)}$。否则，转3。
3. 解上面的线性方程组得方向$d^{(k)}$。
4. 确定步长${\alpha }_k$为集合 { ρ i   ∣   i = 0 , 1 , 2 , ⋅ ⋅ ⋅ } \{ \rho ^i \ | \ i=0,1,2,··· \} {ρi ∣ i=0,1,2,⋅⋅⋅}中使得下面的不等式成立的最大者：
   $$f(x^{(k)}+{\alpha }_k{d}^{(k)})\le f(x^{(k)})+{\sigma }_1{\alpha }_k\nabla f(x^{(k)}{)}^T{d}^{(k)}$$
5. 令$x^{(k=1)}=x^{(k)}+{\alpha }_k{d}^{(k)}，k=k+1$。转2。

② Levenberg-Marquardt算法

Levenberg-Marquardt算法是求解非线性最小二乘问题的另一种常用算法，该算法中$d^{(k)}$是下面的线性方程组的解：
$$[F^{\prime }(x^{(k)}{)}^T{F}^{\prime }(x^{(k)})+{\mu }_{k}I]d+F^{\prime }(x^{(k)}{)}^{T}F(x^{(k)})=0$$
其中，$I\in R^{n\times n}$表示单位矩阵，${\mu }_k>0$。若${\mu }_k=0$，则Levenberg-Marquardt算法还原为Gauss-Newton法。从这种意义上来看，Levenberg-Marquardt算法是一种修正的Gauss-Newton法。算法如下：

1. 给定初始点$x^{(0)}\in R^n$，常数$\rho \in (0,1),{\sigma }_1\in (0,1/2)$，精度$ϵ>0$，令$k=0$。
2. 若$||\nabla f(x^{(k)})||\le ϵ$，则终止算法，得解$x^{(x)}$。否则，转3。
3. 解上面的线性方程组得方向$d^{(k)}$。
4. 确定步长${\alpha }_k$为集合 { ρ i   ∣   i = 0 , 1 , 2 , ⋅ ⋅ ⋅ } \{ \rho ^i \ | \ i=0,1,2,··· \} {ρi ∣ i=0,1,2,⋅⋅⋅}中使得下面的不等式成立的最大者：
   $$f(x^{(k)}+{\alpha }_k{d}^{(k)})\le f(x^{(k)})+{\sigma }_1{\alpha }_k\nabla f(x^{(k)}{)}^T{d}^{(k)}$$
5. 令$x^{(k=1)}=x^{(k)}+{\alpha }_k{d}^{(k)}，k=k+1$。转2。

注：Levenberg-Marquardt算法与Gauss-Newton算法唯一的区别在于确定下降方向的子问题。