Numerical Optimization Ch6. Quasi-Newton Methods

第六章: 拟牛顿法

拟牛顿法由物理学家W.C.Davidson于上世纪五十年代中期首创, 并由Fletcher和Powell证实有效性, 从而一夜之间转变了非线性优化的格局. 在接下来的二十年里, 遂出现了大量的变体.
第二章指出, 拟牛顿法实质上就是改变目标函数在局部的二次近似模型, 具体地, 就是以对称正定阵 $B_k$代替原有的Hessian, 并在迭代中 不断更新. 第三章则指出拟牛顿法的 超线性收敛速度. 由于它的实施过程并不涉及Hessian, 因此在某些问题上拟牛顿法比牛顿法要更加方便有效. 本章主要讨论拟牛顿法在 中小型问题上的应用, 对于大型问题我们将讨论放在下一章.
值得一提的是, 如今Hessian并不需要显式地获取. 我们可以用后面将要介绍的自动微分技术 (automatic differentiaton techniques), 仅用函数值或一阶信息就能在一定精度下计算之. 不过 自动微分在许多场景下仍然是无法应用的, 且二阶信息在实际计算上并不比一阶信息优越. 因此拟牛顿法仍然广受欢迎.

1. BFGS算法

1.1 算法介绍

最受欢迎的拟牛顿法要数BFGS算法, 其命名为发明者的首字母缩写: Broyden-Fletcher-Goldfarb-Shanno. 我们在第二章中曾经提及DFP算法 (另一种拟牛顿法)并简要概述了BFGS算法. 下面我们来详细讨论后者的构成、理论性质和实施要点.
首先是拟牛顿法对目标函数的二次近似: 设当前迭代点为$x_k$, 设二次近似模型为$$m_k(p)=f_k+\nabla f_k^{T}p+\frac{1}{2}{p}^T{B}_{k}p.$$这里$B_k$为$n\times n$的对称正定矩阵, 它在每一步迭代都会被更新. 对这一凸二次模型, 极小点$p_k$具有显式表达式:$$p_k=-B_k^{-1}\nabla f_k,$$从而新的迭代点为$$x_{k+1}=x_k+{\alpha }_k{p}_k,$$其中步长${\alpha }_k$的选取满足Wolfe条件.
下面我们来具体讨论$B_k$的更新方法. 假设在$x_{k+1}$处有新的二次近似模型:$$m_{k+1}(p)=f_{k+1}+\nabla f_{k+1}^{T}p+\frac{1}{2}{p}^T{B}_{k+1}p.$$对于$B_{k+1}$我们应当加什么样的约束呢? 比如:

• $B_{k+1}$的对称性.
• $m_{k+1}$的梯度应当等于目标函数$f$在最近的两个迭代点$x_k,x_{k+1}$上的梯度. 由于$m_{k+1}(0)$已经是$\nabla f_{k+1}$, 因此只需$$\nabla m_{k+1}(-{\alpha }_k{p}_k)=\nabla f_{k+1}-{\alpha }_k{B}_{k+1}{p}_k=\nabla f_k$$$$\Rightarrow B_{k+1}{\alpha }_k{p}_k=\nabla f_{k+1}-\nabla f_k.$$记$$s_k=x_{k+1}-x_k={\alpha }_k{p}_k,y_k=\nabla f_{k+1}-\nabla f_k,$$则之前的条件变成$$B_{k+1}{s}_k=y_k.$$我们将此称为割线方程 (secant equation), 也有称为拟牛顿方程. 注意到$B_{k+1}$是对称正定阵, 因此割线方程隐含条件$s_k^T{y}_k>0$. 事实上当$f$是强凸函数时, 这一隐含条件是自动成立的. 而对非凸函数这一点并不总是对的. 因此我们有必要通过线搜索的方式强加之隐含条件. 不过这在(强)Wolfe条件下已能保证: 由曲率条件, 我们有$$\nabla f_{k+1}^T{s}_k\ge c_2\nabla f_k^T{s}_k$$$$\Rightarrow y_k^T{s}_k\ge (c_2-1){\alpha }_k\nabla f_k^T{p}_k.$$由于$c_2<1$, $p_k$为下降方向 (这是由于$B_k$正定, 后面会说到这可以通过要求初始$B_0$正定获得) , 因此右端为正数, 隐含条件成立.
  注意隐含条件成立的情况下, 割线方程的解仍然可能有无穷个 (我们后面会证明满足隐含条件的正定矩阵是存在的). 这是因为对称$B$的自由度为$n(n+1)/2$, 而隐含条件至多可以去掉$n$个自由度. 我们当然可以通过$B$的正定性得到其$n$个主子式为正, 从而加上$n$个不等式约束, 但这仍然不能保证唯一解. 因此我们还要求
• $B_{k+1}$在某种意义下与当前的$B_k$是最接近的. 这也称为是最小修正原理.

综合以上, 我们就是要求解下面的问题:$$\begin{array}{ll}{\min }_B& \Vert B-B_k\Vert \\ \mathrm{s}.\mathrm{t}.& B=B^T,B{s}_k=y_k,\end{array}$$其中$\Vert \cdot \Vert$为某种范数, 而不同的选择会导致不同的你牛顿方法. 其中最易于求解的且可以导出尺度不变方法的范数为加权Frobenius范数:$$\Vert A{\Vert }_W=\Vert W^{1/2}A{W}^{1/2}{\Vert }_F,$$其中$\Vert \cdot {\Vert }_F$定义为$\Vert C{\Vert }_F^2={\sum }_{i=1}^n{\sum }_{j=1}^n{c}_{ij}^2$, $W$选取为满足关系式$W{y}_k=s_k$的加权矩阵 (这一点保证了尺度不变性). 事实上, $W$有一个天然的选取方式, 即$W={\bar{G}}_k^{-1}$, 其中${\bar{G}}_k$为平均Hessian, 定义为$${\bar{G}}_k={\int }_0^1{\nabla }^{2}f(x_k+\tau {\alpha }_k{p}_k)\mathrm{d}\tau .$$则由Taylor定理, 就有$$y_k={\bar{G}}_k{\alpha }_k{p}_k={\bar{G}}_k{s}_k.$$以上优化问题具有唯一解为$$(\mathrm{D}\mathrm{F}\mathrm{P})B_{k+1}=(I-{\rho }_k{y}_k{s}_k^T)B_k(I-{\rho }_k{s}_k{y}_k^T)+{\rho }_k{y}_k{y}_k^T,$$其中$${\rho }_k=\frac{1}{y_k^T{s}_k}.$$以上更新$B_k$的公式称为DFP公式, 其推导过程在此略去. 之后我们还会介绍BFGS公式并推导之. DFP公式的推导则是类似的. 利用Sherman-Morrison-Woodbury公式 (Low Rank Update), 我们还可以直接更新$B_k$的逆, 从而省去了计算$B_k^{-1}{p}_k$的步骤而只需进行矩阵-向量乘积. 将$B_k$的逆记为$$H_k=B_k^{-1}.$$由DFP公式和Sherman-Morrison-Woodbury公式, 我们就有$$(\mathrm{D}\mathrm{F}\mathrm{P})H_{k+1}=H_k-\frac{H_k{y}_k{y}_k^T{H}_k}{y_k^T{H}_k{y}_k}+\frac{s_k{s}_k^T}{y_k^T{s}_k}.$$注意到DFP公式的后两项都是秩为1的矩阵, 因此DFP公式可以看做是对$H_k$的秩2修正 (类似地, 对$B_k$也是如此). 这就反映了拟牛顿法的本质思想: 与其在每一步耗时耗力地计算$H_k$($B_k$), 我们直接用秩为2的项对其进行修正, 从而大幅度减少计算量. 与此同时, 关于目标函数最新的信息也就集中在这一秩为2的项中.

DFP公式无疑是高效的, 不过它很快就被BFGS公式所超越. 后者的想法与前者的唯一不同在于: 后者是基于$H_{k+1}$与$H_k$的差距不大而导出的, 而前者是基于$B_{k+1}$与$B_k$的差距不大而导出. 尽管句法结构上极其相似, 但实际上确是可能做着完全不同的事情——就算$B_{k+1}$与$B_k$及其接近, 也不免$H_{k+1}$与$H_k$相距甚远. 这与矩阵的条件数相关. 事实证明, 这样的想法能够带来更佳的效果, 从而也使BFGS公式成为最有效的拟牛顿更新公式. 下面我们来推导这一公式:
我们现在的问题转为$$\begin{array}{ll}{\min }_H& \Vert H-H_k{\Vert }_W\\ \mathrm{s}.\mathrm{t}.& H=H^T,H{y}_k=s_k.\end{array}$$注意此时的$H_{k+1}$应当满足的割线方程是$H_{k+1}{y}_k=s_k$; 目标函数中的加权矩阵$W$则应满足$W{s}_k=y_k$ (这取之前的${\bar{G}}_k$即可). 由于皆为等式约束, 因此可用Lagrange乘子法求其稳定点. 为方便起见, 这里短暂记$H_k=H^{\ast }=(h_1^{\ast },h_2^{\ast },\dots ,h_n^{\ast })$. 设$H=(H_{ij})=(h_1,h_2,\dots ,h_n),W=(w_{ij})=(w_1,w_2,\dots ,w_n)$. 则$$\begin{array}{rl}\Lambda & =\sum _{ij}(w_i\cdot (h_j-h_j^{\ast }))(w_j\cdot (h_i-h_i^{\ast }))-\sum _{i<j}{\alpha }_{ij}(H_{ij}-H_{ji})-{\lambda }^T(H{y}_k-s_k)\\ & =\sum _{ij}(w_i\cdot (h_j-h_j^{\ast }))(w_j\cdot (h_i-h_i^{\ast }))-\sum _{ij}{\alpha }_{ij}{H}_{ij}-{\lambda }^T(H{y}_k-s_k)(定义{\alpha }_{ij}=-{\alpha }_{ji}).\end{array}$$这里$A=({\alpha }_{ij}),\lambda$均为Lagrange乘子. 求偏导并置0可得 (以下记$y_k=y,s_k=s$)$$\begin{array}{rl}0=\frac{\partial \Lambda }{\partial H_{ij}}& =\sum _{l}2w_{il}(w_j\cdot (h_l-h_l^{\ast }))-{\alpha }_{ij}-{\lambda }_i{y}_j\\ & =2\sum _l{w}_{il}(W^T(H-H^{\ast }){)}_{jl}-{\alpha }_{ij}-{\lambda }_i{y}_j\\ & =2\sum _l(W^T(H-H^{\ast }){)}_{jl}{w}_{li}-{\alpha }_{ij}-{\lambda }_i{y}_j(\because W的对称性)\\ & =2(W^T(H-H^{\ast })W{)}_{ji}-{\alpha }_{ij}-{\lambda }_i{y}_j\\ & =2(W(H-H^{\ast })W{)}_{ij}-{\alpha }_{ij}-{\lambda }_i{y}_j(\because W和H的对称性).\end{array}$$写成矩阵形式就得到以下式子:$$\begin{array}{rl}0& =2W(H-H^{\ast })W-A-\lambda y^T,\\ A^T& =-A,\\ W^T& =W,\\ H^T& =H,\\ (H^{\ast }{)}^T& =H^{\ast },\\ Hy& =s,\\ Ws& =y.\end{array}$$我们可以利用转置和$H,W$的对称性和$A$的反对称性得到以下关系式:$$\begin{array}{rl}0& =2W(H-H^{\ast })W-A-\lambda y^T,\\ 0& =2W(H-H^{\ast })W+A-y{\lambda }^T,\\ \Rightarrow 0& =4W(H-H^{\ast })W-\lambda y^T-y{\lambda }^T.\end{array}$$我们将利用此式计算$H$, 为此先计算$\frac{\lambda y^T+y{\lambda }^T}{4}$. 上式左右内积$s$得到:$$0=4(y-W{H}^{\ast }y)-\lambda (y\cdot s)-y(\lambda \cdot s).$$再左乘$s^T$得到:$$0=4(y\cdot s)-4(y^T{H}^{\ast }y)-2(y\cdot s)(\lambda \cdot s).$$这表明$$\lambda \cdot s=2\rho y^T(s-H^{\ast }y),其中\rho =1/y\cdot s.$$再把这个代入我们的向量等式:$$\begin{array}{rl}0& =4(y-W{H}^{\ast }y)-\lambda (y\cdot s)-y(\lambda \cdot s)\\ & =4(y-W{H}^{\ast }y)-\lambda (y\cdot s)-y[2\rho y^T(s-H^{\ast }y)]\\ \Rightarrow \lambda & =4\rho (y-W{H}^{\ast }y)-2{\rho }^2[y^T(s-H^{\ast }y)]y.\end{array}$$右乘$y^T$得:$$\lambda y^T=4\rho (y-W{H}^{\ast }y)y^T-2{\rho }^2[y^T(s-H^{\ast }y)]y{y}^T.$$转置:$$y{\lambda }^T=4\rho y(y^T-y^T{H}^{\ast }W)-2{\rho }^2[y^T(s-H^{\ast }y)]y{y}^T.$$结合上面两式:$$\frac{1}{4}(\lambda y^T+y{\lambda }^T)=\rho (2y{y}^T-W{H}^{\ast }y{y}^T-y{y}^T{H}^{\ast }W)-{\rho }^2[y^T(s-H^{\ast }y)]y{y}^T.$$对上式左乘和右乘$W^{-1}$. 注意$Ws=y\Rightarrow s=W^{-1}y\Rightarrow y^T{W}^{-1}=s$. 因此$$\begin{array}{rl}\frac{1}{4}{W}^{-1}(\lambda y^T+y{\lambda }^T)W^{-1}& =2\rho s{s}^T-\rho H^{\ast }y{s}^T-\rho s{y}^T{H}^{\ast }-{\rho }^2(y^{T}s)s{s}^T+{\rho }^2(y^T{H}^{\ast }y)s{s}^T\\ & =2\rho s{s}^T-\rho H^{\ast }y{s}^T-\rho s{y}^T{H}^{\ast }-\rho s{s}^T+s{\rho }^2(y^T{H}^{\ast }y)s^T(由\rho 的定义)\\ & =\rho s{s}^T-\rho H^{\ast }y{s}^T-\rho s{y}^T{H}^{\ast }+s{\rho }^2(y^T{H}^{\ast }y)s^T.\end{array}$$最后得到BFGS公式:$$\begin{array}{rl}0& =4W(H-H^{\ast })W-\lambda y^T-y{\lambda }^T\\ \Rightarrow H& =\frac{1}{4}{W}^{-1}(\lambda y^T+y{\lambda }^T)W^{-1}+H^{\ast }\\ & =\rho s{s}^T-\rho H^{\ast }y{s}^T-\rho s{y}^T{H}^{\ast }+s{\rho }^2(y^T{H}^{\ast }y)s^T+H^{\ast }(由上一段)\\ & =H^{\ast }(I-\rho y{s}^T)+\rho s{s}^T-\rho s{y}^T{H}^{\ast }+(\rho s{y}^T)H^{\ast }(\rho y{s}^T)\\ & =H^{\ast }(I-\rho y{s}^T)+\rho s{s}^T-\rho s{y}^T{H}^{\ast }(I-\rho y{s}^T)\\ & =\rho s{s}^T+(I-\rho s{y}^T)H^{\ast }(I-\rho y{s}^T).\end{array}$$用一个式子表示:$$(\mathrm{B}\mathrm{F}\mathrm{G}\mathrm{S})H_{k+1}=(I-{\rho }_k{s}_k{y}_k^T)H_k(I-{\rho }_k{y}_k{s}_k^T)+{\rho }_k{s}_k{s}_k^T.$$类似地, 结合Sherman-Morrison-Woodbury公式, 我们就有BFGS公式的$B_k$版本:$$(\mathrm{B}\mathrm{F}\mathrm{G}\mathrm{S})B_{k+1}=B_k-\frac{B_k{s}_k{s}_k^T{B}_k}{s_k^T{B}_k{s}_k}+\frac{y_k{y}_k^T}{y_k^T{s}_k}.$$尽管DFP公式与BFGS公式都有两种形式的表示, 但我们一般默认DFP指的是对$H$的公式, BFGS指的是对$B$的公式. 当然这里遗留下来的问题就包括初始矩阵$H_0$如何选取. 我们将在后面讨论之. 结合上述, 我们得到算法1.

算法1 (BFGS Method)
Given starting point $x_0$, convergence tolerance $ϵ>0$, inverse Hessian approximation $H_0$;
$k\leftarrow 0$;
while $\Vert \nabla f_k\Vert >ϵ$
\quad\quad Compute search direction $p_k=-H_k\nabla f_k$;
\quad\quad Set $x_{k+1}=x_k+{\alpha }_k{p}_k$ where ${\alpha }_k$ is computed from a line search procedure to satisfy the Wolfe conditions;
\quad\quad Define $s_k=x_{k+1}-x_k$ and $y_k=\nabla f_{k+1}-\nabla f_k$;
\quad\quad Compute $H_{k+1}$ by means of BFGS formula;
$k\leftarrow k+1$;
end (while)

不难注意到, 算法的每一步迭代只需要$O(n^2)$计算量. 且由之前第三章的分析, 我们已经知道BFGS算法具有超线性收敛速度. 相较牛顿法而言, 尽管速度上由劣势, 但其计算耗费足以相抵.
需要指出的是, 尽管直接应用BFGS公式的$B_k$版本会导致需要求解$B_k{p}_k=-\nabla f_k$, 从而有$O(n^3)$的计算量, 但后面我们会提到如何利用更新$B_k$的Cholesky因子来在$O(n^2)$内求解之.
需要额外提及的是, 若在最小修正中使用$F$-范数且不要求对称, 则会得到另一种更新公式:$${\bar{B}}_{+}=B+\frac{(y-Bs)s^T}{s^{T}s}.$$注意到这是一种秩1修正, 其特点在于最小限度地改变了原矩阵的特征子空间. 我们很容易就能验证上式的正确性:$$\Vert {\bar{B}}_{+}-B{\Vert }_F={\Vert \frac{(y-Bs)s^T}{s^{T}s}\Vert }_F={\Vert (B_{+}-B)\frac{s{s}^T}{s^{T}s}\Vert }_F\le \Vert B_{+}-B{\Vert }_F,\forall B_{+}:B_{+}s=y.$$不保持对称性的更新有用吗? 如果仍然将视野局限在单个目标函数的优化问题上, 似乎没什么用处. 但若放眼于多目标函数优化或者是非线性方程组求解, 由于此时Jacobi矩阵未必对称, 上述公式便有了实践意义.

1.2 算法性质

之前已经指出拟牛顿法的超线性收敛性质. 下面用一个具体的实例展示其效用. 对Rosenbrock函数$$f(x)=100(x_2-x_1^2{)}^2+(1-x_1{)}^2,$$以$(-1.2,1)$为初始点, 实施最速下降法、BFGS算法、非精确牛顿法(后面在大规模环境下会提到), 其中步长选取满足Wolfe条件. 实验表明: 为将梯度范数减小到$10^{-5}$以下, 最速下降法需要5264步迭代, BFGS算法和牛顿法则分别只需要34步和21步迭代. 以下是三者的(部分)收敛表格:

最速下降法	BFGS	牛顿法
1.827e-04	1.70e-03	3.48e-02
1.826e-04	1.17e-03	1.44e-02
1.824e-04	1.34e-04	1.82e-04
1.823e-04	1.01e-06	1.17e-08

下面再对BFGS和DFP算法的推导过程作进一步阐释.

1. BFGS公式是保正定的. 注意到推导出BFGS公式的极小化问题并没有显式地要求$B_k$是正定的. 不过我们可以证明: 如果$H_k$是正定的, 则更新得到的$H_{k+1}$也必定是正定的. 注意到${\rho }_k>0$, 因此BFGS公式是良定的. 对于$\forall z̸=0$, 我们有$$z^T{H}_{k+1}z=w^T{H}_{k}w+{\rho }_k(z^T{s}_k{)}^2\ge 0,$$其中我们定义$w=z-{\rho }_k{y}_k(s_k^{T}z)$. 由$H_k$的正定性, 等号成立仅当$z^T{s}_k=0$. 而这时$w=z̸=0$, 从而第一项为正. 总之, $H_{k+1}$是正定的. 进一步地, 我们说只要初始近似矩阵$B_0(H_0)$是正定阵, 则后面所有的近似矩阵都是正定的.
2. BFGS公式是尺度不变的. 构建BFGS公式时的极小化问题中还需要使得目标函数是尺度不变的. 为此我们构建了矩阵$W$并令目标在加权矩阵范数下极小. 不过在最后的公式中我们并没有发现任何$W$的身影. 这也是这一方法的独到之处. 我们也可以选取其他的$W$诱导加权范数 (但未必保证尺度不变性), 但一般不会显著地优于BFGS算法.
3. BFGS算法在应用到二次函数时具有许多有趣的性质. 我们将这一部分的讨论放在后面进行, 到时会涉及更加广泛的一类更新公式——Broyden族, 而BFGS公式仅仅是它的一种特殊情形.
4. BFGS公式具有自修正性. 若$H_k$错误地反映了目标函数的曲率信息且若这一近似减缓了迭代的进程, 则BFGS公式往往会在后面几步修正近似. 相比之下, DFP算法就没有这么好的自修正性, 这也导致了它在实际应用中表现不佳. 不过, BFGS的自修正性仅仅在有合适的线搜索为前提时才会满足, 例如Wolfe条件.
5. BFGS公式和DFP公式是对偶存在的——$s\leftrightarrow y,B\leftrightarrow H$. 这一对称性并不奇怪, 毕竟我们的算法构成和推导方式就是如此的类似.

1.3 实施细节

算法1还需要增添许多细节才能成为实用的算法.

1.3.1 初始步长的选取

我们提到线搜索需满足(强)Wolfe条件, 应当指明搜索应当以${\alpha }_k=1$为开端. 这是因为单位步长最后在一定条件下总会被接收, 从而达成超线性收敛. 系数$c_1=10^{-4},c_2=0.9$是较为常用的选择.

1.3.2 初始近似矩阵的选取

最简单的选取方式, 就是$H_0\leftarrow \beta I$, 其中$\beta$为一实数. 但这样一来$\beta$的选取就成了问题: 如果$\beta$过大, 则第一步$p_0=-\beta g_0$就会过长, 从而需要更多次获取函数值来确定合适的步长${\alpha }_0$. 一些软件会要求用户预先给定$\delta$作为第一步的长度, 然后设$H_0=\delta \Vert g_0{\Vert }^{-1}I$.
一种高效的启发式方法是, 在计算完第一步后但在第一次BFGS更新之前, 变换初始矩阵的尺度. 我们将临时的$H_0=I$改成$$H_0\leftarrow \frac{y_k^T{s}_k}{y_k^T{y}_k}I.$$注意这里的$s_k,y_k$是第一步尝试后得到的参量. 这使得$H_0$更加接近${\nabla }^{2}f(x_0{)}^{-1}$: 假设平均Hessian是正定矩阵, 则存在${\bar{G}}_k^{1/2}$使得${\bar{G}}_k={\bar{G}}_k^{1/2}{\bar{G}}_k^{1/2}$. 定义$z_k={\bar{G}}_k^{1/2}{s}_k$并由$y_k={\bar{G}}_k{s}_k$, 得到$$\frac{y_k^T{s}_k}{y_k^T{y}_k}=\frac{({\bar{G}}_k^{1/2}{s}_k{)}^T{\bar{G}}_k^{1/2}{s}_k}{({\bar{G}}_k^{1/2}{s}_k{)}^T{\bar{G}}_k{\bar{G}}_k^{1/2}{s}_k}=\frac{z_k^T{z}_k}{z_k^T{\bar{G}}_k{z}_k}.$$上式的倒数便是${\bar{G}}_k$的一个特征值的近似, 从而也是${\nabla }^{2}f(x_0)$的一个特征值的近似. 从而其本身是${\nabla }^{2}f(x_0{)}^{-1}$的一个特征值的近似.

1.3.3 BFGS算法之$B_k$版本的应用

之前我们利用Sherman-Morrison-Woodbury公式得到了BFGS公式的$B_k$更新版本, 并指出直接使用会徒增计算量. 因此, 我们不应当显式地储存$B_k$, 而应当储存这一矩阵的Cholesky因子$L_k{D}_k{L}_k^T$. 我们可以从$B_k$的更新公式推出$L_k,D_k$的更新公式, 而这一操作仅需$O(n^2)$的计算量. 由于Cholesky分解后$B_k{p}_k=-\nabla f_k$也可以在$O(n^2)$内解决, 因此总的运算量与算法1是相差无几的 (至少在量阶上). 这一方法的潜在优势为: 当$D_k$的对角元不够大时, 我们可以对其修正以保证算法的稳定性. 不过实际应用表明这样的优势并不明显, 我们仍然更加倾向于使用算法1. 下面我们来具体地阐述这一方法. 在这之前先证明之前遗留的一个问题: 满足隐含条件且满足割线方程的对称正定矩阵是存在的. 以下省略下标, 并记$B_{k+1}=B_{+}$.

注意矩阵$B_{+}$对称正定当且仅当存在可逆矩阵$J_{+}$使得$B_{+}=J_{+}{J}_{+}^T$. 因此此时割线方程等价于$J_{+}{J}_{+}^{T}s=y$. 而寻求$B_{+}$的存在也就等价于:

1. 选取$J_{+}$: 存在$v$, $J_{+}v=y$;
2. 选取$v$: $J_{+}^{T}s=v$.

设$B$有Cholesky分解$B=L{L}^T$. 事实上, 令$$J_{+}=L+\frac{(y-Lv)v^T}{v^{T}v}$$即可满足1. 进一步推算, $$v=L^{T}s+v\frac{(y-Lv{)}^{T}s}{v^{T}v}.$$因此, $v$有形式$v=\alpha L^{T}s$. 代入上式解得$${\alpha }^2=\frac{y^{T}s}{s^{T}L{L}^{T}s}=\frac{y^{T}s}{s^{T}Bs}.$$取$$v=\sqrt{\frac{y^{T}s}{s^{T}Bs}}{L}^{T}s$$即可. 现在我们要更新$B$的Cholesky因子. 我们以普通的Cholesky分解为例, 即$B=L{L}^T$. 注意到$$B_{+}=J_{+}{J}_{+}^T,J_{+}=L+\frac{(y-Lv)v^T}{v^{T}v}$$中, $J_{+}$不见得是下三角阵, 因此我们需要将分解化成$B_{+}=L_{+}{L}_{+}^T$. 这也是主要的计算量所在. 注意到$J_{+}$实际上是由$L$经秩1修正得到的. 我们可以写作$$J_{+}=L+\tilde{u}{\tilde{v}}^T.$$对秩1矩阵作Givens变换是极其方便的, 因为它任意两行皆成比例. 因此可以同时旋转两行. 我们分别对行号为$(1,2),(2,3),\dots ,(n-1,n)$的行做Givens变换 (一共$n-1$次), 最终将得到只有第一行非零的矩阵. 记这一过程的Givens旋转矩阵为$G$, 它是正交矩阵. 若在$L$上同时作用上$G$, 则$G(L+\tilde{u}{\tilde{v}}^T)$将是一个下Hessenberg矩阵. 于是再用$n-1$次Givens变换即可将其转变为下三角矩阵. 记这一阶段的Givens变换矩阵为$\tilde{G}$, 于是我们就有$$\tilde{G}G(L+\tilde{u}{\tilde{v}}^T)=L_{+}.$$易验证$B_{+}=L_{+}{L}_{+}^T$. 以上计算量可验证在$O(n^2)$.

1.3.4 线搜索条件的重要性

若线搜索不当 (例如没有满足Wolfe条件), 则BFGS算法的表现可能会大打折扣. 例如, 一些软件实施的是Armjo回溯线搜索: 首先尝试步长${\alpha }_k=1$, 再连续减少之直至满足充分下降条件. 对于这样的策略, 我们无法保证隐含条件$y_k^T{s}_k>0$是成立的. 为了克服这一缺陷, 有些人提出直接跳过更新直接令$H_{k+1}=H_k$. 不过这一方案并不提倡, 因为这样的跳步可能会过于频繁以至于$H_k$无法捕捉到重要的曲率信息. 在后面的章节中, 我们会讨论阻滞BFGS更新——它可以更高效地处理隐含条件不成立的情形.

1.3.5 关于内存

在使用$H$版本的BFGS公式时, 我们因不用计算逆矩阵而欣喜若狂, 但仔细一想便会发现$H_k$的使用也并不简单. 首先$H_k$的更新公式中有大量的向量, 因此与负梯度相乘相当于做内积. 这样的计算量在$O(n)$. 但要注意到中间的$H_k$. 如果我们不希望存储$H$, 则就必须回溯之前所有的$H_k$来计算这一次的乘积. 如果每一次都回溯之前所有的向量, 真正实施起来倒不如直接存储$H_k$了. 因此有人提出一种想法: 不妨仅用最近几次 (如$m$个) 的向量来更新, 这样再取前$m$步中的$H_k^0$为某一设定的好算的矩阵, 就可以将计算量和存储量严格控制在最低限度. 此时计算过程甚至不会涉及任何的矩阵-向量乘积. 这被称作是有限内存BFGS算法 (LM-BFGS). 详细内容我们将在后面的章节展示.

2. 对称秩1方法

2.1 介绍

在BFGS公式和DFP公式中, 更新矩阵$B_{k+1}$与原矩阵$B_k$都是差一个秩为2的矩阵, 我们也称之为秩2修正. 基于这样的观察, 我们自然会问: 是否存在更加简洁的秩1修正公式能够同时兼顾矩阵对称性与割线方程呢? 下面我们就来介绍一种对称秩1修正方法. 与之前的秩2修正不同的是, 这种方法不再保证更新后的矩阵仍然是正定的.
设有对称秩1更新公式为$$B_{k+1}=B_k+\sigma v{v}^T,$$其中 σ ∈ { − 1 , + 1 } \sigma\in\{-1,+1\} σ∈{−1,+1}, 且$\sigma ,v$的选取应当使得$B_{k+1}$满足割线方程$y_k=B_{k+1}{s}_k$. 代入可得$$y_k=B_k{s}_k+[\sigma v^T{s}_k]v.$$由于方括号内是个标量, 因此我们可以说$v$是$y_k-B_k{s}_k$的一个倍数, 即存在标量$\delta$使得$v=\delta (y_k-B_k{s}_k)$. 将这一表示代回, 就得到$$y_k-B_k{s}_k=\sigma {\delta }^2[s_k^T(y_k-B_k{s}_k)](y_k-B_k{s}_k).$$显然等式成立当且仅当$$\sigma =\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}[s_k^T(y_k-B_k{s}_k)],\delta =\pm |s_k^T(y_k-B_k{s}_k){|}^{-1/2}.$$因此唯一满足要求的对称秩1公式为$$(\mathrm{S}\mathrm{R}1)B_{k+1}=B_k+\frac{(y_k-B_k{s}_k)(y_k-B_k{s}_k{)}^T}{(y_k-B_k{s}_k{)}^T{s}_k}.$$同样地, 利用Sherman-Morrison-Woodbury公式就得到对应的矩阵逆更新公式:$$(\mathrm{S}\mathrm{R}1)H_{k+1}=H_k+\frac{(s_k-H_k{y}_k)(s_k-H_k{y}_k{)}^T}{(s_k-H_k{y}_k{)}^T{y}_k}.$$注意到SR1公式可以保持可逆性, 但不能保持正定性. 这一点在早期只使用线搜索进行优化的日子里, 一度被视为SR1公式的一项主要缺点. 然而, 随着信赖域方法的出现, SR1公式便重新登上舞台且呈现出强大的能力. 事实上, 它能产生不定Hessian近似的性质反而成为了其主要的优势.
不过SR1公式的缺陷也是明显的, 即公式中的分母可能会很小, 从而导致数值不稳定. 这一点即使放在凸二次函数下也可能发生. 这样一来所谓的SR1公式也就不复存在. 通过对$B_k$版本的探究 (对$H_k$版本也类似), 我们发现了三种可能的情形:

1. $(y_k-B_k{s}_k{)}^T{s}_k̸=0$: 这说明存在唯一的SR1公式满足割线方程, SR1公式良定;
2. $y_k=B_k{s}_k$: 此时唯一满足割线方程的SR1公式就是$B_{k+1}=B_k$;
3. $y_k̸=B_k{s}_k,(y_k-B_k{s}_k{)}^T{s}_k=0$: 此时不存在满足割线方程的SR1公式.

最后一种情形会引发数值不稳定. 这从侧面反映, SR1公式不再能满足所有的要求, 此时应当转向秩2公式.
SR1公式还有许多我们可以探究的方面和优点, 如:

1. 我们可以加入一些防护措施 (Safeguard)来避免上述的数值不稳定现象;
2. SR1公式往往能产生优于BFGS公式的Hessian近似;
3. 在带约束问题的拟牛顿法或求解部分可分离问题的方法中, 我们就不再需要强加隐含条件$y_k^T{s}_k>0$, 此时BFGS公式反倒不合理. 而SR1公式能够更好反映目标函数在当前位置的曲率信息.

下面介绍一种防护措施. 我们在实际应用中发现, 若分母很小, 我们直接跳过更新往往能带来较好的表现. 这就是说, 我们仅在$$|s_k^T(y_k-B_k{s}_k)|\ge r\Vert s_k\Vert \Vert y_k-B_k{s}_k\Vert ,r\in (0,1)$$时更新$B_k$. 这里的$r$往往取较小的数, 如$10^{-8}$. 否则, 我们直接令$B_{k+1}=B_k$.

这里出现一个问题: 为什么SR1方法可以跳步, BFGS方法却不建议跳步呢? 事实上这两种情形是截然不同的. $s_k^T(y_k-B_k{s}_k)\approx 0$较难出现, 这是因为这要求某些向量以特定的方式出现. 当这种情况发生, 跳过更新似乎不会对迭代产生什么负面影响. 这并不奇怪, 因为上面的约等式还可以写成$$s_k^T\bar{G}{s}_k\approx s_k^T{B}_k{s}_k,$$其中$\bar{G}$为上一步的平均Hessian. 这就说明当前的$B_k$沿着$s_k$的曲率已经正确反映真实情况了. 相反, BFGS公式要求的隐含条件$s_k^T{y}_k\ge 0$很容易被破坏 (例如线搜索不满足Wolfe条件), 从而跳步会极为频繁地出现, 从而影响了Hessian近似的质量.

下面给出使用信赖域框架的SR1方法的一个正式的表述. 使用信赖域作为框架是因为信赖域更能接受不定矩阵的使用.

算法2 (SR1 Trust-Region Method)
Given starting point $x_0$, initial Hessian approximation $B_0$, trust-region radius ${\Delta }_0$, convergence tolerance $ϵ>0$, parameters $\eta \in (0,10^{-3})$ and $r\in (0,1)$;
$k\leftarrow 0$;
while $\Vert \nabla f_k\Vert >ϵ$
\quad\quad Compute $s_k$ by solving the subproblem$$\underset{s}{\min }\nabla f_k^{T}s+\frac{1}{2}{s}^T{B}_{k}s\mathrm{s}.\mathrm{t}.\Vert s\Vert \le {\Delta }_k;$$
\quad\quad Compute$$\begin{array}{rl}{y}_k& =\nabla f(x_k+s_k)-\nabla f_k,\\ \mathrm{a}\mathrm{r}\mathrm{e}\mathrm{d}& =f_k-f(x_k+s_k)(\mathrm{a}\mathrm{c}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}\mathrm{r}\mathrm{e}\mathrm{d}\mathrm{u}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}),\\ \mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}& =-\left(\nabla f_k^T{s}_k+\frac{1}{2}{s}_k^T{B}_k{s}_k\right)(\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{c}\mathrm{t}\mathrm{e}\mathrm{d}\mathrm{r}\mathrm{e}\mathrm{d}\mathrm{u}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n});\end{array}$$
\quad\quad if ared/pred$>\eta$
$x_{k+1}=x_k+s_k;$
\quad\quad else
$x_{k+1}=x_k;$
\quad\quad end (if)
\quad\quad if ared/pred$>0.75$
\quad\quad\quad\quad if $\Vert s_k\Vert \le 0.8{\Delta }_k$
${\Delta }_{k+1}={\Delta }_k;$
\quad\quad\quad\quad else
${\Delta }_{k+1}=2{\Delta }_k;$
\quad\quad\quad\quad end (if)
\quad\quad else if $0.1\le$ared/pred$\le 0.75$
${\Delta }_{k+1}={\Delta }_k;$
\quad\quad else
${\Delta }_{k+1}=0.5{\Delta }_k;$
\quad\quad end (if)
\quad\quad if $|s_k^T(y_k-B_k{s}_k)|\ge r\Vert s_k\Vert \Vert y_k-B_k{s}_k\Vert$
\quad\quad\quad\quad Use SR1 to compute $B_{k+1}$(even if $x_{k+1}=x_k$);
\quad\quad else
$B_{k+1}\leftarrow B_k;$
\quad\quad end (if)
$k\leftarrow k+1;$
end (while)

以上是信赖域方法的典型形式. 具体说来, 我们使用了一种特定的更新信赖域半径的策略, 不过这也可以用其他的代替.
为获得较快的收敛速度, 我们有时需要沿着坏方向$s_k$更新$B_k$. 坏方向的产生说明$B_k$不能较好地近似当前方向的Hessian. 如果我们不更新近似, 就有可能产生连锁效应, 阻滞超线性收敛.

2.2 性质

SR1公式的一大优势在于, 能够产生对Hessian的良好近似. 我们首先对目标函数为二次函数的情形进行讨论. 对于这种函数, 步长的选取并不影响更新, 因此为检验更新的效果, 我们不妨假设步长一致为1, 即$$p_k=-H_k\nabla f_k,x_{k+1}=x_k+p_k.$$此时$p_k=s_k$.

定理1 设$f:{\mathbb{R}}^n\to \mathbb{R}$为强凸二次函数$f(x)=b^{T}x+\frac{1}{2}{x}^{T}Ax$, 其中$A$对称正定. 则对于任意初始点$x_0$以及任意初始矩阵$H_0$, 若对于所有的$k$都有$(s_k-H_k{y}_k{)}^T{y}_k̸=0$, 则由SR1产生的序列 { x k } \{x_k\} {xk​}至多经$n$步收敛到极小点. 此外, 若进行了$n$步且搜索方向$p_k$是线性无关的, 则$H_n=A^{-1}$.

证明: 由假设$(s_k-H_k{y}_k{)}^T{y}_k̸=0$, SR1更新总是良定的. 我们先用数学归纳法证明$$H_k{y}_j=s_j,j=0,1,\dots ,k-1.$$换句话说, 我们证明割线方程不仅仅沿着最近的搜索方向成立, 对之前所有的方向也都成立.
由定义, SR1公式满足割线方程, 所以$H_1{y}_0=s_0$. 假设上述对某个$k>1$成立, 下证$k+1$的情形. 由归纳假设,$$(s_k-H_k{y}_k{)}^T{y}_j=s_k^T{y}_j-y_k^T(H_k{y}_j)=s_k^T{y}_j-y_k^T{s}_j=0,\forall j<k,$$其中最后一个等式是因为在二次函数下恒有$y_i=A{s}_i$. 因此由SR1公式可得$$H_{k+1}{y}_j=H_k{y}_j=s_j,\forall j<k.$$而$H_{k+1}{y}_k=s_k$是由割线方程满足的, 因此$k+1$也成立. 因此这对任意$k$都成立.
若算法进行了$n$步且所有的 { s j } \{s_j\} {sj​}线性无关, 则$$s_j=H_n{y}_j=H_{n}A{s}_j,j=0,1,\dots ,n-1.$$这说明$H_{n}A=I\Rightarrow H_n=A^{-1}$. 因此, 在$x_n$处将实施牛顿步, 从而下一迭代点$x_{k+1}$就是解, 算法终止.
现在考虑搜索方向线性相关的情形. 设$s_k$为之前搜索方向的线性组合, 即存在标量${\xi }_i$使得$$s_k={\xi }_0{s}_0+\cdots +{\xi }_{k-1}{s}_{k-1}.$$因此$$\begin{array}{rl}{H}_k{y}_k& =H_{k}A{s}_k\\ & ={\xi }_0{H}_{k}A{s}_0+\cdots +{\xi }_{k-1}{H}_{k}A{s}_{k-1}\\ & ={\xi }_0{H}_k{y}_0+\cdots +{\xi }_{k-1}{H}_k{y}_{k-1}\\ & ={\xi }_0{s}_0+\cdots +{\xi }_{k-1}{s}_{k-1}\\ & =s_k.\end{array}$$这就得到$$H_k(\nabla f_{k+1}-\nabla f_k)=-H_k\nabla f_k,$$从中由$H_k$可逆就推出$\nabla f_{k+1}=0$. 因此$x_{k+1}$就是解. 证毕.

以上定理表明: 当目标函数$f$是二次函数时, 不论线搜索如何, 割线方程都沿之前所有的搜索方向成立. BFGS公式也有类似的结果, 只是需要加上线搜索精确的假设. 我们会最下一节证明之.
而对于一般的非线性函数, SR1公式将在一定条件下持续产生Hessian的良好近似. 我们不加证明的叙述之.

定理2 设$f$二阶连续可微, 其Hessian有界且在一点$x^{\ast }$的一个邻域内Lipschitz连续. 设 { x k } \{x_k\} {xk​}为收敛到$x^{\ast }$的任一序列. 设对于所有的$k$都有$$|s_k^T(y_k-B_k{s}_k)|\ge r\Vert s_k\Vert \Vert y_k-B_k{s}_k\Vert ,r\in (0,1),$$且$s_k$一致线性无关. 则由SR1公式产生的$B_k$满足$$\underset{k\to \infty }{\lim }\Vert B_k-{\nabla }^{2}f(x^{\ast })\Vert =0.$$这里"一致线性无关"大致是指, $s_k$不会都落在小于$n$维的子空间中. 这一假设一般 (但也不总是) 是成立的.

3. Broyden族

至今, 我们已讨论了BFGS、DFP、SR1这三种拟牛顿公式, 但事实上还存在许多其他的公式没有讨论. 其中最吸引人的当属Broyden族, 其更新公式均具有以下形式$$B_{k+1}=B_k-\frac{B_k{s}_k{s}_k^T{B}_k}{s_k^T{B}_k{s}_k}+\frac{y_k{y}_k^T}{y_k^T{s}_k}+{\varphi }_k(s_k^T{B}_k{s}_k)v_k{v}_k^T,$$其中${\varphi }_k$为标量参数, $$v_k=\frac{y_k}{y_k^T{s}_k}-\frac{B_k{s}_k}{s_k^T{B}_k{s}_k}.$$易知, BFGS和DFP公式都是Broyden族的一员——BFGS公式对应于${\varphi }_k=0$, DFP公式对应于${\varphi }_k=1$. 因此我们可以将Broyden族写成这两个公式的线性组合:$$B_{k+1}=(1-{\varphi }_k)B_{k+1}^{\mathrm{B}\mathrm{F}\mathrm{G}\mathrm{S}}+{\varphi }_k{B}_{k+1}^{\mathrm{D}\mathrm{F}\mathrm{P}}.$$这一关系式表明, Broyden族中的所有成员均满足割线方程, 且由于BFGS和DFP在满足隐含条件时保持正定性, 因此在Broyden族中$0\le {\varphi }_k\le 1$也同样保持正定性.
至今已有许多关于所谓的限制Broyden族的研究, 即限制Broyden族至${\varphi }_k\in [0,1]$. 当应用于二次函数时, 它也同SR1公式一样具有一些特殊性质. 由于下面的分析独立于步长选取, 我们假设$$p_k=-B_k^{-1}\nabla f_k,x_{k+1}=x_k+p_k.$$

定理3 设$f:{\mathbb{R}}^n\to \mathbb{R}$为强凸二次函数$f(x)=b^{T}x+\frac{1}{2}{x}^{T}Ax$, 其中$A$对称正定. 设$x_0$为任意初始点, $B_0$为任意初始矩阵, 后续$B_k$由限制Broyden族产生. 令${\lambda }_1^k\le {\lambda }_2^k\le \cdots \le {\lambda }_n^k$为矩阵$$A^{\frac{1}{2}}{B}_k^{-1}{A}^{\frac{1}{2}}$$的特征值, 于是对于所有的$k$, 我们有 min ⁡ { λ i k , 1 } ≤ λ i k + 1 ≤ max ⁡ { λ i k , 1 } , i = 1 , 2 , … , n . \min\{\lambda_i^k,1\}\le\lambda_i^{k+1}\le\max\{\lambda_i^k,1\},\quad i=1,2,\ldots,n. min{λik​,1}≤λik+1​≤max{λik​,1},i=1,2,…,n.需指出, 上述不等式对于${\varphi }_k̸\in [0,1]$不成立.

下面我们不加证明地讨论这个结论的意义. 若$A^{\frac{1}{2}}{B}_k^{-1}{A}^{\frac{1}{2}}$的特征值都是1, 则$B_k$就是二次目标函数的Hessian$A$. 这是理想状态, 因此我们希望这些特征值尽可能地靠近1. 事实上, 上面的不等式就告诉我们特征值序列 { λ i k } \{\lambda_i^k\} {λik​}单调向1收敛 (不一定收敛到1). 更重要的是, 定理3的结论对于非精确线搜索同样成立.

尽管定理3可能会让人觉得最好的更新公式属于限制Broyden族中, 但实际情况并非如此. 已有一些分析和计算结果表明允许${\varphi }_k$为负 (当然要在可控范围内) 的算法往往要比BFGS算法更加优越. SR1公式就是个例子: 它也是Broyden族中的一员, 其中$${\varphi }_k=\frac{s_k^T{y}_k}{s_k^T{y}_k-s_k^T{B}_k{s}_k},$$但它却不属于限制Broyden族, 这是因为这样的${\varphi }_k$可能落在$[0,1]$之外.
在本节接下来部分, 我们来更确切地讨论保持正定性的${\varphi }_k$的取值区间. 这要用到下面的定理:

定理4 (Interlacing Eigenvalue Theorem) 设$A\in {\mathbb{R}}^{n\times n}$为对称阵, 具有特征值${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_n$. 设$z\in {\mathbb{R}}^n:\Vert z\Vert =1$, $\alpha \in \mathbb{R}$. 记$A+\alpha z{z}^T$的特征值为${\xi }_1\ge {\xi }_2\ge \cdots \ge {\xi }_n$, 则对于$\alpha$有$${\xi }_1\ge {\lambda }_1\ge {\xi }_2\ge {\lambda }_2\ge {\xi }_3\ge \cdots \ge {\xi }_n\ge {\lambda }_n,$$且$$\sum _{i=1}^n{\xi }_i-{\lambda }_i=\alpha .$$若$\alpha <0$, 则对应$${\lambda }_1\ge {\xi }_1\ge {\lambda }_2\ge {\xi }_2\ge {\lambda }_3\ge \cdots \ge {\lambda }_n\ge {\xi }_n,$$且$$\sum _{i=1}^n{\xi }_i-{\lambda }_i=\alpha .$$
换句话说, 经秩1修正后的矩阵的特征值与原矩阵的特征值呈交错分布, 且当$\alpha >0$时为非负调整, 当$\alpha <0$时为非正调整. 调整的总量为$\alpha$, 这也正好是秩1修正的2-范数$\Vert \alpha z{z}^T{\Vert }_2$.

证明: 不妨假设$A,A+\alpha z{z}^T$均为加边对角阵, 即$$A=\left(\begin{array}{cc}\gamma & b^T\\ b& \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(a_i)\end{array}\right),A+\alpha z{z}^T=\left(\begin{array}{cc}\gamma +\alpha & b^T\\ b& \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(a_i)\end{array}\right).$$事实上, 对于秩1矩阵$\alpha z{z}^T$, 存在正交矩阵$R$, 使得$$R^T(\alpha z{z}^T)R=\left(\begin{array}{cc}\alpha & O\\ O& O\end{array}\right),$$若记$$R^{T}AR=\left(\begin{array}{cc}\gamma & a^T\\ a& A_{n-1}\end{array}\right),$$则存在$n-1$阶正交阵$S$, 使得$$S^T{A}_{n-1}S=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(a_i).$$若我们定义$Q$为$$Q=R\left(\begin{array}{cc}1& O\\ O& S\end{array}\right),$$则$Q$是正交矩阵, 且$$Q^T(A+\alpha z{z}^T)Q=\left(\begin{array}{cc}\gamma & b^T\\ b& \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(a_i)\end{array}\right)+\left(\begin{array}{cc}\alpha & O\\ O& O\end{array}\right),$$其中$b=S^{T}a$. 因此$A$与$A+\alpha z{z}^T$的特征值就等于$$\left(\begin{array}{cc}\gamma & b^T\\ b& \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(a_i)\end{array}\right)与\left(\begin{array}{cc}\gamma +\alpha & b^T\\ b& \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(a_i)\end{array}\right)$$的特征值.
设$b$只有$s$个分量不为0. 若$b_j$为0, 则$a_j$便是$A$的特征值. 适当选取只与后$n-1$行有关的置换阵$P$, 我们可得矩阵$$B=P^{T}AP=\left(\begin{array}{ccc}\gamma & b^T& O\\ b& \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\beta }_i)& O\\ O& O& \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\gamma }_i)\end{array}\right).$$这里仍用$b$表示边上的向量, 只是现在$b$的分量都不是0. $\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\beta }_i)$是$s$阶的, $\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\gamma }_i)$是$n-1-s$阶的, ${\beta }_i,{\gamma }_i$合在一起是$a_i$的一个排列. 注意某些${\gamma }_i$可以是$\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(a_i)$的重特征值. 因此, $A$的特征值为${\gamma }_i$以及矩阵$$C=\left(\begin{array}{cc}\gamma & b^T\\ b& \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\beta }_i)\end{array}\right)$$的特征值. 若$s=0$, 则$C$就是元素$\gamma$, 此时$A$的特征值为$\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(a_i)$的特征值以及$\gamma$. 若不然, 我们考察$C$的特征多项式, 即$$(\gamma -\lambda )\prod _{i=1}^s({\beta }_i-\lambda )-\sum _{j=1}^s{b}_j^2\prod _{i̸=j}({\beta }_i-\lambda )=0.$$假定$\beta )i$中只有$t$个不同的值, 不失一般性, 可令它们为${\beta }_1,{\beta }_2,\dots ,{\beta }_t$, 重数依次为$r_1,r_2,\dots ,r_t$, 于是$$r_1+r_2+\cdots +r_t=s.$$显然特征多项式有因子$$\prod _{i=1}^t({\beta }_i-\lambda {)}^{r_i-1},$$因此${\beta }_i$为$C$的特征值, 重数为$r_i-1$. 将特征多项式除以${\prod }_{i=1}^t){\beta }_i-\lambda {)}^{r_i}$, 我们得到$C$的其余特征值是方程$$0=(\gamma -\lambda )-\sum _{i=1}^t{c}_i^2({\beta }_i-\lambda {)}^{-1}\equiv \gamma -f(\lambda )$$的根, 其中$c_i^2$为与${\beta }_i$有关的$b_j^2$之和, 因而它严格大于0.

上图给出了$f(\lambda )$的图像, 其中不同的${\beta }_i$按降序排列. 显然$\gamma =f(\lambda )$的$t+1$个根——记为${\delta }_1,{\delta }_2,\dots ,{\delta }_{t+1}$——满足关系式$$\infty >{\delta }_1>{\beta }_1;{\beta }_{i-1}>{\delta }_i>{\beta }_i(i=2,3,\dots ,t);{\beta }_t>{\delta }_{t+1}>-\infty .$$因此$A$的$n$个特征值分成三组:

1. 相应于$b_i=0$的特征值${\gamma }_1,{\gamma }_2,\dots ,{\gamma }_{n-1-s}$, 它们等于$a_i$中的$n-1-s$个数;
2. 由$(r_i-1)$个等于${\beta }_i(i=1,2,\dots ,t)$的值所组成的$s-t$歌特征值, 它们等于$a_i$中另外的$s-t$个数;
3. $t+1$个特征值等于${\delta }_i$. 若$t=0$, 则${\delta }_1=\gamma$.

由此一来, 如果$\alpha >0$, 则第三组的特征值将不减; 如果$\alpha <0$, 则第三组的特征值将不增. 且易得交错分布. 证毕.

因此在Broyden族中, 只要${\varphi }_k$为正, 特征值就不会减少, 从而$B_{k+1}$在${\varphi }_k\ge 0$时都是正定的. 反过来, 若${\varphi }_k<0$, 则特征值将不增. 随着${\varphi }_k$的减小, 更新后的矩阵将变得奇异接着不定. 例如当$${\varphi }_k^c=\frac{1}{1-{\mu }_k},{\mu }_k=\frac{(y_k^T{B}_k^{-1}{y}_k)(s_k^T{B}_k{s}_k)}{(y_k^T{s}_k{)}^2}$$时, $B_{k+1}$就是奇异的 (未验证) . 由Cauchy-Schwarz不等式可知${\mu }_k\ge 1$, 从而${\varphi }_k^c\le 0$. 因此对于任意对称正定的初始矩阵$B_0$, 且若对每个$k$都有$s_k^T{y}_k>0$, ${\varphi }_k>{\varphi }_k^c$, 则由Broyden公式产生的$B_k$都保持对称正定性.
当线搜索精确时, 所有${\varphi }_k\ge {\varphi }_k^c$的Broyden公式将产生相同的迭代序列 (未验证) . 这个结论对一般的非线性函数都是正确的, 此时所有Broyden族中的公式产生的搜索方向只有长度的区别.
类似于SR1方法, Broyden族在应用于二次函数时也有许多良好的性质 (以精确线搜索为前提). 我们不加证明的叙述于下面的定理中.

定理5 设Broyden族中的某一公式应用于强凸二次函数$f(x)=b^{T}x+\frac{1}{2}{x}^{T}Ax$, 其中$x_0$为初始点, $B_0$为任意对称正定的初始矩阵. 假设${\alpha }_k$为精确步长, ${\varphi }_k\ge {\varphi }_k^c,\forall k$. 则下面的叙述正确.

1. 迭代点独立于${\varphi }_k$的选取, 且经至多$n$步收敛于解;
2. 割线方程对所有先前的搜索方向都成立, 即$$B_k{s}_j=y_j,j=k-1,k-2,\dots ,1;$$
3. 若初始矩阵$B_0=I$, 则产生的迭代点与共轭梯度法产生的迭代点相同. 特别地, 搜索方向相互共轭, 即$$s_i^{T}A{s}_j=0,i̸=j;$$
4. 若经$n$步迭代, 则$B_n=A$.

注意1,2,4与定理1是类似的.
我们可以稍微推广定理5: 它在Hessian近似保持非奇异但不必正定时仍然成立. (因此, 只要选取的${\varphi }_k$不致产生奇异矩阵, 我们允许${\varphi }_k<{\varphi }_k^c$.) 我们也可以推广第3条为, 若初始矩阵$B_0$不是单位阵, 则Broyden族等价于预处理共轭梯度法, 其中$B_0$作为预处理子.
类似于定理5这样的结论往往招徕研究人员的注意. 这是因为实际应用时的非精确线搜索可能会导致算法表现出现极大不同. 尽管如此, 这样的分析仍然促进了拟牛顿法的发展.

4. 收敛分析

本节我们讨论BFGS算法与SR1方法的全局与局部收敛性. 由于关于BFGS算法的分析结果更具有一般性和启发性, 我们将着重介绍之, 而只简略介绍SR1方法. 与之前的最速下降法与牛顿法不同的是, 拟牛顿法是以更新公式计算Hessian的近似, 这也增加了分析的复杂性.
提前说明, 即使BFGS算法与SR1方法在实际应用中被证实是极为强健的, 但我们不能对一般的非线性目标函数建立起真正的全局收敛性: 即, 我们无法对拟牛顿法证明, 从任意初始点和任意初始矩阵出发, 算法都能收敛到稳定点. 事实上, 这一性质至今还未得到证明. 在我们的分析中, 我们会假设

• 目标函数为凸; 或
• 迭代点满足一定的性质.

不过就局部而言, 我们在一定的条件下是可以证明它们的超线性收敛速度的.
以下以$\Vert \cdot \Vert$表示欧式向量或矩阵范数, 以$G(x)$表示Hessian矩阵${\nabla }^{2}f(x)$.

4.1 BFGS算法的全局收敛性

我们假定目标函数具有如下性质:

1. $f$二阶连续可微;
2. 水平集 L = { x ∈ R n ∣ f ( x ) ≤ f ( x 0 ) } \mathcal{L}=\{x\in\mathbb{R}^n|f(x)\le f(x_0)\} L={x∈Rn∣f(x)≤f(x0​)}是凸集, 且存在正常数$m,M$使得$$m\Vert z{\Vert }^2\le z^{T}G(x)z\le M\Vert z{\Vert }^2,\forall z\in {\mathbb{R}}^n,x\in \mathcal{L}.$$

其中第2条说明$G(x)$在$\mathcal{L}$上是正定的, 从而$f$在$\mathcal{L}$上有唯一极小点$x^{\ast }$.
我们利用假设考察两个常数$$\frac{y_k^T{s}_k}{s_k^T{s}_k},\frac{y_k^T{y}_k}{y_k^T{s}_k}$$的界. 由平均Hessian矩阵的性质, $$\frac{y_k^T{s}_k}{s_k^T{s}_k}=\frac{s_k^T{\bar{G}}_k{s}_k}{s_k^T{s}_k}\ge m.$$由于${\bar{G}}_k$正定, 因此可定义它的平方根矩阵. 定义$z_k={\bar{G}}_k^{1/2}{s}_k$, 从而$$\frac{y_k^T{y}_k}{y_k^T{s}_k}=\frac{s_k^T{\bar{G}}_k^2{s}_k}{s_k^T{\bar{G}}_k{s}_k}=\frac{z_k^T{\bar{G}}_k{z}_k}{z_k^T{z}_k}\le M.$$下面我们证明BFGS算法的全局收敛性. 我们将不再像在第三章第二节中那样建立关于Hessian近似$B_k$条件数的界, 而是引入矩阵的迹与行列式来估计Hessian近似矩阵的最大和最小特征值.

定理6 设$B_0$为任意对称正定初始矩阵, $x_0$为使假设成立的初始点. 则由算法1 ($ϵ=0$) 产生的序列 { x k } \{x_k\} {xk​}将收敛于$f$的极小点$x^{\ast }$.
证明: 定义$$m_k=\frac{y_k^T{s}_k}{s_k^T{s}_k},M_k=\frac{y_k^T{y}_k}{y_k^T{s}_k},$$由之前的讨论可知$$m_k\ge m,M_k\le M.$$计算BFGS更新公式两端的迹可得$$\mathrm{t}\mathrm{r}(B_{k+1})=\mathrm{t}\mathrm{r}(B_k)-\frac{\Vert B_k{s}_k{\Vert }^2}{s_k^T{B}_k{s}_k}+\frac{\Vert y_k{\Vert }^2}{y_k^T{s}_k}.$$下面计算BFGS更新公式两段的行列式. 这里要用到一个秩2修正的行列式引理.

引理 $$\det (I+x{y}^T+u{v}^T)=(1+y^{T}x)(1+v^{T}u)-(y^{T}u)(v^{T}x).$$
引理的证明: 不妨设$x̸\parallel u$, 则以$x,u$为列向量的矩阵可以扩充成可逆矩阵, 即存在$$Q=(x,u,w_1,\dots ,w_{n-2})$$可逆. 因此$$\begin{array}{rl}\det (I+x{y}^T+u{v}^T)& =\det \left(Q^{-1}\left(I+(x,u)\left(\begin{array}{c}{y}^T\\ v^T\end{array}\right)\right)Q\right)\\ & =\det \left(I+(e_1,e_2)\left(\begin{array}{ccc}{y}^{T}x& y^{T}u& \cdots \\ v^{T}x& v^{T}u& \cdots \end{array}\right)\right)\\ & =\det \left(I+\left(\begin{array}{ccc}{y}^{T}x& y^{T}u& \cdots \\ v^{T}x& v^{T}u& \cdots \end{array}\right)(e_1,e_2)\right)\\ & =\det \left(I+\left(\begin{array}{cc}{y}^{T}x& y^{T}u\\ v^{T}x& v^{T}u\end{array}\right)\right)\\ & =(1+y^{T}x)(1+v^{T}u)-(y^{T}u)(v^{T}x).\end{array}$$回到定理6的证明. 因此$$\begin{array}{rl}\det (B_{k+1})& =\det (B_k)\det \left(I-\frac{s_k{s}_k^T{B}_k}{s_k^T{B}_k{s}_k}+\frac{B_k^{-1}{y}_k{y}_k^T}{y_k^T{s}_k}\right)\\ & =\det (B_k)\frac{y_k^T{s}_k}{s_k^T{B}_k{s}_k}.\end{array}$$我们证明全局收敛性的主要工具仍然是Zoutendijk定理. 现定义$$\cos {\theta }_k=\frac{s_k^T{B}_k{s}_k}{\Vert s_k\Vert \Vert B_k{s}_k\Vert },q_k=\frac{s_k^T{B}_k{s}_k}{s_k^T{s}_k}.$$于是$$\frac{\Vert B_k{s}_k{\Vert }^2}{s_k^T{B}_k{s}_k}=\frac{\Vert B_k{s}_k{\Vert }^2\Vert s_k{\Vert }^2}{(s_k^T{B}_k{s}_k{)}^2}\frac{s_k^T{B}_k{s}_k}{\Vert s_k{\Vert }^2}=\frac{q_k}{{\cos }^2{\theta }_k}.$$此外, 还有$$\det (B_{k+1})=\det (B_k)\frac{y_k^T{s}_k}{s_k^T{s}_k}\frac{s_k^T{s}_k}{s_k^T{B}_k{s}_k}=\det (B_k)\frac{m_k}{q_k}.$$下面结合迹和行列式定义正定矩阵$B$的函数$$\psi (B)=\mathrm{t}\mathrm{r}(B)-\ln (\det (B)).$$易知$\psi (B)>0$. 因此$$\begin{array}{rl}\psi (B_{k+1})& =\mathrm{t}\mathrm{r}(B_k)+M_k-\frac{q_k}{{\cos }^2{\theta }_k}-\ln (\det (B_k))-\ln m_k+\ln q_k\\ & =\psi (B_k)+(M_k-\ln m_k-1)+\left[1-\frac{1_k}{{\cos }^2{\theta }_k}+\ln \frac{q_k}{{\cos }^2{\theta }_k}\right]+\ln {\cos }^2{\theta }_k.\end{array}$$由于函数$h(t)=1-t+\ln t$对所有$t>0$均非正, 从而递推可得$$0<\psi (B_{k+1})\le \psi (B_0)+c(k+1)+\sum _{j=0}^k\ln {\cos }^2{\theta }_j,$$其中$c=M-\ln m-1$(不妨设)为正. 注意到拟牛顿法迭代公式$$x_{k+1}=x_k+{\alpha }_k{p}_k,p_k=-B_k^{-1}\nabla f_k\Rightarrow s_k=-{\alpha }_k{B}_k^{-1}\nabla f_k,$$因此${\theta }_k$也是最速下降方向和搜索方向的夹角. 联系Zoutendijk定理, 序列 { ∥ ∇ f k ∥ } \{\Vert \nabla f_k\Vert\} {∥∇fk​∥}在0的邻域外仅当$\cos {\theta }_j\to 0$.
下面用反证法证明, 假设$\cos {\theta }_j\to 0$. 从而存在$k_1>0$使得对于所有的$j>k_1$, 有$$\ln {\cos }^2{\theta }_j<-2c.$$由此对于$k>k_1$均有:$$\begin{array}{rl}0& <\psi (B_0)+c(k+1)+\sum _{j=0}^{k_1}\ln {\cos }^2{\theta }_j+\sum _{j=k_1+1}^k(-2c)\\ & =\psi (B_0)+\sum _{j=0}^{k_1}\ln {\cos }^2{\theta }_j+2c{k}_1+c-ck.\end{array}$$然而右端对于充分大的$k$显然是负的, 从而矛盾. 因此存在指标 { j k } k = 1 , 2 , … \{j_k\}_{k=1,2,\ldots} {jk​}k=1,2,…​使得$\cos {\theta }_{j_k}\ge \delta >0$. 由Zoutendijk定理这说明$\mathrm{liminf}\Vert \nabla f_k\Vert =0$. 由于问题强凸, 这足以说明$x_k\to x^{\ast }$. 证毕.

定理6已被推广至整个限制Broyden族上, 除了DFP方法 (即${\varphi }_k=0$). 这是因为随着${\varphi }_k$向1逼近, 一些自修正的性质将会极大地弱化.
上述分析的一个推广可以推出迭代序列是线性收敛的. 特别地, 我们可证明序列 { ∥ x k − x ∗ ∥ } \{\Vert x_k-x^*\Vert\} {∥xk​−x∗∥}满足$$\sum _{k=1}^{\infty }\Vert x_k-x^{\ast }\Vert <\infty .$$我们不证明之, 而假定它成立, 进而证明超线性收敛性.

4.2 BFGS算法的超线性收敛性

我们利用第三章中拟牛顿法超线性收敛的充分条件$$\underset{k\to \infty }{\lim }\frac{\Vert (B_k-{\nabla }^{2}f(x^{\ast })p_k\Vert }{\Vert p_k\Vert }=0$$来证明超线性收敛性. 注意这一结论对于一般的非线性目标函数仍然成立. 为此我们需要添加额外的假设: Hessian矩阵$G$在$x^{\ast }$处Lipshcitz连续, 即$$\Vert G(x)-G(x^{\ast })\Vert \le L\Vert x-x^{\ast }\Vert ,$$其中$x$在$x^{\ast }$附近, $L$为正常数.
我们引入如下参量:$${\tilde{s}}_k=G_{\ast }^{1/2}{s}_k,{\tilde{y}}_k=G_{\ast }^{-1/2}{y}_k,{\tilde{B}}_k=G_{\ast }^{-1/2}{B}_k{G}_{\ast }^{-1/2},$$其中$G_{\ast }=G(x^{\ast })$, $x^{\ast }$为$f$的一个极小点. 类似于之前定理6的证明, 我们定义$$\cos {\tilde{\theta }}_k=\frac{{\tilde{s}}_k^T{\tilde{B}}_k{\tilde{s}}_k}{\Vert {\tilde{s}}_k\Vert \Vert {\tilde{B}}_k{\tilde{s}}_k\Vert },{\tilde{q}}_k=\frac{{\tilde{s}}_k^T{\tilde{B}}_k{\tilde{s}}_k}{\Vert {\tilde{s}}_k{\Vert }^2},$$$${\tilde{M}}_k=\frac{\Vert {\tilde{y}}_k{\Vert }^2}{{\tilde{y}}_k^T{\tilde{s}}_k},{\tilde{m}}_k=\frac{{\tilde{y}}_k^T{\tilde{s}}_k}{{\tilde{s}}_k^T{\tilde{s}}_k}.$$在BFGS公式左右端左乘于右乘$G_{\ast }^{-1/2}$并合并整理可得$${\tilde{B}}_{k+1}={\tilde{B}}_k-\frac{{\tilde{B}}_k{\tilde{s}}_k{\tilde{s}}_k^T{\tilde{B}}_k}{{\tilde{s}}_k^T{\tilde{B}}_k{\tilde{s}}_k}+\frac{{\tilde{y}}_k{\tilde{y}}_k^T}{{\tilde{y}}_k^T{\tilde{s}}_k}.$$由于这一表达式与BFGS公式具有同样的形式, 因此同样地有$$\begin{array}{rl}\psi ({\tilde{B}}_{k+1})& =\psi ({\tilde{B}}_k)+({\tilde{M}}_k-\ln {\tilde{m}}_k-1)\\ & +\left[1-\frac{{\tilde{q}}_k}{{\cos }^2{\tilde{\theta }}_k}+\ln \frac{{\tilde{q}}_k}{{\cos }^2{\tilde{\theta }}_k}\right]+\ln {\cos }^2{\tilde{\theta }}_k.\end{array}$$由于$$y_k={\bar{G}}_k{s}_k,$$所以$$y_k-G_{\ast }{s}_k=({\bar{G}}_k-G_{\ast })s_k,$$从而$${\tilde{y}}_k-{\tilde{s}}_k=G_{\ast }^{-1/2}({\bar{G}}_k-G_{\ast })G_{\ast }^{-1/2}{\tilde{s}}_k.$$由假设条件, $$\Vert {\tilde{y}}_k-{\tilde{s}}_k\Vert \le \Vert G_{\ast }^{-1/2}{\Vert }^2\Vert {\tilde{s}}_k\Vert \Vert {\bar{G}}_k-G_{\ast }\Vert \le \Vert G_{\ast }^{-1/2}{\Vert }^2\Vert {\tilde{s}}_k\Vert L{ϵ}_k,$$其中 ϵ k = max ⁡ { ∥ x k + 1 − x ∗ ∥ , ∥ x k − x ∗ ∥ } . \epsilon_k=\max\{\Vert x_{k+1}-x^*\Vert,\Vert x_k-x^*\Vert\}. ϵk​=max{∥xk+1​−x∗∥,∥xk​−x∗∥}.因此$$\frac{\Vert {\tilde{y}}_k-{\tilde{s}}_k\Vert }{\Vert {\tilde{s}}_k\Vert }\le \bar{c}{ϵ}_k,$$其中$\bar{c}$为一正常数. 下面我们将利用此式证明拟牛顿法的超线性收敛.

定理7 设$f$二阶连续可微且BFGS算法产生的迭代点收敛到满足假设条件的极小点$x^{\ast }$. 且设$$\sum _{k=1}^{\infty }\Vert x_k-x^{\ast }\Vert <\infty .$$于是$x_k$超线性收敛于$x^{\ast }$.
证明: 由三角不等式, $$\Vert {\tilde{y}}_k\Vert -\Vert {\tilde{s}}_k\Vert \le \bar{c}{ϵ}_k\Vert {\tilde{s}}_k\Vert ,\Vert {\tilde{s}}_k\Vert -\Vert {\tilde{y}}_k\Vert \le \bar{c}{ϵ}_k\Vert {\tilde{s}}_k\Vert ,$$从而$$(1-\bar{c}{ϵ}_k)\Vert {\tilde{s}}_k\Vert \le \Vert {\tilde{y}}_k\Vert \le (1+\bar{c}{ϵ}_k)\Vert {\tilde{s}}_k\Vert .$$平方后可得$$(1-\bar{c}{ϵ}_k{)}^2\Vert {\tilde{s}}_k{\Vert }^2-2{\tilde{y}}_k^T{\tilde{s}}_k+\Vert {\tilde{s}}_k{\Vert }^2\le \Vert {\tilde{y}}_k{\Vert }^2-2{\tilde{y}}_k^T{\tilde{s}}_k+\Vert {\tilde{s}}_k{\Vert }^2\le {\bar{c}}^2{ϵ}_k^2\Vert {\tilde{s}}_k{\Vert }^2,$$因此$$2{\tilde{y}}_k^T{\tilde{s}}_k\ge (1-2\bar{c}{ϵ}_k+{\bar{c}}^2{ϵ}_k^2+1-{\bar{c}}^2{ϵ}_k^2)\Vert {\tilde{s}}_k{\Vert }^2=2(1-\bar{c}{ϵ}_k)\Vert {\tilde{s}}_k{\Vert }^2.$$由${\tilde{m}}_k$的定义可得$${\tilde{m}}_k=\frac{{\tilde{y}}_k^T{\tilde{s}}_k}{\Vert {\tilde{s}}_k{\Vert }^2}\ge 1-\bar{c}{ϵ}_k.$$类似地我们可以得到$${\tilde{M}}_k=\frac{\Vert {\tilde{y}}_k{\Vert }^2}{{\tilde{y}}_k^T{\tilde{s}}_k}\le \frac{1+\bar{c}{ϵ}_k}{1-\bar{c}{ϵ}_k}.$$因$x_k\to x^{\ast }$, 所以${ϵ}_k\to 0$, 从而存在$c>\bar{c}$使得下面不等式对于充分大的$k$成立:$${\tilde{M}}_k\le 1+\frac{2\bar{c}}{1-\bar{c}{ϵ}_k}{ϵ}_k\le 1+c{ϵ}_k.$$再次利用$h(t)=1-t+\ln t$的非正性, $$\frac{-x}{1-x}-\ln (1-x)=h\left(\frac{1}{1-x}\right)\le 0.$$于是对于充分大的$k$我们可以假设$\bar{c}{ϵ}_k<\frac{1}{2}$, 因此$$\ln (1-\bar{c}{ϵ}_k)\ge \frac{-\bar{c}{ϵ}_k}{1-\bar{c}{ϵ}_k}\ge -2\bar{c}{ϵ}_k.$$所以对于充分大的$k$, $$\ln {\tilde{m}}_k\ge \ln (1-\bar{c}{ϵ}_k)\ge -2\bar{c}{ϵ}_k>-2c{ϵ}_k.$$因此$$0<\psi ({\tilde{B}}_{k+1})\le \psi ({\tilde{B}}_k)+3c{ϵ}_k+\ln {\cos }^2{\tilde{\theta }}_k+\left[1-\frac{{\tilde{q}}_k}{{\cos }^2{\tilde{\theta }}_k}+\ln \frac{{\tilde{q}}_k}{{\cos }^2{\tilde{\theta }}_k}\right].$$累和以上, 我们有$$\sum _{j=0}^{\infty }\left(\ln \frac{1}{{\cos }^2{\tilde{\theta }}_j}-\left[1-\frac{{\tilde{q}}_k}{{\cos }^2{\tilde{\theta }}_k}+\ln \frac{{\tilde{q}}_k}{{\cos }^2{\tilde{\theta }}_k}\right]\right)\le \psi ({\tilde{B}}_0)+3c\sum _{j=0}^{\infty }{ϵ}_j<+\infty .$$由于中括号内非正, 而$\ln \left(1/{\cos }^2{\tilde{\theta }}_j\right)\ge 0$, 因此$$\underset{j\to \infty }{\lim }\ln \frac{1}{{\cos }^2{\tilde{\theta }}_j}=0,\underset{j\to \infty }{\lim }\left[1-\frac{{\tilde{q}}_k}{{\cos }^2{\tilde{\theta }}_k}+\ln \frac{{\tilde{q}}_k}{{\cos }^2{\tilde{\theta }}_k}\right]=0,$$$$\Rightarrow \underset{j\to \infty }{\lim }\cos {\tilde{\theta }}_j=1,\underset{j\to \infty }{\lim }{\tilde{q}}_j=1.$$我们现在只需套用超线性收敛的充分条件. $$\begin{array}{rl}\frac{\Vert G_{\ast }^{-1/2}(B_k-G_{\ast })s_k{\Vert }^2}{\Vert G_{\ast }^{1/2}{s}_k{\Vert }^2}& =\frac{\Vert ({\tilde{B}}_k-I){\tilde{s}}_k{\Vert }^2}{\Vert {\tilde{s}}_k{\Vert }^2}\\ & =\frac{\Vert {\tilde{B}}_k{\tilde{s}}_k{\Vert }^2-2{\tilde{s}}_k^T{\tilde{B}}_k{\tilde{s}}_k+{\tilde{s}}_k^T{\tilde{s}}_k}{{\tilde{s}}_k^T{\tilde{s}}_k}\\ & =\frac{{\tilde{q}}_k^2}{{\cos }^2{\tilde{\theta }}_k}-2{\tilde{q}}_k+1.\end{array}$$由于右端收敛到0, 因此$$\underset{k\to \infty }{\lim }\frac{\Vert (B_k-G_{\ast })s_k\Vert }{\Vert s_k\Vert }=0.$$这说明单位步长${\alpha }_k=1$将最终在解附近满足Wolfe条件, 从而达成超线性收敛. 证毕

4.3 SR1方法的收敛性分析

SR1方法的收敛性质并不如BFGS算法那么好理解. 事实上, 除了之前叙述的关于二次函数的结论, 我们至今没有获得像定理6 (全局收敛性) 和定理7 (局部超线性收敛) 的结论. 然而关于信赖域SR1方法却有一个很有意思的结果: 若目标函数有唯一稳定点, SR1从不跳步且$B_k$上有界, 则$x_k$以$(n+1)-$步超线性收敛速度收敛于$x^{\ast }$. 这一结论不需要信赖域子问题的精确求解. 我们不加证明地叙述之:

定理8 设算法2产生迭代点$x_k$. 假设下面的条件满足:

1. 序列不终止, 但始终落在一个有界闭凸集$D$中, 其上$f$二阶连续可微且具有唯一稳定点$x^{\ast }$.;
2. Hessian${\nabla }^{2}f(x^{\ast })$正定, 且${\nabla }^{2}f(x)$在$x^{\ast }$的一个邻域内Lipschitz连续;
3. 范数序列 { ∥ B k ∥ } \{\Vert B_k\Vert\} {∥Bk​∥}有界;
4. SR1不跳步.

则${\lim }_{k\to \infty }{x}_k=x^{\ast }$, 且我们有$$\underset{k\to \infty }{\lim }\frac{\Vert x_{k+n+1}-x^{\ast }\Vert }{\Vert x_k-x^{\ast }\Vert }=0.$$
相比之下, BFGS算法无需第3条假设. 之前提到, SR1不保证正定性. 实际中, $B_k$可能在任一步迭代不定, 从而信赖域的边界约束会在任意大的$k$处起作用. 不过有意思的是, 可以证明在定理8的条件下, SR1在大多数时间内还是能够产生正定矩阵的. 具体说就是 lim ⁡ k → ∞ { B j } 正 定 的 指 标 数 k = 1. \lim_{k\to\infty}\frac{\{B_j\}正定的指标数}{k}=1. k→∞lim​k{Bj​}正定的指标数​=1.这点与初始Hessian近似正定与否无关.