第 06 章实用牛顿法

学习目标

读完本章，你应当能够：

理解非精确牛顿法：用残差相对大小 $\|r_k\|\le\eta_k\|\nabla f_k\|$ 控制内层求解精度，知道强制序列 $\eta_k$ 如何决定线性、超线性、二次收敛。
实现线搜索 Newton–CG（截断牛顿法），包括 CG 起点为零、负曲率检验和 Hessian-free 的 Hessian–向量积。
说明修正牛顿法的全局收敛条件（有界修正分解性质），比较特征值修正、加单位阵倍数、修正 Cholesky、修正对称不定分解等 Hessian 修正策略。
了解信赖域牛顿法的几种实现（dogleg、二维子空间、精确解、信赖域 Newton–CG）各自的优劣，理解预条件需要与椭球信赖域配合的原因。
理解为什么信赖域牛顿法在解附近会退化为纯牛顿法，从而保持快速局部收敛。
用本章方法修复非半正定的相关系数矩阵，并用 Hessian-free Newton–CG 求解数千资产的风险预算组合。

读前导读

这一章在解决什么问题

前几章已经说明，牛顿法（同时用斜率和弯曲程度）在最优点附近快得惊人，Excel Solver 的"牛顿"选项、计量软件做极大似然估计，最后几步都是这样"瞬间"收敛的。但纯牛顿法有两个实际麻烦，本章逐一解决。

第一个麻烦：Hessian 不正定。牛顿法的前提是"局部看起来像一个碗"。离最优点较远时，目标函数在某些方向可能向下弯（像马鞍），这时牛顿步可能指向上坡，或者长得离谱。这和风控里一个你很熟悉的问题是同一件事：人工设定的压力情景相关矩阵、用成对缺失数据估出来的相关矩阵，常常有负特征值，做不了 Cholesky 分解，也就生成不了相关随机数。两者的修复办法如出一辙：把矩阵"修正"成正定的，同时尽量少改动原有信息。本章介绍几种修正方法（特征值截断、整体加 $\tau I$、修正 Cholesky、对称不定分解），量化实战里直接拿它们修复相关矩阵。

第二个麻烦：规模。每一步牛顿法要解一个 $n\times n$ 的线性方程组。$n$ 是几千时，构造并分解 Hessian 太贵。解决办法是用第 5 章的 CG 近似地解，并且"离最优点远时粗算、近时细算"（非精确牛顿）。更进一步，CG 只需要"Hessian 乘向量"，而这可以不构造 Hessian 直接算出来（Hessian-free），就像用"上下各 bump 1bp 重新定价"算有效久期，而不必推导解析公式。

最后，本章把这些技术放进线搜索和信赖域两个框架，并证明它们在最优点附近都会退化回纯牛顿法，从而保住快速收敛。实战部分用这套方法求解了 2000 只股票的风险预算（风险平价的推广）组合。

需要先想起来的数学

正定、特征值与 Cholesky 分解：对称矩阵正定 ⇔ 所有特征值 $>0$ ⇔ Cholesky 分解 $A=LL^T$ 能顺利完成。所以"尝试做 Cholesky，失败了就说明不正定"是实用的正定检验。见第 00 册第 06 章线性代数速成。
矩阵范数与条件数：$\|A\|$ 衡量矩阵"最多能把向量放大多少倍"，对称矩阵的 2-范数就是最大的 $|\lambda_i|$；Frobenius 范数 $\|A\|_F=\sqrt{\sum_{ij}a_{ij}^2}$ 是把所有元素当成一个长向量算长度。条件数 $\|A\|\|A^{-1}\|$ 是最大与最小 $|\lambda|$ 之比，衡量"最敏感方向与最不敏感方向"差多少。见第 00 册第 06 章。
有限差分：$f'(x)\approx[f(x+h)-f(x)]/h$，误差 $O(h)$。多元时"梯度在方向 $p$ 上的变化率"就是 Hessian 乘 $p$。这和用 bump-and-reprice 估久期、凸性、Greeks 完全同理。见第 00 册第 02 章导数与泰勒展开。
$\limsup$、$O(\cdot)$ 与 $o(\cdot)$：$\limsup a_k\le\eta$ 读作"$a_k$ 最终不会持续超过 $\eta$"（上极限）。$O(\|g\|^2)$ 表示"与 $\|g\|^2$ 同阶或更小"，$o(\|g\|)$ 表示"比 $\|g\|$ 小得多"。见第 00 册第 07 章概率中的分析工具。
对数函数的导数：$(\log w)'=1/w$，$(1/w)'=-1/w^2$。风险预算问题的梯度与 Hessian 只用到这两个。见第 00 册第 02 章。

怎么读这一章

必读：6.1、6.2（定理 6.1、6.2 的结论和强制序列的常用选择）、6.3.1（Newton–CG 的三条规则与 Hessian-free）、6.3.2（修正牛顿法与"有界修正分解"）、6.4.1（特征值修正的例子）、6.6 量化实战。6.4.2–6.4.5 的各种分解算法属于实现细节，第一次读每种方法记住"做什么、优缺点"即可，Algorithm 6.5 的逐行细节可以跳过。6.5 信赖域牛顿法读 6.5.3 的优点与局限；6.5.4 预条件和 6.5.5 局部收敛证明可以留待以后。

6.1 引言：让牛顿法既快又稳

纯牛顿法（单位步长）在接近 $x^*$ 时收敛极快，但从远处出发可能不收敛，在非凸区域的行为也不稳定。本章的目标是让牛顿法在所有情形下都稳健，同时保持效率。

牛顿步是对称线性方程组

\nabla^2f(x_k)p_k^N=-\nabla f(x_k)\tag{6.1}

的解。Hessian 正定时 $p_k^N$ 是下降方向；不正定或接近奇异时，它可能是上升方向，或者长得离谱。保证步长质量有两大策略：

用 CG 解 (6.1)，一旦遇到负曲率就停止——Newton–CG，有线搜索和信赖域两种实现；
在解 (6.1) 之前或过程中修改 Hessian，使它充分正定——修正牛顿法（modified Newton）。

控制计算量也有两条路：Newton–CG 在得到精确解之前就停止 CG，即非精确牛顿（inexact Newton）；直接法则利用 Hessian 的稀疏结构做稀疏高斯消元。

Hessian 的计算通常是主要工作量，没有解析式时可用自动微分或有限差分（第 7 章）。把修正技术、稀疏性和微分技术结合起来，本章的牛顿法是求解中小规模乃至大规模无约束问题最可靠、最强大的方法之一。

6.2 非精确牛顿步

大规模问题中，直接分解 Hessian 的代价很高；而且远离解时二次模型本来就不准，精确求解 (6.1) 并不值得。因此用迭代法求近似解。终止规则基于残差

r_k=\nabla^2f(x_k)p_k+\nabla f(x_k).\tag{6.2}

残差的绝对大小会随目标函数的数乘缩放而变，所以要用相对于右端的大小：

\|r_k\|\le\eta_k\|\nabla f(x_k)\|,\tag{6.3}

其中 $\{\eta_k\}$（$0<\eta_k<1$）称为强制序列（forcing sequence）。

定理 6.1　设 $\nabla f$ 在极小点 $x^*$ 的邻域内连续可微，$\nabla^2f(x^*)$ 正定。考虑迭代 $x_{k+1}=x_k+p_k$，$p_k$ 满足 (6.3)，且 $\eta_k\le\eta$ 对某个常数 $\eta\in[0,1)$ 成立。若 $x_0$ 足够接近 $x^*$，则 $\{x_k\}$ 线性收敛：$\|x_{k+1}-x^*\|\le c\|x_k-x^*\|$，$0<c<1$。

这个条件并不苛刻：如果允许 $\eta_k\ge1$，$p_k=0$ 就总满足 (6.3)，迭代会停滞。

非正式推导：在 $x^*$ 附近 $\|\nabla^2f(x_k)^{-1}\|\le L$，于是 $\|p_k\|\le L(\|\nabla f_k\|+\|r_k\|)\le2L\|\nabla f_k\|$。由 Taylor 定理，

\nabla f(x_{k+1})=\nabla f_k+\nabla^2f_kp_k+O(\|p_k\|^2)=r_k+O(\|\nabla f_k\|^2),\tag{6.4}

\|\nabla f(x_{k+1})\|\le\eta_k\|\nabla f_k\|+O(\|\nabla f_k\|^2).\tag{6.5}

所以梯度每步大约缩小为 $\eta_k$ 倍，$\limsup\frac{\|\nabla f_{k+1}\|}{\|\nabla f_k\|}\le\eta<1$。进一步：若 $r_k=o(\|\nabla f_k\|)$，就是超线性；若 $r_k=O(\|\nabla f_k\|^2)$，则 $\limsup\frac{\|\nabla f_{k+1}\|}{\|\nabla f_k\|^2}$ 有界，恢复二次收敛。迭代点与梯度以相同的速度收敛。

定理 6.2　在定理 6.1 的条件下，若 $x_k\to x^*$，则：$\eta_k\to0$ 时超线性收敛；$\eta_k=O(\|\nabla f(x_k)\|)$ 时二次收敛。

白话解释：(6.5) 把下一步的梯度拆成两部分：$\eta_k\|\nabla f_k\|$ 来自"内层方程没解准"，$O(\|\nabla f_k\|^2)$ 来自"二次模型本身不准"。后者天然是平方级的（牛顿法的优势），所以整体速度由前者决定。数值例：设当前 $\|\nabla f_k\|=10^{-2}$。取常数 $\eta=0.5$：下一步梯度约 $0.5\times10^{-2}$，只是减半——线性收敛，浪费了牛顿法的潜力。取 $\eta_k=\sqrt{\|\nabla f_k\|}=0.1$：下一步约 $10^{-3}$；梯度越小，$\eta_k$ 也越小——超线性。取 $\eta_k=\|\nabla f_k\|=10^{-2}$：下一步约 $10^{-4}$，与二次项同阶——恢复二次收敛。代价是 $\eta_k$ 越小，内层 CG 要迭代越多次。这就是"远处粗算、近处细算"。

常用选择：$\eta_k=\min(0.5,\sqrt{\|\nabla f_k\|})$ 得到超线性收敛；$\eta_k=\min(0.5,\|\nabla f_k\|)$ 得到二次收敛。这些结果（Dembo–Eisenstat–Steihaug）都是局部的，假设迭代已进入解附近并取单位步长。下面几节说明非精确牛顿法可以嵌入实用的线搜索和信赖域框架，而全局化策略不会妨碍最终的快速收敛。

量化启发：强制序列体现了一种"计算经济学"——远离最优时不必精确求解子问题，越接近最优越要精确。日内再平衡等有时间预算的实时优化可以照此设计：前几次迭代粗解保证下降，接近最优时再提高内层精度。

6.3 线搜索牛顿法

迭代 $x_{k+1}=x_k+\alpha_kp_k$，$p_k$ 是牛顿方向或其近似，$\alpha_k$ 满足 Wolfe、Goldstein 或 Armijo 回溯条件，并且总是先试单位步长。

6.3.1 线搜索 Newton–CG

又称截断牛顿法（truncated Newton）：用 CG 解 (6.1)，尝试满足 (6.3)。问题在于 CG 是为正定系统设计的，而远离解时 Hessian 可能有负特征值。对策是：一旦 CG 生成负曲率方向就终止。这样既保证 $p_k$ 是下降方向，又保留了牛顿法的快速收敛。内层 CG（$A=\nabla^2f_k$，$b=-\nabla f_k$，上标 $(i)$ 表示 CG 内部迭代）有三个要求：

(a) CG 的起点 $x^{(0)}=0$；
(b) 负曲率检验：若 CG 的方向满足 $(p^{(i)})^TAp^{(i)}\le0,\tag{6.6}$ 那么：如果这是第一次 CG 迭代（$i=0$），完成这一步（得到 $x^{(1)}$，即沿 $-\nabla f_k$ 的方向）后停止；如果 $i>0$，立即停止并返回 $x^{(i)}$；
(c) 牛顿步 $p_k$ 取最终的 CG 迭代点。

为什么起点取 0？这样如果第一步就遇到负曲率，$p_k$ 就是最速下降方向 $-\nabla f_k$；若 CG 走了不止一步，$p_k$ 也一定是下降方向（练习 2）。CG 中也可以使用预条件。

白话解释：内层 CG 在做的事是"极小化二次模型 $q(p)=\nabla f_k^Tp+\tfrac12p^T\nabla^2f_kp$"。$(p^{(i)})^TAp^{(i)}\le0$ 说明沿这个 CG 方向，模型不是向上弯的碗，而是平的或向下弯的——模型在这个方向上没有最低点，继续按 CG 公式算步长会除以一个非正数，结果没有意义。所以此时停下，返回到目前为止的结果。三条规则合起来保证了一件事：返回的 $p_k$ 至少和最速下降一样"往下走"，并且在 Hessian 正定、CG 跑足时就是牛顿步。这是"稳健"和"快"的折中：最坏退化为最速下降，最好就是牛顿法。 Hessian-free　Newton–CG 不需要显式的 Hessian，只需要 Hessian–向量积 $\nabla^2f(x_k)p$。可以用有限差分

\nabla^2f(x_k)p\approx\frac{\nabla f(x_k+hp)-\nabla f(x_k)}{h},\tag{6.7}

精度 $O(h)$，每次 CG 迭代多算一次梯度；也可以用自动微分精确计算（第 7 章）。这类方法统称 Hessian-free 牛顿法，在深度学习和大规模统计估计中很常用。

金融直觉：(6.7) 就是 bump-and-reprice。算有效凸性时，你把收益率上下 bump 一下，看久期（价格对收益率的一阶导）变了多少；这里是把位置沿方向 $p$ bump 一个小量 $h$，看梯度变了多少，变化率就是 Hessian 在 $p$ 方向上的作用 $\nabla^2f\,p$。好处是只需要一次额外的梯度计算，不必写出 $n\times n$ 的 Hessian。 $h$ 的选择和 bump 大小的取舍一样：太大，截断误差 $O(h)$ 大；太小，两个几乎相等的梯度相减，浮点舍入误差被放大。第 7 章会讨论如何选 $h$。 Algorithm 6.1（线搜索 Newton–CG）

给定 x_0
for k = 0,1,2,...
  对 ∇²f(x_k) p = −∇f_k 从 x^(0)=0 开始做 CG；
  当 ‖r‖ ≤ min(0.5, √‖∇f_k‖)·‖∇f_k‖ 或遇到负曲率（按规则 (b)）时终止
  x_{k+1} = x_k + α_k p_k，α_k 满足 Wolfe、Goldstein 或 Armijo 回溯（先试 α=1）
end

（原书正文中强制序列有一处印作 $\min(0.5,\sqrt{\|\nabla^2f_k\|})$，算法框中为 $\sqrt{\|\nabla f_k\|}$，后者正确。）

弱点（尤其无预条件时）：Hessian 接近奇异时，Newton–CG 方向可能非常长，线搜索要做很多次函数求值，而下降量却很小。规范化牛顿步的好规则很难定（可能破坏良好尺度问题上的快速收敛）；在 (6.6) 中引入阈值更好，但阈值也难选。6.5 节的信赖域 Newton–CG 处理得更好，原书作者略偏好后者。

6.3.2 修正牛顿法

如果想用直接法（如 Cholesky 分解）解 (6.1)，可以在 Hessian 不正定或接近奇异时，在求解前或求解过程中修改它（加一个正对角阵或满矩阵），得到一个正定近似。

Algorithm 6.2（带修正的线搜索牛顿法）

给定 x_0
for k = 0,1,2,...
  分解 B_k = ∇²f(x_k) + E_k：∇²f(x_k) 充分正定时 E_k = 0，否则选 E_k 使 B_k 充分正定
  解 B_k p_k = −∇f(x_k)
  x_{k+1} = x_k + α_k p_k，α_k 满足 Wolfe、Goldstein 或 Armijo 回溯
end

$E_k$ 的选择至关重要；有些方法并不显式计算 $E_k$，而是在标准分解过程中"边分解边修改"。

有界修正分解性质（bounded modified factorization property）：只要 Hessian 序列有界，修正后的矩阵条件数就一致有界：

\mathrm{cond}(B_k)=\|B_k\|\|B_k^{-1}\|\le C.\tag{6.8}

定理 6.3　设 $f$ 在开集 $\mathcal{D}$ 上二阶连续可微，水平集 $\mathcal{L}=\{x\in\mathcal{D}:f(x)\le f(x_0)\}$ 紧，且有界修正分解性质成立。则 Algorithm 6.2 满足 $\lim\nabla f(x_k)=0$。

证明：线搜索保证迭代点留在水平集内；Hessian 连续、水平集紧 ⇒ Hessian 有界 ⇒ (6.8) 成立；由 (3.19)，$\cos\theta_k\ge1/C$，再由 Zoutendijk 定理即得。$\square$

白话解释：为什么只要求"正定"还不够，还要"条件数有界"？方向 $p_k=-B_k^{-1}\nabla f_k$ 与最陡下坡方向的夹角取决于 $B_k$ 的条件数：条件数越大，$B_k^{-1}$ 在不同方向上的放大倍数差别越大，$p_k$ 就越可能被拉得几乎与梯度垂直（$\cos\theta_k\to0$），沿它走几乎不下降。条件数有上界 $C$，就保证 $\cos\theta_k\ge1/C$，方向始终"有一定的下坡分量"，Zoutendijk 定理（第 3 章）随即给出 $\nabla f_k\to0$。下面 6.4.1 的例子 (6.11) 就是反面教材：修正后的矩阵正定，但条件数达 $10^9$，方向几乎完全被一个分量主导。 收敛速度：若 $x_k\to x^*$ 且 $\nabla^2f(x^*)$ 充分正定，使得修正最终为零，那么由定理 3.5 最终 $\alpha_k=1$，算法退化为纯牛顿法，二次收敛。若 $\nabla^2f(x^*)$ 接近奇异，修正可能永远不消失，只能线性收敛。

量化提示：GARCH、Copula、随机波动率模型的极大似然估计中，远离最优点时 Hessian 常常不定，线搜索牛顿法必须修正 Hessian，否则可能沿上升方向走。如果在最优点附近修正仍不消失、收敛变慢，往往说明 Hessian 接近奇异——参数弱识别（例如 GARCH 中 $\alpha+\beta$ 接近 1 时 $\omega$ 难以识别）。这时用 Hessian 的逆估计标准误也不可靠，应先检查模型设定。

6.4 Hessian 修正

目标：$B_k=\nabla^2f(x_k)+E_k$ 充分正定且条件良好（使定理 6.3 成立）；修正尽可能小，以保留二阶信息；分解代价适中。先讲基于特征分解的"理想"策略，再讲实用的分解型方法。

6.4.1 特征值修正

例　设 $\nabla f(x_k)=(1,-3,2)^T$，$\nabla^2f(x_k)=\mathrm{diag}(10,3,-1)$（不定），即

\nabla^2f(x_k)=Q\Lambda Q^T=\sum_{i=1}^n\lambda_iq_iq_i^T,\qquad Q=I.\tag{6.9}

纯牛顿步 $p_k^N=(-0.1,1,2)^T$，而 $\nabla f_k^Tp_k^N=0.9>0$——不是下降方向。

把负特征值换成一个小正数 $\delta=\sqrt u$（$u\approx10^{-16}$ 为机器精度）： $B_k=\sum_{i=1}^2\lambda_iq_iq_i^T+\delta q_3q_3^T=\mathrm{diag}(10,3,10^{-8}),\tag{6.10}$ 它保留了 $q_1,q_2$ 方向的曲率信息，但 $p_k=-B_k^{-1}\nabla f_k=-\sum_{i=1}^2\frac1{\lambda_i}q_i(q_i^T\nabla f_k)-\frac1\delta q_3(q_3^T\nabla f_k)\approx-(2\times10^8)q_3,\tag{6.11}$ 几乎平行于 $q_3$ 且极长。虽然 $f$ 沿这个方向下降，但如此极端的长度违背了牛顿法依赖局部二次模型的精神。
其他选择：翻转负特征值的符号（相当于 $\delta=|\lambda_i|$）；令 (6.11) 的最后一项为零（完全不用负曲率分量）；按"步长不能过长"自适应地选 $\delta$（这已有信赖域的味道）。哪种最理想，尚无共识。

两个"最优修正"的概念：

最小 Frobenius 范数修正：对称矩阵 $A=Q\Lambda Q^T$，使 $\lambda_{\min}(A+\Delta A)\ge\delta$ 且 $\|\Delta A\|_F$ 最小的修正是 $\Delta A=Q\,\mathrm{diag}(\tau_i)Q^T,\qquad\tau_i=\begin{cases}0,&\lambda_i\ge\delta,\\\delta-\lambda_i,&\lambda_i<\delta.\end{cases}\tag{6.12}$ （$\|A\|_F^2=\sum a_{ij}^2$。）它一般不是对角阵，修正后 $A+\Delta A=Q(\Lambda+\mathrm{diag}(\tau_i))Q^T$。(6.10) 就是这种修正。
最小 2-范数修正： $\Delta A=\tau I,\qquad\tau=\max(0,\delta-\lambda_{\min}(A)),\tag{6.13}$ 修正后为 $A+\tau I$（6.14），所有特征值整体平移到 $\ge\delta$。注意它与信赖域子问题中的 $B+\lambda I$ 形式相同。

实用软件不会做特征分解（太贵），而是用高斯消元间接选择修正；数值经验表明这常常（但不总是）给出好的方向。

金融直觉：两种"最优修正"在风险管理里都有对应做法。(6.12) 就是相关矩阵修复中最常用的特征值截断：做 PCA，把负的（或过小的）主成分方差抬到一个小正数，其他主成分原封不动，因此对原矩阵的改动（按所有元素的平方和计）最小。(6.13) 相当于给所有资产的方差统一加一个数，类似协方差收缩估计里向单位阵方向收缩：结构简单，但会"误伤"本来就健康的主成分。 (6.10)(6.11) 的教训也很实际：把一个负特征值改成 $10^{-8}$ 虽然"修好了"，但等于宣称"存在一个几乎零风险的组合"，任何优化器都会在这个方向上加杠杆到 $10^8$ 量级。在组合优化里，这就是用"勉强正定"的协方差矩阵做最小方差组合时出现极端权重的原因。

6.4.2 加单位阵倍数

找 $\tau>0$ 使 $\nabla^2f(x_k)+\tau I$ 充分正定。由 (6.13) 需要知道最小特征值，但通常没有好的估计。可以利用"最大绝对特征值不超过 $\|A\|_F$"来试探：

Algorithm 6.3（加单位阵倍数的 Cholesky）

β ← ‖A‖_F
if min_i a_ii > 0:  τ_0 ← 0   else  τ_0 ← β/2
for k = 0,1,2,...
  尝试 Cholesky 分解 L Lᵀ = A + τ_k I
  if 成功: 返回 L
  else: τ_{k+1} ← max(2τ_k, β/2)
end

（原文此处写作 "incomplete Cholesky"，在这里应理解为对 $A+\tau I$ 做普通 Cholesky 分解。）方法简单，可能比后面的修正分解更可取；缺点是 $\tau$ 可能不必要地大，使方向过度偏向最速下降；而且每个 $\tau_k$ 都要重做一次数值分解（符号分解只需做一次），试很多次时代价高。

6.4.3 修正 Cholesky 分解

思路：对 Hessian 做 Cholesky 分解，在过程中必要时增大对角元。两个目标：保证修正因子存在且其大小相对于 Hessian 的范数有界；Hessian 充分正定时不做任何修改。

先回顾 $LDL^T$ 形式：对称正定矩阵 $A$ 可写成

A=LDL^T,\tag{6.15}

$L$ 单位下三角，$D$ 对角且元素为正。例 6.1（$n=3$）：比较两边各列元素可得 $d_1=a_{11}$，$l_{21}=a_{21}/d_1$，$l_{31}=a_{31}/d_1$；$d_2=a_{22}-d_1l_{21}^2$，$l_{32}=(a_{32}-d_1l_{31}l_{21})/d_2$；$d_3=a_{33}-d_1l_{31}^2-d_2l_{32}^2$。一般算法如下：

Algorithm 6.4（Cholesky，$LDL^T$ 形式）

for j = 1..n
  c_jj ← a_jj − Σ_{s<j} d_s l_js²；  d_j ← c_jj
  for i = j+1..n
    c_ij ← a_ij − Σ_{s<j} d_s l_is l_js；  l_ij ← c_ij / d_j
end

$A$ 正定时所有 $d_j>0$。与标准形式 $A=MM^T$（6.16）的关系是 $M=LD^{1/2}$。$A$ 不定时 $LDL^T$ 可能不存在；即使存在也数值不稳定（元素可以任意大），所以"先分解、再修改对角元"的做法可能失败，或得到与 $A$ 相差极大的矩阵。

改为在分解过程中修改：选参数 $\delta,\beta>0$，要求计算第 $j$ 列时

d_j\ge\delta,\qquad|m_{ij}|\le\beta,\quad i=j+1,\dots,n,\tag{6.17}

其中 $m_{ij}=l_{ij}\sqrt{d_j}$ 是 $M$ 的元素。只需把 $d_j$ 的计算改为

d_j=\max\left(|c_{jj}|,\ \left(\frac{\theta_j}{\beta}\right)^2,\ \delta\right),\qquad\theta_j=\max_{j<i\le n}|c_{ij}|.\tag{6.18}

验证：$c_{ij}=l_{ij}d_j$，所以 $|m_{ij}|=|l_{ij}|\sqrt{d_j}=\dfrac{|c_{ij}|}{\sqrt{d_j}}\le\dfrac{|c_{ij}|\beta}{\theta_j}\le\beta$。由于 $c_{ij}$ 不依赖 $d_j$，$\theta_j$ 可以先于 $d_j$ 算出——这正是引入中间量 $c_{ij}$ 的原因。为了减小修正量，还要做对称的行列交换，使第 $j$ 步的主元是剩余对角元中绝对值最大者。

白话解释：(6.18) 里三个候选值各管一件事。$|c_{jj}|$：主元本来是多少就用多少，若为负则取绝对值（翻转符号）；$\delta$：主元不能小于一个正的下限，保证正定；$(\theta_j/\beta)^2$：主元不能太小，否则 $l_{ij}=c_{ij}/d_j$ 会变得很大，因子里出现巨大元素，修正后的矩阵会面目全非。取三者中最大的，三件事同时满足。如果 $A$ 本来就充分正定，$|c_{jj}|$ 自然是最大的那个，不做任何修改。对照例 6.1 的 $n=3$ 公式看：$d_2=a_{22}-d_1l_{21}^2$ 是"第 2 个变量扣掉能被第 1 个变量解释的部分后剩下的方差"，就像回归里的残差方差。修正 Cholesky 只是在每一步检查这个"剩余方差"是否为正、是否太小，不行就把它抬高。 Algorithm 6.5（修正 Cholesky，Gill–Murray–Wright）

给定 δ > 0, β > 0
for k = 1..n: c_kk ← a_kk                    （初始化对角元）
for j = 1..n
  找 q 使 |c_qq| ≥ |c_ii|, i = j..n；交换第 j 与 q 行、列
  for s = 1..j−1: l_js ← c_js / d_s
  for i = j+1..n: c_ij ← a_ij − Σ_{s<j} l_js c_is
  θ_j ← max_{j<i≤n} |c_ij|（j = n 时为 0）
  d_j ← max{ |c_jj|, (θ_j/β)², δ }
  for i = j+1..n: c_ii ← c_ii − c_ij² / d_j
end

（精读笔记所据文本把选主元一步写在循环外；按原书文字说明，每个 $j$ 步都要选主元，此处已更正。）运算量约 $n^3/6$，与标准 Cholesky 相当，但行列交换的数据移动在大问题上可能代价可观；不需要额外存储（$L,D,c_{ij}$ 可以覆盖 $A$）。

结果可写成：设 $P$ 为置换矩阵，

PAP^T+E=LDL^T=MM^T,\tag{6.19}

$E$ 是非负对角矩阵，$A$ 充分正定时 $E=0$；$e_j=d_j-c_{jj}$，增大 $c_{jj}$ 等价于增大原数据 $a_{jj}$。参数选择：$\delta=u\max(\gamma(A)+\xi(A),1)$，其中 $\gamma=\max_i|a_{ii}|$，$\xi=\max_{i\ne j}|a_{ij}|$；Gill–Murray–Wright 建议

\beta=\max\left(\gamma(A),\ \frac{\xi(A)}{\sqrt{n^2-1}},\ u\right)^{1/2},

目的是让 $\|E\|_\infty$ 尽量小。

例 6.2　$A=\begin{bmatrix}4&2&1\\2&6&3\\1&3&-0.004\end{bmatrix}$，特征值约为 $-1.25,\ 2.87,\ 8.38$。Algorithm 6.5 给出 $E$ 的唯一非零元 $3.008$（加在第三个对角元上），修正后

A+E=\begin{bmatrix}4&2&1\\2&6&3\\1&3&3.004\end{bmatrix},

特征值为 $1.13,\ 3.00,\ 8.87$，条件数约 7.8，相当温和。本章量化实战的第 1 部分用自己实现的 Algorithm 6.5 复现了这个例子。Moré 和 Sorensen 证明：对精确 Hessian 使用 Algorithm 6.5，得到的 $B_k$ 条件数有界，即 (6.8) 成立。

6.4.4 Gershgorin 修正

先对 $A$ 用 Algorithm 6.5 得到 (6.19)；若 $E=0$ 就结束。否则计算 $\lambda_{\min}(A)$ 的两个上界：$b_1$ 由 Gershgorin 圆盘定理给出（保证 $A+b_1I$ 严格对角占优），$b_2=\max_ie_{ii}$。取 $\mu=\min(b_1,b_2)$，对 $A+\mu I$ 做 Cholesky 分解作为修正因子（Schnabel–Eskow）。$b_2$ 的作用是约束过松的 $b_1$。这是否优于单用 Algorithm 6.5 尚无定论。两种方法都不修改"充分正定"的矩阵，但"充分"很难用 $\lambda_{\min}$ 量化，所以它们有可能修改 $\lambda_{\min}>\delta$ 的矩阵。

6.4.5 修正对称不定分解

任何对称矩阵 $A$ 都可以写成

PAP^T=LBL^T,\tag{6.20}

$L$ 单位下三角，$B$ 是由 $1\times1$ 和 $2\times2$ 块组成的块对角矩阵，$P$ 为置换矩阵。允许 $2\times2$ 块，保证了分解总是存在且数值稳定（对不定矩阵算 $LDL^T$ 是不明智的，因为元素可能放大舍入误差）。

例 6.3　$A=\begin{bmatrix}0&1&2&3\\1&2&2&2\\2&2&3&3\\3&2&3&4\end{bmatrix}$，取 $P=[e_1,e_4,e_3,e_2]$，有

L=\begin{bmatrix}1&0&0&0\\0&1&0&0\\\frac19&\frac23&1&0\\\frac29&\frac13&0&1\end{bmatrix},\qquad B=\begin{bmatrix}0&3&0&0\\3&4&0&0\\0&0&\frac79&\frac59\\0&0&\frac59&\frac{10}9\end{bmatrix},\tag{6.21}

两个对角块都是 $2\times2$ 的。

惯性（inertia，正/零/负特征值的个数）：$B$ 与 $A$ 惯性相同（Sylvester 惯性定律）。算法中构造的 $2\times2$ 块总是一正一负两个特征值，所以 $A$ 的正特征值个数 = 正的 $1\times1$ 块数 + $2\times2$ 块数。

白话解释（对例 6.3 的补充说明）：上面"$2\times2$ 块总是一正一负"指的是 Bunch–Parlett、Bunch–Kaufman 等选主元规则实际选出的 $2\times2$ 主元——它们只在对角元相对非对角元很小时才被选中，此时行列式 $<0$，必为一正一负。例 6.3 的分解（经核算 $L B L^T$ 确实等于 $PAP^T$）只是展示分解的形式，并不是由这些规则产生的：第一个块 $\begin{bmatrix}0&3\\3&4\end{bmatrix}$ 行列式 $-9<0$，一正一负；第二个块 $\begin{bmatrix}7/9&5/9\\5/9&10/9\end{bmatrix}$ 行列式 $45/81>0$、迹 $>0$，两个特征值都为正。所以这个例子不能套用上面的计数公式，$A$ 的惯性应按各块实际特征值数：3 个正、1 个负、0 个零。练习 4 中"验证各为一正一负"对第二个块不成立，做题时请按此理解。 Sylvester 惯性定律的直观含义：做"换坐标"式的变换 $A\to LAL^T$（$L$ 可逆）不会改变"有几个方向向上弯、几个方向向下弯"。所以想知道一个大矩阵有几个负特征值，不必求特征值，看分解出的小块就行。- 过程：选一个非奇异的主元块 $E$（单个对角元，或两个对角元连同对应的非对角元），置换使之成为左上角子阵：$\Pi A\Pi^T=\begin{bmatrix}E&C^T\\C&H\end{bmatrix}$（6.22），做块分解
$\Pi A\Pi^T=\begin{bmatrix}I&0\\CE^{-1}&I\end{bmatrix}\begin{bmatrix}E&0\\0&H-CE^{-1}C^T\end{bmatrix}\begin{bmatrix}I&E^{-1}C^T\\0&I\end{bmatrix},$ 然后对 Schur 补 $H-CE^{-1}C^T$ 递归。

运算量约 $n^3/3$，外加选主元和置换的代价（可能可观）。理想的选主元策略应该便宜、使剩余矩阵的元素增长有限、并且不产生过多填充。
Bunch–Parlett：搜索整个工作矩阵，找最大对角元 $\xi_{dia}$ 与最大非对角元 $\xi_{off}$；若以对角元为 $1\times1$ 主元的元素增长（受 $\xi_{dia}/\xi_{off}$ 控制）可以接受就用它，否则取含 $\xi_{off}$ 的 $2\times2$ 块。数值稳定，$L$ 的最大元不超过 2.781；缺点是总共需要 $O(n^3)$ 次比较，总时间可能是 Algorithm 6.5 的两倍。
Bunch–Kaufman：每步至多搜索两列，总共 $O(n^2)$ 次比较；但 $L$ 的元素可能任意大，不适合修正 Cholesky 那样的策略。
有界 Bunch–Kaufman：折中方案——监控 $L$ 元素的大小，增长温和时用便宜的 BK 选择，过大时再进一步搜索；通常代价接近 BK，最坏接近 BP。
大型稀疏矩阵还须考虑主元选择对 $L$ 稀疏性的影响。

修正步骤：算出 (6.20) 后，对块对角矩阵 $B$ 做谱分解 $B=Q\Lambda Q^T$（非常便宜，可以逐块做），构造

F=Q\,\mathrm{diag}(\tau_i)Q^T,\qquad\tau_i=\begin{cases}0,&\lambda_i\ge\delta,\\\delta-\lambda_i,&\lambda_i<\delta,\end{cases}\tag{6.23}

即让 $B+F$ 所有特征值不小于 $\delta$ 的最小 Frobenius 修正。于是 $P(A+E)P^T=L(B+F)L^T$，$E=P^TLFL^TP$（一般非对角）。与修正 Cholesky 不同，这种方法改变的是整个 $A$ 而不只是对角线；它的目标是在 $\lambda_{\min}(A)<\delta$ 时使 $\lambda_{\min}(A+E)\approx\delta$，但并不总能做到（Cheng–Higham）。

量化联系：带等式约束的组合优化，其 KKT 矩阵 $\begin{bmatrix}\Sigma&A^T\\A&0\end{bmatrix}$ 是对称不定的。内点法等求解器内部正是用 Bunch–Kaufman 型的 $LBL^T$ 分解来解这类系统，并通过惯性检查判断当前点是否满足二阶条件（本册第 16a 章 16.3 节）。

6.5 信赖域牛顿法

信赖域方法不要求模型 Hessian 正定，可以直接取 $B_k=\nabla^2f(x_k)$：

\min_p m_k(p)=f_k+\nabla f_k^Tp+\tfrac12p^TB_kp\quad\text{s.t.}\quad\|p\|\le\Delta_k.\tag{6.24}

下面关注大规模实现，讨论四种技术。

6.5.1 Newton–Dogleg 与子空间极小化

Hessian 不定时 dogleg 不能直接用。可以先用 6.4 节的方法得到修正 Hessian 作为 $B_k$，保证模型是凸二次的，再沿 dogleg 路径

\tilde p(\tau)=\begin{cases}\tau p^U,&0\le\tau\le1,\\p^U+(\tau-1)(p^B-p^U),&1\le\tau\le2\end{cases}\tag{6.25}

极小化；或者做二维子空间极小化

\min_p m_k(p)\quad\text{s.t.}\quad\|p\|\le\Delta_k,\ p\in\mathrm{span}\{\nabla f_k,p^B\}.\tag{6.26}

优点：线性代数全部可以用直接法；全局收敛；修正分解在 Hessian 充分正定时不修改，可以保持局部快速收敛。不足：修正分解以一种"近乎随机"的方式扰动 Hessian，偏重某些方向，信赖域方法的好处可能丧失；而且这种修正其实是多余的——信赖域本身就引入了修正（解信赖域问题相当于分解 $\nabla^2f+\lambda I$，$\lambda$ 由半径决定）。结论：dogleg 最适合凸目标（Hessian 处处半正定）；一般情形宜用下面的方法。

6.5.2 信赖域问题的精确解

用 Algorithm 4.4 重复分解 $B_k+\lambda I$。经验上每次外层迭代平均要解 1–3 个线性系统，代价不算过高。得到的算法非常稳健——可以期望收敛到极小点，而不只是驻点（第 4 章定理 4.9）。大规模时每次迭代解多个线性系统可能负担过重。

6.5.3 信赖域 Newton–CG

在 Algorithm 4.3（CG–Steihaug）中取 $B_k=\nabla^2f(x_k)$，对

B_kp_k=-\nabla f_k\tag{6.27}

做 CG，在以下三种情况停止：(i) 近似解越出信赖域；(ii) 达到所需精度（按强制序列 (6.3)）；(iii) 遇到负曲率（此时沿负曲率方向走到边界）。它是 Algorithm 6.1 的信赖域对应物。

控制内层 CG 的精度是降低代价的关键。在良态的解附近，信赖域约束不起作用，算法退化为 6.2 节的非精确牛顿法，强制序列决定后期的收敛速度。
同样可以 Hessian-free（自动微分或 (6.7)）。
优点：全局收敛（第一步沿 $-\nabla f_k$，即 Cauchy 点，之后 CG 只会改进模型值）；不需要矩阵分解，可以利用稀疏性而不担心填充；CG 的核心是矩阵–向量乘，容易并行；Hessian 正定时越接近解越逼近纯牛顿步，收敛快。
相比线搜索 Newton–CG：步长受信赖域控制，不会因 Hessian 近奇异而过长；会探索负曲率方向——经验上有益，有时能让迭代离开非极小的驻点。
局限：它接受遇到的任何负曲率方向，即使由此带来的模型下降微不足道。原书的例子可以写成如下自洽的形式：$m(p)=10^{-3}p_1-10^{-4}p_1^2-p_2^2$，$\|p\|\le1$。在 $p=0$ 处最速下降方向是 $(-10^{-3},0)^T$，它恰好是一个负曲率方向，Algorithm 4.3 第一步就沿它走到边界，模型只下降约 $10^{-3}$；而沿 $e_2$（同样是负曲率方向）走到边界可以下降 1。（PDF 文本此处的符号与作者的论述不一致，这里按论述做了调整。）
补救：Hessian 有负特征值时，步应该在最负特征值对应的特征向量上有显著分量，以便迅速离开非极小驻点。可以用 Lanczos 方法代替 CG：遇到第一个负曲率方向后不终止，继续寻找"足够负"的曲率方向（更稳健，但子问题更贵）。SciPy 的 trust-krylov 就是这一思路。

6.5.4 Newton–CG 的预条件

Hessian 病态时，无预条件的 CG 可能很慢，甚至达不到所需精度，需要预条件：找非奇异的 $D$，使 $D^{-T}B_kD^{-1}$ 的特征值分布更好。

难点在于：预条件 CG 的迭代点不再保持欧氏范数单调增（第 4 章定理 4.2 只对无预条件成立），所以不能一碰到球形边界就停（后续迭代可能回到区域内）。解决办法：存在一个依赖于预条件子的加权范数，使迭代点在这个范数下单调增。考虑椭球信赖域

\min_p m_k(p)\quad\text{s.t.}\quad\|Dp\|\le\Delta_k,\tag{6.28}

令 $\hat p=Dp$，$\hat g_k=D^{-T}\nabla f_k$，$\hat B_k=D^{-T}B_kD^{-1}$，就化为标准形式 (6.24)，对它直接用 CG；Algorithm 4.3 监控的是 $\|\hat p\|=\|Dp\|$（单调增），超过 $\Delta_k$ 就终止。一句话：用椭球信赖域与预条件子相匹配。

常用的预条件子是不完全 Cholesky：$B=LL^T-R$，$L$ 的填充受限（如与 $B$ 的下三角部分同稀疏结构），$R$ 是不精确的部分。还需要处理 Hessian 可能不定的问题：

Algorithm 6.6（非精确修正 Cholesky）

（尺度化）T = diag(‖B e_i‖)；B̄ ← T^{−1/2} B T^{−1/2}；β ← ‖B̄‖
（平移以保证正定）if min_i b̄_ii > 0: α_0 ← 0  else α_0 ← β/2
for k = 0,1,2,...
  尝试不完全 Cholesky 分解：L Lᵀ = B̄ + α_k I
  if 成功: 返回 L   else α_{k+1} ← max(2α_k, β/2)
end

取预条件子 $D=L$。MINPACK-2 中的 NMTR 实现了这种信赖域 Newton–CG；LANCELOT 中也有使用略不同预条件子的 Newton–CG。

6.5.5 信赖域牛顿法的局部收敛

关键是证明：接近解时信赖域约束最终不起作用，近似解在区域内部，并且越来越接近纯牛顿步。满足后一条的步称为渐近精确（asymptotically exact）。

定理 6.4　设 $f$ 二阶 Lipschitz 连续可微，$x_k\to x^*$（$x^*$ 满足二阶充分条件）。设对充分大的 $k$，算法取 $B_k=\nabla^2f(x_k)$，步 $p_k$ 至少达到 Cauchy 下降（$m_k(p_k)\le m_k(p_k^C)$），且当 $\|p_k^N\|\le\frac12\Delta_k$ 时是渐近精确的：

\|p_k-p_k^N\|=o(\|p_k^N\|).\tag{6.29}

则对充分大的 $k$，信赖域约束不起作用。

证明要点：

无论 $\|p_k^N\|\le\frac12\Delta_k$ 与否，都有 $\|p_k\|\le2\|p_k^N\|\le2\|\nabla^2f_k^{-1}\|\|\nabla f_k\|$，即 $\|\nabla f_k\|\ge\frac12\|p_k\|/\|\nabla^2f_k^{-1}\|$。
由 Cauchy 下降估计 (4.34)： $m_k(0)-m_k(p_k)\ge c_1\|\nabla f_k\|\min\left(\Delta_k,\frac{\|\nabla f_k\|}{\|B_k\|}\right)\ge c_1\frac{\|p_k\|^2}{4\|\nabla^2f_k^{-1}\|^2\|\nabla^2f_k\|},$ 由连续性，充分大的 $k$ 时 $\ge c_3\|p_k\|^2$（6.30），$c_3=\dfrac{c_1}{8\|\nabla^2f(x^*)^{-1}\|^2\|\nabla^2f(x^*)\|}$。
由 Hessian 的 Lipschitz 连续性，$|(f(x_k)-f(x_k+p_k))-(m_k(0)-m_k(p_k))|\le\frac L2\|p_k\|^3$，所以 $|\rho_k-1|\le\frac{L}{2c_3}\|p_k\|\le\frac{L}{2c_3}\Delta_k.\tag{6.31}$
半径只在 $\rho_k<\frac14$ 时缩小，而 $\Delta_k$ 小于某个阈值 $\tilde\Delta$ 时 $\rho_k$ 必然接近 1，所以 $\{\Delta_k\}$ 有正下界；另一方面 $\|p_k^N\|\to0$，由 (6.29) $\|p_k\|\to0$，所以约束最终不起作用。$\square$

引理 6.5　设 $x_k\to x^*$（满足二阶充分条件）。用 $B_k=\nabla^2f(x_k)$ 的 dogleg (6.25) 或二维子空间法 (6.26)，充分大的 $k$ 时模型的无约束极小就是 $p_k^N$，所以满足定理 6.4 的条件。Newton–CG 若用 (6.3) 且 $\eta_k\to0$（加上越界/负曲率终止），也满足。

证明：dogleg 和子空间法：$p_k^N$ 在区域内、在 dogleg 路径上、在子空间内，又是区域内模型的极小点，所以 $p_k=p_k^N$。Newton–CG：充分大的 $k$ 时 Hessian 正定，不会因负曲率而停；CG 迭代点的范数递增但不超过 $\|p_k^N\|$，留在区域内；所以只会因 (6.3) 而停，于是

\|p_k-p_k^N\|\le\|\nabla^2f_k^{-1}\|\|r_k\|\le\eta_k\|\nabla^2f_k^{-1}\|\|\nabla f_k\|\le\eta_k\|\nabla^2f_k^{-1}\|\|\nabla^2f_k\|\|p_k^N\|,

满足 (6.29)。$\square$ Algorithm 4.4 也满足定理条件（牛顿步在区域内时 $\lambda=0$，$p_k=p_k^N$）。

结论：最终取精确牛顿步的方法二次收敛；Newton–CG 的渐近速度由 $\eta_k\to0$ 的速度决定（定理 6.2）。

白话解释：定理 6.4 的证明是一个"两头夹"的论证。一头：离解越近，牛顿步越短（$\|p_k^N\|\to0$，因为梯度趋于零）。另一头：步子越短，二次模型越准（误差是 $\|p\|^3$ 量级，而预测下降是 $\|p\|^2$ 量级），$\rho_k$ 越接近 1，半径就不会再被缩小，有一个正的下限。步子趋于零而半径有下限，所以最终步子一定落在圈内，信赖域约束"形同虚设"，算法就是纯牛顿法。这和第 3 章"线搜索牛顿法最终总接受单位步长"是同一个故事的两个版本：全局化机制只在远处起作用，到了近处会自动"让路"。

6.6 量化实战：协方差修复与大规模风险预算

本节有四个部分：

修正 Cholesky：实现 Algorithm 6.5，复现原书例 6.2。
特征值修正：在原书 (6.9) 的例子上比较几种修正得到的方向。
相关系数矩阵修复：风险管理中，人工设定的压力情景相关性、成对缺失数据估计的相关矩阵常常不是半正定的，无法做 Cholesky 分解来生成相关随机数。比较特征值截断 (6.12)、整体平移 (6.13) 和修正 Cholesky 三种修复方式。
Hessian-free Newton–CG：求解 2000 只股票的风险预算（risk budgeting）组合。采用 Spinu 的凸表述 $\min_{w>0}\ f(w)=\tfrac12w^T\Sigma w-\sum_ib_i\log w_i,$ 一阶条件 $w_i(\Sigma w)_i=b_i$ 说明解的风险贡献正比于预算 $b_i$（归一化后即风险预算组合；$b_i$ 相等时就是风险平价）。Hessian 为 $\Sigma+\mathrm{diag}(b_i/w_i^2)$，在 $w>0$ 上正定。$\Sigma$ 采用因子结构，Hessian–向量积只需 $O(nK)$，整个过程不构造任何 $n\times n$ 矩阵。我们手写 Algorithm 6.1（强制序列 $\min(0.5,\sqrt{\|\nabla f\|})$、负曲率检验、先试单位步的回溯，回溯同时保证 $w>0$），并与 SciPy 的 trust-ncg（信赖域 Newton–CG）和 L-BFGS-B 对照。

推导拆解：Spinu 表述为什么给出风险预算组合？第 1 步（梯度）：$\tfrac12w^T\Sigma w$ 的梯度是 $\Sigma w$（第 1 章推导过）；$-b_i\log w_i$ 对 $w_i$ 的导数是 $-b_i/w_i$。所以 $\nabla f=\Sigma w-b/w$（逐元素相除）。第 2 步（一阶条件）：令第 $i$ 个分量为零，$(\Sigma w)_i=b_i/w_i$，两边乘 $w_i$ 得 $w_i(\Sigma w)_i=b_i$。第 3 步（与风险贡献对应）：组合波动率 $\sigma_p=\sqrt{w^T\Sigma w}$，资产 $i$ 的风险贡献 $RC_i=w_i\,\partial\sigma_p/\partial w_i=w_i(\Sigma w)_i/\sigma_p$，所有 $RC_i$ 之和正好是 $\sigma_p$（欧拉分解，CFA 里"边际风险贡献 × 权重"）。一阶条件说 $w_i(\Sigma w)_i\propto b_i$，即 $RC_i\propto b_i$。第 4 步（归一化不影响）：把 $w$ 乘一个正数 $c$，$w_i(\Sigma w)_i$ 整体乘 $c^2$，比例关系不变，所以除以 $\sum w_i$ 后仍是风险预算组合。第 5 步（Hessian 与凸性）：对 $\Sigma w$ 再求导得 $\Sigma$；对 $-b_i/w_i$ 再求导得 $b_i/w_i^2>0$。所以 Hessian $=\Sigma+\mathrm{diag}(b_i/w_i^2)$，半正定加正定，处处正定，问题严格凸。这就是为什么这里可以放心用牛顿法，而直接最小化"风险贡献之差的平方和"（第 1 章提到的非凸表述）则不行。

import numpy as np
from scipy.optimize import minimize
np.set_printoptions(suppress=True)

# ---------- 1) Algorithm 6.5：修正 Cholesky（Gill–Murray–Wright），复现原书例 6.2 ----------
def modified_cholesky(A, delta=None, beta=None):
    """返回置换 P、单位下三角 L、对角 d、对角修正 e，使 P A P' + diag(e) = L diag(d) L'"""
    A = np.array(A, float); n = len(A); u = np.finfo(float).eps
    gamma = np.abs(np.diag(A)).max(); xi = np.abs(A - np.diag(np.diag(A))).max()
    if delta is None: delta = u * max(gamma + xi, 1.0)
    if beta is None:  beta = np.sqrt(max(gamma, xi / np.sqrt(n * n - 1), u))
    C = A.copy(); perm = np.arange(n); L = np.eye(n); d = np.zeros(n); e = np.zeros(n)
    for j in range(n):
        q = j + np.argmax(np.abs(np.diag(C)[j:]))          # 对称选主元：剩余对角元绝对值最大者
        C[[j, q]] = C[[q, j]]; C[:, [j, q]] = C[:, [q, j]]
        L[[j, q], :j] = L[[q, j], :j]; perm[[j, q]] = perm[[q, j]]
        theta = np.abs(C[j + 1:, j]).max() if j < n - 1 else 0.0
        d[j] = max(abs(C[j, j]), (theta / beta)**2, delta)   # (6.18)
        e[j] = d[j] - C[j, j]
        L[j + 1:, j] = C[j + 1:, j] / d[j]
        C[j + 1:, j + 1:] -= np.outer(C[j + 1:, j], C[j + 1:, j]) / d[j]   # 更新 Schur 补
    return perm, L, d, e

A = np.array([[4, 2, 1], [2, 6, 3], [1, 3, -0.004]])
perm, L, d, e = modified_cholesky(A)
P = np.eye(3)[perm]
M = L * np.sqrt(d)
print("例 6.2  A 的特征值:", np.round(np.linalg.eigvalsh(A), 3))
print("置换顺序:", perm, " 对角修正 E =", np.round(e, 3))
print("M =\n", np.round(M, 4))
A_mod = P.T @ (M @ M.T) @ P                           # 回到原顺序的 A + E
print("修正后矩阵:\n", np.round(A_mod, 3))
ev = np.linalg.eigvalsh(A_mod); print("修正后特征值:", np.round(ev, 2), " 条件数 %.1f" % (ev[-1] / ev[0]))

# ---------- 2) 特征值修正：原书 (6.9)-(6.13) 的例子 ----------
g = np.array([1.0, -3.0, 2.0]); H = np.diag([10.0, 3.0, -1.0])
pN = -np.linalg.solve(H, g)
print("\n纯牛顿步 %s，g'p = %.2f（>0，不是下降方向）" % (pN, g @ pN))
for name, Bk in [("负特征值换成 δ=1e-8 (6.10)", np.diag([10, 3, 1e-8])),
                 ("翻转负特征值符号", np.diag([10, 3, 1.0])),
                 ("整体平移 τI, δ=1 (6.13)", H + 2.0 * np.eye(3))]:
    p = -np.linalg.solve(Bk, g)
    print("%-28s p = %s  g'p = %.3g  ||p|| = %.3g" % (name, np.round(p, 3), g @ p, np.linalg.norm(p)))

# ---------- 3) 量化：修复非半正定的相关系数矩阵 ----------
C = np.array([[1.0, 0.9, 0.7, 0.2],      # 压力情景下人工设定的相关系数，互相不一致
              [0.9, 1.0, -0.4, 0.3],
              [0.7, -0.4, 1.0, 0.5],
              [0.2, 0.3, 0.5, 1.0]])
print("\n设定的相关矩阵特征值:", np.round(np.linalg.eigvalsh(C), 4))
lam, Q = np.linalg.eigh(C); delta = 1e-4
C_clip = Q @ np.diag(np.maximum(lam, delta)) @ Q.T          # (6.12) 最小 Frobenius 修正
s = 1 / np.sqrt(np.diag(C_clip)); C_clip = s[:, None] * C_clip * s[None, :]   # 再把对角线拉回 1
tau = max(0.0, delta - lam[0]); C_shift = (C + tau * np.eye(4)) / (1 + tau)   # (6.13) 平移后再归一
_, Lm, dm, em = modified_cholesky(C)
for name, X in [("特征值截断+归一", C_clip), ("平移 τI+归一", C_shift)]:
    print("%-16s 最小特征值 %.4f  与原矩阵 Frobenius 距离 %.4f" %
          (name, np.linalg.eigvalsh(X)[0], np.linalg.norm(X - C)))
print("修正 Cholesky   对角修正 E = %s，Frobenius 距离 %.4f（对角线不再为 1）" % (np.round(em, 4), np.linalg.norm(em)))
Lc = np.linalg.cholesky(C_clip)
Z = np.random.default_rng(0).standard_normal((200000, 4)) @ Lc.T
print("用修复后矩阵生成相关正态样本，样本相关系数:\n", np.round(np.corrcoef(Z.T), 3))

# ---------- 4) Hessian-free 线搜索 Newton–CG（Algorithm 6.1）：大规模风险预算组合 ----------
# Spinu 凸表述：min f(w) = 1/2 w'Σw - Σ b_i log w_i，w>0；解按比例缩放后风险贡献 w_i(Σw)_i ∝ b_i
rng = np.random.default_rng(6)
n, K = 2000, 8
Bx = rng.standard_normal((n, K)) * 0.3; Fc = np.diag(np.r_[0.04, 0.01 * np.ones(K - 1)])
dv = rng.uniform(0.15, 0.45, n)**2
Sv = lambda v: Bx @ (Fc @ (Bx.T @ v)) + dv * v              # Σv，O(nK)，从不构造 Σ
b = rng.uniform(0.5, 1.5, n); b /= b.sum()                   # 风险预算
f    = lambda w: 0.5 * w @ Sv(w) - b @ np.log(w) if (w > 0).all() else np.inf
grad = lambda w: Sv(w) - b / w
hessp = lambda w, v: Sv(v) + (b / w**2) * v                  # Hessian–向量积，无需 Hessian 矩阵

def newton_cg(w, tol=1e-10, maxit=100):
    total_cg = 0
    for k in range(maxit):
        g = grad(w); gn = np.linalg.norm(g)
        if gn < tol: return w, k, total_cg
        eta = min(0.5, np.sqrt(gn))                          # 强制序列：超线性
        z = np.zeros(n); r = g.copy(); d = -r                # 内层 CG 解 H p = -g，起点 0
        for j in range(200):
            Hd = hessp(w, d); dHd = d @ Hd
            if dHd <= 0:                                     # 负曲率检验 (6.6)
                if j == 0: z = d                             # 首步即负曲率：取最速下降
                break
            a = (r @ r) / dHd; z = z + a * d; r_new = r + a * Hd
            total_cg += 1
            if np.linalg.norm(r_new) <= eta * gn: break      # (6.3)
            d = -r_new + (r_new @ r_new) / (r @ r) * d; r = r_new
        alpha = 1.0                                          # 总先试单位步；回溯同时保证 w>0
        while f(w + alpha * z) > f(w) + 1e-4 * alpha * (g @ z): alpha *= 0.5
        w = w + alpha * z
    return w, maxit, total_cg

w0 = np.full(n, 1.0)
w, iters, ncg = newton_cg(w0)
x = w / w.sum(); rc = x * Sv(x); rc /= rc.sum()
print("\nNewton–CG：外层 %d 次，内层 CG 共 %d 次 Hessian–向量积；风险贡献与预算最大相对偏差 %.1e" %
      (iters, ncg, np.abs(rc / b - 1).max()))
res = minimize(f, w0, jac=grad, hessp=hessp, method="trust-ncg", options={"gtol": 1e-10})
print("SciPy trust-ncg：外层 %d 次，梯度范数 %.1e，与上面解的差 %.1e" %
      (res.nit, np.linalg.norm(res.jac), np.abs(res.x - w).max()))
res2 = minimize(f, w0, jac=grad, method="L-BFGS-B", bounds=[(1e-8, None)] * n,
                options={"gtol": 1e-10, "ftol": 1e-15, "maxiter": 5000})
print("对照 L-BFGS-B：%d 次迭代，与 Newton–CG 解的差 %.1e" % (res2.nit, np.abs(res2.x - w).max()))

运行输出：

例 6.2  A 的特征值: [-1.251  2.869  8.379]
置换顺序: [1 0 2]  对角修正 E = [0.    0.    3.008]
M =
 [[2.4495 0.     0.    ]
 [0.8165 1.8257 0.    ]
 [1.2247 0.     1.2264]]
修正后矩阵:
 [[4.    2.    1.   ]
 [2.    6.    3.   ]
 [1.    3.    3.004]]
修正后特征值: [1.13 3.   8.87]  条件数 7.9

纯牛顿步 [-0.1  1.   2. ]，g'p = 0.90（>0，不是下降方向）
负特征值换成 δ=1e-8 (6.10)         p = [-1.e-01  1.e+00 -2.e+08]  g'p = -4e+08  ||p|| = 2e+08
翻转负特征值符号                     p = [-0.1  1.  -2. ]  g'p = -7.1  ||p|| = 2.24
整体平移 τI, δ=1 (6.13)          p = [-0.083  0.6   -2.   ]  g'p = -5.88  ||p|| = 2.09

设定的相关矩阵特征值: [-0.4215  0.7799  1.4578  2.1838]
特征值截断+归一         最小特征值 0.0001  与原矩阵 Frobenius 距离 0.5134
平移 τI+归一         最小特征值 0.0001  与原矩阵 Frobenius 距离 0.5689
修正 Cholesky   对角修正 E = [0.     0.     0.7806 1.65  ]，Frobenius 距离 1.8253（对角线不再为 1）
用修复后矩阵生成相关正态样本，样本相关系数:
 [[ 1.     0.665  0.503  0.241]
 [ 0.665  1.    -0.242  0.231]
 [ 0.503 -0.242  1.     0.424]
 [ 0.241  0.231  0.424  1.   ]]

Newton–CG：外层 16 次，内层 CG 共 137 次 Hessian–向量积；风险贡献与预算最大相对偏差 1.6e-09
SciPy trust-ncg：外层 33 次，梯度范数 1.0e-09，与上面解的差 2.3e-09
对照 L-BFGS-B：91 次迭代，与 Newton–CG 解的差 1.7e-07

解读

例 6.2 复现：选主元后的顺序是 $(2,1,3)$，修正量 $E$ 只加在最后一个主元（即原矩阵第 3 个对角元）上，大小 3.008，与原书一致；修正后的特征值 $1.13,3.00,8.87$、条件数约 7.8（输出四舍五入为 7.9）也一致。$M$ 的行顺序与原书排列方式不同，但对应同一个分解。
特征值修正：纯牛顿步 $g^Tp=0.9>0$，是上升方向。把负特征值换成 $10^{-8}$ 得到的方向虽然下降，长度却是 $2\times10^8$，完全不可用；翻转符号或整体平移 $\tau I$ 都给出长度约 2 的合理方向。这说明"正定"只是必要条件，修正后的矩阵还必须条件良好——正是有界修正分解性质 (6.8) 的意义。
相关矩阵修复：设定的相关矩阵最小特征值为 $-0.42$。特征值截断（再把对角线归一）与原矩阵的 Frobenius 距离最小（0.51），平移法次之（0.57）。修正 Cholesky 虽然直接给出了可用的因子，但修正全加在对角线上（输出中按选主元后的顺序列出），距离 1.83，而且对角线不再为 1——它是为优化算法设计的（要求快、条件数有界），不是为"最近相关矩阵"设计的。实务中修复相关矩阵应首选特征值截断，或更精确的 Higham 最近相关矩阵算法（交替投影）。修复后矩阵生成的模拟样本相关性，在原来"不一致"的位置（如资产 1–2 与 2–3）都做了折中。
风险预算：手写的 Newton–CG 外层 16 次、共 137 次 Hessian–向量积就收敛到梯度 $10^{-10}$，风险贡献与预算的最大相对偏差 $1.6\times10^{-9}$。SciPy 的 trust-ncg 得到同一个解；L-BFGS-B 需要 91 次迭代且精度略低。计算量方面，每次 Hessian–向量积只需约 $2nK$ 次乘法，整个求解不到一秒；若构造稠密的 $2000\times2000$ Hessian 并做 Cholesky 分解，每次外层迭代就要约 $n^3/3\approx2.7\times10^9$ 次运算。

本章小结

实用牛顿法要解决两个问题：Hessian 不正定时如何保证下降与全局收敛，大规模时如何控制计算量。非精确牛顿法用强制序列控制内层求解精度：$\eta_k\le\eta<1$ 线性收敛，$\eta_k\to0$ 超线性，$\eta_k=O(\|\nabla f_k\|)$ 二次。线搜索 Newton–CG 从零出发做 CG、遇负曲率即停，只需 Hessian–向量积；修正牛顿法把 Hessian 修正为条件数有界的正定矩阵，修正方法包括特征值修正、加 $\tau I$、修正 Cholesky、Gershgorin 修正和修正对称不定分解。信赖域牛顿法不需要修正 Hessian：dogleg 适合凸问题，精确解法最稳健，信赖域 Newton–CG 最适合大规模问题并能利用负曲率；预条件要与椭球信赖域配合。在解附近，信赖域约束失效，各方法都恢复牛顿法的快速局部收敛。

概念	公式/要点
非精确牛顿	$\Vert \nabla^2f_kp_k+\nabla f_k\Vert \le\eta_k\Vert \nabla f_k\Vert $
强制序列与速度	$\eta_k\le\eta<1$ 线性；$\eta_k\to0$ 超线性；$\eta_k=O(\Vert \nabla f_k\Vert )$ 二次
Newton–CG 规则	CG 起点 0；$d^TAd\le0$ 即停；首步即负曲率则取 $-\nabla f_k$
Hessian–向量积	$\nabla^2f(x)p\approx[\nabla f(x+hp)-\nabla f(x)]/h$，或自动微分
有界修正分解	$\Vert B_k\Vert \Vert B_k^{-1}\Vert \le C$ ⇒ $\lim\nabla f_k=0$
Frobenius 最优修正	把小于 $\delta$ 的特征值抬到 $\delta$
2-范数最优修正	$A+\tau I$，$\tau=\max(0,\delta-\lambda_{\min})$
修正 Cholesky	$d_j=\max(\vert c_{jj}\vert ,(\theta_j/\beta)^2,\delta)$；$PAP^T+E=LDL^T$
对称不定分解	$PAP^T=LBL^T$，$B$ 含 $1\times1$/$2\times2$ 块，惯性同 $A$
信赖域牛顿局部收敛	渐近精确步 ⇒ 约束最终失效 ⇒ 牛顿速度

练习

基础

对 $f(x)=\frac12x^Tx+\frac14\sigma(x^TAx)^2$ 写出梯度、Hessian 和 Hessian–向量积的表达式，说明 Hessian–向量积可以不构造 Hessian 而以 $O(n^2)$（$A$ 稠密时）或更低代价计算。提示：$\nabla f=x+\sigma(x^TAx)Ax$，$\nabla^2f=I+\sigma(x^TAx)A+2\sigma Axx^TA$。
（原书 6.2）证明线搜索 Newton–CG 产生的方向总是下降方向。提示：CG 从 0 出发，每一步都降低二次模型 $q(p)=\nabla f_k^Tp+\frac12p^T\nabla^2f_kp$，且第一步就是负梯度方向的正倍数。
（原书 6.4）对 $A=\mathrm{diag}(-2,12,4)$ 手算 Algorithm 6.5 的结果（取 $\delta$ 很小，$\beta$ 按建议公式），说明修正后矩阵是什么。提示：选主元顺序为 12、4、−2；对角阵的 $\theta_j=0$，于是 $d_j=\max(|c_{jj}|,\delta)$，$-2$ 被修正为 $2$。
（原书 6.3）计算 (6.21) 中两个 $2\times2$ 块的特征值，验证各为一正一负，并据此确定 $A$ 的惯性。
说明为什么在 Algorithm 6.1 中强制序列取 $\min(0.5,\|\nabla f_k\|)$ 而不是常数 $0.5$，会让后期收敛更快但每步内层 CG 更贵。

进阶

（原书 6.6）证明定理 6.2。
（原书 6.7）证明块对角矩阵的特征分解可以逐块计算，并据此说明 (6.23) 的修正为什么便宜。
（原书 6.1）实现无线搜索的纯牛顿迭代，用 CG 求方向，分别选择强制序列使收敛为线性、超线性、二次。在 $f(x)=\frac12x^Tx+0.25\sigma(x^TAx)^2$ 上测试，其中 $A=\begin{bmatrix}5&1&0&0.5\\1&4&0.5&0\\0&0.5&3&0\\0.5&0&0&2\end{bmatrix}$，$x_1=(\cos70°,\sin70°,\cos70°,\sin70°)^T$，$\sigma=1$ 或更大。观察 $\sigma$ 增大（偏离二次）时的影响。
在本章实战第 4 部分，把强制序列改为 $\min(0.5,\|\nabla f_k\|)$ 和常数 $0.5$，记录外层迭代次数和 Hessian–向量积总数，解释结果。
用 Higham 的交替投影算法（在"半正定矩阵集合"与"单位对角矩阵集合"之间交替投影）求本章实战第 3 部分相关矩阵的最近相关矩阵，与特征值截断法比较 Frobenius 距离。

原书推荐习题：6.1（强制序列与收敛速度的数值验证）、6.2（Newton–CG 方向的下降性证明）、6.4（手算修正 Cholesky）、6.5（无预条件 CG–Steihaug 信赖域 Newton–CG，并构造在原点附近有负曲率的函数观察负曲率步）、6.6（证明定理 6.2）。

原书对照

本章内容	原书位置	PDF 页码
6.1 引言	Ch.6 开头	PDF p.155–156
6.2 非精确牛顿，定理 6.1、6.2	§6.1 Inexact Newton Steps	PDF p.156–158
6.3.1 线搜索 Newton–CG，Algorithm 6.1	§6.2 Line Search Newton–CG	PDF p.159–161
6.3.2 修正牛顿法，定理 6.3	Modified Newton's Method，Algorithm 6.2	PDF p.161–162
6.4.1 特征值修正	§6.3 Eigenvalue Modification	PDF p.163–164
6.4.2 加单位阵倍数，Algorithm 6.3	Adding a Multiple of the Identity	PDF p.164–165
6.4.3 修正 Cholesky，Algorithm 6.4、6.5，例 6.1、6.2	Modified Cholesky Factorization	PDF p.165–170
6.4.4 Gershgorin 修正	Gershgorin Modification	PDF p.170
6.4.5 修正对称不定分解，例 6.3	Modified Symmetric Indefinite Factorization	PDF p.171–174
6.5.1–6.5.3 Newton–dogleg、精确解、信赖域 Newton–CG	§6.4 Trust-Region Newton Methods	PDF p.174–177
6.5.4 Newton–CG 的预条件，Algorithm 6.6	Preconditioning the Newton–CG Method	PDF p.177–179
6.5.5 局部收敛，定理 6.4，引理 6.5	Local Convergence of Trust-Region Newton Methods	PDF p.179–182
注释与参考、习题	Notes and References，Exercises 6.1–6.7	PDF p.182–183

页码换算：原书正文页码 = PDF 页码减 20（已用 PDF 页眉核对；第 1–2 章减 21，第 3–11 章减 20，第 12–16 章减 19，第 17 章至附录减 18）。

概念	公式/要点
非精确牛顿	\(\Vert \nabla^2f_kp_k+\nabla f_k\Vert \le\eta_k\Vert \nabla f_k\Vert \)
强制序列与速度	\(\eta_k\le\eta<1\) 线性；\(\eta_k\to0\) 超线性；\(\eta_k=O(\Vert \nabla f_k\Vert )\) 二次
Newton–CG 规则	CG 起点 0；\(d^TAd\le0\) 即停；首步即负曲率则取 \(-\nabla f_k\)
Hessian–向量积	\(\nabla^2f(x)p\approx[\nabla f(x+hp)-\nabla f(x)]/h\)，或自动微分
有界修正分解	\(\Vert B_k\Vert \Vert B_k^{-1}\Vert \le C\) ⇒ \(\lim\nabla f_k=0\)
Frobenius 最优修正	把小于 \(\delta\) 的特征值抬到 \(\delta\)
2-范数最优修正	\(A+\tau I\)，\(\tau=\max(0,\delta-\lambda_{\min})\)
修正 Cholesky	\(d_j=\max(\vert c_{jj}\vert ,(\theta_j/\beta)^2,\delta)\)；\(PAP^T+E=LDL^T\)
对称不定分解	\(PAP^T=LBL^T\)，\(B\) 含 \(1\times1\)/\(2\times2\) 块，惯性同 \(A\)
信赖域牛顿局部收敛	渐近精确步 ⇒ 约束最终失效 ⇒ 牛顿速度

第 06 章 实用牛顿法