第 05 章 共轭梯度法

学习目标

读完本章，你应当能够：

1. 理解"求解对称正定方程组 \(Ax=b\)"与"极小化凸二次函数 \(\frac12x^TAx-b^Tx\)"的等价性，以及残差就是梯度。
2. 说清共轭方向法为什么至多 \(n\) 步收敛（定理 5.1）、扩展子空间极小化性质（定理 5.2），以及 CG 为什么只需前一个方向就能自动保持共轭（定理 5.3）。
3. 写出并实现标准 CG（Algorithm 5.2）与预条件 CG（Algorithm 5.3）。
4. 用特征值分布解释 CG 的收敛速度：\(r\) 个不同特征值 \(r\) 步终止、特征值成簇时极快、条件数界 \(\left(\frac{\sqrt\kappa-1}{\sqrt\kappa+1}\right)^k\) 优于最速下降。
5. 掌握非线性 CG（FR、PR、PR+、HS）的公式、对线搜索的要求及各自的收敛性质，知道实践中为什么推荐 PR/PR+。
6. 在"因子 + 特质"结构的大规模协方差矩阵上，用矩阵无关（matrix-free）的预条件 CG 求解组合问题。

读前导读

这一章在解决什么问题

先回到 Excel。没有约束的均值–方差问题，最优权重是 \(w^*=\gamma^{-1}\Sigma^{-1}\mu\)，在 Excel 里就是 MMULT(MINVERSE(Σ), μ)，Solver 都不用开。资产只有 10 只时这毫无问题；但做全市场选股、\(n=3000\) 时，协方差矩阵有 900 万个元素，求逆的计算量与 \(n^3\) 成正比，内存和时间都吃紧。更麻烦的是，很多大型问题根本不会把 \(\Sigma\) 完整存下来，只存因子暴露和特质方差。

本章的共轭梯度法（CG）提供了另一条路：不求逆，也不需要把矩阵完整写出来，只要你能算"矩阵乘以一个向量"，就能一步步逼近 \(\Sigma w=\mu\) 的解。它的核心想法有一个很好的金融类比：如果你能找到一组两两不相关的组合方向（关于 \(\Sigma\) 共轭），那么沿每个方向分别优化一次，前面已经调好的方向不会被后面的调整破坏，\(n\) 步就能做完。CG 的巧妙之处在于，它边走边造出这样一组方向，每次只需要记住上一个方向。

更实际的结论是：CG 的速度取决于矩阵特征值的分布，而不只是条件数。"因子 + 特质"风险模型的特征值天然只有少数几个"簇"，再配合合适的预条件（相当于先按特质波动率给股票做标准化），3000 只股票的组合十几步就能解完。本章后半部分把 CG 推广到一般的非线性函数（非线性 CG），它是第 9 章 L-BFGS 出现以前大规模优化的主力。

需要先想起来的数学

• 线性方程组与二次函数的关系：\(\phi(x)=\tfrac12x^TAx-b^Tx\) 的梯度是 \(Ax-b\)，令梯度为零就是解 \(Ax=b\)。所以"解方程"和"求二次函数最低点"是一回事。一元时：\(\tfrac12ax^2-bx\) 的最低点是 \(x=b/a\)。见 第 00 册第 05 章 多元微积分与优化。
• 线性无关、span 与子空间：\(\mathrm{span}\{v_1,\dots,v_k\}\) 是 \(v_1,\dots,v_k\) 所有线性组合的集合，类比"用这 \(k\) 个基础组合能复制出来的所有组合"。\(k\) 个线性无关的向量张成一个 \(k\) 维子空间。见 第 00 册第 06 章 线性代数速成。
• 特征值与特征向量：\(Av=\lambda v\)。对协方差矩阵，特征向量是主成分组合，特征值是该主成分的方差。"\(A\) 有 \(r\) 个不同特征值"就是"只有 \(r\) 种不同大小的主成分方差"。见第 00 册第 06 章。
• 多项式：\(P_k(\lambda)=c_0+c_1\lambda+\cdots+c_k\lambda^k\)。对矩阵也可以代入：\(P_k(A)=c_0I+c_1A+\cdots+c_kA^k\)。本章用"多项式在特征值处取值大小"来衡量 CG 的误差。
• 加权范数：\(\|z\|_A=\sqrt{z^TAz}\)。\(A\) 取协方差矩阵时，\(\|w-w^*\|_\Sigma\) 就是组合 \(w\) 相对目标组合 \(w^*\) 的跟踪误差。

怎么读这一章

必读：5.2.1、5.2.2（共轭的定义、定理 5.1、定理 5.2 的结论）、5.2.4 的 Algorithm 5.2、5.2.5 的定理 5.4 和条件数界 (5.35) 的结论、5.2.6 预条件的思想、5.4 量化实战。定理 5.3 的证明、定理 5.5 的构造第一次可以跳过。5.3 非线性 CG 部分，读 5.3.1、5.3.2 的公式和"推荐 PR/PR+"的结论即可；5.3.5–5.3.7 的收敛分析可以留到需要时再看。

5.1 引言

共轭梯度法（conjugate gradient, CG）有双重身份：

• 线性 CG：求解大型对称正定线性方程组的最有用的方法之一。Hestenes 与 Stiefel 在 1950 年代提出，原本作为高斯消元的替代（一种精确解法）；多年以后人们才认识到它是一种迭代法，往往在远少于 \(n\) 步时就给出很好的近似解，这被视为稀疏线性代数最重要的进展之一。它的表现与系数矩阵的特征值分布密切相关，通过预条件（preconditioning）改善特征值分布可以显著加速。
• 非线性 CG：Fletcher 与 Reeves 在 1960 年代把它推广到一般非线性优化，这是最早的大规模非线性优化技术之一。它的关键特点是不需要存储矩阵，又比最速下降快得多。

对量化来说，线性 CG 是求解大规模组合与风险模型方程（\(\Sigma w=\mu\)、广义最小二乘、岭回归正规方程）的利器，也是第 4、6 章信赖域/截断牛顿法的内核；非线性 CG 则是内存受限时的大规模优化工具（现在更多被 L-BFGS 取代，见本册第 9 章）。

5.2 线性共轭梯度法

5.2.1 问题的两种形式

求解

等价于极小化

二者有相同的唯一解。\(\phi\) 的梯度就是方程组的残差：

组合优化里，\(\min_w\frac\gamma2w^T\Sigma w-\mu^Tw\) 与 \(\gamma\Sigma w=\mu\) 正是这样一对。

5.2.2 共轭方向法

定义（共轭）　非零向量组 \(\{p_0,\dots,p_l\}\) 称为关于对称正定矩阵 \(A\) 共轭（conjugate），如果

共轭向量组必定线性无关（练习 1），所以 \(A\) 至多有 \(n\) 个共轭方向。

金融直觉：把 \(A\) 换成协方差矩阵 \(\Sigma\)，\(p_i^T\Sigma p_j\) 就是组合 \(p_i\) 与组合 \(p_j\) 收益的协方差。所以"关于 \(\Sigma\) 共轭"就是"两个组合不相关"。普通的"正交"（\(p_i^Tp_j=0\)）只是说权重向量在几何上垂直，并不意味着收益不相关；共轭才是金融上有意义的"互不干扰"。 为什么不相关这么重要？组合方差 \(\mathrm{Var}(\sum_i c_ip_i)=\sum_i c_i^2\,p_i^T\Sigma p_i\)，交叉项全消失了。于是在每个方向上选多少仓位（\(c_i\)）可以各自独立地决定，调整第 3 个方向的仓位不会改变前两个方向的最优仓位。这就是共轭方向法能"逐个方向一次做对"的原因。主成分组合是一组特殊的共轭方向（既正交又不相关），但求它们需要做完整的特征分解，太贵。 共轭方向法：给定 \(x_0\) 和一组共轭方向 \(\{p_0,\dots,p_{n-1}\}\)，令

\(\alpha_k\) 是沿 \(p_k\) 方向精确极小化 \(\phi\) 的步长（即第 3 章的 (3.39)）。

定理 5.1　对任意 \(x_0\)，共轭方向法至多 \(n\) 步收敛到 (5.1) 的解 \(x^*\)。

证明：方向线性无关，张成 \(\mathbb{R}^n\)，可写 \(x^*-x_0=\sum_k\sigma_kp_k\)。左乘 \(p_k^TA\) 并用共轭性：

又 \(x_k=x_0+\sum_{i<k}\alpha_ip_i\)，由共轭性 \(p_k^TA(x_k-x_0)=0\)，于是

所以 \(\sigma_k=\alpha_k\)：第 \(k\) 步恰好消掉误差在 \(p_k\) 方向上的分量，\(n\) 步后误差为零。\(\square\)

推导拆解：证明的逻辑是"先算出理想答案，再证明算法恰好给出它"。 第 1 步：\(n\) 个线性无关的方向构成一组"基"，所以从起点到终点的总位移 \(x^*-x_0\) 一定能写成 \(\sum_k\sigma_kp_k\)，就像任何组合都能用 \(n\) 只基础资产复制。 第 2 步：求系数 \(\sigma_k\)。两边左乘 \(p_k^TA\)，右边求和中 \(i\ne k\) 的项 \(p_k^TAp_i=0\)（共轭性），只剩 \(\sigma_kp_k^TAp_k\)，于是得到 (5.7)。这一步和"用不相关的因子做回归，每个系数可以单独算"是同一个道理。 第 3 步：证明算法的步长 \(\alpha_k\) 等于 \(\sigma_k\)。\(x_k-x_0\) 是 \(p_0,\dots,p_{k-1}\) 的组合，与 \(p_k\) 共轭，所以 \(p_k^TA(x_k-x_0)=0\)，可以把 (5.7) 分子里的 \(x_0\) 换成 \(x_k\)。再用 \(Ax^*=b\) 和 \(r_k=Ax_k-b\)：\(A(x^*-x_k)=b-Ax_k=-r_k\)。代回就得到 \(\sigma_k=-p_k^Tr_k/p_k^TAp_k=\alpha_k\)。 几何解释　如果 \(A\) 是对角阵，等高线椭圆的轴与坐标轴对齐，依次沿 \(e_1,\dots,e_n\) 做一维极小化，\(n\) 步就到达解（原书图 5.1）；\(A\) 非对角时，坐标轮换不再有限步收敛（图 5.2）。做变量变换 \(\hat x=S^{-1}x\)，\(S=[p_0\ p_1\cdots p_{n-1}]\)（5.8），则

由共轭性 \(S^TAS\) 是对角阵。所以共轭方向法就是"在使 Hessian 对角化的坐标系里做坐标下降"——这把第 3 章的坐标下降与本章联系了起来。

对角情形还有一个性质：每次坐标极小化都把解的一个分量定准了，\(k\) 步后已经在 \(e_1,\dots,e_k\) 张成的子空间上极小化了 \(\phi\)。一般情形由下面的定理给出。先注意残差的递推：

定理 5.2（扩展子空间极小化）　对任意 \(x_0\)，共轭方向法生成的 \(x_k\) 满足

且 \(x_k\) 是 \(\phi\) 在仿射集

上的极小点。

证明：\(\tilde x\) 是 \(\phi\) 在 (5.11) 上的极小点，当且仅当 \(r(\tilde x)^Tp_i=0\)，\(i<k\)（因为 \(h(\sigma)=\phi(x_0+\sum\sigma_ip_i)\) 是严格凸二次函数，令偏导数为零并用链式法则）。下面归纳证明 (5.10)。\(k=1\)：\(r_1^Tp_0=0\)，因为 \(\alpha_0\) 是一维极小步。设 \(r_{k-1}^Tp_i=0\)（\(i\le k-2\)）。由 (5.9)，

"当前残差与之前所有搜索方向正交"这一性质 (5.10) 在本章中被反复使用。

剩下的问题是：共轭方向从哪里来？\(A\) 的特征向量既正交又共轭，但大规模问题求全部特征向量太贵；用修改的 Gram–Schmidt 过程也能生成共轭方向，但要存储整个方向组，同样昂贵。CG 给出了一个巧妙的答案。

5.2.3 CG 的基本性质

CG 是一种特殊的共轭方向法：生成新方向 \(p_k\) 只用到前一个方向 \(p_{k-1}\)，而它会自动与所有更早的方向共轭——存储和计算都极少。令

由要求 \(p_{k-1}^TAp_k=0\) 得 \(\beta_k=\dfrac{r_k^TAp_{k-1}}{p_{k-1}^TAp_{k-1}}\)；第一个方向 \(p_0\) 取最速下降方向 \(-r_0\)。

Algorithm 5.1（CG 初步版）　\(r_0=Ax_0-b\)，\(p_0=-r_0\)；当 \(r_k\ne0\) 时：

定义 Krylov 子空间（Krylov subspace）

定理 5.3　若第 \(k\) 个迭代点不是解，则

因此 \(\{x_k\}\) 至多 \(n\) 步收敛到 \(x^*\)。

证明（对 (5.16)–(5.18) 归纳）：

• (5.16)：由归纳假设 \(r_k,p_k\in\mathcal{K}(r_0;k)\)，所以 \(Ap_k\in\mathrm{span}\{Ar_0,\dots,A^{k+1}r_0\}\)（5.19）；由 (5.9)，\(r_{k+1}\in\mathcal{K}(r_0;k+1)\)，得"⊂"。反过来，\(A^{k+1}r_0=A(A^kr_0)\in\mathrm{span}\{Ap_0,\dots,Ap_k\}\)，而 \(Ap_i=(r_{i+1}-r_i)/\alpha_i\)，所以 \(A^{k+1}r_0\in\mathrm{span}\{r_0,\dots,r_{k+1}\}\)。
• (5.17)：由 (5.13e)，\(\mathrm{span}\{p_0,\dots,p_{k+1}\}=\mathrm{span}\{p_0,\dots,p_k,r_{k+1}\}=\mathrm{span}\{r_0,\dots,A^kr_0,r_{k+1}\}=\mathcal{K}(r_0;k+1)\)。
• (5.18)：\(p_{k+1}^TAp_i=-r_{k+1}^TAp_i+\beta_{k+1}p_k^TAp_i\)（5.20）。\(i=k\) 时由 \(\beta_{k+1}\) 的定义为零。\(i\le k-1\) 时：由于 \(p_0,\dots,p_k\) 已是共轭的，定理 5.2 给出 \(r_{k+1}^Tp_i=0\)（\(i\le k\)）（5.21）；又 \(Ap_i\in\mathrm{span}\{Ar_0,\dots,A^{i+1}r_0\}\subset\mathrm{span}\{p_0,\dots,p_{i+1}\}\)（5.22），所以 \(r_{k+1}^TAp_i=0\)；第二项由归纳假设为零。
• (5.15)：由 (5.10) \(r_k^Tp_i=0\)，而 \(p_i=-r_i+\beta_ip_{i-1}\) 说明 \(r_i\in\mathrm{span}\{p_i,p_{i-1}\}\)，所以 \(r_k^Tr_i=0\)。\(\square\)

两点注意：证明依赖于 \(p_0=-r_0\)，换别的初始方向结论就不成立；另外，CG 中是梯度（残差）相互正交、搜索方向关于 \(A\) 共轭，所以"共轭梯度"这个名字其实不太准确。

白话解释：定理 5.3 的四条结论可以这样串起来理解。Krylov 子空间 \(\mathcal{K}(r_0;k)\) 是"从初始残差出发，反复乘 \(A\) 能到达的方向"。CG 每步只做一次矩阵–向量乘法，所以 \(k\) 步后它能"看到"的信息恰好就是这个子空间，(5.16)(5.17) 说的就是这一点。在这个前提下，CG 做到了最好：新残差与所有旧残差正交 (5.15)、新方向与所有旧方向共轭 (5.18)，而这一切只靠和上一个方向做一次组合就自动成立——这是 CG 能用 \(O(n)\) 内存解大问题的关键。 用 \(\Sigma\) 来想：\(\Sigma r_0\) 可以理解为"把初始误差信号通过协方差结构传播一次"。每多乘一次 \(\Sigma\)，就多吸收一层资产间的相关信息。

5.2.4 实用形式

利用 (5.13e) 和 (5.10)，\(\alpha_k\) 可化为 \(\dfrac{r_k^Tr_k}{p_k^TAp_k}\)；利用 \(\alpha_kAp_k=r_{k+1}-r_k\) 和残差正交性，\(\beta_{k+1}\) 可化为 \(\dfrac{r_{k+1}^Tr_{k+1}}{r_k^Tr_k}\)。得到标准形式：

推导拆解：补上两处化简的中间式。 \(\alpha_k\)：把 \(p_k=-r_k+\beta_kp_{k-1}\) 代入分子 \(-r_k^Tp_k\)，得 \(r_k^Tr_k-\beta_kr_k^Tp_{k-1}\)；由 (5.10)，当前残差与旧方向正交，\(r_k^Tp_{k-1}=0\)，所以分子就是 \(r_k^Tr_k\)。 \(\beta_{k+1}\)：由 (5.9)，\(Ap_k=(r_{k+1}-r_k)/\alpha_k\)。分子 \(r_{k+1}^TAp_k=(r_{k+1}^Tr_{k+1}-r_{k+1}^Tr_k)/\alpha_k=r_{k+1}^Tr_{k+1}/\alpha_k\)（残差两两正交）。分母 \(p_k^TAp_k=r_k^Tr_k/\alpha_k\)（就是刚才 \(\alpha_k\) 公式的变形）。两者相除，\(\alpha_k\) 约掉。 新形式的好处：\(r_k^Tr_k\) 上一步已经算过，可以直接复用，每步只需一次矩阵–向量乘法 \(Ap_k\)。

Algorithm 5.2（CG）

给定 x_0；r_0 ← A x_0 − b；p_0 ← −r_0；k ← 0
while r_k ≠ 0
  α_k ← r_kᵀr_k / p_kᵀA p_k              (5.23a)
  x_{k+1} ← x_k + α_k p_k                  (5.23b)
  r_{k+1} ← r_k + α_k A p_k                (5.23c)
  β_{k+1} ← r_{k+1}ᵀr_{k+1} / r_kᵀr_k      (5.23d)
  p_{k+1} ← −r_{k+1} + β_{k+1} p_k         (5.23e)
  k ← k+1
end

• 只需保存最近一步的 \(x,r,p\)（新值覆盖旧值）。每步主要计算是一次矩阵–向量乘 \(Ap_k\)，外加两个内积、三个向量加法（\(O(n)\) 次运算）。
• 原书强调：CG 只推荐用于大规模问题；小问题应该用高斯消元、Cholesky 或 SVD 等分解法，它们对舍入误差更不敏感。大规模时 CG 的优点是：不改变系数矩阵，不产生填充（fill-in）；只要能算 \(Av\)，甚至不需要显式存储 \(A\)（matrix-free）；而且有时收敛非常快。

5.2.5 收敛速度

精确算术下 CG 至多 \(n\) 步终止；更重要的是，特征值分布有利时，远少于 \(n\) 步。

由 (5.23b) 和 (5.17)，

\(P_k^*\) 是某个 \(k\) 次多项式。用 \(A\)-范数 \(\|z\|_A^2=z^TAz\)（5.26），有 \(\frac12\|x-x^*\|_A^2=\phi(x)-\phi(x^*)\)（5.27）。由定理 5.2，\(P_k^*\) 恰是

的解。也就是说，在所有前 \(k\) 步迭代限制在 Krylov 子空间中的方法里，CG 使 \(A\)-范数误差最小——这是 CG 的最优性。

因为 \(r_0=A(x_0-x^*)\)，有 \(x_{k+1}-x^*=[I+P_k^*(A)A](x_0-x^*)\)（5.29）。设 \(A\) 的特征值 \(0<\lambda_1\le\cdots\le\lambda_n\) 对应正交特征向量 \(v_i\)，\(x_0-x^*=\sum_i\xi_iv_i\)（5.30），则

所以收敛速度由

刻画：找一个在 0 处取值 1、在所有特征值处都尽量小的 \(k+1\) 次多项式。特征值越少、越集中，这样的多项式越容易找。

白话解释：(5.31) 的意思是：把初始误差分解到各个主轴（特征向量）上，第 \(i\) 个分量经过 \(k+1\) 步后被乘上了 \(Q(\lambda_i)=1+\lambda_iP_k(\lambda_i)\) 这个"衰减系数"。CG 自动选出最好的多项式 \(Q\)，让所有特征值处的衰减系数都尽量小。约束是 \(Q(0)=1\)（多项式的常数项被固定）。 小例子：特征值只有 2 和 5 两个。取 \(Q(\lambda)=(1-\lambda/2)(1-\lambda/5)\)，它在 0 处为 1，在 2 和 5 处都为 0，是 2 次多项式，所以 CG 2 步就把误差清零——这就是定理 5.4。如果特征值是 2、2.01、5、5.02，同一个 \(Q\) 在 2.01 和 5.02 处的值也非常接近 0，所以 2 步后误差已经很小，这就是"成簇"的威力。 金融含义：特征值"成簇"相当于很多主成分的方差几乎一样大。CG 不在乎有多少只股票，只在乎有多少种"不同大小的风险"。 定理 5.4　若 \(A\) 只有 \(r\) 个不同的特征值，则 CG 至多 \(r\) 步终止于解。

证明：设不同特征值为 \(\tau_1<\cdots<\tau_r\)，令

则 \(Q_r(\lambda_i)=0\) 对所有特征值成立，且 \(Q_r(0)=1\)。于是 \(\bar P_{r-1}(\lambda)=(Q_r(\lambda)-1)/\lambda\) 是 \(r-1\) 次多项式，代入 (5.33)（\(k=r-1\)）得到 0，所以 \(x_r=x^*\)。\(\square\)

定理 5.5（Luenberger）

思路：取 \(Q_{k+1}(\lambda)=1+\lambda\bar P_k(\lambda)\)，让它的根落在最大的 \(k\) 个特征值 \(\lambda_n,\dots,\lambda_{n-k+1}\) 以及 \(\lambda_1\) 与 \(\lambda_{n-k}\) 的中点；可以验证它在其余特征值上的最大绝对值恰为 \((\lambda_{n-k}-\lambda_1)/(\lambda_{n-k}+\lambda_1)\)。

应用：特征值成簇　设 \(A\) 有 \(m\) 个大特征值，其余 \(n-m\) 个都聚集在 1 附近（原书图 5.3）。令 \(\epsilon=\lambda_{n-m}-\lambda_1\)，由定理 5.5，\(m+1\) 步后 \(\|x_{m+1}-x^*\|_A\approx\epsilon\|x_0-x^*\|_A\)。\(\epsilon\) 小时，\(m+1\) 步就得到很好的近似——不管那 \(m\) 个大特征值有多大。

• 原书图 5.4 的数值例：5 个大特征值，其余在 \([0.95,1.05]\) 内。定理 5.5 预测第 6 步误差陡降，实际第 5 步就降了（定理只是上界）；第 7 步再次陡降。作为对比，特征值随机均匀分布的问题收敛更慢也更均匀。
• 一般规律：特征值分成 \(r\) 个簇时，CG 大约 \(r\) 步就近似解出（构造在每簇内都有零点的多项式）。图 5.5：\(n=14\)，四簇（140、120 各一个，10 附近 10 个，其余在 \([0.95,1.05]\)），4 步后误差已很小，6 步后精确。

基于条件数的界　记 \(\kappa(A)=\|A\|_2\|A^{-1}\|_2=\lambda_n/\lambda_1\)，有

（精读所据 PDF 文本中此式写作指数 \(2k\)、无系数 2，条件数也误印为 \(\lambda_1/\lambda_n\)；这里采用数值线性代数的标准形式。）这个界常常严重高估误差，但在只知道极端特征值时很有用。把它与最速下降的 (3.29) 对比：形式相同，但依赖 \(\sqrt\kappa\) 而不是 \(\kappa\)——这是 CG 比最速下降快得多的根本原因。例如 \(\kappa=10^4\) 时，最速下降每步误差的缩减因子约为 \(1-2\times10^{-4}\)，CG 约为 \(1-0.02\)。

白话解释：把这两个数字换算成"误差缩小 10 倍需要几步"：最速下降约 \(\ln 10/(2\times10^{-4})\approx 11500\) 步，CG 约 \(\ln10/0.02\approx 115\) 步，差 100 倍，正好是 \(\sqrt{\kappa}\)。这两个缩减因子来自近似 \(\frac{\kappa-1}{\kappa+1}\approx1-\frac2\kappa\) 和 \(\frac{\sqrt\kappa-1}{\sqrt\kappa+1}\approx1-\frac2{\sqrt\kappa}\)（\(\kappa\) 很大时成立）。 为什么 CG 能做到 \(\sqrt\kappa\)？最速下降每一步只"记得"当前梯度，相当于每步只用 1 次多项式；CG 第 \(k\) 步用的是 \(k\) 次多项式，而在区间 \([\lambda_1,\lambda_n]\) 上，高次的切比雪夫多项式能把衰减系数压得比 \(\left(\frac{\kappa-1}{\kappa+1}\right)^k\) 小得多。

5.2.6 预条件

既然收敛速度取决于特征值分布，就可以通过变量变换来改善它。令

\(C\) 非奇异，则

即求解 \((C^{-T}AC^{-1})\hat x=C^{-T}b\)，收敛速度取决于 \(C^{-T}AC^{-1}\) 的特征值。目标：选 \(C\) 使 \(C^{-T}AC^{-1}\) 的条件数远小于 \(\kappa(A)\)，或者使其特征值成簇。推导之后会发现算法不需要显式做变换，只用到 \(M=C^TC\)（对称正定），称为预条件子（preconditioner）。

金融直觉：预条件就是"换单位"。取 \(C=\mathrm{diag}(\sigma_1,\dots,\sigma_n)\)（各资产波动率），\(\hat x=Cx\) 就是把"权重"换成"风险单位"（权重 × 波动率），\(C^{-T}\Sigma C^{-1}\) 变成相关系数矩阵——这正是第 2 章量化实战里条件数从 16442 降到 2.5 的做法。预条件子 \(M\) 要像 \(A\)（这样 \(M^{-1}A\approx I\)），但又要让 \(My=r\) 容易解。对角矩阵是最极端的"容易解"：除一下就行。 理想情况 \(M=A\) 时一步收敛，但解 \(My=r\) 就等于解原问题，没有意义。实用的预条件子是两者之间的折中。 Algorithm 5.3（预条件 CG）

给定 x_0，预条件子 M；r_0 ← A x_0 − b；解 M y_0 = r_0；p_0 ← −y_0；k ← 0
while r_k ≠ 0
  α_k ← r_kᵀy_k / p_kᵀA p_k               (5.38a)
  x_{k+1} ← x_k + α_k p_k                  (5.38b)
  r_{k+1} ← r_k + α_k A p_k                (5.38c)
  解 M y_{k+1} = r_{k+1}                    (5.38d)
  β_{k+1} ← r_{k+1}ᵀy_{k+1} / r_kᵀy_k      (5.38e)
  p_{k+1} ← −y_{k+1} + β_{k+1} p_k         (5.38f)
end

（PDF 文本中初始化印作 \(p_0=-r_0\)，应为 \(p_0=-y_0\)，此处已更正。）\(M=I\) 时退化为标准 CG。残差的正交性变为

与无预条件相比，主要额外代价是每步求解一次 \(My=r\)。

实用预条件子

• 没有对所有矩阵都"最好"的预条件子，需要在三者间权衡：\(M\) 的有效性、构造和存储 \(M\) 的代价、求解 \(My=r\) 的代价。
• 特定类型的矩阵有好策略。例如偏微分方程离散化产生的矩阵，\(My=r\) 常取原系统的简化版本（如更粗网格上的离散）。了解问题的结构和来源是设计有效预条件子的关键。
• 通用预条件子：对称逐次超松弛（SSOR）、不完全 Cholesky（incomplete Cholesky）、带状预条件子。不完全 Cholesky 通常最有效：按 Cholesky 过程计算，但只保留比真实因子 \(L\) 更稀疏的近似 \(\tilde L\)（通常不比 \(A\) 的下三角部分更稠密），\(A\approx\tilde L\tilde L^T\)。取 \(C=\tilde L^T\)，则 \(M=\tilde L\tilde L^T\)，\(C^{-T}AC^{-1}=\tilde L^{-1}A\tilde L^{-T}\approx I\)。求解 \(My=r\) 只需两次三角回代，代价与一次 \(Ap\) 相当。
• 陷阱：不完全分解的结果可能不（够）正定，需要增大对角元；为保稀疏而丢弃元素可能导致数值不稳定甚至分解中断，可以允许更多填充，但代价更高。

最简单的预条件子是对角（Jacobi）预条件 \(M=\mathrm{diag}(A)\)。在量化的因子风险模型里，它恰好是"最懂问题结构"的选择，见本章量化实战。

5.3 非线性共轭梯度法

5.3.1 Fletcher–Reeves 方法

对 Algorithm 5.2 做两处改动，就得到一般非线性函数 \(f\) 的 CG：(1) 步长 \(\alpha_k\) 改为沿 \(p_k\) 做线搜索得到的近似极小；(2) 残差 \(r\) 换成 \(f\) 的梯度。

Algorithm 5.4（FR-CG）

给定 x_0；计算 f_0, ∇f_0；p_0 ← −∇f_0；k ← 0
while ∇f_k ≠ 0
  线搜索得 α_k，x_{k+1} = x_k + α_k p_k；计算 ∇f_{k+1}
  β^{FR}_{k+1} ← ∇f_{k+1}ᵀ∇f_{k+1} / ∇f_kᵀ∇f_k       (5.40a)
  p_{k+1} ← −∇f_{k+1} + β^{FR}_{k+1} p_k             (5.40b)
end

• \(f\) 为强凸二次函数且用精确线搜索时，它就是线性 CG。每步只需函数值和梯度，没有矩阵运算，只存几个向量，适合大规模问题。
• 下降性依赖线搜索：\(\nabla f_k^Tp_k=-\|\nabla f_k\|^2+\beta_k^{FR}\nabla f_k^Tp_{k-1}\)（5.41）。精确线搜索时 \(\nabla f_k^Tp_{k-1}=0\)，\(p_k\) 必为下降方向；非精确线搜索时第二项可能占优，使 \(p_k\) 变成上升方向。解决办法是要求步长满足强 Wolfe 条件
  \[f(x_k+\alpha_kp_k)\le f(x_k)+c_1\alpha_k\nabla f_k^Tp_k,\qquad|\nabla f(x_k+\alpha_kp_k)^Tp_k|\le-c_2\nabla f_k^Tp_k,\tag{5.42}\]
  并且 \(0<c_1<c_2<\frac12\)（注意比第 3 章的 \(c_2<1\) 更严）。下面的引理 5.6 说明这样就保证了所有方向都是下降方向。

5.3.2 Polak–Ribière 方法及其他

• 强凸二次 + 精确线搜索时梯度相互正交 (5.15)，\(\beta^{PR}=\beta^{FR}\)；但对一般非线性函数和非精确线搜索，两者表现差异显著，数值经验表明 PR-CG 更稳健、更高效。
• 一个意外的事实：强 Wolfe 条件不能保证 PR 方向是下降方向。若定义
  \[\beta_{k+1}^+=\max\{\beta_{k+1}^{PR},0\},\tag{5.44}\]
  得到 PR+ 算法，对强 Wolfe 条件稍作修改即可保证下降性。
• Hestenes–Stiefel 公式：
  \[\beta_{k+1}^{HS}=\frac{\nabla f_{k+1}^T(\nabla f_{k+1}-\nabla f_k)}{(\nabla f_{k+1}-\nabla f_k)^Tp_k},\tag{5.45}\]
  理论和实践都与 PR 相似。推导：要求相邻方向关于线段 \([x_k,x_{k+1}]\) 上的平均 Hessian \(\bar G_k=\int_0^1\nabla^2f(x_k+\tau\alpha_kp_k)d\tau\) 共轭；由 \(\nabla f_{k+1}=\nabla f_k+\alpha_k\bar G_kp_k\)，令 \(p_{k+1}^T\bar G_kp_k=0\) 即得 (5.45)。
• 其他各种 \(\beta\) 的选择都没有显著优于 PR。

5.3.3 二次终止与重启

• 实现中通常在线搜索里包含沿 \(p_k\) 的二次或三次插值，保证当 \(f\) 是严格凸二次函数时步长是精确的，于是算法退化为线性 CG。
• 重启：每 \(n\) 步令 \(\beta_k=0\)（走一次最速下降），周期性地丢掉可能无益的旧信息。理论上可得 \(n\) 步二次收敛：
  \[\|x_{k+n}-x^*\|=O(\|x_k-x^*\|^2).\tag{5.46}\]
  直观理解：若 \(f\) 在解附近是强凸二次函数，迭代进入该区域后某次重启，之后就是线性 CG，\(n\) 步内终止。重启之所以重要，是因为线性 CG 的有限终止性要求 \(p_0\) 是负梯度。
• 但实践意义有限：非线性 CG 只推荐用于大 \(n\)，常常在少于 \(n\) 步时就得到近似解，重启可能根本不发生。所以实践中要么不重启，要么用别的准则。最流行的准则基于二次函数梯度正交性 (5.15)：当相邻梯度远非正交时重启：
  \[\frac{|\nabla f_k^T\nabla f_{k-1}|}{\|\nabla f_k\|^2}\ge\nu,\qquad\nu\approx0.1.\tag{5.47}\]
• (5.44) 也可以看成一种重启（\(\beta^{PR}<0\) 时回到最速下降），但 \(\beta^{PR}\) 多数时候为正，这种重启很少发生。

5.3.4 数值表现

原书表 5.1 比较了不重启的 FR、PR、PR+（强 Wolfe 参数 \(c_1=10^{-4}\)、\(c_2=0.1\)；终止条件 \(\|\nabla f_k\|_\infty<10^{-5}(1+|f_k|)\)，或 10000 次迭代，记为 *）。表中为"迭代次数/函数–梯度求值次数"，mod 列是 PR+ 中 (5.44) 实际起作用的次数。

FR 在 GENROS 上远离解时步长极短、几乎没有改进。PR/PR+ 并非总比 FR 好，而且多需要存一个向量，但作者推荐尽量使用 PR 或 PR+。

5.3.5 Fletcher–Reeves 方法的行为

假设水平集 \(\mathcal{L}\) 有界、\(f\) 二阶连续可微，由引理 3.1 存在满足强 Wolfe 条件的步长。

引理 5.6　设 Algorithm 5.4 的步长满足强 Wolfe 条件 (5.42) 且 \(0<c_2<\frac12\)，则所有 \(p_k\) 都是下降方向，且

证明：函数 \(t(\xi)=(2\xi-1)/(1-\xi)\) 在 \([0,\frac12]\) 上单调增，\(t(0)=-1\)，\(t(\frac12)=0\)，所以 \(c_2\in(0,\frac12)\) 时

因此 (5.48) 一旦成立就说明是下降方向。归纳：\(k=0\) 时中间项为 \(-1\)，满足。由 (5.40)，

由强 Wolfe 第二条件 \(|\nabla f_{k+1}^Tp_k|\le-c_2\nabla f_k^Tp_k\)，得

代入归纳假设的左端 \(\frac{\nabla f_k^Tp_k}{\|\nabla f_k\|^2}\ge-\frac1{1-c_2}\)，得 \(-1-\frac{c_2}{1-c_2}\le\cdot\le-1+\frac{c_2}{1-c_2}\)，即 (5.48) 对 \(k+1\) 成立。\(\square\)

证明只用到了强 Wolfe 的第二个条件。(5.48) 限制了 \(\|p_k\|\) 增长的速度，在收敛分析中起关键作用。

FR 的弱点　由 (5.48)，存在常数 \(\chi_1,\chi_2>0\) 使

若某一步 \(p_k\) 很差（\(\cos\theta_k\approx0\)），就意味着 \(\|\nabla f_k\|\ll\|p_k\|\)；步长可能极小，\(x_{k+1}\approx x_k\)，\(\nabla f_{k+1}\approx\nabla f_k\)，于是 \(\beta_{k+1}^{FR}\approx1\)（5.53），\(p_{k+1}\approx p_k\)——新方向几乎没有改进，接下来是一长串无效迭代。

PR 在同样情形下：\(\nabla f_{k+1}\approx\nabla f_k\) 使 \(\beta^{PR}_{k+1}\approx0\)，\(p_{k+1}\) 接近最速下降方向，\(\cos\theta_{k+1}\approx1\)——PR 遇到坏方向后会自动重启。PR+ 和 HS 同理。原书引用 Gilbert–Nocedal 的实例：在 \(n=100\) 的问题上，\(\cos\theta_k\sim10^{-2}\) 持续数百步、步长 \(\sim10^{-2}\)，FR 需要数千步，PR 只要 37 步；FR 若周期性地沿最速下降重启会好得多。结论：FR 不应在没有重启策略的情况下使用。

5.3.6 全局收敛

与线性 CG 不同，非线性 CG 的收敛性质令人意外，有时甚至古怪，理论至今不完整。

假设 5.1　(i) 水平集 \(\mathcal{L}=\{x:f(x)\le f(x_0)\}\) 有界；(ii) 在 \(\mathcal{L}\) 的某个邻域 \(\mathcal{N}\) 内梯度 Lipschitz 连续（5.54）。由此存在 \(\bar\gamma\) 使 \(\|\nabla f(x)\|\le\bar\gamma\)，\(x\in\mathcal{L}\)（5.55）。

定理 5.7（Zoutendijk 定理的重述）　在假设 5.1 下，若每个 \(p_k\) 是下降方向、步长满足 Wolfe 条件，则 \(\sum_{k\ge1}\cos^2\theta_k\|\nabla f_k\|^2<\infty\)（5.56）。

• 周期重启版本：重启步 \(\cos\theta=1\)，所以 \(\sum_{\text{重启步}}\|\nabla f_k\|^2<\infty\)（5.57）。若两次重启间隔不超过 \(\bar n\)，重启步有无穷多，得 \(\liminf\|\nabla f_k\|=0\)（5.58）。这对本章所有算法的重启版都成立。
• 更有意义的是不重启的版本，因为大规模问题（\(n\ge1000\)）通常在远少于 \(n\) 步内收敛，重启从不发生。

定理 5.8（FR 的全局收敛）　设假设 5.1 成立，Algorithm 5.4 的步长满足强 Wolfe 条件且 \(0<c_1<c_2<\frac12\)，则

证明（反证，设 \(\|\nabla f_k\|\ge\gamma>0\) 对所有 \(k\) 成立（5.60））：

1. 由引理 5.6，\(\cos\theta_k\ge\dfrac{1-2c_2}{1-c_2}\dfrac{\|\nabla f_k\|}{\|p_k\|}\)（5.61），代入 Zoutendijk 条件得
   \[\sum_k\frac{\|\nabla f_k\|^4}{\|p_k\|^2}<\infty.\tag{5.62}\]
   （PDF 文本中 (5.61) 的系数印作 \(\frac1{1-c_2}\)；按 (5.48) 的右端应为 \(\frac{1-2c_2}{1-c_2}\)，两者都是正常数，不影响结论。）
2. 由强 Wolfe 条件和引理 5.6，\(|\nabla f_k^Tp_{k-1}|\le-c_2\nabla f_{k-1}^Tp_{k-1}\le\frac{c_2}{1-c_2}\|\nabla f_{k-1}\|^2\)（5.63），于是
   \[\|p_k\|^2\le\|\nabla f_k\|^2+2\beta_k^{FR}|\nabla f_k^Tp_{k-1}|+(\beta_k^{FR})^2\|p_{k-1}\|^2\le c_3\|\nabla f_k\|^2+(\beta_k^{FR})^2\|p_{k-1}\|^2,\quad c_3=\frac{1+c_2}{1-c_2}.\]
   反复递推，并利用 \((\beta_k^{FR})^2(\beta_{k-1}^{FR})^2\cdots(\beta_{k-i}^{FR})^2=\|\nabla f_k\|^4/\|\nabla f_{k-i-1}\|^4\)，得
   \[\|p_k\|^2\le c_3\|\nabla f_k\|^4\sum_{j=1}^k\|\nabla f_j\|^{-2}.\tag{5.64}\]
3. 由 (5.55)(5.60)，\(\|p_k\|^2\le\dfrac{c_3\bar\gamma^4}{\gamma^2}k\)（5.65），所以 \(\sum_k\dfrac1{\|p_k\|^2}\ge\gamma_4\sum_k\dfrac1k=\infty\)（5.66）。
4. 另一方面由 (5.62) 和 (5.60)，\(\gamma^4\sum_k\frac1{\|p_k\|^2}\le\sum_k\frac{\|\nabla f_k\|^4}{\|p_k\|^2}<\infty\)（5.67），矛盾。\(\square\)

这个结果适用于 FR 的实用实现和一般非线性函数（不要求凸），比只对凸函数成立的结论更令人满意。一般地，若存在常数 \(c_4,c_5>0\) 使 \(\cos\theta_k\ge c_4\frac{\|\nabla f_k\|}{\|p_k\|}\) 且 \(\frac{\|\nabla f_k\|}{\|p_k\|}\ge c_5\)，则由定理 5.7 得 \(\lim\|\nabla f_k\|=0\)。对 PR 方法，这可以在 \(f\) 强凸 + 精确线搜索时证明。

但对一般非凸函数，无法对 PR 证明类似定理 5.8 的结论（尽管实践中 PR 更好）：

定理 5.9（Powell 反例）　考虑 PR 方法加理想线搜索（取 \(t(\alpha)=f(x_k+\alpha p_k)\) 的第一个驻点）。存在二阶连续可微的 \(f:\mathbb{R}^3\to\mathbb{R}\) 和初始点，使 \(\|\nabla f_k\|\) 始终远离零。

也就是说 PR 可能无限循环而不接近解。这里假设的步长（第一个驻点）可以被任何实用线搜索接受。证明需要相邻搜索方向几乎互为相反，而在理想线搜索下这只有 \(\beta_k<0\) 时才可能——这启发了 \(\beta_k^+=\max\{\beta_k^{PR},0\}\)（5.68），即 PR+。配合一个保证下降的修改版 Wolfe 线搜索，PR+ 对一般函数全局收敛（Gilbert–Nocedal）。

5.3.7 历史与延伸

• 收敛速度方面（多假设精确线搜索）：Crowder–Wolfe 证明了线性收敛，并构造例子说明不能 Q-超线性；Cohen、Burmeister 证明一般函数上 \(n\) 步二次收敛，Ritter 进一步证明是 \(n\) 步超二次 \(o(\|x_k-x^*\|^2)\)。
• Powell 分析了 FR + 精确线搜索的低效：若迭代进入 \(f=\frac12x^Tx\) 的某个二维区域，梯度与搜索方向的夹角保持不变，若接近 90°，收敛会极慢，甚至比最速下降还慢。
• Nemirovsky–Yudin 指出，在强凸问题上 FR、PR 达不到最优的复杂度界。Nesterov 提出了达到最优界的算法，原书（1999 年）作者认为它"不太可能实用"。后来的发展证明恰恰相反：Nesterov 加速梯度法及其近端版本（FISTA）成了大规模凸优化（如 Lasso 型稀疏组合、带 \(\ell_1\) 惩罚的因子选择）的主力算法之一。这是原书第 1 版成书后的发展，读者阅读时应注意。

5.4 量化实战：大规模因子协方差下的组合求解

场景　风险模型常采用"因子 + 特质"结构

\(B\) 是 \(n\times K\) 的因子暴露矩阵，\(F\) 是 \(K\times K\) 的因子协方差，\(D\) 是对角的特质方差矩阵。全市场选股时 \(n\) 可达数千，\(K\) 通常为 10–50。无约束的均值–方差组合需要求解 \(\Sigma w=\mu\)（最小方差、最大夏普组合的核心计算也都是这个方程）。

两个观察：

1. 矩阵无关：\(\Sigma v=B(F(B^Tv))+Dv\) 只需 \(O(nK)\) 次运算，CG 根本不需要构造 \(n\times n\) 的 \(\Sigma\)（\(n=3000\) 时稠密矩阵就占 72 MB，求解需要 \(O(n^3)\)）。
2. 特征值天然成簇：\(\Sigma\) 有 \(K\) 个大特征值（因子）和大量较小特征值（特质方差，但它们分布在 \([\min d_i,\max d_i]\) 上，未必集中）。用 \(M=D\) 做对角预条件后，
   \[D^{-1/2}\Sigma D^{-1/2}=I+D^{-1/2}BFB^TD^{-1/2}=I+\text{秩 }K\text{ 矩阵},\]
   它的特征值恰好是：\(n-K\) 个等于 1，外加 \(K\) 个大于 1 的特征值——只有 \(K+1\) 个"簇"。由定理 5.4，PCG 至多 \(K+1\) 步就收敛，与 \(n\) 无关，也与条件数无关。

推导拆解：为什么 \(I+U F U^T\)（\(U=D^{-1/2}B\) 是 \(n\times K\)）有 \(n-K\) 个特征值等于 1？任取一个与 \(U\) 的 \(K\) 列都垂直的向量 \(v\)，即 \(U^Tv=0\)，则 \((I+UFU^T)v=v+UF(U^Tv)=v\)，所以 \(v\) 是特征值为 1 的特征向量。与 \(K\) 个向量都垂直的方向构成 \(n-K\) 维子空间，于是特征值 1 至少有 \(n-K\) 重。剩下 \(K\) 个特征值来自 \(U\) 的列空间，等于 1 加上 \(FU^TU\) 的特征值，都大于 1。 金融含义：\(v\) 是"对所有因子暴露为零"的组合方向（按特质波动率标准化后）。这类组合只有特质风险，标准化后每个方向的风险都是 1，所以它们在 PCG 眼里"一模一样"，一步就能统一处理。真正需要逐个处理的只有 \(K\) 个因子方向。 实测迭代次数为 12，比理论的 \(K+1=11\) 多一步，原因是浮点舍入使共轭性不再严格成立，且停止准则要求相对残差 \(\le10^{-8}\)。 我们还顺便验证定理 5.4：只有 \(r\) 个不同特征值时，CG 恰好 \(r\) 步收敛。

import numpy as np, time
rng = np.random.default_rng(5)

def pcg(matvec, b, Minv=None, tol=1e-8, maxit=5000):
    """Algorithm 5.3（预条件 CG）；Minv=None 时即 Algorithm 5.2。返回解与迭代次数"""
    x = np.zeros_like(b); r = -b.copy()            # r = A x - b
    y = r if Minv is None else Minv(r); p = -y; ry = r @ y
    for k in range(maxit):
        if np.linalg.norm(r) <= tol * np.linalg.norm(b): return x, k
        Ap = matvec(p); alpha = ry / (p @ Ap)
        x += alpha * p; r += alpha * Ap
        y = r if Minv is None else Minv(r)
        ry_new = r @ y; p = -y + (ry_new / ry) * p; ry = ry_new
    return x, maxit

# ---------- 1) 定理 5.4：只有 r 个不同特征值 => 至多 r 步 ----------
n = 500
Qm, _ = np.linalg.qr(rng.standard_normal((n, n)))
for r_ in [3, 7, 20]:
    lam = rng.choice(np.linspace(1, 100, r_), size=n)          # 只有 r_ 个不同特征值
    A = (Qm * lam) @ Qm.T
    _, k = pcg(lambda v: A @ v, rng.standard_normal(n), tol=1e-10)
    print("不同特征值个数 %2d：CG 迭代 %d 次" % (r_, k))

# ---------- 2) 因子结构协方差 Σ = B F B' + D：矩阵无关（matrix-free）的 CG ----------
n, K = 3000, 10
Bx = rng.standard_normal((n, K)) * np.r_[1.0, 0.5 * np.ones(K - 1)] * 0.8   # 因子暴露
Fc = np.diag(np.r_[0.04, 0.01 * np.ones(K - 1)])                           # 因子协方差（年化）
d = rng.uniform(0.02, 0.25, n)**2                                          # 特质方差，波动率 2%~25%
mu = rng.normal(0.05, 0.03, n)
matvec = lambda v: Bx @ (Fc @ (Bx.T @ v)) + d * v                          # O(nK) 次运算

Sigma = Bx @ Fc @ Bx.T + np.diag(d)                                         # 仅用于对照
ev = np.linalg.eigvalsh(Sigma)
print("\nn=%d, K=%d：κ(Σ)=%.0f；最大 %d 个特征值 %s" % (n, K, ev[-1] / ev[0], K, np.round(ev[-K:][::-1], 1)))
Dh = 1 / np.sqrt(d)
ev_p = np.linalg.eigvalsh(Dh[:, None] * Sigma * Dh[None, :])
print("预条件后 D^{-1/2} Σ D^{-1/2}：特征值等于 1 的个数 %d，其余 %d 个范围 [%.1f, %.1f]，κ=%.0f" %
      ((np.abs(ev_p - 1) < 1e-8).sum(), (np.abs(ev_p - 1) >= 1e-8).sum(),
       ev_p[np.abs(ev_p - 1) >= 1e-8].min(), ev_p.max(), ev_p.max() / ev_p.min()))

t = time.perf_counter(); w_cg, k1 = pcg(matvec, mu); t1 = time.perf_counter() - t
t = time.perf_counter(); w_pcg, k2 = pcg(matvec, mu, Minv=lambda r: r / d); t2 = time.perf_counter() - t
t = time.perf_counter(); w_chol = np.linalg.solve(Sigma, mu); t3 = time.perf_counter() - t
print("CG（无预条件）   : %4d 次迭代, %6.1f ms" % (k1, 1e3 * t1))
print("PCG（M = D）     : %4d 次迭代, %6.1f ms" % (k2, 1e3 * t2))
print("稠密直接法 solve : %6.1f ms（还不含构造 n×n 矩阵的时间与 %.0f MB 内存）" % (1e3 * t3, Sigma.nbytes / 1e6))
print("与直接法的相对误差: CG %.1e, PCG %.1e" %
      (np.linalg.norm(w_cg - w_chol) / np.linalg.norm(w_chol), np.linalg.norm(w_pcg - w_chol) / np.linalg.norm(w_chol)))
print("理论界 (5.35) 需要的迭代数（无预条件，误差 1e-8）≈ %.0f" %
      (np.log(2e8) / -np.log((np.sqrt(ev[-1] / ev[0]) - 1) / (np.sqrt(ev[-1] / ev[0]) + 1))))

运行输出（计时因机器而异）：

不同特征值个数  3：CG 迭代 3 次
不同特征值个数  7：CG 迭代 7 次
不同特征值个数 20：CG 迭代 20 次

n=3000, K=10：κ(Σ)=180637；最大 10 个特征值 [73.4  5.2  5.1  4.9  4.8  4.7  4.6  4.6  4.5  4.3]
预条件后 D^{-1/2} Σ D^{-1/2}：特征值等于 1 的个数 2990，其余 10 个范围 [760.8, 13522.3]，κ=13522
CG（无预条件）   :  254 次迭代,    6.2 ms
PCG（M = D）     :   12 次迭代,    0.3 ms
稠密直接法 solve :   61.7 ms（还不含构造 n×n 矩阵的时间与 72 MB 内存）
与直接法的相对误差: CG 9.3e-09, PCG 8.2e-14
理论界 (5.35) 需要的迭代数（无预条件，误差 1e-8）≈ 4062

解读

• 定理 5.4 得到精确验证：3、7、20 个不同特征值分别对应 3、7、20 步。
• 协方差矩阵条件数约 \(1.8\times10^5\)。条件数界 (5.35) 预测无预条件 CG 需要约 4000 步，实际只用了 254 步——条件数界只用到最大、最小两个特征值，忽略了分布信息：这里最大的那个特征值（市场因子）是孤立的，按定理 5.5 的思路，CG 花一两步就能"消掉"它，此后有效条件数要小得多。
• 对角预条件后条件数仍高达 13522，但特征值只有 11 个簇（2990 个精确等于 1，另 10 个分散），PCG 只用 12 步就达到 \(10^{-13}\) 的精度。这说明：决定 CG 速度的是特征值分布，而不是条件数；好的预条件子未必要降低条件数，让特征值成簇同样有效。
• 计算量对比：PCG 每步 \(O(nK)\)，总量约 \(12\times2\times3000\times10\) 次乘法，比稠密直接法的 \(O(n^3)\) 少几个数量级，内存也只需 \(O(nK)\)。

实务讨论

• 对这个特定的无约束问题，用 Woodbury 恒等式 \((D+BFB^T)^{-1}=D^{-1}-D^{-1}B(F^{-1}+B^TD^{-1}B)^{-1}B^TD^{-1}\) 也能 \(O(nK^2)\) 地精确求解。CG 的优势在于通用性：一旦目标中加入其他结构（例如交易成本的二次项 \(\frac12(w-w_0)^T\Lambda(w-w_0)\)、多期优化的块结构、或者作为第 6 章截断牛顿法的内层求解器），Woodbury 不再直接适用，而 CG 只需要你能计算矩阵–向量乘积。
• 小规模问题（几十到几百个资产）直接用 Cholesky 分解更稳健，这与原书"CG 只推荐用于大规模问题"的建议一致。
• 非线性 CG 在量化中的典型用法是内存受限的大规模校准或机器学习；SciPy 的 minimize(method='CG') 实现的是 PR 型非线性 CG。不过现在 L-BFGS（本册第 9 章）通常是更好的默认选择。

本章小结

线性 CG 是求解对称正定方程组（等价于极小化凸二次函数）的迭代法：它是一种共轭方向法，每一步只用前一个方向就自动保持与所有旧方向共轭，残差两两正交，迭代在逐步扩大的 Krylov 子空间上使 \(A\)-范数误差最小。收敛速度由特征值分布决定——\(r\) 个不同特征值 \(r\) 步终止，特征值成簇时几步即可，最坏情况由 \(\sqrt\kappa\) 控制而非最速下降的 \(\kappa\)。预条件通过改善特征值分布来加速，它的设计依赖于对问题结构的理解。非线性 CG 用线搜索和梯度替换线性 CG 的步长与残差：FR 需要强 Wolfe 且 \(c_2<\frac12\) 保证下降和 \(\liminf\|\nabla f_k\|=0\)，但遇到坏方向会陷入低效；PR 会自动重启，实践更好，但非凸时可能不收敛，PR+ 修正了这一点。

练习

基础

1. （原书 5.2）证明关于对称正定矩阵 \(A\) 共轭的非零向量组线性无关。 提示：若 \(\sum\sigma_ip_i=0\)，左乘 \(p_j^TA\)。
2. （原书 5.3）验证 (5.6) 是沿 \(p_k\) 精确极小化 \(\phi\) 的步长。
3. （原书 5.6）证明 (5.23d) 与 (5.13d) 等价。
4. 用 CG 手算求解 \(\begin{bmatrix}4&1\\1&3\end{bmatrix}x=\begin{bmatrix}1\\2\end{bmatrix}\)，\(x_0=0\)，验证两步得到精确解，并验证 \(r_1^Tr_0=0\)、\(p_1^TAp_0=0\)。 提示：解为 \(x^*=(1/11,7/11)\)。
5. （原书 5.11）证明：对二次函数和精确线搜索，PR 和 HS 公式都退化为 FR 公式。

进阶

6. （原书 5.1）实现 CG 求解 Hilbert 矩阵 \(A_{ij}=1/(i+j-1)\) 的方程组，\(b=(1,\dots,1)^T\)，\(x_0=0\)，\(n=5,8,12,20\)，报告残差降到 \(10^{-6}\) 所需的迭代次数，并与理论"至多 \(n\) 步"对比，解释差异。 提示：Hilbert 矩阵极度病态，舍入误差破坏了共轭性，迭代次数可能超过 \(n\)。
7. （原书 5.9、5.10）从变量变换 (5.36)(5.37) 出发推导预条件 CG（Algorithm 5.3），并验证 (5.39)。
8. （原书 5.8 改编）在本章实战代码中，把特质方差改为全部相等（\(d_i\equiv0.01\)），不加预条件时 CG 需要几步？解释原因。再把因子数 \(K\) 从 10 改为 50，观察 PCG 迭代次数的变化。 提示：此时 \(\Sigma=0.01I+\) 秩 \(K\) 矩阵，至多 \(K+1\) 个不同特征值。
9. （原书 5.12）证明引理 5.6 对任何满足 \(|\beta_k|\le\beta_k^{FR}\) 的方法都成立。由此说明为什么把 PR 的 \(\beta\) 截断到 \([-\beta^{FR},\beta^{FR}]\) 可以继承 FR 的下降性。
10. 用 scipy.optimize.minimize 分别以 method='CG' 和 method='L-BFGS-B' 求解扩展 Rosenbrock 函数（\(n=1000\)），比较函数–梯度求值次数。

原书推荐习题：5.1（Hilbert 矩阵上的 CG，体会病态与舍入误差）、5.8（构造不同特征值分布验证聚类效应，可直接用因子协方差做实验）、5.9 与 5.10（推导预条件 CG）、5.11（FR/PR/HS 的联系）、5.12（引理 5.6 的推广）。

原书对照

页码换算：原书正文页码 = PDF 页码减 20（已用 PDF 页眉核对；第 1–2 章减 21，第 3–11 章减 20，第 12–16 章减 19，第 17 章至附录减 18）。精读笔记中"预条件"两小节标注的页码（PDF p.118–119、119–120）有误，已对照 PDF 更正为 Preconditioning 在 PDF p.138–139、Practical Preconditioners 在 PDF p.139–140。