第 08 章 拟牛顿法

学习目标

读完本章，你应当能够：

1. 从"让模型在最近两个迭代点上梯度都对"出发推导割线方程 \(B_{k+1}s_k=y_k\)，并说明为什么必须配合 Wolfe 线搜索才能保证曲率条件 \(s_k^Ty_k>0\)。
2. 写出 BFGS 的逆 Hessian 更新公式和 Hessian 更新公式，证明它保持正定，并能按 Algorithm 8.1 自己实现一个可用的 BFGS（含初始矩阵缩放）。
3. 推导 SR1 更新，说明它何时会失效、跳过规则为什么合理，以及它为什么适合放在信赖域框架里。
4. 理解 Broyden 族把 BFGS、DFP、SR1 串成一条线，以及二次函数上"\(n\) 步终止"和"与共轭梯度等价"等性质。
5. 读懂 BFGS 全局收敛和超线性收敛证明的主线（势函数 \(\psi(B)=\operatorname{tr}B-\ln\det B\) 与 Dennis–Moré 条件），并据此判断：BFGS 末步的 \(H_k\) 不能直接当作 MLE 的协方差矩阵。

读前导读

这一章在解决什么问题

牛顿法很快，但每一步都要算 Hessian（全部二阶导数），这在 GARCH 似然、Heston 校准这类问题里又贵又麻烦。拟牛顿法的办法是：不算 Hessian，而是边走边"猜"Hessian。每走一步，你都能看到梯度变了多少；梯度的变化除以位置的变化，就是这个方向上的曲率。把这些观测一点点累积进一个矩阵 \(B_k\)，它就越来越像真正的 Hessian。

用你熟悉的语言说：牛顿法相当于每次都用解析公式算 Gamma；拟牛顿法相当于只有 Delta 的报价，用"两次 Delta 之差 ÷ 标的价格之差"去估 Gamma，并且不断用新的观测修正这个估计。BFGS 是这类方法中最好用的一个，scipy.optimize.minimize 不给 Hessian 时默认用的就是它，你做 MLE 时几乎肯定用过。

本章还有一个对实务很要紧的结论：BFGS 结束时给出的逆 Hessian 近似 hess_inv 只在"走过的方向"上准确，不能直接拿来算 MLE 的标准误。8.8.3 节用 GARCH 的例子说明了偏差可以有多大。

需要先想起来的数学

1. 矩阵与向量的乘法、转置、外积。 \(s^Ty\) 是内积，结果是一个数（对应分量相乘再求和）；\(sy^T\) 是外积，结果是 \(n\times n\) 矩阵，第 \((i,j)\) 元是 \(s_iy_j\)，它的秩为 1（所有列都是 \(s\) 的倍数）。例：\(s=(1,2)^T\)、\(y=(3,4)^T\)，\(s^Ty=11\)，\(sy^T=\begin{bmatrix}3&4\\6&8\end{bmatrix}\)。本章的更新公式全是"旧矩阵 + 一两个外积"。见 第 00 册第 06 章 线性代数速成。

2. 正定矩阵。 对称矩阵 \(B\) 正定，指对任何非零向量 \(z\) 都有 \(z^TBz>0\)，等价于所有特征值为正。直观上，正定矩阵对应"向任何方向走都是往上弯的碗"。协方差矩阵就是半正定的：\(w^T\Sigma w\) 是组合方差，不可能为负。记号 \(A\succeq mI\) 表示 \(A-mI\) 半正定，即 \(A\) 的所有特征值都不小于 \(m\)。见 第 00 册第 06 章。

3. 梯度、Hessian 与二次模型。 二次函数 \(m(p)=f+g^Tp+\tfrac12p^TBp\) 的梯度是 \(g+Bp\)，令其为零得极小点 \(p=-B^{-1}g\)（\(B\) 正定时）。这就是"牛顿步"。一元类比：\(f+gp+\tfrac12bp^2\) 的极小点 \(p=-g/b\)。见 第 00 册第 05 章 多元微积分与优化。

4. 迹、行列式与特征值。 迹 \(\operatorname{tr}B\) 是对角元之和，等于全部特征值之和；行列式 \(\det B\) 等于全部特征值之积。所以 \(\ln\det B=\sum\ln\lambda_i\)。8.7 节的收敛证明就是靠迹盯住"最大特征值不能太大"、靠行列式盯住"最小特征值不能太小"。例：\(\operatorname{diag}(2,3)\) 的迹为 5、行列式为 6。见 第 00 册第 06 章。

5. 收敛速度的说法。 线性收敛：误差每步乘一个固定比例（如 0.5）；超线性：这个比例本身趋于 0；二次：新误差约为旧误差的平方（\(10^{-3}\to10^{-6}\to10^{-12}\)）。在第 02、03 章已经出现过，本章 8.3.5 的表格是很好的例子。

怎么读这一章

核心必读：8.2（割线方程与曲率条件）、8.3.1–8.3.4（BFGS 公式与正定性）、8.4（实现要点）、8.8（实战）。这几节读懂，你就能看懂优化器的报错、知道该怎么调。8.5 的 SR1 读到 8.5.2 即可掌握要点（何时用、为什么能跳过），8.5.3–8.5.4 可以只看结论。8.6 Broyden 族和 8.7 收敛性分析第一次可以只读每小节开头的结论和 8.7.2 末尾的"要点"，证明细节留到需要时再看。

建议顺序：8.1 → 8.2 → 8.3 → 8.4 → 8.8 → 8.5 → 8.7.2 的"要点" → 其余。

8.1 动机：只用梯度，却想要牛顿法的速度

先说结论：拟牛顿法（quasi-Newton methods）每步只需要梯度，却能达到超线性收敛，在不想或不能计算 Hessian 的场合，它是最常用的无约束优化算法。scipy.optimize.minimize 在无约束且没给 Hessian 时默认用的就是 BFGS；带简单上下界时常用的 L-BFGS-B 是它的大规模版本（第 09 章）。

回顾两类基本方法：

• 最速下降法只用梯度，每步便宜，但只有线性收敛，在病态问题上极慢；
• 牛顿法用 Hessian，局部二次收敛，但每步要算二阶导数并解一个 \(n\times n\) 线性方程组。

拟牛顿法的想法是：每走一步，梯度的变化 \(\nabla f_{k+1}-\nabla f_k\) 本身就包含了沿这一步方向的曲率信息。把这些信息一点一点累积进一个近似矩阵，就能逐步"学会"Hessian，而不必显式计算它。

这个想法有一段有趣的历史。1950 年代中期，Argonne 国家实验室的物理学家 W. C. Davidon 用坐标下降法做长时间的优化计算，当时计算机很不稳定，总在算完之前崩溃。为了加快迭代，他发明了第一个拟牛顿算法。Fletcher 与 Powell 很快证明它比当时已有的方法快且可靠得多，非线性优化领域"一夜之间"被改变，此后二十年出现了大量变体和数百篇论文。讽刺的是，Davidon 的论文当年没有被接受发表，作为技术报告存在了三十多年，直到 1991 年才刊登在 SIAM Journal on Optimization 的创刊号上。

自动微分（原书第 7 章）出现后，手工推导二阶导数的麻烦和出错风险大大降低，拟牛顿法的吸引力有所下降，但只是有限程度：在许多问题上它仍然很有竞争力，尤其是函数本身计算昂贵、Hessian 稠密或难以得到的时候。本章讨论中小规模问题，第 09 章讨论大规模推广。割线方程、SR1 与 BFGS 公式已在第 02 章 2.3.2 节初次出现，Wolfe 条件与 Dennis–Moré 条件见第 03 章；本章给出它们的推导与收敛理论，记号（\(B_k\) 近似 Hessian、\(H_k\) 近似逆 Hessian、\(s_k\)、\(y_k\)）与前几章一致。

8.2 割线方程与曲率条件

8.2.1 二次模型

在当前迭代点 \(x_k\) 处建立二次模型

其中 \(B_k\) 是 \(n\times n\) 对称正定矩阵，每步更新。模型在 \(p=0\) 处的函数值和梯度都与 \(f\) 一致。凸二次模型的极小点

作为搜索方向，迭代为

步长 \(\alpha_k\) 满足 Wolfe 条件（原书第 3 章 (3.6)）：

形式上这就是线搜索牛顿法，只不过用近似矩阵 \(B_k\) 代替了真 Hessian。

8.2.2 割线方程

Davidon 的关键想法是：不要每步从头计算 \(B_k\)，而是更新它，把最近一步测到的曲率加进去。走完一步后得到新模型

一个合理的要求是：\(m_{k+1}\) 的梯度在最近的两个迭代点 \(x_k\)、\(x_{k+1}\) 处都与 \(f\) 的梯度一致。在 \(x_{k+1}\)（即 \(p=0\)）处自动满足；在 \(x_k\)（即 \(p=-\alpha_kp_k\)）处要求

记

就得到割线方程（secant equation）

直观理解：一维时它就是割线法 \(f''\approx\dfrac{f'(x_{k+1})-f'(x_k)}{x_{k+1}-x_k}\)；高维时它要求新矩阵把"刚走的一步"映射成"刚观测到的梯度变化"。由 Taylor 定理，\(y_k=\bar G_ks_k\)，其中

是沿这一步的平均 Hessian。所以割线方程要求 \(B_{k+1}\) 沿 \(s_k\) 方向与平均 Hessian 一致。

推导拆解：从模型到割线方程。

1. 对 \(m_{k+1}(p)\) 求梯度：常数项 \(f_{k+1}\) 求导为零；\(\nabla f_{k+1}^Tp\) 的梯度是 \(\nabla f_{k+1}\)；\(\tfrac12p^TB_{k+1}p\) 的梯度是 \(B_{k+1}p\)（\(B\) 对称，类比一元的 \(\tfrac12bp^2\) 求导得 \(bp\)）。所以 \(\nabla m_{k+1}(p)=\nabla f_{k+1}+B_{k+1}p\)。
2. 旧点 \(x_k\) 相对于新点 \(x_{k+1}\) 的位移是 \(p=-\alpha_kp_k=-s_k\)。代入得 \(\nabla f_{k+1}-B_{k+1}s_k=\nabla f_k\)。
3. 移项：\(B_{k+1}s_k=\nabla f_{k+1}-\nabla f_k=y_k\)。 (8.11) 中的积分是"沿着线段把 Hessian 取平均"：梯度的变化等于一路上 Hessian 乘以位移的累积，和"价格变化 = 沿路径对 Delta 积分"是同一个道理。

金融直觉：一维时割线方程就是 \(b_{k+1}=\dfrac{\Delta_{\text{新}}-\Delta_{\text{旧}}}{S_{\text{新}}-S_{\text{旧}}}\)，也就是用两次观测到的 Delta 估 Gamma。高维时一次观测只能告诉你"沿刚走的那个方向"的曲率，其他方向一无所知。所以 \(B_{k+1}\) 只被割线方程确定了一部分，剩下的自由度要靠"别离旧矩阵太远"来定（8.3 节）。

8.2.3 曲率条件与 Wolfe 线搜索

对称正定的 \(B_{k+1}\) 要把 \(s_k\) 映射为 \(y_k\)，必须满足（在 (8.6) 两边左乘 \(s_k^T\) 即见）

这叫曲率条件（curvature condition）。\(f\) 强凸时它对任意两点都成立；非凸时不一定成立，需要线搜索来强制。Wolfe 条件的第二条给出 \(\nabla f_{k+1}^Ts_k\ge c_2\nabla f_k^Ts_k\)，于是

因为 \(c_2<1\) 且 \(p_k\) 是下降方向。这就是拟牛顿法必须用 Wolfe 线搜索、而不能只用 Armijo 回溯的根本原因：只检查充分下降的回溯法可能接受一个太短的步长，使得 \(s_k^Ty_k\le0\)，更新就无法保持正定。

推导拆解：(8.7) 和 (8.8) 各一步。

• (8.7)：在 \(B_{k+1}s_k=y_k\) 两边左乘 \(s_k^T\)，得 \(s_k^TB_{k+1}s_k=s_k^Ty_k\)。左边是正定矩阵的二次型，必须为正，所以右边也必须为正。
• (8.8)：\(y_k^Ts_k=\nabla f_{k+1}^Ts_k-\nabla f_k^Ts_k\)。Wolfe 第二条给出 \(\nabla f_{k+1}^Ts_k\ge c_2\nabla f_k^Ts_k\)（\(s_k=\alpha_kp_k\)，两边同乘 \(\alpha_k>0\)）。于是 \(y_k^Ts_k\ge(c_2-1)\nabla f_k^Ts_k=(c_2-1)\alpha_k\nabla f_k^Tp_k\)。\(c_2-1<0\)，下降方向又使 \(\nabla f_k^Tp_k<0\)，负负得正。 白话解释：\(s_k^Ty_k>0\) 的意思是"沿着刚走的方向，坡度在变缓"，也就是函数在这个方向上是往上弯的（正曲率）。Wolfe 第二条要求新点的斜率"没那么陡了"，正好保证了这一点。只检查"函数值降够了"的回溯法不管斜率，可能停在斜率依然很陡、甚至更陡的地方。

8.3 DFP 与 BFGS：最近矩阵问题

8.3.1 割线方程的解不唯一

曲率条件成立时，割线方程有无穷多个对称正定解：对称矩阵有 \(n(n+1)/2\) 个自由度，割线方程只提供 \(n\) 个条件；正定性给出 \(n\) 个不等式（所有顺序主子式为正），仍不足以确定剩余的自由度。于是自然地要求 \(B_{k+1}\) 在某种意义下最接近 \(B_k\)：

不同的范数给出不同的拟牛顿法。原书采用加权 Frobenius 范数

权矩阵 \(W\) 取任一满足 \(Wy_k=s_k\) 的矩阵，例如平均 Hessian 的逆 \(W=\bar G_k^{-1}\)。这个权使范数无量纲，因而所得方法具有尺度不变性：对变量做线性变换不改变算法的本质行为。

8.3.2 DFP 公式

(8.9) 在上述范数下的唯一解是

它由 Davidon 在 1959 年提出，Fletcher 与 Powell 研究、实现并推广，故称 DFP。计算方向时更方便的是逆近似 \(H_k=B_k^{-1}\)，由 Sherman–Morrison–Woodbury 公式得

后两项都是秩一矩阵，所以这是对 \(H_k\) 的秩二修正。这就是拟牛顿更新的基本思想：不从头重算，而是用一个简单修正把最新观测到的信息与已有近似结合起来。

白话解释：几个新名词。

• Frobenius 范数 \(\|C\|_F\)：把矩阵所有元素平方求和再开方，就是把矩阵当成一个长向量算长度。"离 \(B_k\) 最近"就是"改动的元素平方和最小"。加权版本先用 \(W^{1/2}\) 把坐标"拉成同一尺度"再算，避免单位不同的变量（比如 GARCH 里的 \(\omega\) 和 \(\beta\)）主导距离。
• 秩一矩阵 \(uv^T\)：所有列都是 \(u\) 的倍数，只往一个方向"加东西"。秩二修正就是两个这样的矩阵之和，改动非常"轻"，只碰两个方向。
• Sherman–Morrison–Woodbury 公式：给出"矩阵加低秩修正后的逆"的显式表达。最简单的版本是 \((A+uv^T)^{-1}=A^{-1}-\dfrac{A^{-1}uv^TA^{-1}}{1+v^TA^{-1}u}\)。它的意义是：\(B_k\) 做低秩修正，\(H_k=B_k^{-1}\) 也只需做低秩修正，不必重新求逆（\(O(n^3)\)）。 这和更新协方差矩阵的思路类似：新来一个观测，样本协方差加上一个外积项，而不是重新算全部数据。

8.3.3 BFGS 公式

DFP 很有效，但很快被 BFGS（Broyden、Fletcher、Goldfarb、Shanno 四人的姓氏首字母）取代，后者至今被公认为最有效的拟牛顿公式。推导方式是把同样的思路用在逆矩阵上：要求 \(H_{k+1}\) 对称正定且满足 \(H_{k+1}y_k=s_k\)，

加权 Frobenius 范数的权 \(W\) 满足 \(Ws_k=y_k\)（如 \(W=\bar G_k\)）。唯一解为

对比 (8.13) 与 (8.16)：互换 \(B\leftrightarrow H\)、\(s\leftrightarrow y\)，DFP 就变成 BFGS。两者互为对偶。

用 Sherman–Morrison–Woodbury 公式可得 BFGS 的 Hessian 形式

它的结构很好记：减去 \(B_k\) 沿 \(s_k\) 的"旧曲率"，加上观测到的"新曲率"。

推导拆解：验证 (8.16) 确实满足 \(H_{k+1}y_k=s_k\)。把 \(y_k\) 从右边乘进去：

1. 先看最右边的括号：\((I-\rho_ky_ks_k^T)y_k=y_k-\rho_ky_k(s_k^Ty_k)=y_k-y_k=0\)，因为 \(\rho_k(s_k^Ty_k)=1\)。
2. 所以第一大项整体为零，只剩 \(\rho_ks_ks_k^Ty_k=\rho_k(s_k^Ty_k)s_k=s_k\)。 同理验证 (8.19)：\(B_{k+1}s_k=B_ks_k-B_ks_k\cdot\dfrac{s_k^TB_ks_k}{s_k^TB_ks_k}+y_k\cdot\dfrac{y_k^Ts_k}{y_k^Ts_k}=y_k\)。前两项抵消，正好说明"减去旧曲率、加上新曲率"：沿 \(s_k\) 方向，旧的认识被完全替换成新观测。 还可以看出 \(\rho_k=1/(y_k^Ts_k)\) 出现在分母里，曲率条件一旦失败（\(y_k^Ts_k\le0\)），整个公式就失去意义。

Algorithm 8.1（BFGS 方法）

给定 x_0，容差 ε > 0，逆 Hessian 近似 H_0；k ← 0
while ‖∇f_k‖ > ε
    p_k = −H_k ∇f_k                                  (8.18)
    x_{k+1} = x_k + α_k p_k，α_k 由满足 Wolfe 条件的线搜索给出
    s_k = x_{k+1} − x_k，y_k = ∇f_{k+1} − ∇f_k
    按 (8.16) 计算 H_{k+1}
    k ← k+1
end

每步只有 \(O(n^2)\) 运算（外加函数与梯度求值），不需要解线性方程组或做矩阵–矩阵乘法这类 \(O(n^3)\) 运算。它稳健、超线性收敛，对多数实际用途已经足够快。牛顿法虽然是二次收敛，但每步代价更高；BFGS 更重要的优势是不需要二阶导数。

8.3.4 正定性自动保持

(8.15) 并没有显式要求正定，但只要 \(H_k\) 正定且 \(y_k^Ts_k>0\)，\(H_{k+1}\) 就正定。证明：对任意 \(z\ne0\)，

右端为零要求两项都为零。第二项为零意味着 \(s_k^Tz=0\)，此时 \(w=z\ne0\)，第一项严格为正。所以 \(H_{k+1}\) 正定。这个证明再次说明曲率条件的地位：\(\rho_k>0\) 是一切的前提。

推导拆解：等式 \(z^TH_{k+1}z=w^TH_kw+\rho_k(z^Ts_k)^2\) 从哪里来。

1. 把 (8.16) 夹在 \(z^T\) 和 \(z\) 之间：\(z^TH_{k+1}z=z^T(I-\rho_ks_ky_k^T)H_k(I-\rho_ky_ks_k^T)z+\rho_kz^Ts_ks_k^Tz\)。
2. 令 \(w=(I-\rho_ky_ks_k^T)z=z-\rho_ky_k(s_k^Tz)\)（\(s_k^Tz\) 是一个数，可以挪到后面）。转置一下，\(w^T=z^T(I-\rho_ks_ky_k^T)\)，正好是左边那个括号。于是第一项就是 \(w^TH_kw\)。
3. 第二项 \(z^Ts_ks_k^Tz=(s_k^Tz)^2\)。 第一项因 \(H_k\) 正定而 \(\ge0\)，第二项因 \(\rho_k>0\) 而 \(\ge0\)。这种"凑成平方和"的手法，与证明组合方差 \(w^T\Sigma w\ge0\) 时把它写成 \(\operatorname{Var}(w^TR)\) 是一个思路。

8.3.5 数值对比：Rosenbrock 函数

原书在 Rosenbrock 函数 \(f(x)=100(x_2-x_1^2)^2+(1-x_1)^2\) 上，从 \((-1.2,1)\) 出发，三种方法都用 Wolfe 步长，把梯度范数降到 \(10^{-5}\)：最速下降 5264 次迭代，BFGS 34 次，非精确牛顿 21 次。最后几步的误差 \(\|x_k-x^*\|\)：

最速下降的误差比值接近 1（线性且极慢），BFGS 的比值趋于 0（超线性），牛顿法每步误差指数翻倍（二次）。8.8 节会用自己实现的代码复现这个对比。

8.3.6 自校正性质

如果某一步 \(y_k^Ts_k\) 很小但为正，由 (8.16) 知 \(H_{k+1}\) 会变得很大。一个坏的近似会不会一直坏下去？分析和大量实验给出的答案是：BFGS 有很强的自校正性质（self-correcting properties）——即使 \(H_k\) 错误估计了曲率并拖慢了迭代，几步之内近似就会自行修正。DFP 的自校正能力弱得多，这被认为是它实际表现较差的原因。8.7 节的势函数证明会给出这一性质的数学形式。前提仍然是合格的线搜索：Wolfe 条件保证梯度在"能让模型捕获恰当曲率"的点上被采样。

尺度不变性方面，权矩阵 \(W\) 的选择保证 (8.9)、(8.15) 的目标对变量线性变换不变。换别的 \(W\) 会得到别的公式，但人们搜寻了很多年，没有找到显著优于 BFGS 的公式。

8.4 BFGS 的实现要点

线搜索。 应满足 Wolfe 或强 Wolfe 条件，并且总是先试 \(\alpha_k=1\)：在一定条件下单位步长最终总会被接受，这正是超线性收敛的来源。经验表明较不精确的线搜索更省函数求值，常用 \(c_1=10^{-4}\)、\(c_2=0.9\)。

初始矩阵。 \(H_0=\beta I\) 没有通用的好选法。\(\beta\) 太大，首步 \(p_0=-\beta g_0\) 过长，线搜索要很多次函数求值。有的软件让用户给定首步长度 \(\delta\)，取 \(H_0=\delta\|g_0\|^{-1}I\)。原书推荐一个很有效的启发式：先用 \(H_0=I\) 算出第一步，但在第一次 BFGS 更新之前把它重设为

理由：设平均 Hessian \(\bar G_k\) 正定，令 \(z_k=\bar G_k^{1/2}s_k\)，则

右端的倒数是 \(\bar G_k\) 的 Rayleigh 商，近似 Hessian 的某个特征值；因此 (8.20) 让 \(H_0\) 的尺度与真实逆 Hessian 相当。其他缩放因子也可以用，但这一个在实践中最成功。

白话解释：Rayleigh 商 \(\dfrac{z^TAz}{z^Tz}\) 是"矩阵 \(A\) 沿方向 \(z\) 的平均放大倍数"，它总落在 \(A\) 的最小特征值和最大特征值之间。(8.21) 的推导只需把 \(y_k=\bar G_ks_k\)、\(s_k=\bar G_k^{-1/2}z_k\) 代入：分子 \(y_k^Ts_k=z_k^Tz_k\)，分母 \(y_k^Ty_k=z_k^T\bar G_kz_k\)。 一维看最清楚：\(f''=b\) 时 \(y=bs\)，\(\dfrac{ys}{y^2}=\dfrac1b\)，正好是逆 Hessian。所以 (8.20) 等于"用第一步测到的曲率，给逆 Hessian 定一个合适的量级"。 量级为什么重要：BFGS 第一步是 \(-H_0\nabla f_0\)。若 \(H_0=I\) 而真实曲率是 \(10^4\)，首步就长了一万倍，线搜索要回退很多次。

用 \(B_k\) 的实现。 也可以不存 \(H_k\)，而存 \(B_k\) 的 Cholesky 型分解 \(L_kD_kL_k^T\)，由 (8.19) 推出 \(O(n^2)\) 的因子更新公式，解 \(B_kp_k=-\nabla f_k\) 也只需 \(O(n^2)\)（两次三角回代）。好处是能在 \(D_k\) 的对角元过小时加以修正，但经验上没有实质优势，作者更偏好直接更新 \(H_k\) 的 Algorithm 8.1。

曲率条件失败时怎么办。 若只用 Armijo 回溯（先试 \(\alpha=1\) 再缩小直到充分下降），不能保证 \(y_k^Ts_k>0\)，因为满足曲率条件可能需要大于 1 的步长。有些实现在 \(y_k^Ts_k\le0\) 或接近零时跳过更新（令 \(H_{k+1}=H_k\)），原书不推荐：跳过可能发生得太频繁，使 \(H_k\) 捕获不到重要的曲率信息。更好的对策是原书第 18 章的阻尼 BFGS（damped BFGS），它把 \(y_k\) 换成 \(y_k\) 与 \(B_ks_k\) 的凸组合，使修正后的曲率条件总成立。

8.5 SR1 方法

8.5.1 推导

BFGS 和 DFP 都是秩二修正。能否用更简单的对称秩一（symmetric-rank-1, SR1）修正 \(B_{k+1}=B_k+\sigma vv^T\)（\(\sigma=\pm1\)）满足割线方程？代入 (8.6)：

方括号内是标量，所以 \(v\) 必须是 \(y_k-B_ks_k\) 的倍数，\(v=\delta(y_k-B_ks_k)\)。代回得

当且仅当 \(\sigma=\operatorname{sign}[s_k^T(y_k-B_ks_k)]\)、\(\delta=\pm|s_k^T(y_k-B_ks_k)|^{-1/2}\) 时成立。于是满足割线方程的唯一对称秩一更新是

由 Sherman–Morrison 公式，其逆形式为

注意 SR1 是自对偶的：(8.24) 换 \(B\to H\)、\(s\leftrightarrow y\) 就得到 (8.25)。推导如此简单，以致 SR1 被多次"重新发现"。

8.5.2 优点与缺点

不保证正定。 即使 \(B_k\) 正定，\(B_{k+1}\) 也可能不定。在只用线搜索的年代这被看作大缺陷；但在信赖域方法中，能产生不定近似恰恰是 SR1 的主要优点之一——真实 Hessian 本来就可能不定。

分母可能为零。 即使目标是凸二次函数，也可能某一步不存在满足割线方程的对称秩一更新。分三种情形：

1. \((y_k-B_ks_k)^Ts_k\ne0\)：唯一的秩一更新就是 (8.24)；
2. \(y_k=B_ks_k\)：唯一的"更新"是 \(B_{k+1}=B_k\)；
3. \(y_k\ne B_ks_k\) 但 \((y_k-B_ks_k)^Ts_k=0\)：由 (8.23)，不存在满足割线方程的对称秩一更新。

情形 3 说明秩一的自由度不够，可能导致数值不稳定甚至崩溃。

仍然重视 SR1 的理由。 (i) 一个简单的保护措施就足以防止崩溃；(ii) SR1 产生的矩阵往往是很好的 Hessian 近似，常常优于 BFGS；(iii) 在约束优化和部分可分优化（原书第 18、9 章）中，曲率条件 \(y_k^Ts_k>0\) 未必能保证，BFGS 不宜使用，而不定近似正好能反映真实 Hessian 的不定性。

跳过规则。 只有当

时才更新（\(r\in(0,1)\) 很小，如 \(10^{-8}\)），否则令 \(B_{k+1}=B_k\)。大多数 SR1 实现都用这类规则。

为什么 SR1 可以跳过而 BFGS 不宜跳过？ \(s_k^T(y_k-B_ks_k)\approx0\) 要求向量以特定方式对齐，很少发生；而且一旦发生，意味着 \(s_k^T\bar G_ks_k\approx s_k^TB_ks_k\)，即 \(B_k\) 沿 \(s_k\) 的曲率已经正确，跳过无害。BFGS 的曲率条件则在线搜索不施加 Wolfe 条件时很容易失败，跳过会频繁发生，损害近似质量。

8.5.3 SR1 信赖域方法

Algorithm 8.2（SR1 信赖域方法）：给定 \(x_0\)、\(B_0\)、半径 \(\Delta_0\)、容差 \(\epsilon>0\)、\(\eta\in(0,10^{-3})\)、\(r\in(0,1)\)。当 \(\|\nabla f_k\|>\epsilon\) 时：

1. 求解信赖域子问题 \(\min_s\nabla f_k^Ts+\tfrac12s^TB_ks\) s.t. \(\|s\|\le\Delta_k\)，得 \(s_k\)；
2. 计算 \(y_k=\nabla f(x_k+s_k)-\nabla f_k\)，实际下降 ared \(=f_k-f(x_k+s_k)\)，预测下降 pred \(=-(\nabla f_k^Ts_k+\tfrac12s_k^TB_ks_k)\)；
3. 若 ared/pred \(>\eta\)，接受 \(x_{k+1}=x_k+s_k\)，否则 \(x_{k+1}=x_k\)；
4. 半径更新：若 ared/pred \(>0.75\)，则 \(\|s_k\|\le0.8\Delta_k\) 时保持、否则加倍；若 \(0.1\le\) ared/pred \(\le0.75\)，保持；否则减半；
5. 若 (8.26) 成立，用 (8.24) 更新 \(B_{k+1}\)——即使这一步被拒绝（\(x_{k+1}=x_k\)）也更新；否则 \(B_{k+1}=B_k\)。

选信赖域框架，是因为它不需要把 Hessian 近似改成充分正定。第 5 步的要点值得记住：一步失败说明 \(B_k\) 在该方向近似得差；如果不加以修正，后面会反复生成类似的方向又反复被拒绝，阻碍超线性收敛。

8.5.4 SR1 的理论性质

定理 8.1（二次函数上的有限终止）：设 \(f(x)=b^Tx+\tfrac12x^TAx\)，\(A\) 对称正定。对任意 \(x_0\) 和任意对称（不必正定）的 \(H_0\)，只要对所有 \(k\) 有 \((s_k-H_ky_k)^Ty_k\ne0\)，以单位步长 \(p_k=-H_k\nabla f_k\)、\(x_{k+1}=x_k+p_k\) 进行的 SR1 迭代至多 \(n\) 步收敛到极小点；若执行了 \(n\) 步且方向线性无关，则 \(H_n=A^{-1}\)。

证明思路是归纳证明遗传性（hereditary property）：

即割线方程沿所有历史方向都成立。归纳步的关键计算是：对 \(j<k\)，由归纳假设和 \(y_i=As_i\)，

所以 (8.25) 中的修正项作用在 \(y_j\) 上为零，\(H_{k+1}y_j=H_ky_j=s_j\)；再加上新的割线方程 \(H_{k+1}y_k=s_k\)。若 \(n\) 个方向线性无关，\(H_nAs_j=s_j\) 对 \(n\) 个无关向量成立，故 \(H_nA=I\)，第 \(n+1\) 步就是牛顿步，直接到达解。若某一步方向与之前的方向线性相关，可推出 \(H_k\nabla f_{k+1}=0\)，由 \(H_k\) 非奇异得 \(\nabla f_{k+1}=0\)，已到达解。

意义：二次函数上 SR1 的遗传性与线搜索方式无关；BFGS 的类似结论只在精确线搜索下成立（见 8.6 节）。

白话解释：遗传性说的是"学过的不会忘"。每一步 SR1 都会学到一个方向上的正确曲率 \(H y_j=s_j\)，而且以后的更新不会破坏已经学到的方向。二次函数的 Hessian 是常数矩阵 \(A\)，\(n\) 维空间只要学满 \(n\) 个线性无关的方向，就把 \(A^{-1}\) 完全学会了，下一步就是精确的牛顿步。 关键计算 \((s_k-H_ky_k)^Ty_j=0\) 的含义是：SR1 的修正方向 \(s_k-H_ky_k\) 与所有旧的 \(y_j\) "互不干扰"，所以修正只改变新方向上的信息。用到的事实只有两个：归纳假设 \(H_ky_j=s_j\)，以及二次函数上 \(y_i=As_i\)（梯度变化严格等于 Hessian 乘位移）。

定理 8.2（一般非线性函数）：设 \(f\) 二阶连续可微，Hessian 在 \(x^*\) 附近有界且 Lipschitz 连续，\(\{x_k\}\) 是任意收敛到 \(x^*\) 的序列，(8.26) 对所有 \(k\) 成立，且步 \(s_k\) 一致线性无关（大致指步长不会趋于落在维数小于 \(n\) 的子空间里），则

这是 BFGS 不具备的性质：BFGS 矩阵一般不收敛到真 Hessian，只在步方向上收敛（8.7 节）。"一致线性无关"在实践中通常（但不总是）成立。

8.6 Broyden 族

把 BFGS 和 DFP 放进一个单参数族：

\(\phi_k=0\) 是 BFGS，\(\phi_k=1\) 是 DFP，一般地 \(B_{k+1}=(1-\phi_k)B_{k+1}^{\text{BFGS}}+\phi_kB_{k+1}^{\text{DFP}}\)。所以族内所有成员都满足割线方程；当 \(s_k^Ty_k>0\)、\(0\le\phi_k\le1\) 时保持正定。\(\phi_k\in[0,1]\) 的子族称为受限 Broyden 族（restricted Broyden class）。

定理 8.3：\(f\) 为强凸二次函数（Hessian \(A\)），用单位步长 \(p_k=-B_k^{-1}\nabla f_k\)，\(B_0\) 对称正定，按 \(\phi_k\in[0,1]\) 的 Broyden 公式更新。记 \(A^{1/2}B_k^{-1}A^{1/2}\) 的特征值为 \(\lambda_1^k\le\dots\le\lambda_n^k\)，则对所有 \(k\)，

若 \(\phi_k\notin[0,1]\)，此性质不再成立。

白话解释：\(A^{1/2}\) 是满足 \(A^{1/2}A^{1/2}=A\) 的对称正定矩阵（类比"标准差是方差的平方根"）。\(A^{1/2}B_k^{-1}A^{1/2}\) 衡量的是"近似 \(B_k\) 与真 Hessian \(A\) 的比值"：一维时就是 \(a/b_k\)。比值等于 1 说明估得正好，大于 1 说明低估了曲率，小于 1 说明高估了曲率。用 \(A^{1/2}\) 两边夹，是为了让这个"比值"仍是对称矩阵、有实特征值。

解读：特征值全为 1 等价于 \(B_k=A\)，即理想情形。(8.36) 说明特征值（非严格地）单调趋向 1。比如某步最小特征值 \(\lambda_1^k=0.7\)，下一步 \(\lambda_1^{k+1}\in[0.7,1]\)；若 \(\phi_k\) 在 \([0,1]\) 之外，它可能掉到 0.7 以下。这个结论即使线搜索不精确也成立。

不过最好的公式未必在受限族内：有分析和实验表明，在严格控制下允许 \(\phi_k\) 取负值可能优于 BFGS。SR1 也属于 Broyden 族，对应

它可能落在 \([0,1]\) 之外，所以不属于受限族。

奇异点。 (8.32) 的最后一项是秩一修正，由特征值交错定理，\(\phi_k<0\) 时特征值减小；\(\phi_k\) 继续减小，矩阵先奇异后不定。使 \(B_{k+1}\) 奇异的值为

由 Cauchy–Schwarz 不等式 \(\mu_k\ge1\)，所以 \(\phi_k^c\le0\)。若 \(B_0\) 对称正定、\(s_k^Ty_k>0\) 且 \(\phi_k>\phi_k^c\)，所有 \(B_k\) 都保持对称正定。

精确线搜索下的等价性。 精确线搜索时，所有 \(\phi_k\ge\phi_k^c\) 的 Broyden 族方法产生相同的迭代序列（对一般非线性函数也成立）：各方法的方向只差一个长度因子，精确线搜索找到同一个一维极小点。

定理 8.4（精确线搜索下的二次函数）：Broyden 族方法用于强凸二次函数，\(B_0\) 对称正定，\(\alpha_k\) 为精确步长，\(\phi_k\ge\phi_k^c\)，则

1. 至多 \(n\) 步收敛到解；
2. 割线方程对所有历史方向成立：\(B_ks_j=y_j\)，\(j=k-1,\dots,1\)；
3. 若 \(B_0=I\)，迭代与共轭梯度法完全相同，搜索方向两两共轭：\(s_i^TAs_j=0\)（\(i\ne j\)）；若 \(B_0\ne I\)，等价于以 \(B_0\) 为预条件子的预条件共轭梯度法；
4. 若执行 \(n\) 次迭代，则 \(B_n=A\)。

这类结果主要有理论意义——实际的非精确线搜索使各方法表现差异很大——但这类分析指导了拟牛顿法的大部分发展。

8.7 收敛性分析

总体结论先放在前面：BFGS 和 SR1 实践中都非常稳健，但对一般非凸目标，至今无法证明从任意起点出发的全局收敛。理论结果要么假设目标凸，要么假设迭代满足某些性质。局部超线性收敛则在合理假设下成立。本节记 \(G(x)=\nabla^2f(x)\)。

8.7.1 BFGS 的全局收敛（凸情形）

假设 8.1：(1) \(f\) 二阶连续可微；(2) 水平集 \(\mathcal L=\{x:f(x)\le f(x_0)\}\) 凸，且存在 \(0<m\le M\) 使

于是 \(f\) 在 \(\mathcal L\) 上强凸，有唯一极小点 \(x^*\)。利用 \(y_k=\bar G_ks_k\)，

定理 8.5：\(B_0\) 为任意对称正定矩阵，\(x_0\) 满足假设 8.1，则 BFGS（Algorithm 8.1）产生的 \(\{x_k\}\) 收敛到 \(x^*\)。

证明主线。 难点在于无法直接控制 \(B_k\) 的条件数。原书的巧妙之处是同时跟踪迹（控制最大特征值）和行列式（控制最小特征值）。对 (8.19) 取迹和行列式：

定义 \(s_k\) 与 \(B_ks_k\) 的夹角余弦和 Rayleigh 商

由于 \(s_k=-\alpha_kB_k^{-1}\nabla f_k\)，\(B_ks_k\) 与 \(-\nabla f_k\) 同向，所以 \(\theta_k\) 正是搜索方向与最速下降方向的夹角——这是线搜索全局收敛理论（Zoutendijk 定理）的核心量。

引入正定矩阵的势函数

因为 \(t-\ln t\ge1\)，\(\psi\) 同时惩罚过大的特征值（迹项）和过小的特征值（\(-\ln\det\) 项）。把 (8.44)(8.45) 代入并整理：

函数 \(h(t)=1-t+\ln t\le0\)（\(t>0\)），所以方括号非正；\(M_k-\ln m_k-1\le M-\ln m-1\equiv c\)。累加得

现在反证：若 \(\cos\theta_j\to0\)，则 \(\ln\cos^2\theta_j\to-\infty\)，后期每项都小于 \(-2c\)，右端最终变为负数，与 \(\psi>0\) 矛盾。所以存在子列 \(\cos\theta_{j_k}\ge\delta>0\)，由 Zoutendijk 结果得 \(\liminf\|\nabla f_k\|=0\)，强凸性再保证整个序列收敛到 \(x^*\)。

直观地说：\(\psi\) 有界就意味着 \(B_k\) 的特征值不能同时失控，于是方向不会长期与梯度接近正交——这就是"自校正"的数学体现。

推导拆解：这段证明的骨架只有四步，其余是代数整理。

1. 为什么要管 \(\cos\theta_k\)：\(\theta_k\) 是搜索方向与"最陡下坡方向"的夹角。Zoutendijk 定理（第 03 章）说，只要夹角不长期趋于 90°（\(\cos\theta_k\) 不趋于 0），配合 Wolfe 线搜索，梯度就会被压到零。所以全局收敛归结为证明 \(\cos\theta_k\) 不会趋于 0。
2. 为什么用 \(\psi\)：\(\psi=\sum_i(\lambda_i-\ln\lambda_i)\)。函数 \(t-\ln t\) 在 \(t=1\) 处最小（值为 1），\(t\to\infty\) 或 \(t\to0\) 时都趋于无穷。所以 \(\psi\) 不爆炸 ⇔ 没有特征值跑到无穷大或塌到零。\(\operatorname{tr}\) 对应 \(\sum\lambda_i\)，\(\ln\det\) 对应 \(\sum\ln\lambda_i\)，(8.44)(8.45) 正好给出了这两者每步怎么变。
3. 每步的增量：(8.50) 把 \(\psi\) 的增量拆成三块：第一块被常数 \(c\) 封顶；方括号形如 \(h(t)=1-t+\ln t\le0\)；最后一块 \(\ln\cos^2\theta_k\le0\)，而且 \(\cos\theta_k\) 越小，这一项越负。
4. 反证：如果 \(\cos\theta_k\) 一直趋于 0，最后一项会负得越来越厉害，压过每步至多 \(+c\) 的增长，累加到 \(k\) 步后 \(\psi\) 会变成负数，但 \(\psi\) 永远为正，矛盾。 一句话：每步"坏方向"都会让 \(\psi\) 付出代价，而 \(\psi\) 的预算是有限的，所以坏方向不可能一直出现。这有点像风控里的"亏损预算"：每次越界都扣额度，额度有限，越界就不能无限持续。

推广：定理 8.5 对整个受限 Broyden 族除 DFP 外都成立（\(\phi_k\in[0,1)\)）；\(\phi_k\to1\) 时论证失效，因为自校正性质被大大削弱。进一步可证收敛至少是线性的，且

8.7.2 BFGS 的超线性收敛

回忆 Dennis–Moré 条件（原书第 3 章 (3.32)）：拟牛顿方向的迭代超线性收敛，当且仅当

注意它只要求 \(B_k\) 沿搜索方向逼近真 Hessian。

假设 8.2：Hessian 在 \(x^*\) 处 Lipschitz 连续，\(\|G(x)-G(x^*)\|\le L\|x-x^*\|\)。

定理 8.6：\(f\) 二阶连续可微，BFGS 迭代收敛到满足假设 8.2 的极小点 \(x^*\)，且 (8.52) 成立，则 \(x_k\) 超线性收敛到 \(x^*\)。

证明主线。 做坐标变换 \(G_*=G(x^*)\)，\(\tilde s_k=G_*^{1/2}s_k\)，\(\tilde y_k=G_*^{-1/2}y_k\)，\(\tilde B_k=G_*^{-1/2}B_kG_*^{-1/2}\)。BFGS 公式在这一变换下形式不变（尺度不变性的又一次体现），所以 (8.50) 在新坐标下照样成立。在新坐标下理想矩阵是单位阵。由 Lipschitz 条件，

由此 \(\tilde M_k\le1+c\epsilon_k\)、\(\ln\tilde m_k\ge-2c\epsilon_k\)，代入 (8.50) 的类比式：

求和，并利用 (8.52) 得 \(\sum\epsilon_k<\infty\)：非负的两项 \(-\ln\cos^2\tilde\theta_j\) 与 \(-[\cdots]\) 的累加和有限，所以各自趋于 0，即 \(\cos\tilde\theta_j\to1\)、\(\tilde q_j\to1\)。而

换回原坐标就是 Dennis–Moré 条件。再结合"单位步长在解附近满足 Wolfe 条件"，得到超线性收敛。

要点：超线性收敛只需要 \(B_k\) 在步方向上逼近真 Hessian，不需要也不保证 \(B_k\to\nabla^2f(x^*)\)。8.8 节会看到这一点在统计推断中的后果。

金融直觉：把 \(B_k\) 想成用历史数据估出来的协方差矩阵。如果过去只观测到某几个因子组合的波动，那么沿这些组合方向的方差估计会很准，其他方向则基本是先验（初始值 \(H_0\) 留下的痕迹）。优化器只需要"走的方向"准就能快速收敛，所以它没有动力去修正其他方向；但算标准误时需要的是每一个参数方向的曲率，正好撞上了没被校准过的那些方向。 实务结论：res.hess_inv 可以当作"粗略看一眼"的参考，正式报告的标准误、t 统计量和置信区间，要用在最优点重新计算的 Hessian（解析、AD 或 7.2.3 节的有限差分）。

8.7.3 SR1 的收敛

SR1 的理论不如 BFGS 完善：除二次函数外，没有类似定理 8.5 的全局结果，也没有类似定理 8.6 的局部结果。但对信赖域 SR1 有：

定理 8.7：设 \(x_k\) 由 Algorithm 8.2 产生，且 (c1) 迭代不终止，停留在闭、有界、凸集 \(D\) 中，\(f\) 在 \(D\) 上二阶连续可微且有唯一驻点 \(x^*\)；(c2) \(\nabla^2f(x^*)\) 正定，\(\nabla^2f\) 在 \(x^*\) 附近 Lipschitz 连续；(c3) \(\{B_k\}\) 范数有界；(c4) 每步 (8.26) 成立。则 \(x_k\to x^*\) 且

即 \((n+1)\) 步超线性收敛。不要求精确求解子问题。BFGS 的分析不需要有界性假设 (c3)。还可以证明，在这些假设下 SR1 矩阵"大部分时间"是半正定的：半正定的迭代比例趋于 1，与 \(B_0\) 是否正定无关。

8.8 量化实战

8.8.1 在哪里用

• 参数估计与模型校准。 GARCH 族、随机波动率模型、利率期限结构模型的最大似然估计，Heston/SABR 等模型对市场价格的校准，绝大多数都用 BFGS 或 L-BFGS（scipy.optimize.minimize(method='BFGS')、R 的 optim、MATLAB 的 fminunc）。优化器报 "Desired error not necessarily achieved due to precision loss"、"line search failed" 时，本章的知识就是排查工具：检查梯度是否准确（有限差分误差）、变量尺度是否悬殊、线搜索是否满足 Wolfe 条件。
• 变量尺度。 GARCH 中 \(\omega\sim10^{-6}\)，而 \(\alpha,\beta\sim0.1\)。\(H_0\) 取单位阵时，初始几步在 \(\omega\) 方向上的步长完全不合尺度。(8.20) 的标量缩放只能修正整体量级，修不了各变量之间的相对量级，因此先把变量缩放到相近量级（例如用百分比收益率，或优化 \(\omega\times10^6\)）几乎总是值得的。
• BFGS 逆 Hessian 的误用。 常见做法是把 BFGS 结束时的 \(H_k\)（res.hess_inv）直接当作 MLE 的渐近协方差矩阵来算标准误。定理 8.6 告诉我们 \(H_k\) 只在搜索方向上逼近 \(\nabla^2f(x^*)^{-1}\)，其他方向上可能差得很远。正确做法是在最优点用解析 Hessian 或有限差分 Hessian 重新计算（统计背景见第 03 册第 09 章「参数推断」中 Fisher 信息与 MLE 渐近正态性的内容）。
• 非凸问题。 带非凸交易成本的组合优化、神经网络因子模型训练中常遇到不定 Hessian，此时 SR1 + 信赖域（scipy.optimize.minimize(method='trust-constr') 默认用 BFGS，也可指定 hess=SR1()）比线搜索 BFGS 更合适。

8.8.2 代码一：自己实现 BFGS，复现原书对比

下面按 Algorithm 8.1 实现 BFGS（含 (8.20) 缩放），用 SciPy 的强 Wolfe 线搜索（\(c_1=10^{-4}\)，\(c_2=0.9\)），与最速下降、线搜索牛顿法在 Rosenbrock 函数上对比。

import numpy as np
from scipy.optimize import line_search

def f(x):  return 100*(x[1]-x[0]**2)**2 + (1-x[0])**2
def g(x):  return np.array([-400*x[0]*(x[1]-x[0]**2) - 2*(1-x[0]), 200*(x[1]-x[0]**2)])
def h(x):  return np.array([[1200*x[0]**2-400*x[1]+2, -400*x[0]], [-400*x[0], 200.0]])

xstar = np.ones(2)

def run(method, x0, tol=1e-5, maxit=20000):
    x = x0.astype(float); n = len(x); H = np.eye(n); errs = []; k = 0
    gk = g(x)
    while np.linalg.norm(gk) > tol and k < maxit:
        if method == 'sd':
            p = -gk
        elif method == 'newton':
            p = np.linalg.solve(h(x), -gk)
            if gk @ p >= 0: p = -gk            # 非下降时退回最速下降
        else:
            p = -H @ gk                        # (8.18)
        a = line_search(f, g, x, p, gfk=gk, c1=1e-4, c2=0.9)[0]
        if a is None: a = 1e-3                 # 线搜索失败时的保护
        xn = x + a*p; gn = g(xn)
        s, y = xn - x, gn - gk
        if method == 'bfgs':
            if k == 0: H = (y @ s)/(y @ y) * np.eye(n)   # (8.20) 首次更新前缩放 H0
            rho = 1.0/(y @ s)
            V = np.eye(n) - rho*np.outer(y, s)
            H = V.T @ H @ V + rho*np.outer(s, s)          # (8.16)
        x, gk = xn, gn; k += 1
        errs.append(np.linalg.norm(x - xstar))
    return k, errs

x0 = np.array([-1.2, 1.0])
for m in ['sd', 'bfgs', 'newton']:
    k, e = run(m, x0)
    print(f"{m:7s} 迭代 {k:5d} 次，最后 4 步误差:", " ".join(f"{v:.2e}" for v in e[-4:]))

输出：

sd      迭代  7265 次，最后 4 步误差: 2.17e-05 2.17e-05 2.17e-05 2.17e-05
bfgs    迭代    40 次，最后 4 步误差: 2.57e-04 3.49e-05 8.64e-07 1.69e-07
newton  迭代    21 次，最后 4 步误差: 1.84e-02 1.95e-03 2.01e-05 2.85e-09

迭代次数与原书（5264 / 34 / 21）量级一致，差别来自线搜索实现细节。三种收敛速率的特征清楚可见：最速下降误差几乎不动；BFGS 相邻误差比为 0.14、0.025、0.20，整体趋于 0；牛顿法最后两步误差从 \(2\times10^{-5}\) 直接降到 \(3\times10^{-9}\)。

8.8.3 代码二：GARCH(1,1) 最大似然——变量缩放与标准误

模拟 3000 天的 GARCH(1,1) 收益率（\(\omega=2\times10^{-6}\)，\(\alpha=0.08\)，\(\beta=0.90\)），用 BFGS 最大化高斯似然。对比：(1) 不缩放 \(\omega\)；(2) 优化 \(\omega\times10^6\)。然后对比 res.hess_inv 与最优点处数值 Hessian 给出的标准误。

import numpy as np
from scipy.optimize import minimize

rng = np.random.default_rng(42)
T, omega, alpha, beta = 3000, 2e-6, 0.08, 0.90
r = np.empty(T); s2 = omega/(1-alpha-beta)
for t in range(T):
    r[t] = np.sqrt(s2)*rng.standard_normal()
    s2 = omega + alpha*r[t]**2 + beta*s2

def nll(theta, scale):
    w, a, b = theta[0]*scale, theta[1], theta[2]
    if w <= 0 or a < 0 or b < 0 or a + b >= 1:
        return 1e10                       # 不可行点给大值，线搜索会自动缩步
    s2 = np.empty(T); s2[0] = r.var()
    for t in range(1, T):
        s2[t] = w + a*r[t-1]**2 + b*s2[t-1]
    return 0.5*np.sum(np.log(s2) + r**2/s2)

def num_hess(fun, x, h=1e-4):
    n = len(x); H = np.zeros((n, n)); E = np.eye(n)*h
    for i in range(n):
        for j in range(n):
            H[i, j] = (fun(x+E[i]+E[j]) - fun(x+E[i]-E[j]) - fun(x-E[i]+E[j]) + fun(x-E[i]-E[j]))/(4*h*h)
    return H

# (1) 不缩放：ω 以原始量级出现
res0 = minimize(nll, [1e-6, 0.05, 0.9], args=(1.0,), method='BFGS')
print("不缩放 :", res0.nit, "次迭代, 参数", np.round(res0.x, 8), res0.message)
# (2) 缩放：优化变量取 ω×1e6，使三个参数量级相近
res = minimize(nll, [1.0, 0.05, 0.9], args=(1e-6,), method='BFGS')
print("缩放后 :", res.nit, "次迭代, ω=%.3e α=%.4f β=%.4f" % (res.x[0]*1e-6, res.x[1], res.x[2]))

# 标准误：BFGS 末步 H_k  vs  在最优点重新数值计算的 Hessian
se_bfgs = np.sqrt(np.diag(res.hess_inv))
Hn = num_hess(lambda th: nll(th, 1e-6), res.x)
se_num = np.sqrt(np.diag(np.linalg.inv(Hn)))
print("标准误（BFGS 的 H_k）:", np.round(se_bfgs, 4))
print("标准误（数值 Hessian）:", np.round(se_num, 4))

输出：

不缩放 : 13 次迭代, 参数 [1.8300000e-06 7.0327390e-02 9.1223426e-01] Desired error not necessarily achieved due to precision loss.
缩放后 : 17 次迭代, ω=1.921e-06 α=0.0711 β=0.9105
标准误（BFGS 的 H_k）: [0.9986 0.0117 0.0198]
标准误（数值 Hessian）: [0.6444 0.0112 0.0151]

两点结论。第一，不缩放时优化器以 "precision loss" 提前停下，\(\omega\) 停在 \(1.83\times10^{-6}\) 而非真正的极大似然点 \(1.92\times10^{-6}\)：\(\omega\) 方向的梯度量级与 \(\alpha,\beta\) 方向差了好几个数量级，线搜索无法同时照顾。第二，即使收敛良好，用 hess_inv 算出的 \(\omega\) 标准误（单位 \(10^{-6}\)）是 1.00，而正确的 0.64，高估约 55%；\(\beta\) 的标准误也高估约 30%。这正是定理 8.6 的实际后果：\(H_k\) 只在走过的方向上准确。做统计推断时，务必在最优点重新计算 Hessian。

本章小结

拟牛顿法通过割线方程 \(B_{k+1}s_k=y_k\) 把每步观测到的梯度变化累积为 Hessian 近似，只用梯度就得到超线性收敛。BFGS 是"在加权 Frobenius 范数下离旧逆矩阵最近、满足割线方程的对称矩阵"，它保持正定、自校正能力强，是最常用的通用算法，但必须配合 Wolfe 线搜索以保证曲率条件。SR1 是唯一满足割线方程的对称秩一更新，可产生不定近似，配合信赖域和跳过规则效果很好，并且常常是更准的 Hessian 近似。Broyden 族把这些方法统一起来。收敛理论的核心是势函数 \(\psi(B)=\operatorname{tr}B-\ln\det B\)（全局收敛）和 Dennis–Moré 条件（超线性收敛）；后者也告诉我们 BFGS 矩阵不收敛到真 Hessian，不能直接用来算标准误。

练习

基础

1. 证明：若 \(f\) 强凸（存在 \(m>0\) 使 \(\nabla^2f\succeq mI\)），则对任意两点 \(x_k\ne x_{k+1}\) 有 \(s_k^Ty_k\ge m\|s_k\|^2>0\)。 提示：\(y_k=\bar G_ks_k\)，\(\bar G_k\succeq mI\)。
2. 设一元函数 \(f\) 满足 \(f'(0)=-1\)，从 \(x_0=0\) 以步长 \(\alpha=1\) 走到 \(x_1=1\)。若 \(f'(1)=-\tfrac14\)，曲率条件是否成立？若 \(f'(1)=-2\) 呢？给出一个满足后一种情形的具体函数，并说明 Wolfe 线搜索会如何处理这一步。 提示：\(s=1\)，\(y=f'(1)-f'(0)\)；\(y=0.75>0\) 成立，\(y=-1<0\) 不成立。例如 \(f(x)=-x-x^2/2\)。此时 \(|f'(1)|>c_2|f'(0)|\)，Wolfe 第二条件不满足，线搜索会继续延长步长（本例函数无下界，实际问题中会在曲率变正处停下）。
3. 证明强 Wolfe 条件 \(|\nabla f(x_k+\alpha p_k)^Tp_k|\le c_2|\nabla f_k^Tp_k|\) 蕴含曲率条件 \(s_k^Ty_k>0\)。
4. 用 Sherman–Morrison–Woodbury 公式验证 (8.16) 与 (8.19) 互逆；用 Sherman–Morrison 公式验证 (8.24) 与 (8.25) 互逆。
5. 在 8.8.2 的代码中把 (8.20) 的缩放去掉（始终用 \(H_0=I\)），比较迭代次数和函数求值次数。再把 line_search 换成只满足 Armijo 条件的简单回溯，观察是否出现 \(y_k^Ts_k\le0\)。 提示：去掉缩放后首步线搜索需要更多次函数求值；纯回溯在 Rosenbrock 的弯曲谷底处容易出现曲率条件失败。

进阶

6. 证明 \(h(t)=1-t+\ln t\le0\) 对 \(t>0\) 成立，等号仅在 \(t=1\)；进而证明 \(\psi(B)=\sum(\lambda_i-\ln\lambda_i)>0\)。
7. 证明 \(\det(I+xy^T+uv^T)=(1+y^Tx)(1+v^Tu)-(x^Tv)(y^Tu)\)，并由此推出 (8.45)。 提示：把 \(B_{k+1}\) 写成 \(B_k(I-\cdots)\) 的形式，先用 \(\det(I+xy^T)=1+y^Tx\) 热身。
8. 用 SR1 逆形式 (8.25) 在三维强凸二次函数上做单位步长迭代，数值验证定理 8.1：3 步后 \(H_3=A^{-1}\)。
9. 在 8.8.3 中改用百分比收益率（\(r\times100\)）而不是缩放 \(\omega\)，重做估计。为什么这与缩放 \(\omega\) 等价？此时还应对哪个量的标准误做换算？ 提示：\(r\to100r\) 时 \(\omega\to10^4\omega\)，\(\alpha,\beta\) 不变；似然只差常数。
10. 对 8.8.3 的估计结果，分别用 res.hess_inv 与数值 Hessian 构造 \(\alpha+\beta\) 的 95% 置信区间，比较两者对"波动率持续性是否接近 1"这一判断的影响。

原书推荐习题：8.2（强 Wolfe ⇒ 曲率条件）、8.4（SR1 逆公式）、8.7、8.8、8.9、8.10（全局收敛证明所需的迹、行列式与 \(\psi\) 函数工具）。原书第 8 章共有 8.1–8.11 题。

原书对照

页码换算：本段原书正文页码约等于 PDF 页码减 20。