Numerical Optimization Ch7. Large-Scale Unconstrained Optimization

第七章: 大规模无约束优化

许多现实的问题都需要对 成千上万个变量进行优化. 此时优化算法的 存储和计算量就成了主要的问题. 至今已有许多适用于大规模场景下的无约束优化算法, 其中有一些是第 三、 四、 六章中方法的直接推广, 而有一些则是基于这些方法的改动. 第 五章中的非线性共轭梯度法可以直接应用于大规模场景, 这是由于其在存储和计算上的双重俱佳表现.
第三、四章的算法 往往需要分解Hessian矩阵. 在大规模场景下, 这些分解可用 稀疏消元技术达成. 若这样的稀疏分解方法的计算量和存储量在允许的范围内, 且若能够解析地构造Hessian矩阵, 此时 牛顿法将成为解决问题的 不二之选.
然而一般说来, 分解Hessian矩阵并不可取. 在第1节中, 我们讨论 线搜索和信赖域框架下的非精确牛顿法. 我们指出, 最终得到的算法具有 喜人的全局收敛性, 基于某些参数的选取 甚至是超线性收敛的. 这些算法即使在真实的Hessian不定时仍然能求得高效的搜索方向, 且不涉及Hessian的解析计算与存储.
第2节主要处理 拟牛顿法在大规模场景下的推广. 事实上, 当真实的Hessian足够稀疏时, 拟牛顿法产生的Hessian近似仍然可能是稠密的, 从而不适用于大规模场景. 我们将介绍拟牛顿法的 有限内存变体, 其只需少量的向量计算与存储. 它们都十分 强健、经济且易于实施, 但它们的 收敛速度并不快. 我们还将在第3节简要介绍另一种方法, 它对于拟牛顿近似 $B_k$的定义能够保持稀疏性.
第4节主要讨论大型问题中目标函数的常见结构性质—— 部分可分性. 对此可采取高效的牛顿和拟牛顿法, 从而 快速收敛且过程强健. 不过有时目标函数的详细信息并不容易获取, 从而无法直接应用.

1. 非精确牛顿法

回忆原始的牛顿步是通过求解系数矩阵为$n\times n$对称矩阵的线性系统$${\nabla }^2{f}_k{p}_k^N=-\nabla f_k$$得到的. 本节我们介绍如何经济地近似求解$p_k^N$, 且近似的方向会是好的搜索方向. 这些方法主要基于共轭梯度法或Lanczos方法求解牛顿方程, 其中会对Hessian中负曲率的出现进行修正. 以下将对线搜索和信赖域两种框架分别讨论. 这一类方法我们统称为非精确牛顿法.
注意使用迭代法求解牛顿方程时, 我们

• 无需考虑分解Hessian的花销以及此过程中产生的填充 (fill-in)
• 能保证在非精确求解下仍不损失快速收敛的性质
• 计算过程不会涉及Hessian, 从而无需显式地计算或储存Hessian.

我们首先考察非精确步的计算将如何决定非精确牛顿法的局部收敛性质. 然后再基于共轭梯度法 (可能带有预处理) 考虑线搜索和信赖域的框架. 最后讨论使用Lanczos方法求解牛顿方程.

1.1 非精确牛顿法的局部收敛性

大多数迭代求解牛顿方程的算法, 其收敛均基于残差$$r_k={\nabla }^2{f}_k{p}_k+\nabla f_k,$$其中$p_k$就是非精确牛顿步. 通常若$$\Vert r_k\Vert \le {\eta }_k\Vert \nabla f_k\Vert ,$$我们就终止正在进行的迭代. 其中序列 { η k : 0 &lt; η k &lt; 1 } \{\eta_k:0&lt;\eta_k&lt;1\} {ηk​:0<ηk​<1}称为强迫序列(forcing sequence). 下面我们讨论基于以上判定方式的非精确牛顿法, 其收敛速度受强迫序列怎样的影响. 以下的两个定理不仅适用于牛顿—共轭梯度法(Newton-CG), 同样也适用于所有基于以上判别方式的非精确牛顿法.

第一个定理仅需${\eta }_k$与1保持一定的距离.
定理1 设${\nabla }^{2}f(x)$在极小点$x^{\ast }$的一个邻域内存在且连续, ${\nabla }^{2}f(x^{\ast })$正定. 考虑迭代$x_{k+1}=x_k+p_k$, 其中$p_k$满足上面的收敛判定, 且假设${\eta }_k\le \eta \in [0,1)$. 则若初始点$x_0$距离$x^{\ast }$充分近, 就有序列 { x k } \{x_k\} {xk​}收敛到$x^{\ast }$且满足$$\Vert {\nabla }^{2}f(x^{\ast })(x_{k+1}-x^{\ast })\Vert \le \hat{\eta }\Vert {\nabla }^{2}f(x^{\ast })(x_k-x^{\ast })\Vert ,$$其中$\eta <\hat{\eta }<1$.
证明: 由于${\nabla }^{2}f$在$x^{\ast }$处正定且在其附近连续, 从而存在正常数$L$对于任意充分接近$x^{\ast }$的$x_k$均有$\Vert ({\nabla }^{2}f(x_k{)}^{-1}\Vert \le L$. 从而由$p_k$满足$$\Vert p_k\Vert \le L(\Vert \nabla f_k\Vert +\Vert r_k\Vert )\le 2L\Vert \nabla f_k\Vert .$$利用Taylor定理以及${\nabla }^{2}f(x)$的连续性, 我们有$$\begin{array}{rl}\nabla f_{k+1}& =\nabla f_k+{\nabla }^2{f}_k{p}_k+{\int }_0^1[{\nabla }^{2}f(x_k+t{p}_k)-{\nabla }^{2}f(x_k)]p_k\mathrm{d}t\\ & =\nabla f_k+{\nabla }^2{f}_k{p}_k+o(\Vert p_k\Vert )\\ & =\nabla f_k+(r_k-\nabla f_k)+o(\Vert \nabla f_k\Vert )\\ & =r_k+o(\Vert \nabla f_k\Vert ).\end{array}$$在两边取范数可得,$$\Vert \nabla f_{k+1}\Vert \le {\eta }_k\Vert \nabla f_k\Vert +o(\Vert \nabla f_k\Vert )\le ({\eta }_k+o(1))\Vert \nabla f_k\Vert .$$当$x_k$距离$x^{\ast }$充分近时, $o(1)$可以用$(1-\eta )/2$控制, 从而$$\Vert \nabla f_{k+1}\Vert \le ({\eta }_k+(1-\eta )/2)\Vert \nabla f_k\Vert \le \frac{1+\eta }{2}\Vert \nabla f_k\Vert .$$这说明梯度范数序列至少以$(1+\eta )/2$的速度下降. 由光滑性的假设, $$\nabla f_k={\nabla }^{2}f(x^{\ast })(x_k-x^{\ast })+o(\Vert x_k-x^{\ast }\Vert ).$$因此, 当$x_k$离$x^{\ast }$很近时, 梯度与"缩放过"的误差${\nabla }^{2}f(x^{\ast })(x_k-x^{\ast })$相差是很小的. 因此类似的估计对 { ∥ x k ∥ } \{\Vert x_k\Vert\} {∥xk​∥}也成立. 证毕.

由于$$\frac{\Vert \nabla f_{k+1}\Vert }{\Vert \nabla f_k\Vert }\le {\eta }_k+o(1),$$因此若${\eta }_k\to 0(k\to \infty )$, 则$$\frac{\Vert \nabla f_{k+1}\Vert }{\Vert \nabla f_k\Vert }\to 0(k\to \infty ).$$这说明梯度范数序列具有Q-超线性收敛性. 对于 { ∥ x k ∥ } \{\Vert x_k\Vert\} {∥xk​∥}也有类似结论:由于$$\nabla f_{k+1}-\nabla f^{\ast }={\nabla }^{2}f({\gamma }_{k+1})(x^{k+1}-x^{\ast }),\exists {\gamma }_{k+1},$$所以$$\begin{array}{rl}\Vert x^{k+1}-x^{\ast }\Vert & =\Vert {\nabla }^{2}f({\gamma }_{k+1}{)}^{-1}\nabla f_{k+1}\Vert \\ & \le \Vert \nabla f_{k+1}\Vert =o(\Vert \nabla f_k\Vert )\\ & =o(\Vert {\nabla }^{2}f({\gamma }_k)(x^k-x^{\ast })\Vert )\\ & \le \tilde{L}o(\Vert x^k-x^{\ast }\Vert ),\exists \tilde{L}.\end{array}$$若对${\nabla }^{2}f(x)$附加上在$x^{\ast }$附近Lipschitz连续的条件, 我们就能得到二次收敛性. 此时$$\nabla f_{k+1}=r_k+O(\Vert \nabla f_k{\Vert }^2).$$选取强迫序列${\eta }_k=O(\Vert \nabla f_k\Vert )$, 我们就有$$\Vert \nabla f_{k+1}\Vert =O(\Vert \nabla f_k{\Vert }^2).$$从而有梯度范数序列的二次收敛性 ( { ∥ x k ∥ } \{\Vert x_k\Vert\} {∥xk​∥}的二次收敛性也类似). 以下第二个定理就总结了这些讨论.
定理2 设定理1的条件成立, 且设由非精确牛顿法产生的序列 { x k } \{x_k\} {xk​}收敛到$x^{\ast }$. 则当${\eta }_k\to 0$, 就有超线性收敛; 若${\nabla }^{2}f(x)$在$x^{\ast }$附近Lipschitz连续, ${\eta }_k=O(\Vert \nabla f_k\Vert )$, 则就有二次收敛.

对于超线性收敛, 我们可设置${\eta }_k=\min (0.5,\sqrt{\Vert \nabla f_k\Vert })$; 对于二次收敛, 我们可设置${\eta }_k=\min (0.5,\Vert \nabla f_k\Vert )$.

上面的结论均是局部意义上的. 它们假设序列 { x k } \{x_k\} {xk​}会最终进入$x^{\ast }$的一个邻域内, 并假定单位步长${\alpha }_k=1$, 使得全局的策略不会影响收敛速度. 下面我们说明, 非精确牛顿法可被包含在实用线搜索和信赖域的框架内, 从而有良好的局部和全局收敛性质. 我们先从线搜索算法开始.

1.2 线搜索Newton-CG算法

线搜索Newton-CG算法也被称为截断牛顿法(truncated Newton method), 其中我们使用共轭梯度法求解牛顿方程, 并以前述的非精确的牛顿法的收敛条件作为算法终止的判据. 然而共轭梯度法仅适用于正定系统, 但当$x_k$接近$x^{\ast }$时${\nabla }^2{f}_k$很可能拥有负特征值. 因此这样一来, 当算法产生负曲率的方向时, 我们就应当立即终止CG迭代. 这一改动使得CG能总是产生下降方向. 进一步地, 这也能将原始牛顿法的快速收敛性质嫁接到此处, 当然这需要步长选取满足一定的条件 (例如${\alpha }_k=1$总是头一个被检测).
下面我们给出将修正共轭梯度迭代作为内迭代的线搜索算法. 为记号方便统一, 我们将牛顿方程写成$$B_{k}p=-\nabla f_k,$$其中$B_k$在这里就是${\nabla }^2{f}_k$. 而对于CG内迭代, 则以$d_j$表示搜索方向, 以$z_j$表示产生的迭代项. 当$B_k$正定, 内迭代将收敛到牛顿步$p_k^N$, 从而求解了牛顿方程. 在每一步迭代中, 我们都定义容忍限${ϵ}_k$以保证计算出来的解的精确程度. 同时选取强迫序列为${\eta }_k=\min (0.5,\sqrt{\Vert \nabla f_k\Vert })$以得到超线性收敛. 当然这里强迫序列的选取不唯一.

算法1 (Line Search Newton-CG)
Given initial point $x_0$;
for $k=0,1,2,...$
\quad\quad Define tolerance ${ϵ}_k=\min (0.5,\sqrt{\Vert \nabla f_k\Vert })\Vert \nabla f_k\Vert$;
\quad\quad Set $z_0=0,r_0=\nabla f_k,d_0=-r_0=-\nabla f_k$;
\quad\quad for $j=0,1,2,...$
\quad\quad\quad\quad if $d_j^T{B}_k{d}_j\le 0$
\quad\quad\quad\quad\quad\quad if $j=0$
\quad\quad\quad\quad\quad\quad\quad\quad return $p_k=-\nabla f_k$;
\quad\quad\quad\quad\quad\quad else
\quad\quad\quad\quad\quad\quad\quad\quad return $p_k=z_j$;
\quad\quad\quad\quad Set ${\alpha }_j=r_j^T{r}_j/d_j^T{B}_k{d}_j$;
\quad\quad\quad\quad Set $z_{j+1}=z_j+{\alpha }_j{d}_j$;
\quad\quad\quad\quad Set $r_{j+1}=r_j+{\alpha }_j{B}_k{d}_j$;
\quad\quad\quad\quad if $\Vert r_{j+1}\Vert <{ϵ}_k$
\quad\quad\quad\quad\quad\quad return $p_k=z_{j+1}$;
\quad\quad\quad\quad Set ${\beta }_{j+1}=r_{j+1}^T{r}_{j+1}/r_j^T{r}_j$;
\quad\quad\quad\quad Set $d_{j+1}=-r_{j+1}+{\beta }_{j+1}{d}_j$;
\quad\quad end (for)
\quad\quad Set $x_{k+1}=x_k+{\alpha }_k{p}_k$, where ${\alpha }_k$ satisfies the Wolfe, Goldstein, or Armijo backtracking conditions (using ${\alpha }_k=1$ if possible)
end (for)

算法1的内循环与第五章的实用CG算法相比有以下几点不同:

1. 算法1中指定了初始点为$z_0=0$;
2. 容忍限${ϵ}_k$的使用使得CG可在非精确解上终止;
3. 负曲率检测$d_j^T{B}_k{d}_j\le 0$保证了$p_k$一定是$f$在$x_k$的下降方向. 事实上, 若此前均有$d_j^T{B}_k{d}_j>0$, 则$$p_k^T\nabla f_k=\sum _{i=0}^{j-1}\frac{\Vert r_i{\Vert }^2}{d_i^T{B}_k{d}_i}{d}_i^T\nabla f_k.$$而$$r_i=B_k{z}_i+\nabla f_k=B_k\sum _{m=0}^{i-1}{\alpha }_m{d}_m+\nabla f_k,$$因此$$d_i^T\nabla f_k=d_i^T{r}_i-\sum _{m=0}^{i-1}{\alpha }_m{d}_i^T{B}_k{d}_m=d_i^T{r}_i(\because conjugacy).$$可以验证$d_i^T{r}_i$均负. 从而$p_k$是个下降方向(Con’t d). 若在$j=0$就检测到了负曲率方向, 则$p_k=-\nabla f_k$也显然是一个下降方向(尽管曲率非正).

当然, 我们也可以用类似于第五章中的方法对算法1中的CG迭代预处理.
算法1相当适用于处理大规模问题, 但它有一个缺陷: 当Hessian矩阵${\nabla }^2{f}_k$接近奇异时, 线搜索Newton-CG方向可能会很长且质量相当不好, 从而之后尽管在线搜索中计算了很多次函数值也仅带来很小的下降. 为缓解这一问题, 我们可以标准化牛顿步, 但这样一来何谓"好的"标准化又成了问题. 我们需要考虑这样做的风险: 可能会使算法丧失原有的收敛速度.
另一个细节是, 在负曲率检测中我们希望使用的是阈值而不是单纯地判定符号, 否则无法保证舍入误差带来的后果. 同样地, 这样一来阈值的选取也成了问题. 之后将介绍的信赖域框架下的Newton-CG迭代则会更有效地处理这一问题. 从这个角度说, 我们更倾向于使用信赖域Newton-CG.
需要指明的是, 线搜索Newton-CG算法并不需要显式的Hessian矩阵, 而仅需要我们提供矩阵-向量的乘积${\nabla }^2{f}_{k}d$. 这向来可以看做CG的优点. 特别是, 当我们无法计算二阶导数或者Hessian需要巨大的存储量时, 这一特点显得尤为实用. 我们可以用第八章中的自动微分和有限差分技巧计算矩阵-向量乘积. 这样的方法被称作无Hessian牛顿法.

1.3 信赖域Newton-CG算法

在第四章中, 我们讨论了两种改进Cauchy点、近似求解信赖域子问题的方法. 这里我们定义一种修正的CG算法用于求解子问题, 可见下面的算法2. 优化$f$的完整算法则是将第四章中的信赖域框架和算法2结合起来得到, 其中每步迭代${ϵ}_k$需要适当选取.
我们使用与之前类似的记号定义信赖域子问题:$$\underset{p\in {\mathbb{R}}^n}{\min }{m}_k(p)\stackrel{def}{=}{f}_k+(\nabla f_k{)}^{T}p+\frac{1}{2}{p}^T{B}_{k}p\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}\Vert p\Vert \le {\Delta }_k,$$其中$B_k={\nabla }^2{f}_k$. 同算法1, 我们使用$d_j$表示每步CG的搜索方向, $z_j$表示其产生的迭代项.

算法2 (CG-Steihaug)
Given tolerance ${ϵ}_k>0$;
Set $z_0=0,r_0=\nabla f_k,d_0=-r_0=-\nabla f_k$;
if $\Vert r_0\Vert <{ϵ}_k$
\quad\quad return $p_k=z_0=0$;
for $j=0,1,2,...$
\quad\quad if $d_j^T{B}_k{d}_j\le 0$
\quad\quad\quad\quad Find $\tau$ such that $p_k=z_j+\tau d_j$ minimizes $m_k(p_k)$ and satisfies $\Vert p_k\Vert ={\Delta }_k$;
\quad\quad\quad\quad return $p_k$;
\quad\quad Set ${\alpha }_j=r_j^T{r}_j/d_j^T{B}_k{d}_j$;
\quad\quad Set $z_{j+1}=z_j+{\alpha }_j{d}_j$;
\quad\quad if $\Vert z_{j+1}\Vert \ge {\Delta }_k$
\quad\quad\quad\quad Find $\tau \ge 0$ such that $p_k=z_j+\tau d_j$ satisfies $\Vert p_k\Vert ={\Delta }_k$;
\quad\quad\quad\quad return $p_k$;
\quad\quad Set $r_{j+1}=r_j+{\alpha }_j{B}_k{d}_j$;
\quad\quad if $\Vert r_{j+1}\Vert <{ϵ}_k$
\quad\quad\quad\quad return $p_k=z_{j+1}$;
\quad\quad Set ${\beta }_{j+1}=r_{j+1}^T{r}_{j+1}/r_j^T{r}_j$;
\quad\quad Set $d_{j+1}=-r_{j+1}+{\beta }_{j+1}{d}_j$;
end (for)

算法中, 循环里第一个if将在算法遇到负曲率方向时终止算法, 而第二个if将在$z_{j+1}$越过信赖域边界时终止算法. 两种情况下, 算法都将返回当前方向与信赖域边界相交得到的$p_k$.
算法2中容忍限${ϵ}_k$的选取在控制整个算法的花销上至关重要. 在接近$x^{\ast }$时, 信赖域边界约束将不再积极, 从而约减为前述的非精确牛顿法. 此时我们可以用算法1中类似的方式选取${ϵ}_k$.
算法2内循环与第五章实用CG迭代的本质不同在于: 算法2终止于

1. 信赖域边界被穿越时;
2. 算法遇到负曲率方向时;
3. 满足收敛条件时.

从这个角度看, 算法2与算法1二者的内循环是很相像的.

需要说明, 算法2中指定$z_0=0$至关重要, 它的效果有二:

1. 若$\Vert \nabla f_k\Vert \ge {ϵ}_k$, 则算法2将得到$p_k$满足$m_k(p_k)\le m_k(p_k^C)$. 我们分几种情形来证明这一论述.

   • 若$d_0^T{B}_k{d}_0=(\nabla f_k{)}^T{B}_k\nabla f_k\le 0$. 则算法将返回Cauchy点.
   • 否则, 算法将定义$$z_1={\alpha }_0{d}_0=\frac{r_0^T{r}_0}{d_0^T{B}_k{d}_0}{d}_0=-\frac{(\nabla f_k{)}^T\nabla f_k}{(\nabla f_k{)}^T{B}_k\nabla f_k}\nabla f_k.$$若$\Vert z_1\Vert <{\Delta }_k$, 则$z_1$就是Cauchy点. 后续产生的迭代步必定满足$m_k(p_k)\le m_k(z_1)$. 若$\Vert z_1\Vert \ge {\Delta }_k$, 则算法也在Cauchy点终止.
     以上, 我们就能说算法2是全局收敛的.

2. 迭代项$z_j$的范数是严格递增的. 这就暗示: 只要一越过信赖域边界, 算法就可以终止, 毕竟在信赖域内再没有迭代项会给出更小的模型函数值$m_k$了. 我们利用共轭梯度法的扩张子空间性质证明这一论断.

定理3 算法2产生的序列 { z j } \{z_j\} {zj​}满足$$0=\Vert z_0\Vert <\cdots <\Vert z_j\Vert <\Vert z_{j+1}\Vert <\cdots <\Vert p_k\Vert \le {\Delta }_k.$$
证明: 我们首先证明算法2产生的序列满足$z_j^T{r}_j=0,j\ge 0;z_j^T{d}_j>0,j\ge 1$. 我们将计算$z_j$的算式显式地写出来,$$z_j=z_0+\sum _{i=0}^{j-1}{\alpha }_i{d}_i=\sum _{i=0}^{j-1}{\alpha }_i{d}_i.$$两端内积$r_j$并利用共轭梯度法的扩张子空间性质, 我们有$$z_j^T{r}_j=\sum _{i=0}^{j-1}{\alpha }_i{d}_i^T{r}_j=0.$$我们使用数学归纳法证明第二条. 由扩张子空间性质, 有$$z_1^T{d}_1=({\alpha }_0{d}_0{)}^T(-r_1+{\beta }_1{d}_0)={\alpha }_0{\beta }_1{d}_0^T{d}_0>0.$$假设$z_j^T{d}_j>0$, 下证$z_{j+1}^T{d}_{j+1}>0$. 由前面证明, $z_{j+1}^T{r}_{j+1}=0$, 从而$$\begin{array}{rl}{z}_{j+1}^T{d}_{j+1}& =z_{j+1}^T(-r_{j+1}+{\beta }_{j+1}{d}_j)\\ & ={\beta }_{j+1}{z}_{j+1}^T{d}_j\\ & ={\beta }_{j+1}(z_j+{\alpha }_j{d}_j{)}^T{d}_j\\ & ={\beta }_{j+1}{z}_j^T{d}_j+{\alpha }_j{\beta }_{j+1}{d}_j^T{d}_j.\end{array}$$由归纳假设以及${\beta }_{j+1},{\alpha }_j$为正, 得最后一个表达式为正.
下面证明定理. 若算法因$d_j^T{B}_k{d}_j\le 0$或$\Vert z_{j+1}\Vert \ge {\Delta }_k$终止, 则最后的$p_k$必定满足$\Vert p_k\Vert ={\Delta }_k$. 为涵盖所有的可能性, 我们需证明$\Vert z_j\Vert <\Vert z_{j+1}\Vert$, 其中$z_{j+1}=z_j+{\alpha }_j{d}_j,j\ge 1$. 注意到$$\Vert z_{j+1}{\Vert }^2=(z_j+{\alpha }_j{d}_j{)}^T(z_j+{\alpha }_j{d}_j)=\Vert z_j{\Vert }^2+2{\alpha }_j{z}_j^T{d}_j+{\alpha }_j^2\Vert d_j{\Vert }^2>\Vert z_j{\Vert }^2.$$证毕.

从这一定理我们发现, 算法2实际上就是在沿着从$z_1$到最终$p_k$的一条插值路径移动, 且每一步都会增加其余初始点的距离. 当$B_k={\nabla }^2{f}_k$正定时, 我们可以拿这条路径与折线算法中的路径作比较: 二者均以$-\nabla f_k$为初始方向, 均最终指向$p_k^N$直到打破信赖域边界约束.

1.4 预处理信赖域Newton-CG算法

如同第五章讨论的那般, 预处理可用于加速CG迭代. 主要就是要找到一个非奇异的矩阵$D$使得$D^{-T}{\nabla }^2{f}_k{D}^{-1}$的特征值具有更好的分布. 类似地推广定理3, 我们就可证明迭代项$z_j$的加权范数$\Vert D\cdot \Vert$是单调递增的. 为保持一致性, 我们需要在原来的范数(即2-范数)下重新叙述信赖域子问题 (注意信赖域的边界约束改变了):$$\underset{p\in {\mathbb{R}}^n}{\min }{m}_k(p)\stackrel{def}{=}{f}_k+\nabla f_k^{T}p+\frac{1}{2}{p}^T{B}_{k}p,\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}\Vert Dp\Vert \le {\Delta }_k.$$作变量替换$\hat{p}=Dp$, 并定义$${\hat{g}}_k=D^{-T}\nabla f_k,{\hat{B}}_k=D^{-T}({\nabla }^2{f}_k)D^{-1},$$我们可把问题重写为$$\underset{\hat{p}\in {\mathbb{R}}^n}{\min }{f}_k+{\hat{g}}_k^T\hat{p}+\frac{1}{2}{\hat{p}}^T{\hat{B}}_k\hat{p},\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}\Vert \hat{p}\Vert \le \Delta .$$这就与原来的问题形式相同. 我们可将算法2原样照搬用于上述子问题, 这等价于用带预处理的算法2求解改变边界约束的信赖域子问题.
在第五章中我们说到没有普适的预处理子, 但仍然有一些选项可供选择, 起哄较为有效的一种, 是不完全Cholesky分解. 正定矩阵$B$的不完全Cholesky分解是指寻找下三角矩阵$L$使得$$B=L{L}^T-R,$$其中我们需要限制$L$中元素的填充 (例如我们要求$L$具有与$B$下三角部分同样的稀疏性, 或者要求与$B$具有相似的非零元个数). 而$R$则反映了近似分解的非精确程度. 当${\nabla }^2{f}_k$不定时, 可能会使情形复杂化; 次卧我们必须同时处理不定性与稀疏性. 下面的算法结合了不完全Cholesky分解以及一种修正的Cholesky分解, 定义了一种对于信赖域Newton-CG算法的预处理子.

算法3 (Inexact Modified Cholesky)
Compute $T=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\Vert B{e}_1\Vert ,\Vert B{e}_2\Vert ,\dots ,\Vert B{e}_n\Vert )$, where $e_i$ is the $i$th coordinate vector;
Set $\bar{B}\leftarrow T^{-1/2}B{T}^{-1/2}$; Set $\beta \leftarrow \Vert \bar{B}\Vert$;
(compute a shift to ensure positive definiteness)
if ${\min }_i{b}_{ii}>0$
\quad\quad ${\alpha }_0\leftarrow 0$;
else
\quad\quad ${\alpha }_0\leftarrow \beta /2$;
for $k=0,1,2,...$
\quad\quad Attempt to apply incomplete Cholesky algorithm to obtain$$L{L}^T=\bar{B}+{\alpha }_{k}I;$$ \quad\quad if the factorization is completed successfully
\quad\quad\quad\quad stop and return $L$
\quad\quad else
\quad\quad\quad\quad ${\alpha }_{k+1}\leftarrow \max (2{\alpha }_k,\beta /2)$;
end (for)

之后我们就可设$D=L^T$.

1.5信赖域Newton-Lanczos算法(Cont’d)

算法2的一个缺陷在于, 它对负曲率方向不具有选择性, 即它会接受任何负曲率的方向, 即使这个方向能给出的函数值下降并不显著. 例如, 考虑下面的子问题$$\underset{p}{\min }m(p)=10^{-3}{p}_1-10^{-4}{p}_1^2-p_2^2,\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}\Vert p\Vert \le 1.$$在$p=0$处的最速下降方向是$(-10^{-3},0{)}^T$, 这是一个负曲率方向. 算法2便会沿着这一方向直至信赖域的边界, 给出函数值大约$10^{-3}$的下降. 然而沿着另一负曲率方向$e_2$搜索便会得到大小约为1的下降.
对于这一现象至今已提出了许多方案. 在第四章我们了解到, 当Hessian矩阵${\nabla }^2{f}_k$具有负特征值, 对应的搜索方向就会在${\nabla }^2{f}_k$的最小特征值对应的特征向量上具有相当大的分量, 这是因为$$({\nabla }^2{f}_k+\lambda I)=-\nabla f_k.$$详细可见第四章定理1. 这一特征可使得算法迅速远离不是极小点的稳定点. 一种实现的方式, 是使用第四章第三节中的迭代方法近似精确地计算信赖域子问题. 这一实现方式需要每次求搜索方向时, 求解几个具有形如$B_k+\lambda I$的系数矩阵的线性系统. 尽管这样在大规模情境下需要相当的计算量, 但它的确总是能产生良好的搜索方向.
另一种较为实用的方式是使用Lanczos算法而不是CG求解线性系统$B_{k}p=-\nabla f_k$. Lanczos算法可以看做是CG在不定系统下的推广, 有了它我们可以在收集负曲率信息时继续CG过程.
在$j$步后, Lanczos算法产生了一个$n\times j$矩阵$Q_j$, 其列向量相互正交且张成由此算法产生的Krylov子空间. 这一矩阵满足$Q_j^{T}B{Q}_j=T_j$, 其中$T_j$是三对角矩阵. 我们可以利用三对角的结构, 在$Q_j$的值域中找信赖域子问题的一个近似解. 为此, 我们需求解$$\underset{w\in {\mathbb{R}}^j}{\min }{f}_k+e_1^T{Q}_j(\nabla f_k)e_1^{T}w+\frac{1}{2}{w}^T{T}_{j}w,\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}\Vert w\Vert \le {\Delta }_k,$$再定义子问题的近似解为$p_k=Q_{j}w$. 由于$T_j$三对角, 因此上述问题可以通过分解$T_j+\lambda I$并利用第四章第三节中的迭代方法求解.
如同Newton-CG中一样, Lanczos迭代也可能会由类似于$$\Vert r_k\Vert \le {\eta }_k\Vert \nabla f_k\Vert$$的判据所终止. 我们也可以在Lanczos迭代中加入预处理以加速收敛. 要注意: 与Newton-CG算法相比, 这一求解信赖域的算法的强健性是以高昂地求解子问题为代价的. 至今为止, 我们依然不能说明Newton-Lanczos方法显著优于Steihaug算法.

2. 有限内存的拟牛顿法

我们曾经在第六章提到拟牛顿法在大规模情景下实施可能存在的问题. 有限内存的拟牛顿法对于求解Hessian难以获得或者不稀疏的大型问题相当有效. 这些算法并不储存完整的$n\times n$近似矩阵, 而是仅仅储存能够隐式地代表近似矩阵的几个长度为$n$的向量. 尽管储存上需求小, 但它们往往能带来良好(尽管是线性)的收敛性. 至今已有许多有限内存的算法, 我们仅讨论称为L-BFGS的算法. 它也正如它的名字一样, 与BFGS公式紧密相关. 这一算法的主要想法是, 只使用最近几次迭代的曲率信息(如向量$s_k,y_k$)"组建"Hessian的近似. 而更早的曲率信息则被认为对当前Hessian的近似没有什么显著的贡献, 从而被弃之不顾以节省内存.
这一节我们将首先介绍L-BFGS算法及其收敛性质, 之后说明它与非线性共轭梯度法的关联. 接着我们将介绍利用近似Hessian信息的所谓"紧表示(compact representation)"实施有限内存算法. 这一技巧不仅适用于L-BFGS而且也适用于其他的拟牛顿法的有限内存版本, 例如SR1. 最后我们讨论能够带来稀疏近似的拟牛顿更新方法.

2.1 有限内存BFGS

我们先回忆一下BFGS. BFGS算法的每一步迭代皆形如$$x_{k+1}=x_k-{\alpha }_k{H}_k\nabla f_k,$$其中${\alpha }_k$为步长, $H_k$在每步迭代中以公式$$H_{k+1}=V_k^T{H}_k{V}_k+{\rho }_k{s}_k{s}_k^T$$更新, 其中$${\rho }_k=\frac{1}{y_k^T{s}_k},V_k=I-{\rho }_k{y}_k{s}_k^T,$$$$s_k=x_{k+1}-x_k,y_k=\nabla f_{k+1}-\nabla f_k.$$由于Hessian近似的逆$H_k$中的元素分布可能会逐渐稠密, 当变量数很大时显式地对其储存和运算都是不允许的. 为规避这一问题, 我们隐式地储存$H_k$的一种修正形式, 即储存一定数目(如$m$)的校正对 { s i , y i } \{s_i,y_i\} {si​,yi​}. 此时乘积$H_k\nabla f_k$就可通过一系列的内积和求和得到. 在新的迭代项产生后, 校正对集中最"老"的的一对 { s i , y i } \{s_i,y_i\} {si​,yi​}就被新的校正对 { s k , y k } \{s_k,y_k\} {sk​,yk​}所替代. 这样, 校正对集中永远只存储最近$m$步迭代的曲率信息. 实际操作中不太大的$m$(如3~20)就可以带来较好的效果. 这就是有限内存的BFGS算法.
下面对算法进行详述. 在第$k$步迭代, 当前迭代点为$x_k$, 校正对集中有校正对 { s i , y i } , i = k − m , … , k − 1 \{s_i,y_i\},i=k-m,\ldots,k-1 {si​,yi​},i=k−m,…,k−1. 我们首先选取某个初始Hessian的近似$H_k^0$(这与标准的BFGS迭代不同, L-BFGS的初始近似是允许随着迭代而改变的) . 反复使用更新公式, 我们发现L-BFGS中的近似$H_k$满足以下关系式:$$\begin{array}{rl}{H}_k=& (V_{k-1}^T\cdots V_{k-m}^T)H_k^0(V_{k-m}\cdots V_{k-1})\\ & +{\rho }_{k-m}(V_{k-1}^T\cdots V_{k-m+1}^T)s_{k-m}{s}_{k-m}^T(V_{k-m+1}\cdots V_{k-1})\\ & +{\rho }_{k-m+1}(V_{k-1}^T\cdots V_{k-m+2}^T)s_{k-m+1}{s}_{k-m+1}^T(V_{k-m+2}\cdots V_{k-1})\\ & +\cdots \\ & +{\rho }_{k-1}{s}_{k-1}{s}_{k-1}^T.\end{array}$$由这一表示, 我们可以得到一种高效计算乘积$H_k\nabla f_k$的递归程序:

算法4 (L-BFGS two-loop recursion)
$q\leftarrow \nabla f_k$;
for $i=k-1,k-2,\dots ,k-m$
\quad\quad ${\alpha }_i\leftarrow {\rho }_i{s}_i^{T}q$;
\quad\quad $q\leftarrow q-{\alpha }_i{y}_i$;
end (for)
$r\leftarrow H_k^{0}q$;
for $i=k-m,k-m+1,\dots ,k-1$
\quad\quad $\beta \leftarrow {\rho }_i{y}_i^{T}r$;
\quad\quad $r\leftarrow r+s_i({\alpha }_i-\beta )$;
end (for)
stop with result $H_k\nabla f_k=r$.

说明:

1. 第一个循环计算完成后, $$q=(V_{k-m}\cdots V_{k-1})\nabla f_k,r=H_k^{0}q=H_k^0(V_{k-m}\cdots V_{k-1})\nabla f_k.$$此时$r$就是$H_k$表达式的第一行的右半部分与$\nabla f_k$乘积的结果. 计算过程中的${\alpha }_i$则可以显式的写为$${\alpha }_i={\rho }_i{s}_i(V_{i+1}\cdots V_{k-1})\nabla f_k.$$
2. 第二个循环乍一看有些晦涩, 我们从第一步开始解读. 当$i=k-m$时, 此时$r=H_k^0(V_{k-m}\cdots V_{k-1})\nabla f_k$, 第一步循环后$$\beta ={\rho }_{k-m}{y}_{k-m}^T{H}_k^0(V_{k-m}\cdots V_{k-1})\nabla f_k,$$$$\begin{array}{rl}r=& H_k^0(V_{k-m}\cdots V_{k-1})\nabla f_k\\ & +{\rho }_{k-m}{s}_{k-m}{s}_{k-m}^T(V_{k-m+1}\cdots V_{k-1})\nabla f_k\\ & -{\rho }_{k-m}{s}_{k-m}{y}_{k-m}^T{H}_k^0(V_{k-m}\cdots V_{k-1})\nabla f_k\end{array}$$$$\begin{array}{rl}& =V_{k-m}^T{H}_k^0(V_{k-m}\cdots V_{k-1})\nabla f_k\\ & +{\rho }_{k-m}{s}_{k-m}{s}_{k-m}^T(V_{k-m+1}\cdots V_{k-1})\nabla f_k\end{array}.$$于是我们可以看出, 第二个循环所做的, 就是计算$H_k^0$表达式中其他部分与$q$的乘积. 其中${\alpha }_i$为与第一个循环共有的部分. 如果将第二个循环中的$r\leftarrow r+s_i({\alpha }_i-\beta )$写作$r\leftarrow (r-s_i\beta )+s_i{\alpha }_i$可能会更佳易于理解.
3. 这样的表达式还能保证$s_k^T{y}_k>0$、满足割线方程吗? 答案是肯定的. 事实上由于$$V_{k-1}{y}_{k-1}=(I-{\rho }_{k-1}{y}_{k-1}{s}_{k-1}^T)y_{k-1}=0,$$因此$$H_k{y}_{k-1}={\rho }_{k-1}{s}_{k-1}{s}_{k-1}^T{y}_{k-1}=s_{k-1},$$即满足割线方程. 若$H_k^0$是正定的, 我们可以类似得到$H_k$是正定的, 从而必定有$s_k^T{y}_k>0$.
4. 考虑计算量. 在不考虑乘积$H_k^{0}q$时, 二重循环递归需要$4mn$次乘积; 若$H_k^0$是对角阵, 则再需要额外的$n$次乘积. 除了计算量友好以外, 这一递归格式还有另一个优势: 初始矩阵$H_k^0$与剩下的计算无关, 因此可以相对自由地选择且每步迭代可以不相同. 我们甚至可以隐式地选取$H_k^0$: 定义对于Hessian的某个初始近似矩阵$B_k^0$, 之后求解线性系统$B_k^{0}r=q$.
5. 这种递归形式在Broyden族的其他成员身上还未得到发展, 如SR1, DFP.

下面讨论算法细节. 首先是对$H_k^0$的选取. 一种较为高效的设定方案是令$H_k^0={\gamma }_{k}I$, 其中$${\gamma }_k=\frac{s_{k-1}^T{y}_{k-1}}{y_{k-1}^T{y}_{k-1}}.$$正如第六章所说明的, ${\gamma }_k$是利用最近的搜索方向对真实Hessian特征值的一种估计. 这一选取方式使得算出的搜索方向$p_k$是尺度良好的, 从而在大多数迭代中步长${\alpha }_k=1$都将被采用, 实现较快的收敛速度. 同时, 满足Wolfe条件或强Wolfe条件的线搜索使得BFGS更新更加稳定.

有限内存的BFGS算法叙述如下.
算法5 (L-BFGS)
Choose starting point $x_0$, interger $m>0$;
$k\leftarrow 0$;
repeat
\quad\quad Choose $H_k^0$;
\quad\quad Compute $p_k\leftarrow -H_k\nabla f_k$ from Algorithm 4;
\quad\quad Compute $x_{k+1}\leftarrow x_k+{\alpha }_k{p}_k$, where ${\alpha }_k$ is chosen to satisfy the Wolfe conditions;
\quad\quad if $k>m$
\quad\quad\quad\quad Discard the vector pair { s k − m , y k − m } \{s_{k-m},y_{k-m}\} {sk−m​,yk−m​} from storage;
\quad\quad Compute and save $s_k\leftarrow x_{k+1}-x_k,y_k=\nabla f_{k+1}-\nabla f_k$;
\quad\quad $k\leftarrow k+1$;
until convergence.

可以发现, 在L-BFGS的头$m-1$步迭代中, 如果$H_0=H_k^0$, 算法5是等价于标准的BFGS算法的.
下表展现了算法5在不同$m$下的实验结果, 其中给出了函数值和梯度的计算次数(nfg)以及总的CPU时间. 测试问题均来自CUTE问题集, 变量数用$n$表示, 算法终止标准为$\Vert \nabla f_k\Vert \le 10^{-5}$.

Problem	L-BFGS $m=3$	L-BFGS $m=5$	L-BFGS $m=17$	L-BFGS $m=29$
$n$	nfg time	nfg time	nfg time	nfg time
DIXMAANL 1500	146 16.5	134 17.4	120 28.2	125 44.4
EIGENALS 110	821 21.5	569 15.7	363 16.2	168 12.5
FREUROTH 1000	>999 —	>999 —	69 8.1	38 6.3
TRIDIA 1000	876 46.6	611 41.4	531 84.6	462 127.1

上表说明, 储存的信息越多, 所需要的函数值就越少; 但与此同时储存的耗费也在增加, 且最佳的CPU时间均在较小的$m$上获得. 显然, 最优的$m$是依赖于问题本身的.
由于其他竞争的算法并不高效, 算法5往往是处理大规模问题(且真实Hessian并不稀疏)的不二之选. 特别地, 在这样的情况下计算精确的Hessian并使用牛顿法求解非常不现实. L-BFGS算法可能要比无Hessian的牛顿法(如Newton-CG)更好, 后者的Hessian-向量乘积利用有限差分或自动微分计算.
L-BFGS的主要缺点在于, 它在病态问题上的收敛相当慢, 也就是说它的表现依赖于问题的条件数. 从这个角度说, 非线性共轭梯度法要比有限内存的拟牛顿法更好.

2.2 与共轭梯度法的关系

有限内存的算法的提出是希望能够改进非线性共轭梯度法, 早期的工作中共轭梯度法比拟牛顿法更加常用. 二者的关联成为了无记忆BFGS(memoryless BFGS)迭代的基础.
我们首先考虑HS形式的非线性共轭梯度法, 其中$s_k={\alpha }_k{p}_k,y_k=\nabla f_{k+1}-\nabla f_k,{\beta }_{k+1}=\frac{\nabla f_{k+1}^T{y}_k}{y_k^T{p}_k}$. 所以有$$p_{k+1}=-\nabla f_{k+1}+\frac{\nabla f_{k+1}^T{y}_k}{y_k^T{p}_k}{p}_k=-\left(I-\frac{s_k{y}_k^T}{y_k^T{s}_k}\right)\nabla f_{k+1}\equiv -{\hat{H}}_{k+1}\nabla f_{k+1}.$$上述式子仿照了拟牛顿迭代的写法, 但这里的${\hat{H}}_{k+1}$既不是对称矩阵也不是正定矩阵. 我们自然可以将其对称化为${\hat{H}}_{k+1}^T{\hat{H}}_{k+1}$, 但这样一来又不再满足割线方程${\hat{H}}_{k+1}{y}_k=s_k$. 一种对称、正定且满足割线方程的迭代矩阵是$$H_{k+1}=\left(I-\frac{s_k{y}_k^T}{y_k^T{s}_k}\right)\left(I-\frac{y_k{s}_k^T}{y_k^T{s}_k}\right)+\frac{s_k{s}_k^T}{y_k^T{s}_k}.$$这一矩阵恰巧就是对单位阵做一次BFGS更新所得到的矩阵. 因此, 产生搜索方向为$p_{k+1}=-H_{k+1}\nabla f_{k+1}$的算法可以看做是一种"无记忆"的BFGS算法, 即前一次的Hessian近似总是设置成单位阵, 且仅用最近的一次曲率信息$(s_k,y_k)$更新它. 换句话说, 这就是算法5中$m=1,H_k^0=I$的特例.
如果我们考虑带精确线搜索的无记忆BFGS算法, 则这种联系将变得更加直接. 这是因为对$\forall k$, 我们有$\nabla f_{k+1}^T{p}_k=0$. 从而$$p_{k+1}=-H_{k+1}\nabla f_{k+1}=-\nabla f_{k+1}+\frac{\nabla f_{k+1}^T{y}_k}{y_k^T{p}_k}{p}_k,$$这不就是HS共轭梯度法吗? 进一步地, 当$\nabla f_{k+1}^T{p}_k=0$时, HS公式就变成了PR公式. 即使精确线搜索的假设不切实际, 但BFGS公式以这种方式和共轭梯度法联系在一起还是非常有趣的一件事.

2.3 一般的有限内存更新

有限内存的拟牛顿近似在许多优化问题中都很有用处. 算法5即L-BFGS算法是一种求解无约束优化问题的线搜索算法, 其每步迭代将更新Hessian矩阵的逆$H_k$; 而信赖域方法则需要对于Hessian本身的近似$B_k$. 我们也希望能够对于SR1公式建立起有限内存的版本, 这将成为BFGS的一种有效的替代.
本节中我们考虑一般情形下的有限内存更新. 我们将拟牛顿矩阵以紧(或外积)的形式表示, 从而可以将有限内存的方法方便地推广到所有(广泛使用的)拟牛顿更新公式及它们的逆形式. 紧表示在设计有约束优化的有限内存方法时同样有用, 当然这就涉及到对于Lagrange函数的Hessian或者既约Hessian的近似. 我们将在较后的内容中介绍.
以下我们仅考虑会连续更新校正对$(s_k,y_k)$的有限内存方法. 与之相比有, 我们也可以一直储存校正对直到内存炸裂, 之后再释放所有的内存(或者保留一个校正对), 重新开始迭代. 实际的数值实验告诉我们, 这样的方式效率并不高.
本节中我们以$B_k$表示对于Hessian的近似, 以$H_k$表示对于Hessian的逆的近似. 我们总有$B_k^{-1}=H_k$.

2.3.1 BFGS更新的紧表示

下面我们介绍一种基于拟牛顿紧表示的有限内存更新. 我们以BFGS的$B_k$形式为例.

定理4 设$B_0$对称正定, 并假设校正对 { s i , y i } i = 0 k − 1 \{s_i,y_i\}_{i=0}^{k-1} {si​,yi​}i=0k−1​满足$s_i^T{y}_i>0$. 设$B_k$为对$B_0$做$k$次BFGS更新得到的矩阵. 于是我们有$$B_k=B_0-\left[\begin{array}{cc}{B}_0{S}_k& Y_k\end{array}\right]{\left[\begin{array}{cc}{S}_k^T{B}_0{S}_k& L_k\\ L_k^T& -D_k\end{array}\right]}^{-1}\left[\begin{array}{c}{S}_k^T{B}_0\\ Y_k^T\end{array}\right],$$其中$S_k,Y_k$为$n\times k$矩阵, 定义为$$S_k=[s_0,\dots ,s_{k-1}],Y_k=[y_0,\dots ,y_{k-1}],$$$L_k,D_k$为$k\times k$矩阵$$(L_k{)}_{i,j}=\{\begin{array}{ll}{s}_{i-1}^T{y}_{j-1}& i>j,\\ 0& otherwise,\end{array}$$$$D_k=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}[s_0^T{y}_0,\dots ,s_{k-1}^T{y}_{k-1}].$$

与之平行的有如下对$H_k$的定理.

定理4’ 设$H_0$对称正定, 并假设校正对 { s i , y i } i = 0 k − 1 \{s_i,y_i\}_{i=0}^{k-1} {si​,yi​}i=0k−1​满足$s_i^T{y}_i>0$. 设$H_k$为对$H_0$做$k$次BFGS更新得到的矩阵. 于是我们有$$H_k=H_0+\left[\begin{array}{cc}{S}_k& H_0{Y}_k\end{array}\right]\left[\begin{array}{cc}{R}_k^{-T}(D_k+Y_k^T{H}_0{Y}_k)R_k^{-1}& -R_k^{-T}\\ -R_k^{-1}& O\end{array}\right]\left[\begin{array}{c}{S}_k^T\\ Y_k^T{H}_0\end{array}\right],$$其中$S_k,Y_k,D_k$如之前定理4中定义, $R_k$为$k\times k$矩阵$$(R_k{)}_{i,j}=\{\begin{array}{ll}{s}_{i-1}^T{y}_{j-1}& i\le j,\\ 0& otherwise.\end{array}$$
我们从证明定理4’出发, 最后用Sherman-Morrison-Woodbury公式即可得到定理4. 为证明定理4’, 我们先介绍一个刻画多个投影矩阵乘积的引理.

引理 $V_0\cdots V_{k-1}=I-Y_k{R}_k^{-1}{S}_k^T,k\ge 1$.
证明: 使用数学归纳法证明.
i. 当$k=1$时, 显然公式成立;
ii. 假设对于$k$, 上述公式成立, 即有$V_0\cdots V_{k-1}=I-Y_k{R}_k^{-1}{S}_k^T$, 下面证明$k+1$的情形. 从右端开始, 首先考察$R_{k+1}$. 从$R_k$的定义可得$$R_{k+1}=\left[\begin{array}{cc}{R}_k& S_k^T{y}_k\\ & {\rho }_k^{-1}\end{array}\right]\Rightarrow R_{k+1}^{-1}=\left[\begin{array}{cc}{R}_k^{-1}& -{\rho }_k{R}_k^{-1}{S}_k^T{y}_k\\ & {\rho }_k\end{array}\right].$$因此$$\begin{array}{rl}I-Y_{k+1}{R}_{k+1}^{-1}{S}_{k+1}^T& =I-\left[\begin{array}{cc}{Y}_k& y_k\end{array}\right]\left[\begin{array}{cc}{R}_k^{-1}& -{\rho }_k{R}_k^{-1}{S}_k^T{y}_k\\ & {\rho }_k\end{array}\right]\left[\begin{array}{c}{S}_k^T\\ s_k^T\end{array}\right]\\ & =I-(Y_k{R}_k^{-1}{S}_k^T-{\rho }_k{Y}_k{R}_k^{-1}{S}_k^T{y}_k{s}_k^T+{\rho }_k{y}_k{s}_k^T)\\ & =(I-Y_k{R}_k^{-1}{S}_k^T)(I-{\rho }_k{y}_k{s}_k^T)\\ & =V_0\cdots V_{k-1}{V}_k.\end{array}$$最后一步使用归纳假设. 证毕.

定理4’的证明: 我们将$H_k$的更新公式简记为$$H_k=M_k+N_k,$$其中$$M_0=H_0,M_{k+1}=V_k^T{M}_k{V}_k,$$$$N_1={\rho }_0{s}_0{s}_0^T,N_{k+1}=V_k^T{N}_k{V}_k+{\rho }_k{s}_k{s}_k^T.$$由引理, $$\begin{array}{rl}{M}_k& =V_{k-1}^T{M}_{k-1}{V}_{k-1}\\ & =\cdots \\ & =V_{k-1}^T{V}_{k-2}^T\cdots V_0^T{H}_0{V}_0\cdots V_{k-2}{V}_{k-1}\\ & =(V_0\cdots V_{k-1}{)}^T{H}_0(V_0\cdots V_{k-1})\\ & =(I-S_k{R}_k^{-T}{Y}_k^T)H_0(I-Y_k{R}_k^{-1}{S}_k^T)\\ & =H_0-S_k{R}_k^{-T}{Y}_k^T{H}_0-H_0{Y}_k{R}_k^{-1}{S}_k^T+S_k{R}_k^{-T}{Y}_k^T{H}_0{Y}_k{R}_k^{-1}{S}_k^T.\end{array}$$对比结果, 我们只需要证明$$N_k=S_k{R}_k^{-T}{D}_k{R}_k^{-1}{S}_k^T.$$同样用数学归纳法证明.
i. 当$k=1$, 显然有$N_1=S_1{R}_1^{-T}{D}_1{R}_1^{-1}{S}_1^T={\rho }_0{s}_0{s}_0^T$.
ii. 假设结论对于$k$成立. 下面对$k+1$证明. 由$N_k$的定义,$$\begin{array}{rl}{N}_{k+1}& =V_k^T{N}_k{V}_k+{\rho }_k{s}_k{s}_k^T\\ & =V_k^T{S}_k{R}_k^{-T}{D}_k{R}_k^{-1}{S}_k^T{V}_k+{\rho }_k{s}_k{s}_k^T.\end{array}$$其中$$\begin{array}{rl}{R}_k^{-1}{S}_k^T{V}_k& =R_k^{-1}{S}_k^T(I-{\rho }_k{y}_k{s}_k^T)\\ & =R_k^{-1}{S}_k^T-{\rho }_k{R}_k^{-1}{S}_k^T{y}_k{s}_k^T\\ & =\left[\begin{array}{cc}I& 0\end{array}\right]R_{k+1}^{-1}\left[\begin{array}{c}{S}_k^T\\ s_k^T\end{array}\right]\\ & =\left[\begin{array}{cc}I& 0\end{array}\right]R_{k+1}^{-1}{S}_{k+1}^T.\end{array}$$而利用$R_{k+1}$的表达式, 我们又可以将$s_k$写成$$s_k=S_{k+1}{R}_{k+1}^{-T}{e}_{k+1}/{\rho }_k.$$因此$$\begin{array}{rl}{N}_{k+1}=& S_{k+1}{R}_{k+1}^T\left[\begin{array}{c}I\\ 0\end{array}\right]D_k\left[\begin{array}{cc}I& 0\end{array}\right]R_{k+1}^{-1}{S}_{k+1}^T\\ & +S_{k+1}{R}_{k+1}^{-T}\left[\begin{array}{cccc}0& & & \\ & \ddots & & \\ & & 0& \\ & & & {\rho }_k^{-1}\end{array}\right]R_{k+1}^{-1}{S}_{k+1}^T\end{array}$$$$=S_{k+1}{R}_{k+1}^{-T}{D}_{k+1}{R}_{k+1}^{-1}{S}_{k+1}^T.$$证毕.

定理4的证明: 将定理4’中$H_k$的表达式记作$$H_k=H_0+U_k{C}_k{U}_k^T,$$其中$$U_k=\left[\begin{array}{cc}{S}_k& H_0{Y}_k\end{array}\right],$$$$C_k=\left[\begin{array}{cc}{R}_k^{-T}(D_k+Y_k^T{H}_0{Y}_k)R_k^{-1}& -R_k^{-T}\\ -R_k^{-1}& O\end{array}\right].$$从而易得$$C_k^{-1}=\left[\begin{array}{cc}O& -R_k\\ -R_k^T& -(D_k+Y_k^T{H}_0{Y}_k)\end{array}\right].$$由Sheman-Morrison-Woodbury公式可得$$\begin{array}{rl}{B}_k& =B_0-B_0{U}_k(I+C_k{U}_k^T{B}_0{U}_k{)}^{-1}{C}_k{U}_k^T{B}_0\\ & =B_0-B_0{U}_k(C_k^{-1}+U_k^T{B}_0{U}_k{)}^{-1}{U}_k^T{B}_0,\end{array}$$其中$$\begin{array}{rl}{U}_k^T{B}_0{U}_k& =\left[\begin{array}{c}{S}_k^T\\ Y_k^T{H}_0\end{array}\right]B_0\left[\begin{array}{cc}{S}_k& H_0{Y}_k\end{array}\right]\\ & =\left[\begin{array}{cc}{S}_k^T{B}_0{S}_k& S_k^T{Y}_k\\ Y_k^T{S}_k& Y_k^T{H}_0{Y}_k\end{array}\right](\because B_0{H}_0=I).\end{array}$$所以$$C_k^{-1}+U_k^T{B}_0{U}_k=\left[\begin{array}{cc}{S}_k^T{B}_0{S}_k& S_k^T{Y}_k-R_k\\ Y_k^T{S}_k-R_k^T& -D_k\end{array}\right].$$由$L_k,R_k$的定义可得$R_k+L_k=S_k^T{Y}_k$, 代入即可得证.

注意$s_i^T{y}_i>0,i=0,1,\dots ,k-1$可以保证$B_k$表达式中中间矩阵的可逆性, 从而这一表达式是良定的. 事实上$$\left[\begin{array}{cc}{S}_k^T{B}_0{S}_k& L_k\\ L_k^T& -D_k\end{array}\right]\sim \left[\begin{array}{cc}{S}_k^T{B}_0{S}_k+L_k{D}_k^{-1}{L}_k^T& O\\ L_k^T& -D_k\end{array}\right],$$因此我们仅需关注$S_k^T{B}_0{S}_k+L_k{D}_k^{-1}{L}_k^T$的可逆性. 显然这一矩阵是半正定的($D_k$的对角元均正). 对$\forall u$, 若$$u^T(S_k^T{B}_0{S}_k+L_k{D}_k^{-1}{L}_k^T)u=0,$$则必定有$$L_k^{T}u=0,S_{k}u=0.$$因为$L_k^T=Y_k^T{S}_k-R_k^T$, 所以这又得到$$R_k^{T}u=0.$$但$R_k$是可逆阵, 从而$u=0$. 因此$S_k^T{B}_0{S}_k+L_k{D}_k^{-1}{L}_k^T$正定从而可逆.
这一表示表面上仅仅简洁而毫无实质益处, 实际上在考察有限内存更新时它将起重要的作用. 下面对这一点进行说明.
L-BFGS算法中, 我们存储最近$m$个校正对 { s i , y i } \{s_i,y_i\} {si​,yi​}并在每步迭代结束后剔除最老的和增加最新的校正对.
在头$m$步迭代中, 定理4的更新策略可以不加改动地应用于上, 除了通常我们选取初始矩阵为$B_k^0={\delta }_{k}I$, 其中${\delta }_k=1/{\gamma }_k$.
在之后的第$k>m$步迭代中, 定理4的更新策略就需要有所改进以反映校正对中的变化. 定义$n\times m$矩阵$S_k,Y_k$为$$S_k=[s_{k-m},\dots ,s_{k-1}],Y_k=[y_{k-m},\dots ,y_{k-1}],$$从而根据定理4, 对$B_0^{(k)}={\delta }_{k}I$做$m$次更新后得到的$B_k$为$$B_k={\delta }_{k}I-\left[\begin{array}{cc}{\delta }_k{S}_k& Y_k\end{array}\right]{\left[\begin{array}{cc}{\delta }_k{S}_k^T{S}_k& L_k\\ L_k^T& -D_k\end{array}\right]}^{-1}\left[\begin{array}{c}{\delta }_k{S}_k^T\\ Y_k^T\end{array}\right],$$其中$L_k,D_k$为$m\times m$矩阵, 定义为$$(L_k{)}_{i,j}=\{\begin{array}{ll}(s_{k-m-1+i}{)}^T(y_{k-m-1+j})& i>j,\\ 0& otherwise,\end{array}$$$$D_k=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}[s_{k-m}^T{y}_{k-m},\dots ,s_{k-1}^T{y}_{k-1}].$$在产生新的迭代点$x_{k+1}$后, 将从$S_k$中删除$s_{k-m}$加入$s_k$得到$S_{k+1}$, 类似得到$Y_{k+1}$. $L_k,D_k$也做相应变动成为$L_{k+1},D_{k+1}$.
由于中间矩阵的维数$2m$很小, 因此其分解需要的计算量基本可以忽略不计. 紧表示的思想在于, 对于初始矩阵的校正可表示为两个长条形矩阵——$[{\delta }_k{S}_k{Y}_k]$及其转置——的外积, 其中中间以一个$2m\times 2m$的小矩阵"加权". 可见下图.

$B_k$的有限内存的更新策略大概需要$2mn+O(m^3)$的计算量(其中$O(m^3)$表示分解的计算量, $2mn$则表示两长条形矩阵乘积的计算量, 不包括加法). 矩阵-向量乘积$B_{k}v$则需要$(4m+1)n+O(m^2)$次乘法. 这些都表明有限内存更新和计算乘积都是十分经济的.
这一近似$B_k$可用于无约束优化中的信赖域框架中, 或者更甚, 用于带边界约束的优化和一般约束优化的问题中. 前者, 需要重复计算$B_k$向约束梯度张成子空间的投影. 后者, 可能需使用紧表示来近似Lagrange函数的Hessian矩阵.
我们也可以类似地将上述想法用于$H_k$的有限内存更新和乘积计算中. 计算量是类似的. 不过在BFGS的场景下, $B_k$的使用将多一层线性系统的求解, 而$H_k$只需要矩阵-向量乘积.
SR1公式也有相应的紧表示. 对初始对称矩阵$B_0$使用校正对 { s i , y i } i = 0 k − 1 \{s_i,y_i\}_{i=0}^{k-1} {si​,yi​}i=0k−1​做$k$次SR1更新, 得到的矩阵$B_k$可以表示成$$B_k=B_0+(Y_k-B_0{S}_k)(D_k+L_k+L_k^T-S_k^T{B}_0{S}_k{)}^{-1}(Y_k-B_0{S}_k{)}^T.$$由Sherman-Morrison-Woodbury公式可得到对应的$H_k$形式(当然也可以直接利用SR1的自对偶性, 将上式中的$B,s,y$分别换成$H,y,s$即可).

2.3.2 展开更新

我们可能会怀疑有限内存更新能否以更简单的方式实施. 事实上, 我们下面将说明, 有限内存BFGS更新的最显然直接的实施策略要比基于紧表示的策略要昂贵得多.
BFGS公式可以写作$$B_{k+1}=B_k-a_k{a}_k^T+b_k{b}_k^T,$$其中$$a_k=\frac{B_k{s}_k}{(s_k^T{B}_k{s}_k{)}^{\frac{1}{2}}},b_k=\frac{y_k}{(y_k^T{s}_k{)}^{\frac{1}{2}}}.$$我们可以继续储存校正对 { s i , y i } \{s_i,y_i\} {si​,yi​}但用上面的公式计算矩阵-向量乘积. 有限内存的BFGS算法则对应$$B_k=B_k^0+\sum _{i=k-m}^{k-1}[b_i{b}_i^T-a_i{a}_i^T].$$而向量对 { a i , b i } i = k − m k − 1 \{a_i,b_i\}_{i=k-m}^{k-1} {ai​,bi​}i=k−mk−1​将以如下方式从储存的 { s i , y i } i = k − m k − 1 \{s_i,y_i\}_{i=k-m}^{k-1} {si​,yi​}i=k−mk−1​中得到:

算法6 (Unrolling the BFGS formula)
for $i=k-m,k-m+1,\dots ,k-1$
\quad\quad $b_i\leftarrow y_i/(y_i^T{s}_i{)}^{1/2}$;
\quad\quad $a_i\leftarrow B_k^0{s}_i+{\sum }_{j=k-m}^{i-1}\left[(b_j^T{s}_i)b_j-(a_j^T{s}_i)a_j\right]$;
\quad\quad $a_i\leftarrow a_i/(s_i^T{a}_i{)}^{1/2}$;
end (for)

注意向量$a_i$在每步迭代中都得重新计算. 而相比之下, $b_i$和$b_j^T{s}_i$可以储存以备后用, 每步迭代只需要计算最新的$b_{k-1},b_j^T{s}_{k-1}$.
下面来考虑计算量. 假设$B_k^0=I$, 则为了更新有限内存, 大约需要$\frac{3}{2}{m}^{2}n$次运算. 而实际的$B_{m}v$则需要$4mn$次运算. 总的来说, 这样的方式要比紧表示低效.

3. 稀疏拟牛顿更新

下面我们介绍一种大规模情境下较为吸引人的拟牛顿法: 我们要求拟牛顿近似矩阵$B_k$拥有与真实Hessian矩阵相同(或差不多)的稀疏模式. 这一方法可以减少算法所需的存储量并可能给出对于Hessian更加精确的近似.
假设我们已经知道在定义域的某点处Hessian有哪些元素不为零, 即我们已知 Ω = d e f { ( i , j ) ∣ [ ∇ 2 f ( x ) ] i j ≠ 0 &ThinSpace;&ThinSpace; f o r &ThinSpace;&ThinSpace; s o m e &ThinSpace;&ThinSpace; x &ThinSpace;&ThinSpace; i n &ThinSpace;&ThinSpace; t h e &ThinSpace;&ThinSpace; d o m a i n &ThinSpace;&ThinSpace; o f &ThinSpace;&ThinSpace; f } . \Omega\xlongequal{def}\{(i,j)|[\nabla^2f(x)]_{ij}\ne0 \mathrm{\,\,for\,\,some\,\, }x\mathrm{\,\, in\,\, the\,\, domain\,\, of\,\, }f\}. Ωdef {(i,j)∣[∇2f(x)]ij​̸​=0forsomexinthedomainoff}.再假定当前的Hessian近似$B_k$反映了真实Hessian的非零结构, 即$(B_k{)}_{ij}=0,\forall (i,j)\in \Omega$. 在更新$B_k$到$B_{k+1}$时, 我们应当要求$B_{k+1}$

• 满足割线方程;
• 具有相同的稀疏模式;
• 与$B_k$尽量相近.

确切地说, 我们定义$B_{k+1}$为以下二次规划的解:$$\begin{gathered} \underset{B}{\min }\Vert B-B_k{\Vert }_F^2=\sum _{(i,j)\in \Omega }[B_{ij}-(B_k{)}_{ij}{]}^2, \\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}B{s}_k=y_k,B=B^T,B_{ij}=0,\forall (i,j)̸\in \Omega . \end{gathered}$$我们可以证明解$B_{k+1}$可以通过求解一个$n\times n$的稀疏模式为$\Omega$的线性系统得到. 求得$B_{k+1}$后, 我们可将其套进信赖域的框架, 得到下一迭代点$x_{k+1}$. 注意这里用信赖域不用线搜索的原因是, 这里$B_{k+1}$不保证正定.
由于这一方法有诸多弊端, 我们不详细讨论它的细节了. 其一, 它不具有尺度不变性; 其二, 也是最重要的一点, 它的实际表现往往比较糟糕. 主要原因在于, 目标函数可能本身就是个不恰当的模型, 从而会产生较差的Hessian近似.
一种代替方案的想法是, 放松对于割线方程的要求, 仅需它在最近几步上均"近似"满足而不是在最近的一步上严格满足, 即$B_{k+1}$为如下问题的解:$$\begin{gathered} \underset{B}{\min }\Vert B{S}_k-Y_k{\Vert }_F^2 \\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}B=B^T,B_{ij}=0,\forall (i,j)̸\in \Omega . \end{gathered}$$注意算得的$B_{k+1}$仍不保证正定, 从而信赖域的框架是必须的. 这一凸优化问题有解, 不过较难计算. 而且这一方法可能会产生奇异的或者是病态的(即条件数较大)Hessian近似. 即使它要比前一种方法要好, 但其在大型问题上的表现仍然不尽人意.

4. 部分可分函数及其相应算法

4.1 可分与部分可分

在可分的无约束优化问题中, 目标函数可以分解成一些更简单的、可独立优化的函数之和. 例如对于$$f(x)=f_1(x_1,x_3)+f_2(x_2,x_4,x_6)+f_3(x_5),$$我们可以通过独立地极小化$f_i,i=1,2,3$求出最终的最优点$x$和最优值. 求解$m$个低维的优化问题要比求解一个$n$维的优化问题要经济得多.
不过在许多大型问题中, 目标函数$f:{\mathbb{R}}^n\to \mathbb{R}$是不可分的, 但它仍然可以写成一些更简单的函数之和, 这些函数被称为基本函数. 每一基本函数都具有如下性质: 当我们沿着许多线性无关的方向移动时, 函数值不变. 如果这一性质成立, 我们就称$f$是部分可分的. 所有Hessian稀疏的函数均是部分可分的, 但实际上不稀疏的也不少. 部分可分性使我们能够更加经济地表示问题、进行高效地自动微分以及高效的拟牛顿更新.
部分可分性最简单的一种形式体现在目标函数可以写作$$f(x)=\sum _{i=1}^{ne}{f}_i(x)$$时, 其中每一个基本函数$f_i$仅依赖于$x$的少许分量. 从而梯度$\nabla f_i$和Hessian${\nabla }^2{f}_i$也仅有少量的非零元素. 我们有$$\nabla f(x)=\sum _{i=1}^{ne}\nabla f_i(x),{\nabla }^{2}f(x)=\sum _{i=1}^{ne}{\nabla }^2{f}_i(x).$$

4.2 算法讨论

一个自然的问题是: 分别对每个基本函数的Hessian${\nabla }^2{f}_i(x)$做拟牛顿近似是否比直接对整个${\nabla }^{2}f(x)$近似更加高效呢? 我们将说明这个答案是肯定的, 当然前提是拟牛顿近似能够完全挖掘每个基本函数Hessian的结构.
我们从一个简单的例子出发. 考虑目标函数$$\begin{array}{rl}f(x)& =(x_1-x_3^2{)}^2+(x_2-x_4^2{)}^2+(x_3-x_2^2{)}^2+(x_4-x_1^2{)}^2\\ & \triangleq f_1(x)+f_2(x)+f_3(x)+f_4(x).\end{array}$$这里每个基本函数$f_i$的Hessian均是$4\times 4$的稀疏奇异矩阵, 且仅有4个非零元素.
下以$f_1$为例, 其他基本函数则是类似的. 尽管$f_1$看上去依赖于$x$的所有分量, 但实际上它仅依赖于$x_1,x_3$, 此时称这两个为$f_1$的基本变量. 我们将基本变量集成到一个向量$x_{[1]}$中, 即$$x_{[1]}=\left[\begin{array}{c}{x}_1\\ x_3\end{array}\right],$$且它与$x$的关系是$$x_{[1]}=U_{1}x,U_1=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right].$$若定义函数${\varphi }_1$为$${\varphi }_1(z_1,z_2)=(z_1-z_2^2{)}^2,$$则我们有$f_1(x)={\varphi }_1(U_{1}x)$. 应用链式法则, 有$$\nabla f_1(x)=U_1^T\nabla {\varphi }_1(U_{1}x),{\nabla }^2{f}_1(x)=U_1^T{\nabla }^2{\varphi }_1(U_{1}x)U_1.$$具体写出来就是$${\nabla }^2{\varphi }_1(U_{1}x)=\left[\begin{array}{cc}2& -4x_3\\ -4x_3& 12x_3^2-4x_1\end{array}\right],{\nabla }^2{f}_1(x)=\left[\begin{array}{cccc}2& 0& -4x_3& 0\\ 0& 0& 0& 0\\ -4x_3& 0& 12x_3^2-4x_1& 0\\ 0& 0& 0& 0\end{array}\right].$$矩阵$U_1$称为紧致矩阵, 它将低维函数${\varphi }_1$的导数信息映射到基本函数的导数信息.
下面就是关键的想法: 与其对${\nabla }^2{f}_1$做拟牛顿近似, 我们转而近似${\nabla }^2{\varphi }_1$得到$2\times 2$的$B_{[1]}$并使用$${\nabla }^2{f}_1(x)=U_1^T{\nabla }^2{\varphi }_1(U_{1}x)U_1$$将它转换到对于${\nabla }^2{f}_1$的近似. 假设现在从$x$迭代至$x^{+}$, 我们有$$s_{[1]}=x_{[1]}^{+}-x_{[1]},y_{[1]}=\nabla {\varphi }_1(x_{[1]}^{+})-\nabla {\varphi }_1(x_{[1]}).$$接着使用BFGS或SR1公式得到新的近似$B_{[1]}^{+}$. 接着转换$${\nabla }^2{f}_1(x^{+})\approx U_1^T{B}_{[1]}^{+}{U}_1.$$依此对其他基本函数进行操作, 我们有$$f(x)=\sum _{i=1}^{ne}{\varphi }_i(U_{i}x),B=\sum _{i=1}^{ne}{U}_1^T{B}_{[i]}{U}_i.$$我们可以将这一近似Hessian用于信赖域的框架中, 进而得到线性系统$$B_k{p}_k=-\nabla f_k$$的近似解.
注意我们不需要显式地得到$B_k$, 我们可以使用从而梯度法求解, 其中只需要矩阵-向量乘积, 这可以用$U_i$和$B_{[i]}$分步操作.
为阐释这一从基本函数角度出发的更新技术的优势所在, 我们现在考虑具有同样形式的问题, 但此时牵涉1000个变量. 函数${\varphi }_i$依然只依赖于2个变量, 从而每个Hessian近似$B_{[i]}$是$2\times 2$矩阵. 经过几步迭代后, 我们已经有了充足的信息$s_{[i]}$使得$B_{[i]}$能够成为${\nabla }^2{\varphi }_i$的良好近似, 从而类比地可以由${\nabla }^{2}f(x)$的优良近似. 相反, 若对整个Hessian近似而不考虑部分可分性, 我们将需要近似一个$1000\times 1000$的矩阵. 当变量数$n$很大时, 在拟牛顿近似良好之前将需要很多次迭代. 这与基于部分可分性设计的算法就要逊色得多.
当然实际中并不总允许我们使用BFGS公式更新$B_{[i]}$, 这是因为我们不能保证隐含条件$s_{[i]}^T{y}_{[i]}>0$. 这就是说, 即使${\nabla }^{2}f(x)$在解$x^{\ast }$处至少是半正定的, 依然可能有一些${\nabla }^2{\varphi }_i(\cdot )$是不定的. 一种解决方法是使用SR1公式更新.
以上基于部分可分性的拟牛顿方法的主要缺陷在于:

• 我们需要求解$B_k{p}_k=-\nabla f_k$获取步长, 而这一步的计算量相当于计算牛顿步. 这一线性系统的求解是必须的, 因为我们没能得到对于$H_k$的部分可分更新方式.
• 其次, 识别函数结构中的部分可分性也是一个难点. 当我们识别出了问题的最佳部分分解时, 拟牛顿法的效果将会很好.