Gauss-Newton BFGS方法所产生的迭代矩阵序列的收敛性

The Convergence of Iterate Matrices Sequences Generated by Gauss-Newton BFGS Methods

摘要: 收敛速度的快慢是决定一个算法好坏的重要因素。在拟牛顿算法中，算法的收敛性在某种程度上等价于Dennis-Moré条件，但这并不意味着算法所产生的迭代矩阵就会收敛到Hessian矩阵。本文证明了由求解对称非线性方程组的Gauss-Newton BFGS方法所产生的迭代矩阵序列的收敛性，并通过数值实验对结论进行验证。

Abstract: The speed of convergence is an important factor that determines the quality of an algorithm. In the quasi-Newton algorithm, the convergence of the algorithm is equivalent to the Dennis-Moré condi-tion to some extent, but this does not mean that the iterative matrix generated by the algorithm will converge to the Hessian matrix. This paper proves the convergence of the iterative matrix se-quence generated by the Gauss-Newton BFGS method for solving symmetric nonlinear equations, and validates the conclusion by numerical experiments.

1. 引言

我考虑如下非线性方程组：

$g\left(x\right)=0$, $x\in R^n$.

其中g是 $R^n$ 到 $R^n$ 上的可微映射且g的雅可比阵 $\nabla g\left(x\right)$ 是对称的。在这种情况下我把g看作某个函数f从 $R^n$ 到R的梯度映射，从而 $g\left(x\right)=\nabla f\left(x\right)$， $\nabla g\left(x\right)={\nabla }^{2}f\left(x\right)$， $g\left(x\right)=0$ 就是无约束优化问题：

$\underset{x\in R^n}{\min }f\left(x\right)$.

的一阶必要条件。

DONGHUI LI and MASAO FUKUSHIMA在 [1] 中提出了处理此类问题的Gauss-Newton BFGS方法，其迭代格式为：

$x_{k+1}=x_k+{\lambda }_k{p}_k,$

$p_k:=-\frac{B_k^{-1}\left(g\left(x_k+{\lambda }_{k-1}{g}_k\right)-g_k\right)}{{\lambda }_{k-1}}$, $g_k:=g\left(x_k\right)$.

其中 $p_k$ 是搜索方向，步长 ${\lambda }_k>0$ 且使得序列 $\left\{g_k\right\}$ 在总体上具有近似范数下降性。矩阵 $B_k$ 由BFGS公式进行更新来确保它的正定性从而 $p_k$ 就是映射g在迭代点 $x_k$ 处的下降方向，因此 $B_k$ 满足割线方程：

$B_{k+1}{s}_k=y_k$,

其中

$s_k:=x_{k+1}-x_k$, $y_k:=g\left(x_k+{\delta }_k\right)-g_k$, ${\delta }_k=g_{k+1}-g_k$.

和常用的BFGS公式不同的是我经常用 $y_k$ 来表示梯度差而这里用 ${\delta }_k$ 表示，在这种情况下我有以下近似关系：

$y_k\approx \nabla g_{k+1}{\delta }_k\approx \nabla g_{k+1}\nabla g_{k+1}{s}_k$.

又因为 $B_{k+1}$ 满足割线方程且雅可比阵 $\nabla g_k$ 是对称的，我又有以下近似关系：

$B_{k+1}{s}_k\approx \nabla g_{k+1}\nabla g_{k+1}{s}_k\approx \nabla g_{k+1}\nabla g_{k+1}^{\text{T}}{s}_k$.

这就意味着 $B_{k+1}$ 沿 $s_k$ 方向近似于 $\nabla g_{k+1}\nabla g_{k+1}^{\text{T}}$，又因为

${\lambda }_{k-1}^{-1}\left(g\left(x_k+{\lambda }_{k-1}{g}_k\right)-g_k\right)\approx \nabla g_k{g}_k$,

所以我得到

$B_k{p}_k+\nabla g_k{g}_k\approx 0$.

因此我把 $p_k$ 看作Gauss-Newton的近似方向，把这种方法称为Gauss-Newton-based BFGS方法。该方法的优点在于不需要计算雅可比阵并且在合适的条件下具有全局收敛性和超线性收敛性。

我让函数f满足以下假设：

1) 函数f在 $R^n$ 上二阶连续可微。

2) $\exists M>0$，对 $\forall x\in {\Omega }_1$， $\Vert {\nabla }^{2}f\left(x\right)\Vert \le M$。

3) ${\nabla }^{2}f\left(x\right)$ 在 ${\Omega }_1$ 上一致非奇异。

其中 ${\Omega }_1$ 的定义见 [1] 中假设A。

在以上假设条件下，由 [1] 我可以得到以下结论：

4) $\exists x^{*}$ 使得 $g\left(x^{*}\right)=0$。

5) ${\displaystyle \underset{k}{\sum }\Vert x_k-x^{*}\Vert }<\infty$, ${\displaystyle \underset{k}{\sum }\Vert s_k\Vert }<\infty$。

6) 当 $k\to \infty$ 时， $B_{k+1}$ 总是由BFGS公式更新。

我进一步假设：

7) ${\nabla }^{2}f\left(x^{*}\right)$ 是正定的且 ${\nabla }^{2}f\left(x\right)$ 在 $x^{*}$ 处Lipschitz连续，即 $\Vert {\nabla }^{2}f\left(x\right)-{\nabla }^{2}f\left(x^{*}\right)\Vert \le L\Vert x-x^{*}\Vert$，其中x属于 $x^{*}$ 的一个邻域。

同样我从 [1] 中得到以下结论：

8) $\underset{k}{\sup }\Vert B_k\Vert <\infty$, $\underset{k}{\sup }\Vert B_k^{-1}\Vert <\infty$。

9) 当 $k\to \infty$ 时， ${\lambda }_k$ 恒等于1。

类似于常用的BFGS方法，超线性收敛性在某种程度上等价于Dennis-Moré条件 [2]，但由于 $y_k$ 的取不同，条件的表现形式也有所改变，具体如下：

$\underset{k\to \infty }{\lim }\frac{\Vert \left(B_k-\nabla g_{*}{\nabla }^{\text{T}}{g}_{*}\right)s_k\Vert }{\Vert s_k\Vert }=0$.

但这并不意味着 $B_k$ 就一定最终等于 $\nabla g_{*}{\nabla }^{\text{T}}{g}_{*}$ (由 [2] 中的反例我可以得到这个结论)。最早证明相关结论的是Ge Ren-Pu and Powell [3] 中证明了DFP方法和BFGS方法取恒定步长为一所产生的迭代矩阵序列的收敛性，该证明不要求 $B_k$ 最终等于 ${\nabla }^{2}f\left(x^{*}\right)$，之后Stoer [4] 把该结论推广到Broyden族方法上且步长最终收敛到1即可。由于1)到9)的结论满足 [4] 的假设B，所以本文的证明参照 [4] 进行了一点点的改动，以下我统称为假设最后我给出Gauss-Newton-based BFGS算法的框架。

Gauss-Newton-based BFGS算法

Step 0. Choose an initial point $x_0\in R^n$, an initial symmetric positive definite matrix $B_0\in R^{n\text{×}n}$, a positive

sequence $\left\{{\omega }_k\right\}$ satisfying ${\displaystyle \underset{k=0}{\overset{\infty }{\sum }}{\omega }_k}<\infty$, and constants $r,\rho \in \left(0,1\right)$, ${\sigma }_1,{\sigma }_2>0,{\lambda }_{-1}>0$. Let $k:=0$.

Step 1. Stop if $g_k=0$. Otherwise, solve the following linear equation to get $p_k$ :

$B_{k}p+{\lambda }_{k-1}^{-1}\left(g\left(x_k+{\lambda }_{k-1}{g}_k\right)-g_k\right)=0$.

Step 2. If

$\Vert g\left(x_k+p_k\right)\Vert \le \rho \Vert g_k\Vert$,

then take ${\lambda }_k=1$ and go to Step 4. Otherwise go to Step 3.

Step 3. Let $i_k$ be the smallest nonnegative integer i such that:

${\Vert g\left(x_k+\lambda p_k\right)\Vert }^2-{\Vert g_k\Vert }^2\le -{\sigma }_1{\Vert \lambda g_k\Vert }^2-{\sigma }_2{\Vert \lambda p_k\Vert }^2+{\omega }_k{\Vert g_k\Vert }^2$ holds for $\lambda =r^i$.

Let ${\lambda }_k=r^{i_k}$.

Step 4. Let the next iterate be $x_{k+1}=x_k+{\lambda }_k{p}_k$.

Step 5. Put $s_k=x_{k+1}-x_k={\lambda }_k{p}_k$, ${\delta }_k=g_{k+1}-g_k$, and $y_k=g\left(x_k+{\delta }_k\right)-g\left(x_k\right)$. If $y_k^{\text{T}}{s}_k\le 0$, then $B_{k+1}=B_k$ and go to Step 6. Otherwise, update $B_k$ by the BFGS formula:

$B_{k+1}=B_k-\frac{B_k{s}_k{s}_k^{\text{T}}{B}_k}{s_k^{\text{T}}{B}_k{s}_k}+\frac{y_k^{\text{T}}{y}_k}{s_k^{\text{T}}{y}_k}$.

Step 6. Let $k:=k+1$. Go to Step 1.

2. 收敛性的证明

由BFGS方法的不变性，我假设 ${\nabla }^{2}f\left(x^{*}\right){\nabla }^{2}f{\left(x^{*}\right)}^{\text{T}}=I$。如果函数 $f\left(x\right)$ 是二次函数，那么 $g\left(x\right)$ 就是线性的，那么此时

$p_k=-B_k^{-1}g$, $y_k=g_{k+1}-g_k={\delta }_k$,

我直接由 [3] 得到结论，对于一般情况我有

$\begin{array}{c}{y}_k=g\left(x_k+{\delta }_k\right)-g_k={\displaystyle {\int }_0^1{\nabla }^{2}f\left(x_k+\theta {\delta }_k\right)\text{d}\theta {\delta }_k}\\ ={\displaystyle {\int }_0^1{\nabla }^{2}f\left(x_k+\theta {\delta }_k\right)\text{d}\theta }{\displaystyle {\int }_0^1{\nabla }^{2}f\left(x_k+\theta s_k\right)\text{d}\theta s_k}=Q_k{P}_k{s}_k\end{array}$

其中 $Q_k={\displaystyle {\int }_0^1{\nabla }^{2}f\left(x_k+\theta {\delta }_k\right)\text{d}\theta }$, $P_k={\displaystyle {\int }_0^1{\nabla }^{2}f\left(x_k+\theta s_k\right)\text{d}\theta }$。

有假设得知当k充分大时， $P_k$ 是正定矩阵且

$P_k={\rm I}+O\left(\Vert s_k\Vert \right)$,

同理当k充分大时

${\delta }_k=g_{k+1}-g_k={\displaystyle {\int }_0^1{\nabla }^{2}f\left(x_k+t{s}_k\right)\text{d}\theta s_k}=s_k+O\left({\Vert s_k\Vert }^2\right)$,

$y_k=\left[{\rm I}+O\left(\Vert {\delta }_k\Vert \right)\right]\left[{\rm I}+O\left(\Vert s_k\Vert \right)\right]\text{d}\theta s_k=s_k+O\left({\Vert s_k\Vert }^2\right)$.

因此由BFGS公式我有

$\begin{array}{c}{B}_{k+1}=B_k-\frac{B_k{s}_k{s}_k^{\text{T}}{B}_k}{s_k^{\text{T}}{B}_k{s}_k}+\frac{y_k^{\text{T}}{y}_k}{s_k^{\text{T}}{y}_k}\\ =B_k-\frac{B_k{s}_k{s}_k^{\text{T}}{B}_k}{s_k^{\text{T}}{B}_k{s}_k}+\frac{s_k^{\text{T}}{s}_k}{s_k^{\text{T}}{s}_k}+O\left(\Vert s_k\Vert \right)\\ ={B^{\prime }}_{k+1}+O\left(\Vert s_k\Vert \right)\end{array}$ (1)

我定义 $E_k:=B_k-{\rm I}$， ${\lambda }_{ki},i=1,2,\cdots ,n$ 是它的特征值且 $\left|{\lambda }_{k1}\right|\le \left|{\lambda }_{k2}\right|\le \cdots \le \left|{\lambda }_{kn}\right|$， ${\upsilon }_{ki},i=1,2,\cdots ,n$ 是特征值对应的正交特征向量， $E_k{\upsilon }_{ki}={\lambda }_{ki}{\upsilon }_{ki}$。由 [5] 中的结论我得到对于 $E_{k+1}:=B_{k+1}-{\rm I}$ 的特征值 $\left\{{\lambda }_{k+1,i}\right\}$ 按照同样的方式排列 $\left|{\lambda }_{k+1,1}\right|\le \left|{\lambda }_{k+1,2}\right|\le \cdots \le \left|{\lambda }_{k+1,n}\right|$，它满足

$\left|{\lambda }_{k+1,i}\right|\le \left|{\lambda }_{k,i}\right|$,

$sign{\lambda }_{k+1,i}=sign{\lambda }_{k,i}$, $i=1,2,\cdots ,n$.

由(1)以及 ${\displaystyle \underset{k}{\sum }\Vert s_k\Vert }$ 收敛，根据文献 [3] 中引理4可知对每个 $i=1,2,\cdots ,n$，极限 $\underset{k\to \infty }{\lim }{\lambda }_{ki}$ 存在。我假设极限趋于0的特征值个数为m，其余的都大于 $\beta >0$，即

$\underset{k}{\lim }{\lambda }_{ki}=0$, $i=1,2,\cdots ,n$

$\left|{\lambda }_{ki}\right|\ge \beta$, $i=m+1,\cdots ,n$.

又因为

$\frac{s_k^{\text{T}}{E}_{k+1}{s}_k}{s_k^{\text{T}}{s}_k}=\frac{s_k^{\text{T}}\left(y_k-s_k\right)}{s_k^{\text{T}}{s}_k}=\frac{s_k^{\text{T}}\left(Q_k{P}_k-{\rm I}\right)s_k}{s_k^{\text{T}}{s}_k}=O\left(\Vert s_k\Vert \right)$

而 $\underset{k}{\lim }\Vert s_k\Vert =0$，所以 $m\ge 1$。

由 [3] 我对 $E_k$ 进行如下分解

$E_k:={\Delta }_k+H_k$

其中

${\Delta }_k:={\displaystyle \underset{i=1}{\overset{m}{\sum }}{\lambda }_{ki}}{\upsilon }_{ki}{\upsilon }_{ki}^{\text{T}}$, $H_k:={\displaystyle \underset{i=m+1}{\overset{n}{\sum }}{\lambda }_{ki}}{\upsilon }_{ki}{\upsilon }_{ki}^{\text{T}}$,

$S_k:=span\left\{{\upsilon }_{k1},{\upsilon }_{k2},\cdots ,{\upsilon }_{km}\right\}$.

那么 $S_k^{\perp }=span\left\{{\upsilon }_{k,m+1},\cdots ,{\upsilon }_{k,n}\right\}$。

由以上定义我有 ${\Delta }_k{S}_k\subset S_k$， $H_k{S}_k^{\perp }\subset S_k^{\perp }$， ${\Delta }_k{S}_k^{\perp }=H_k{S}_k=0$，根据m的定义有 $\underset{k}{\lim }{\Delta }_k=0$。

通过以上讨论，我将通过证明 ${\displaystyle \underset{k}{\sum }\Vert H_{k+1}-H_k\Vert }<\infty$ 来证明序列 $\left\{B_k\right\}$ 的收敛性。我将证明

${\Vert H_{k+1}-H_k\Vert }_F\le c_1\Vert {\eta }_k\Vert /\Vert s_k\Vert +c_2\Vert s_k\Vert$ (2)

对充分大的k都成立。这里 ${\eta }_k$ 取自 $s_k$ 的正交分解

$s_k={\gamma }_k+{\eta }_k,{\gamma }_k\in S_k,{\eta }_k\in S_k^{\perp }$.

(2)式的证明就是 [3] 中对DFP方法的证明，这里不再赘述。

因为 ${\displaystyle \underset{k}{\sum }\Vert s_k\Vert }<\infty$，由(2)我只需证明 ${\displaystyle \underset{k}{\sum }\Vert {\eta }_k\Vert /\Vert s_k\Vert }<\infty$。首先证明一个重要不等式

${\displaystyle \underset{k}{\sum }\frac{{\Vert E_k{s}_k\Vert }^2}{{\Vert s_k\Vert }^2}}<\infty$ (3)

$\begin{array}{c}{E}_{k+1}=B_{k+1}-{\rm I}={B^{\prime }}_{k+1}-{\rm I}+O\left(\Vert s_k\Vert \right)\\ =E_k-\frac{\left({\rm I}+E_k\right)s_k{s}_k^{\text{T}}\left({\rm I}+E_k\right)}{{\Vert s_k\Vert }^2+s_k^{\text{T}}{E}_k{s}_k}+\frac{s_k{s}_k^{\text{T}}}{{\Vert s_k\Vert }^2}+O\left(\Vert s_k\Vert \right)\end{array}$

由文献 [6] 的定理2我有

${\Vert E_{k+1}\Vert }_F^2\le {\Vert E_k\Vert }_F^2-\frac{{\left(E_k{s}_k\right)}^{\text{T}}{B}_k\left(E_k{s}_k\right)}{s_k^{\text{T}}{B}_k{s}_k}+O\left(\Vert s_k\Vert \right)$,

不等式两边对k进行累加，又因为 $\left\{B_k\right\}$ 和 $\left\{B_k^{-1}\right\}$ 都是正定有界的，(3)式得证。从文献 [4] 我立马又得到 ${\displaystyle \underset{k}{\sum }\Vert B_{k+1}-B_k\Vert }<\infty$，证明没有任何改变。

接下来我证明 ${\displaystyle \underset{k}{\sum }\Vert {\eta }_k\Vert /\Vert s_k\Vert }<\infty$，由 $k\to \infty$ 时 ${\lambda }_k\equiv 1$ 得，存在常数 $k_1$，当 $k\ge k_1$ 时， ${\lambda }_k=1$。下面我证明几个不等式，其中 $\bar{E}=B_k-Q_k{P}_k$。

a) ${\displaystyle \underset{k_1}{\sum }\frac{{\Vert E_k{\bar{E}}_k{s}_k\Vert }^2}{{\Vert {\bar{E}}_k{s}_k\Vert }^2}}<\infty$,

b) ${\displaystyle \underset{k_1}{\sum }\frac{{\Vert {\bar{E}}_k{s}_k\Vert }^2}{{\Vert s_k\Vert }^2}}<\infty$,

c) ${\displaystyle \underset{k_1}{\sum }\frac{\Vert E_k\bar{E}{s}_k\Vert }{\Vert s_k\Vert }}<\infty$,

d) ${\displaystyle \underset{k_1}{\sum }\frac{\Vert E_k^2{s}_k\Vert }{\Vert s_k\Vert }}<\infty$.

证明和 [4] 中类似但是有两个小变化。一个是当 $k\ge k_1$ 时步长 ${\lambda }_k=1$，这时

$p_k=s_k=-B_k^{-1}\left(g\left(x_k+g_k\right)-g_k\right)=-B_k^{-1}{\displaystyle {\int }_0^1{\nabla }^{2}f\left(x_k+l{g}_k\right)\text{d}l{g}_k}=-B_k^{-1}{G}_k{g}_k$

另一个是

$\begin{array}{c}{\bar{E}}_k{s}_k=B_k{s}_k-y_k=-G_k{g}_k-Q_k{\delta }_k\\ =\left(Q_k-G_k\right)g_k-Q_k{g}_{k+1}\\ =\left(Q_k-G_k\right)g_k+Q_k{G}_{k+1}^{-1}{B}_{k+1}{s}_{k+1}\end{array}$

当k充分大时，由于函数f二阶可微可得 $Q_k-G_k\to 0$，所以上式最后一个等式的第一部分趋于0，由假设2)，3)以及 $B_k$ 的有界性得证。

现在我的结论显而易见，由 $E_k:={\Delta }_k+H_k$, $s_k={\gamma }_k+{\eta }_k$ 以及 $\left|{\lambda }_{ki}\right|\ge \beta ,i=m+1,\cdots ,n$ 我有

$\frac{\Vert E_k^2{s}_k\Vert }{\Vert s_k\Vert }=\frac{\Vert {\Delta }_k^2{\gamma }_k+H_k^2{\eta }_k\Vert }{\Vert s_k\Vert }\ge \frac{\Vert H_k^2{\eta }_k\Vert }{\Vert s_k\Vert }\ge {\beta }^2\frac{\Vert {\eta }_k\Vert }{\Vert s_k\Vert }$

所以 ${\displaystyle \underset{k}{\sum }\Vert {\eta }_k\Vert /\Vert s_k\Vert }<\infty$，得证。

3. 数值实验

在这部分我通过画图来观察 $\Vert B_k-\nabla g_{*}{\nabla }^{\text{T}}{g}_{*}\Vert$ 的收敛性，我在三个不同的方面进行比较：初始点、维度、和精度。我把 $\nabla g_{*}$ 记作j，这里所取的问题就是 [1] 中的原问题，是为了验证自己的编程是否正确。

问题1

$g\left(x\right)\triangleq Ax+\frac{1}{{\left(n+1\right)}^2}F\left(x\right)=0$,

其中 $A\in R^{n\text{×}n}$

$A=\left(\begin{array}{ccccc}2& -1& & & \\ -1& 2& -1& & \\ & \ddots & \ddots & \ddots & \\ & & \ddots & \ddots & -1\\ & & & -1& 2\end{array}\right)$,

$F\left(x\right)={\left(F_1\left(x\right),F_2\left(x\right),\cdots ,F_n\left(x\right)\right)}^{\text{T}}$

$F_i\left(x\right)=\sin x_i-1$, $i=1,2,\cdots ,n$.

以下是不同初始点(见表1)的影响(见图1)，维度n = 19。

Table 1. Different initial points and their manifestations

表1. 不同初始点及其表现形式

Figure 1. Performance of different starting points

图1. 不同初始点的表现

以下是不同维度(见表2)的影响(见图2)，初始点为全1向量。

Table 2. Different dimensions and their manifestations

表2. 不同维度及其表现形式

Figure 2. Performance of different dimensions

图2. 不同维度的表现

以下是不同精度(见表3)的影响(见图3)，维度n = 19，初始点设置为全1向量。

Table 3. Precision

表3. 精度

Figure 3. Performance of different precision

图3. 不同精度的表现

在实验中对应的参数为 $r=0.1$, $\rho =\sqrt{0.9}$, ${\sigma }_1={\sigma }_2={10}^{-5}$, ${\lambda }_{-1}=0.01$, and ${\omega }_k=k^{-2}$，初始矩阵 $B_0$ 为单位阵。在前两个实验中精度固定为10⁻⁵，第二个实验我选取全1向量的原因是因为问题是对称问题，在这种情况下初始点的运动轨迹是一样的。而在最后一个实验中由于精度的提高我们将MATLAB程序显示的有效数字调整为15位来得到更好的准确性。从这三张图我们可以清晰的看出当迭代点趋近最小点时函数图像基本保持水平，所以结论是有效的。

4. 总结

本文证明了处理特殊问题的Gauss-Newton-based BFGS方法所产生的迭代矩阵序列的收敛性，由此可以猜测这种收敛性是公式本身的性质，这也启发我如果在应用相应公式做算法时，如果算法不具有超线性收敛性，是否可以通过修正 $B_k$ 使他倾向于Hessian阵从而具有超线性收敛性。