从零发明数值最优化 4：拟牛顿法

发表于 更新于 分类于 ， ， ， ，

牛顿法要存的 Hessian 矩阵实在是太占地方了，解起来也很耗时。而且考虑到复杂地形中没法做到处处正定，牛顿法就不得不在曲率上做出很多妥协，把它硬掰成正定的，也拟合的没有很好。

改良的想法是我们干脆不算 Hessian 矩阵了，而是使用相邻迭代点处的一阶梯度差来积累曲率信息。反正一阶展开还是准的，二阶曲率的准确性要求倒也没那么迫切了。或者说，在工程上，我们愿意多牺牲一点二阶曲率来换取更快的时间和更少的空间。

根据这种思想，我们构造近似模型 f(xk+pk)≈mk(pk)=f(xk)+∇f(xk)⊤pk+12pk⊤Bkpkf(x_k + p_k) \approx m_{k}(p_k) = f(x_k) + \nabla f(x_k)^\top p_k + \frac{1}{2} p_k^\top B_k p_k。接下来的任务就是求出一个合适的 BkB_k。

割线约束

为使模型包含历史观测信息，强制要求该二次模型在 xk−1x_{k-1} 处的预测梯度与当时实际计算出的真实梯度 ∇f(xk−1)\nabla f(x_{k-1}) 一致，则有

∇mk(pk⋆)=∇f(xk−1),pk⋆=xk−1−xk\nabla m_k(p_k^\star) = \nabla f(x_{k-1}), \quad p_k^\star = x_{k-1} - x_k

化简，得到

∇f(xk)+Bk(xk−1−xk)=∇f(xk−1)\nabla f(x_k) + B_k (x_{k-1} - x_k) = \nabla f(x_{k-1})

∇f(xk)−∇f(xk−1)=Bk(xk−xk−1)\nabla f(x_k) - \nabla f(x_{k-1}) = B_k (x_k - x_{k-1})

习惯记法令 sk−1=xk−xk−1s_{k-1} = x_k - x_{k-1}，yk−1=∇f(xk)−∇f(xk−1)y_{k-1} = \nabla f(x_k) - \nabla f(x_{k-1})，故上式也可记为：

yk−1=Bksk−1y_{k-1} = B_k s_{k-1}

为了求出一个唯一的 pkp_k 让 mk(pk)m_k(p_k) 最小，我们肯定要保证 BkB_k 是可逆的。记 Hk=Bk−1H_k = B_k^{-1}，所以有

sk−1=Hkyk−1s_{k-1} = H_k y_{k-1}

这是最通用的一种表示方法，也是 BFGS 使用的表示方法。值得注意的是，这里的 HkH_k 和 Hessian 矩阵一点关系都没有！

最小扰动

然而，yk−1=Bksk−1y_{k-1} = B_k s_{k-1} 只表现了 BkB_k 需要满足什么约束，却无法求出 BkB_k 是什么。因为 n×nn \times n 对称矩阵包含 n(n+1)2\frac{n(n+1)}{2} 个独立自由度，而割线方程仅提供 nn 个独立等式约束，所以解空间是不定的。

于是我们再提出一个假设。我们希望算法在迭代过程中数值保持平稳变化，所以我们希望曲率信息每次不能改变太多，这就是最小扰动假设。

说是最小扰动，可我们要最小化什么扰动呢？有人认为，我们可以让 Bk=Bk−1+σvv⊤B_k = B_{k-1} + \sigma v v^\top，也就是加上一个对称的秩一矩阵。这样未知量就变成了 σ\sigma 和 vv，满足割线方程的解是唯一的。这种方法叫 SR1 法。然而，SR1 法无法保证正定性，还可能出现分母趋近于零的数值崩溃。

后来又有一种叫 DFP 的算法，通过最小化 BkB_k 的 Frobenius 范数来求出 BkB_k。不过我们发现，最后我们是要求方向 pkp_k 的，而 pk=−Bk−1∇f(xk)=−Hk∇f(xk)p_k = - B_k^{-1} \nabla f(x_k) = -H_k \nabla f(x_k)，所以为求平稳我们最小化 ∥H−Hk∥W2\Vert H - H_k \Vert_{W}^2 就好。

BFGS

这就是 BFGS 算法，即

min⁡Hk∥Hk−Hk−1∥W2s.t.Hk=Hk⊤,Hkyk−1=sk−1\min_{H_k} \Vert H_k - H_{k-1} \Vert_{W}^2 \quad \text{s.t.} \quad H_k = H_k^\top, \quad H_k y_{k-1} = s_{k-1}

有趣的是，DFP 算法要求 min⁡Bk∥Bk−Bk−1∥W2s.t.Bk=Bk⊤,Bksk−1=yk−1\min_{B_k} \Vert B_k - B_{k-1} \Vert_{W}^2 \quad \text{s.t.} \quad B_k = B_k^\top, \quad B_k s_{k-1} = y_{k-1}，与 BFGS 完全对偶。

我们跳过复杂的证明过程（见附录），直接给出结论：

Hk=Hk−1−sk−1yk−1⊤Hk−1+Hk−1yk−1sk−1⊤yk−1⊤sk−1+sk−1sk−1⊤yk−1⊤sk−1+yk−1⊤Hk−1yk−1(yk−1⊤sk−1)2sk−1sk−1⊤H_k= H_{k-1} - \frac{s_{k-1} y_{k-1}^\top H_{k-1} + H_{k-1} y_{k-1} s_{k-1}^\top}{y_{k-1}^\top s_{k-1}} + \frac{s_{k-1} s_{k-1}^\top}{y_{k-1}^\top s_{k-1}} + \frac{y_{k-1}^\top H_{k-1} y_{k-1}}{(y_{k-1}^\top s_{k-1})^2} s_{k-1} s_{k-1}^\top

令 ρk−1=1yk−1⊤sk−1\rho_{k-1} = \frac{1}{y_{k-1}^\top s_{k-1}}，将其重组因式分解，即得

Hk=(I−ρk−1sk−1yk−1⊤)Hk−1(I−ρk−1yk−1sk−1⊤)+ρk−1sk−1sk−1⊤H_k = (I - \rho_{k-1} s_{k-1} y_{k-1}^\top) H_{k-1} (I - \rho_{k-1} y_{k-1} s_{k-1}^\top) + \rho_{k-1} s_{k-1} s_{k-1}^\top

Hk≻0H_k \succ 0 的充要条件是 yk−1⊤sk−1>0y_{k-1}^\top s_{k-1} > 0（证明见附录），而这个条件可以由 Wolfe 线搜索保证（证明见附录），所以 BFGS 一般配合 Wolfe 线搜索使用。

收敛速度

根据 Dennis-Moré 准则（内容与证明见附录），一个拟牛顿法在极小值点邻域产生超线性收敛的充分必要条件是搜索方向渐近地逼近牛顿方向，即：

lim⁡k→∞∥(Bk−∇2f(x⋆))sk∥∥sk∥=0\lim_{k \to \infty} \frac{\Vert (B_k - \nabla^2 f(x^\star)) s_k \Vert}{\Vert s_k \Vert} = 0

在严格凸二次展开且步长通过 Wolfe 条件进行线搜索的条件下，BFGS 更新公式构造的对称正定矩阵序列 {Bk}\{B_k\} 严格满足 Dennis-Moré 准则，因此其收敛速度满足：

lim⁡k→∞∥xk+1−x⋆∥∥xk−x⋆∥=0\lim_{k \to \infty} \frac{\Vert x_{k+1} - x^\star \Vert}{\Vert x_k - x^\star \Vert} = 0

因此，BFGS 是超线性收敛的。

算法复杂度分析

时间复杂度

或许你会说，这个地方要用很大的一个矩阵乘法从 Hk−1H_{k-1} 更新到 HkH_k，不还是 O(n3)O(n^3) 的复杂度吗？但这里的低秩结构帮了大忙。回看那个因式分解前的式子，你会发现它只有一个真正的矩阵 Hk−1H_{k-1}。又考虑到矩阵的乘法结合律 ABC=(AB)C=A(BC)ABC = (AB)C = A(BC)，这样的结构当然可以拆成多个矩阵乘向量的形式。也就是说，它实质上是 O(n2)O(n^2) 的！

空间复杂度

要记录的最大的东西只有一个 n×nn \times n 的矩阵 HH，所以空间复杂度是 O(n2)O(n^2)。

L-BFGS

然而，有些时候 nn 实在是太大了。当它到达千万或上亿级别的时候，即使 O(n2)O(n^2) 的时间复杂度与空间复杂度也是不可接受的。

让我们来回顾一下因式分解后的 BFGS 逆 Hessian 更新公式：

Hk=(I−ρk−1sk−1yk−1⊤)Hk−1(I−ρk−1yk−1sk−1⊤)+ρk−1sk−1sk−1⊤H_k = (I - \rho_{k-1} s_{k-1} y_{k-1}^\top) H_{k-1} (I - \rho_{k-1} y_{k-1} s_{k-1}^\top) + \rho_{k-1} s_{k-1} s_{k-1}^\top

引入紧凑算子 Vk−1=I−ρk−1sk−1yk−1⊤V_{k-1} = I - \rho_{k-1} s_{k-1} y_{k-1}^\top 来简化运算，所以 Hk=Vk−1⊤Hk−1Vk−1+ρk−1sk−1sk−1⊤H_k = V_{k-1}^\top H_{k-1} V_{k-1} + \rho_{k-1} s_{k-1} s_{k-1}^\top。

我们可以向前继续追溯，把 Hk−1H_{k-1} 用 Hk−2H_{k-2} 表示，有

Hk=Vk−1⊤(Vk−2⊤Hk−2Vk−2+ρk−2sk−2sk−2⊤)Vk−1+ρk−1sk−1sk−1⊤=Vk−1⊤Vk−2⊤Hk−2Vk−2Vk−1+ρk−2Vk−1⊤sk−2sk−2⊤Vk−1+ρk−1sk−1sk−1⊤\begin{aligned} H_k &= V_{k-1}^\top \left( V_{k-2}^\top H_{k-2} V_{k-2} + \rho_{k-2} s_{k-2} s_{k-2}^\top \right) V_{k-1} + \rho_{k-1} s_{k-1} s_{k-1}^\top \\ &= V_{k-1}^\top V_{k-2}^\top H_{k-2} V_{k-2} V_{k-1} + \rho_{k-2} V_{k-1}^\top s_{k-2} s_{k-2}^\top V_{k-1} + \rho_{k-1} s_{k-1} s_{k-1}^\top \\ \end{aligned}

如果一直往前展开，最终会追溯到初始矩阵 H0H_0。这意味着 HkH_k 中其实“蕴含”着之前所有步的信息。但随着迭代步数增加，历史越久远的曲率信息对当前的参考价值越低，反而占用巨量内存。

L-BFGS 的思想是：只保留最近 mm 步的历史向量对（通常实践中取 m∈[3,20]m \in [3, 20]）。在第 kk 步时，仅利用最近的 mm 个历史位移对 si,yi,i=k−1,k−2,⋯ ,k−m{s_i, y_i},\quad i = k-1, k-2, \cdots, k-m。当回退到 k−mk-m 时，不再回溯。展开式为：

Hk=(Vk−1⊤⋯Vk−m⊤)Hkinit(Vk−m⋯Vk−1)+ρk−m(Vk−1⊤⋯Vk−m+1⊤)sk−msk−m⊤(Vk−m+1⋯Vk−1)+⋯+ρk−1sk−1sk−1⊤\begin{aligned} H_k = &\left( V_{k-1}^\top \cdots V_{k-m}^\top \right) H_k^\text{init} \left( V_{k-m} \cdots V_{k-1} \right) \\ &+ \rho_{k-m} \left( V_{k-1}^\top \cdots V_{k-m+1}^\top \right) s_{k-m} s_{k-m}^\top \left( V_{k-m+1} \cdots V_{k-1} \right) \\ &+ \cdots \\ &+ \rho_{k-1} s_{k-1} s_{k-1}^\top \end{aligned}

复杂度优化

逆序循环

直接计算 HkH_k 仍然需要 O(n2)O(n^2) 的时间与空间，但是我们的目标是 pk=Hk∇f(xk)p_k = H_k \nabla f(x_k)。

观察 ViV_i 作用于任意向量 qq 的效果：

Viq=(I−ρiyisi⊤)q=q−ρiyi(si⊤q)=q−(ρisi⊤q)yiV_i q = (I - \rho_i y_i s_i^\top) q = q - \rho_i y_i (s_i^\top q) = q - (\rho_i s_i^\top q) y_i

这意味着 ViV_i 的作用等效于在 qq 上减去一个沿着 yiy_i 方向的分量，而这一操作是 O(n)O(n) 的。因此，Vk−m⋯Vk−1∇f(xk)V_{k-m} \cdots V_{k-1} \nabla f(x_k) 在时间上是 O(mn)O(mn) 的。而且，因为我们只需要存储最新 mm 步的 si,yi{s_i, y_i}，所以在空间上这也是 O(mn)O(mn) 的。

基底变换

在经过一系列变换后，一个向量终于来到了 HkinitH_k^\text{init} 前。那是不是总要维护一个 O(n2)O(n^2) 的矩阵呢？

不是这样的。L-BFGS 都已经为空间舍弃了这么多信息，也不在乎再多丢一点了。所以 HkinitH_k^\text{init} 通常直接取为当前局部尺度的标量矩阵：

Hkinit=γkI,γk=sk−1⊤yk−1yk−1⊤yk−1H_k^\text{init} = \gamma_k I, \quad \gamma_k = \frac{s_{k-1}^\top y_{k-1}}{y_{k-1}^\top y_{k-1}}

所以这一步还是 O(n)O(n) 的。

正序循环

类似的，Vi⊤V_i^\top 作用于任意向量 qq 也是 O(n)O(n) 的。所以，这里也是一个总计 O(mn)O(mn) 的操作。

综上所述，L-BFGS 在时间和空间上都是 O(mn)O(mn) 的。

收敛速度

L-BFGS 丢弃了 mm 步之前的曲率信息。因此，近似 Hessian 矩阵无法随着 k→∞k \to \infty 充分积累全空间的二阶几何特征，此时 Dennis-Moré 条件不再成立。所以理论上其渐近收敛速度退化为线性收敛，其收缩系数由保存的步数 mm 和函数的条件数决定。

附录

BFGS 的 HkH_k 推导

BFGS 的目标是在所有满足对称性与割线约束

Hk=Hk⊤,Hkyk−1=sk−1H_k = H_k^\top, \qquad H_k y_{k-1} = s_{k-1}

的矩阵中，挑出一个与上一步的 Hk−1H_{k-1} 差别最小的。这里的“差别”用带权 Frobenius 范数度量：

∥A∥W2=tr⁡ ⁣(W1/2AWAW1/2)=∥W1/2AW1/2∥F2\Vert A \Vert_W^2 = \operatorname{tr}\!\left( W^{1/2} A W A W^{1/2} \right) = \Vert W^{1/2} A W^{1/2} \Vert_F^2

权重 WW 取任意满足 Wsk−1=yk−1W s_{k-1} = y_{k-1} 的对称正定矩阵（可以理解为割线意义下曲率的“替身”）。于是优化问题写成

min⁡Hk∥Hk−Hk−1∥W2s.t.Hk=Hk⊤,Hkyk−1=sk−1\min_{H_k} \Vert H_k - H_{k-1} \Vert_W^2 \quad \text{s.t.} \quad H_k = H_k^\top, \quad H_k y_{k-1} = s_{k-1}

坐标变换

记增量矩阵 E=Hk−Hk−1E = H_k - H_{k-1}，并作无量纲化变换

E^=W1/2EW1/2,H^k−1=W1/2Hk−1W1/2\hat E = W^{1/2} E W^{1/2}, \qquad \hat H_{k-1} = W^{1/2} H_{k-1} W^{1/2}

目标函数就化为普通的 Frobenius 范数平方 ∥E^∥F2=tr⁡(E^2)\Vert \hat E \Vert_F^2 = \operatorname{tr}(\hat E^2)。再处理约束：把 Eyk−1=sk−1−Hk−1yk−1E y_{k-1} = s_{k-1} - H_{k-1} y_{k-1} 两边左乘 W1/2W^{1/2}，并在左边插入 I=W−1/2W1/2I = W^{-1/2} W^{1/2}，得

W1/2EW1/2 W−1/2yk−1=W1/2sk−1−W1/2Hk−1yk−1W^{1/2} E W^{1/2} \, W^{-1/2} y_{k-1} = W^{1/2} s_{k-1} - W^{1/2} H_{k-1} y_{k-1}

由 Wsk−1=yk−1W s_{k-1} = y_{k-1} 可知 W1/2sk−1=W−1/2yk−1W^{1/2} s_{k-1} = W^{-1/2} y_{k-1}。令

v=W−1/2yk−1=W1/2sk−1,r=v−H^k−1vv = W^{-1/2} y_{k-1} = W^{1/2} s_{k-1}, \qquad r = v - \hat H_{k-1} v

则约束变为 E^v=r\hat E v = r，原问题简化为对称矩阵空间上的标准问题

min⁡E^=E^⊤12tr⁡ ⁣(E^2)s.t.E^v=r\min_{\hat E = \hat E^\top} \frac{1}{2} \operatorname{tr}\!\left( \hat E^2 \right) \quad \text{s.t.} \quad \hat E v = r

求解子问题

构造拉格朗日函数，其中 λ∈Rn\lambda \in \mathbb{R}^n 为乘子向量：

L(E^,λ)=12tr⁡ ⁣(E^2)−λ⊤(E^v−r)\mathcal{L}(\hat E, \lambda) = \frac{1}{2} \operatorname{tr}\!\left( \hat E^2 \right) - \lambda^\top (\hat E v - r)

利用 λ⊤E^v=tr⁡ ⁣(E^⋅vλ⊤+λv⊤2)\lambda^\top \hat E v = \operatorname{tr}\!\left( \hat E \cdot \frac{v \lambda^\top + \lambda v^\top}{2} \right)，对对称矩阵 E^\hat E 求梯度并令其为零：

∂L∂E^=E^−12(vλ⊤+λv⊤)=0⟹E^=12(vλ⊤+λv⊤)\frac{\partial \mathcal{L}}{\partial \hat E} = \hat E - \frac{1}{2}(v \lambda^\top + \lambda v^\top) = 0 \quad \Longrightarrow \quad \hat E = \frac{1}{2}(v \lambda^\top + \lambda v^\top)

代回约束：

12((λ⊤v)v+∥v∥2λ)=r\frac{1}{2}\left( (\lambda^\top v) v + \Vert v \Vert^2 \lambda \right) = r

两边左乘 v⊤v^\top，得 λ⊤v=v⊤r∥v∥2\lambda^\top v = \dfrac{v^\top r}{\Vert v \Vert^2}，再回代解出

λ=2∥v∥2r−v⊤r∥v∥4v\lambda = \frac{2}{\Vert v \Vert^2} r - \frac{v^\top r}{\Vert v \Vert^4} v

因此

E^=rv⊤+vr⊤∥v∥2−v⊤r∥v∥4vv⊤\hat E = \frac{r v^\top + v r^\top}{\Vert v \Vert^2} - \frac{v^\top r}{\Vert v \Vert^4} v v^\top

将 r=v−H^k−1vr = v - \hat H_{k-1} v 代入并整理，即

E^=vv⊤∥v∥2−H^k−1vv⊤+vv⊤H^k−1∥v∥2+v⊤H^k−1v∥v∥4vv⊤\hat E = \frac{v v^\top}{\Vert v \Vert^2} - \frac{\hat H_{k-1} v v^\top + v v^\top \hat H_{k-1}}{\Vert v \Vert^2} + \frac{v^\top \hat H_{k-1} v}{\Vert v \Vert^4} v v^\top

由 H^k=H^k−1+E^\hat H_k = \hat H_{k-1} + \hat E，得到变换空间中的更新式

H^k=H^k−1−H^k−1vv⊤+vv⊤H^k−1∥v∥2+vv⊤∥v∥2+v⊤H^k−1v∥v∥4vv⊤\hat H_k = \hat H_{k-1} - \frac{\hat H_{k-1} v v^\top + v v^\top \hat H_{k-1}}{\Vert v \Vert^2} + \frac{v v^\top}{\Vert v \Vert^2} + \frac{v^\top \hat H_{k-1} v}{\Vert v \Vert^4} v v^\top

逆变换还原

由 Hk=W−1/2H^kW−1/2H_k = W^{-1/2} \hat H_k W^{-1/2}，逐项还原：

∥v∥2=yk−1⊤sk−1,W−1/2v=sk−1,W−1/2H^k−1v=Hk−1yk−1,v⊤H^k−1v=yk−1⊤Hk−1yk−1\Vert v \Vert^2 = y_{k-1}^\top s_{k-1}, \qquad W^{-1/2} v = s_{k-1}, \qquad W^{-1/2} \hat H_{k-1} v = H_{k-1} y_{k-1}, \qquad v^\top \hat H_{k-1} v = y_{k-1}^\top H_{k-1} y_{k-1}

代回后即得

Hk=Hk−1−sk−1yk−1⊤Hk−1+Hk−1yk−1sk−1⊤yk−1⊤sk−1+sk−1sk−1⊤yk−1⊤sk−1+yk−1⊤Hk−1yk−1(yk−1⊤sk−1)2sk−1sk−1⊤H_k = H_{k-1} - \frac{s_{k-1} y_{k-1}^\top H_{k-1} + H_{k-1} y_{k-1} s_{k-1}^\top}{y_{k-1}^\top s_{k-1}} + \frac{s_{k-1} s_{k-1}^\top}{y_{k-1}^\top s_{k-1}} + \frac{y_{k-1}^\top H_{k-1} y_{k-1}}{(y_{k-1}^\top s_{k-1})^2} s_{k-1} s_{k-1}^\top

整个过程里权重 WW 被完全消去，说明结果与 WW 的选取无关。代入 ρk−1=1yk−1⊤sk−1\rho_{k-1} = \frac{1}{y_{k-1}^\top s_{k-1}} 并重组因式分解，就得到正文中使用的紧凑形式

Hk=(I−ρk−1sk−1yk−1⊤)Hk−1(I−ρk−1yk−1sk−1⊤)+ρk−1sk−1sk−1⊤H_k = (I - \rho_{k-1} s_{k-1} y_{k-1}^\top) H_{k-1} (I - \rho_{k-1} y_{k-1} s_{k-1}^\top) + \rho_{k-1} s_{k-1} s_{k-1}^\top

Hk≻0H_k \succ 0 的充要条件

设 Hk−1≻0H_{k-1} \succ 0，则 Hk≻0H_k \succ 0 当且仅当 yk−1⊤sk−1>0y_{k-1}^\top s_{k-1} > 0。

割线方程恒成立

先验证 Hkyk−1=sk−1H_k y_{k-1} = s_{k-1}。由 ρk−1yk−1sk−1⊤yk−1=ρk−1(yk−1⊤sk−1)yk−1=yk−1\rho_{k-1} y_{k-1} s_{k-1}^\top y_{k-1} = \rho_{k-1} (y_{k-1}^\top s_{k-1}) y_{k-1} = y_{k-1} 可知

(I−ρk−1yk−1sk−1⊤)yk−1=0(I - \rho_{k-1} y_{k-1} s_{k-1}^\top) y_{k-1} = 0

将这一结果代入因式分解式：

Hkyk−1=(I−ρk−1sk−1yk−1⊤)Hk−1⋅0+ρk−1sk−1(sk−1⊤yk−1)=sk−1H_k y_{k-1} = (I - \rho_{k-1} s_{k-1} y_{k-1}^\top) H_{k-1} \cdot 0 + \rho_{k-1} s_{k-1} (s_{k-1}^\top y_{k-1}) = s_{k-1}

可见因式分解式自动满足割线方程。

必要性

若 Hk≻0H_k \succ 0，取 z=yk−1z = y_{k-1}（由 Hkyk−1=sk−1H_k y_{k-1} = s_{k-1} 且 HkH_k 可逆可知 yk−1≠0y_{k-1} \neq 0），有

yk−1⊤Hkyk−1=yk−1⊤sk−1y_{k-1}^\top H_k y_{k-1} = y_{k-1}^\top s_{k-1}

左边是正定矩阵的二次型，恒为正，于是 yk−1⊤sk−1>0y_{k-1}^\top s_{k-1} > 0。

充分性

设 yk−1⊤sk−1>0y_{k-1}^\top s_{k-1} > 0，则 ρk−1>0\rho_{k-1} > 0。对任意非零向量 zz，由因式分解式

z⊤Hkz=w⊤Hk−1w+ρk−1(sk−1⊤z)2,w=(I−ρk−1yk−1sk−1⊤)zz^\top H_k z = w^\top H_{k-1} w + \rho_{k-1} (s_{k-1}^\top z)^2, \qquad w = (I - \rho_{k-1} y_{k-1} s_{k-1}^\top) z

因为 Hk−1≻0H_{k-1} \succ 0 且 ρk−1>0\rho_{k-1} > 0，右边两项均非负。若 z⊤Hkz=0z^\top H_k z = 0，则两项必须同时为零：

  1. ρk−1(sk−1⊤z)2=0⟹sk−1⊤z=0\rho_{k-1} (s_{k-1}^\top z)^2 = 0 \Longrightarrow s_{k-1}^\top z = 0；
  2. w⊤Hk−1w=0⟹w=0w^\top H_{k-1} w = 0 \Longrightarrow w = 0（因 Hk−1H_{k-1} 严格正定）。

把条件 1 代回 ww 的表达式，得 w=z−ρk−1yk−1(sk−1⊤z)=zw = z - \rho_{k-1} y_{k-1}(s_{k-1}^\top z) = z；再结合条件 2 的 w=0w = 0，推出 z=0z = 0，与假设矛盾。因此对任意 z≠0z \neq 0 都有 z⊤Hkz>0z^\top H_k z > 0，即 Hk≻0H_k \succ 0。

Wolfe 条件保证曲率条件

现在说明为什么 yk−1⊤sk−1>0y_{k-1}^\top s_{k-1} > 0 可以由 Wolfe 线搜索保证。第 k−1k-1 步的方向取为 pk−1=−Hk−1∇f(xk−1)p_{k-1} = -H_{k-1} \nabla f(x_{k-1})。因为 Hk−1≻0H_{k-1} \succ 0，当 ∇f(xk−1)≠0\nabla f(x_{k-1}) \neq 0 时

∇f(xk−1)⊤pk−1=−∇f(xk−1)⊤Hk−1∇f(xk−1)<0\nabla f(x_{k-1})^\top p_{k-1} = -\nabla f(x_{k-1})^\top H_{k-1} \nabla f(x_{k-1}) < 0

即 pk−1p_{k-1} 确实是下降方向。此时 xk=xk−1+αk−1pk−1x_k = x_{k-1} + \alpha_{k-1} p_{k-1}，故

sk−1=xk−xk−1=αk−1pk−1s_{k-1} = x_k - x_{k-1} = \alpha_{k-1} p_{k-1}

Wolfe 条件中的曲率准则为 ϕ′(αk−1)≥c2ϕ′(0)\phi'(\alpha_{k-1}) \ge c_2 \phi'(0)，其中 0<c2<10 < c_2 < 1，而 ϕ′(α)=∇f(xk−1+αpk−1)⊤pk−1\phi'(\alpha) = \nabla f(x_{k-1} + \alpha p_{k-1})^\top p_{k-1}，所以

∇f(xk)⊤pk−1≥c2∇f(xk−1)⊤pk−1\nabla f(x_k)^\top p_{k-1} \ge c_2 \nabla f(x_{k-1})^\top p_{k-1}

两边减去 ∇f(xk−1)⊤pk−1\nabla f(x_{k-1})^\top p_{k-1}：

yk−1⊤pk−1≥(c2−1)∇f(xk−1)⊤pk−1y_{k-1}^\top p_{k-1} \ge (c_2 - 1) \nabla f(x_{k-1})^\top p_{k-1}

再乘以正步长 αk−1>0\alpha_{k-1} > 0，并利用 sk−1=αk−1pk−1s_{k-1} = \alpha_{k-1} p_{k-1}：

yk−1⊤sk−1≥(c2−1)αk−1∇f(xk−1)⊤pk−1y_{k-1}^\top s_{k-1} \ge (c_2 - 1) \alpha_{k-1} \nabla f(x_{k-1})^\top p_{k-1}

由于 0<c2<10 < c_2 < 1 有 c2−1<0c_2 - 1 < 0，而 ∇f(xk−1)⊤pk−1<0\nabla f(x_{k-1})^\top p_{k-1} < 0、αk−1>0\alpha_{k-1} > 0，两个负因子相乘得正，因此

yk−1⊤sk−1>0y_{k-1}^\top s_{k-1} > 0

这正好是 Hk≻0H_k \succ 0 所需要的条件，所以 BFGS 一般配合 Wolfe 线搜索使用。强 Wolfe 条件对导数绝对值额外施加了 ∣∇f(xk)⊤pk−1∣≤c2∣∇f(xk−1)⊤pk−1∣\vert \nabla f(x_k)^\top p_{k-1} \vert \le c_2 \vert \nabla f(x_{k-1})^\top p_{k-1} \vert 的约束，其曲率准则与标准 Wolfe 一致，同样保证 yk−1⊤sk−1>0y_{k-1}^\top s_{k-1} > 0。

Dennis–Moré 准则

内容

设 ff 二次连续可微，x⋆x^\star 满足 ∇f(x⋆)=0\nabla f(x^\star) = 0，且 ∇2f(x⋆)\nabla^2 f(x^\star) 对称正定（为书写方便，下文记 G=∇2f(x⋆)G = \nabla^2 f(x^\star)）。沿用正文的记号，第 kk 步的拟牛顿方向为

pk=−Bk−1∇f(xk)=−Hk∇f(xk),p_k = - B_k^{-1} \nabla f(x_k) = - H_k \nabla f(x_k),

迭代点为 xk+1=xk+αkpkx_{k+1} = x_k + \alpha_k p_k，故位移

sk=xk+1−xk=αkpks_k = x_{k+1} - x_k = \alpha_k p_k

（这与正文的 sk−1=xk−xk−1=αk−1pk−1s_{k-1} = x_k - x_{k-1} = \alpha_{k-1} p_{k-1} 是同一套下标约定）。在考察极小值点邻域的超线性收敛时，可设线搜索最终接受单位步长 αk=1\alpha_k = 1，于是

sk=−Bk−1∇f(xk),等价地Bksk=−∇f(xk).(DM-0)s_k = - B_k^{-1} \nabla f(x_k), \qquad\text{等价地}\qquad B_k s_k = - \nabla f(x_k). \tag{DM-0}

假设迭代序列 xk→x⋆x_k \to x^\star。Dennis–Moré 准则指出：

lim⁡k→∞∥xk+1−x⋆∥∥xk−x⋆∥=0⟺lim⁡k→∞∥(Bk−∇2f(x⋆))sk∥∥sk∥=0.(DM)\lim_{k \to \infty} \frac{\Vert x_{k+1} - x^\star \Vert}{\Vert x_k - x^\star \Vert} = 0 \quad\Longleftrightarrow\quad \lim_{k \to \infty} \frac{\Vert (B_k - \nabla^2 f(x^\star)) s_k \Vert}{\Vert s_k \Vert} = 0. \tag{DM}

也就是说，拟牛顿法在极小值点邻域超线性收敛的充分必要条件是 BkB_k 沿迭代方向 sks_k 渐近地趋近于真实 Hessian GG。

直观解释

条件 (DM) 只要求矩阵 BkB_k 在搜索方向 sks_k 上等于 GG 的越来越好的近似，而不要求在整个 Rn\mathbb{R}^n 上接近 GG。由于 sk=−Bk−1∇f(xk)s_k = - B_k^{-1} \nabla f(x_k)，若 Bksk≈GskB_k s_k \approx G s_k，则 pk=−Bk−1∇f(xk)≈−G−1∇f(xk)p_k = -B_k^{-1}\nabla f(x_k) \approx -G^{-1}\nabla f(x_k)，即搜索方向渐近地逼近正文中的牛顿方向

pkNewton=−(∇2f(xk))−1∇f(xk).p_k^{\text{Newton}} = - \left( \nabla^2 f(x_k) \right)^{-1} \nabla f(x_k).

由于 xk→x⋆x_k \to x^\star 且 Hessian 连续，∇2f(xk)→∇2f(x⋆)=G\nabla^2 f(x_k) \to \nabla^2 f(x^\star) = G，因此该方向在极限意义下就是 −G−1∇f(xk)-G^{-1}\nabla f(x_k)，这正是正文“搜索方向渐近地逼近牛顿方向”的含义。

证明

记位移误差 ek=xk−x⋆e_k = x_k - x^\star（以区别于正文中表示函数值误差的 Δk\Delta_k）。因为 ∇f\nabla f 在 x⋆x^\star 处可微、∇f(x⋆)=0\nabla f(x^\star) = 0、∇2f(x⋆)=G\nabla^2 f(x^\star) = G，由导数的定义（Fréchet 导数）有

∇f(xk)=Gek+rk,∥rk∥=o(∥ek∥).(1)\nabla f(x_k) = G e_k + r_k, \qquad \Vert r_k \Vert = o(\Vert e_k \Vert). \tag{1}

又由拟牛顿步关系式 (DM-0) 可得

Gsk=Bksk−(Bk−G)sk=−∇f(xk)−(Bk−G)sk.(2)G s_k = B_k s_k - (B_k - G) s_k = - \nabla f(x_k) - (B_k - G) s_k. \tag{2}

充分性。 设 (DM) 成立，即 (Bk−G)sk=o(∥sk∥)(B_k - G) s_k = o(\Vert s_k \Vert)。把 (1) 代入 (2)：

Gsk=−Gek−rk−(Bk−G)sk.G s_k = - G e_k - r_k - (B_k - G) s_k.

两边左乘 G−1G^{-1} 并取范数，利用 ∥rk∥=o(∥ek∥)\Vert r_k \Vert = o(\Vert e_k \Vert)：

∥sk∥≤∥ek∥+∥G−1∥∥rk∥+∥G−1∥ o(∥sk∥)=O(∥ek∥)+o(∥sk∥),\Vert s_k \Vert \le \Vert e_k \Vert + \Vert G^{-1} \Vert \Vert r_k \Vert + \Vert G^{-1} \Vert \, o(\Vert s_k \Vert) = O(\Vert e_k \Vert) + o(\Vert s_k \Vert),

由此推出 ∥sk∥=O(∥ek∥)\Vert s_k \Vert = O(\Vert e_k \Vert)，从而 o(∥sk∥)=o(∥ek∥)o(\Vert s_k \Vert) = o(\Vert e_k \Vert)。将其代回上式：

sk=−ek+G−1o(∥ek∥)=−ek+o(∥ek∥).s_k = - e_k + G^{-1} o(\Vert e_k \Vert) = - e_k + o(\Vert e_k \Vert).

于是

ek+1=ek+sk=o(∥ek∥),即∥xk+1−x⋆∥∥xk−x⋆∥→0,e_{k+1} = e_k + s_k = o(\Vert e_k \Vert), \qquad\text{即}\qquad \frac{\Vert x_{k+1} - x^\star \Vert}{\Vert x_k - x^\star \Vert} \to 0,

所以迭代超线性收敛。

必要性。 设迭代超线性收敛，即 ek+1=o(∥ek∥)e_{k+1} = o(\Vert e_k \Vert)。由 sk=ek+1−eks_k = e_{k+1} - e_k 得

∥sk∥≥∥ek∥−∥ek+1∥=(1−o(1))∥ek∥,\Vert s_k \Vert \ge \Vert e_k \Vert - \Vert e_{k+1} \Vert = \big( 1 - o(1) \big) \Vert e_k \Vert,

所以 ∥ek∥≤(1+o(1))∥sk∥\Vert e_k \Vert \le \big( 1 + o(1) \big) \Vert s_k \Vert。另一方面，由 (1) 与 (2)：

(Bk−G)sk=−∇f(xk)−Gsk=−Gek−rk−Gsk=−G(ek+sk)−rk=−Gek+1−rk.(B_k - G) s_k = - \nabla f(x_k) - G s_k = - G e_k - r_k - G s_k = - G (e_k + s_k) - r_k = - G e_{k+1} - r_k.

取范数并利用 ek+1=o(∥ek∥)e_{k+1} = o(\Vert e_k \Vert)、∥rk∥=o(∥ek∥)\Vert r_k \Vert = o(\Vert e_k \Vert)：

∥(Bk−G)sk∥≤∥G∥∥ek+1∥+∥rk∥=o(∥ek∥)≤o(∥sk∥),\Vert (B_k - G) s_k \Vert \le \Vert G \Vert \Vert e_{k+1} \Vert + \Vert r_k \Vert = o(\Vert e_k \Vert) \le o(\Vert s_k \Vert),

这正是 (DM)。□\square

几点说明

  • 该定理把“超线性收敛”这一关于迭代点序列的性质，完全转化为关于矩阵 BkB_k 的一个可检验条件：只需 BkB_k 在迭代方向 sks_k 上趋近于 ∇2f(x⋆)\nabla^2 f(x^\star)，而无需在整个空间上趋近。BFGS、DFP、SR1 等拟牛顿法的经典超线性收敛定理，正是在适当线搜索下证明其 {Bk}\{B_k\} 满足 (DM)。
  • 对 L-BFGS 而言，由于只保留最近 mm 步的曲率信息，BkB_k（或 HkH_k）无法在任意方向上恢复全空间二阶信息，(DM) 一般不成立，因此渐近收敛速度退化为线性收敛，收缩率依赖于保存步数 mm 与函数条件数。
  • 与拟牛顿法的联系：BFGS 维护的是逆 Hessian 近似 HkH_k，而定理中的 Bk=Hk−1B_k = H_k^{-1}；由于 (DM) 用 BkB_k 表述，实际验证时也可等价地考察 Hk−1H_k^{-1} 在方向 sks_k 上是否收敛到 ∇2f(x⋆)\nabla^2 f(x^\star)。