第 03b 章 Weyr 标准形与三角分解

对应原书：Horn & Johnson《Matrix Analysis》第 2 版，第 3 章的 3.4.2–3.5 节（书 p.203–223，PDF p.223–243）。本章接第 03a 章（Jordan 标准形）。

本章分两部分，分量很不相同。前半部分是 Weyr 标准形：它和 Jordan 形包含相同的信息，但在研究"与给定矩阵交换的矩阵"时好用得多。它是纯理论工具，量化中没有直接应用，这里讲清楚它是什么、为什么存在、何时需要回原书，并顺带得到"投影矩阵的酉标准形"这个对回归分析有用的推论。后半部分是 LU 类三角分解：LU、LDU、PLU、LPU。它们是数值线性代数的主力，线性方程组、协方差矩阵求逆、有限差分定价都离不开它；对称正定情形的 \(LDL^T\) 与 Cholesky 分解还有清楚的统计含义——主元就是逐次条件方差。

学习目标

读完本章，你应当能够：

1. 说明 Weyr 块与 Weyr 标准形的构造，知道它与 Jordan 形置换相似、包含相同的信息；理解它的核心优点：与 Weyr 矩阵交换的矩阵都是分块上三角的。
2. 了解酉 Weyr 形（Littlewood 定理），并用它得到投影矩阵的酉标准形：正交投影的非零奇异值全为 1，斜投影的"斜度"由大于 1 的奇异值刻画。
3. 掌握 LU 分解的存在性判据（行包含性质；非奇异时等价于所有顺序主子式非零），LDU 分解的唯一性与主元公式 \(d_i=\det A_{[i]}/\det A_{[i-1]}\)。
4. 理解对称正定矩阵的 \(LDL^T\) 分解：主元 \(d_i\) 是第 \(i\) 个变量对前 \(i-1\) 个变量回归的残差方差，\(L^{-1}\) 给出回归系数。
5. 知道 PLU（选主元）与 LPU 分解对任意矩阵存在，理解选主元的数值必要性，能用三对角 LU 实现有限差分期权定价。

读前导读

这一章在解决什么问题。 本章其实是两章拼在一起，对你的价值相差很大。

前半（3.10–3.11）的 Weyr 标准形是 Jordan 形的"换一种排法"，用来研究"哪些矩阵和 \(A\) 可交换"。这是纯理论工具，量化中几乎用不到。第一遍可以整体跳过 3.10 和 3.11.1，只读 3.11.2 的"直观"三条：OLS 的帽子矩阵是"正交投影"，工具变量（IV）估计对应"斜投影"，弱工具让斜投影极度放大噪声——这一点你在 CFA 里没见过，但对理解计量结果为什么不稳定很有用。

后半（3.12）的三角分解是你每天都在间接使用的东西。解线性方程组 \(Ax=b\)（均值方差最优权重、回归正规方程、有限差分定价每一步）时，软件做的就是把 \(A\) 拆成下三角 \(L\) 乘上三角 \(U\)，然后一行一行往回代。对协方差矩阵，这个拆法有一个漂亮的统计意义：\(\Sigma=LDL^T\) 中的 \(d_i\) 就是第 \(i\) 个资产对前面所有资产做回归后剩下的残差方差，\(L\) 里藏着全部"逐步对冲比率"。如果你做过 CFA 里"用股指期货对冲个股、最优对冲比率 \(=\rho\sigma_S/\sigma_F\)、对冲后剩余方差 \(=\sigma_S^2(1-\rho^2)\)"的计算，那就是 \(2\times2\) 的 \(LDL^T\)；本章把它推广到 \(n\) 个资产逐个加入的情形。Cholesky 分解（蒙特卡洛生成相关随机数的标准工具）是它的直接推论。

需要先想起来的数学。

• 三角矩阵与回代：上三角矩阵的主对角线以下全为 0。解 \(\begin{bmatrix}2&1\\0&3\end{bmatrix}x=\begin{bmatrix}5\\6\end{bmatrix}\)：先由第 2 行得 \(x_2=2\)，再代入第 1 行 \(2x_1+2=5\) 得 \(x_1=1.5\)。三角矩阵的行列式等于对角元之积。见 第 00 册第 06 章 线性代数速成。
• 高斯消元：用第 1 行的倍数把第 1 列下面的元素消成 0，再用第 2 行处理第 2 列……这就是 LU 分解的计算过程；每一步除以的那个对角元叫"主元"（pivot）。
• 顺序主子式：左上角 \(k\times k\) 子矩阵 \(A_{[k]}\) 的行列式。例：\(\begin{bmatrix}4&2\\2&5\end{bmatrix}\) 的顺序主子式是 \(4\) 和 \(4\cdot5-2\cdot2=16\)。
• 简单回归与条件方差：\(y\) 对 \(x\) 回归，斜率 \(\beta=\operatorname{Cov}(x,y)/\operatorname{Var}(x)\)，残差方差 \(\operatorname{Var}(y)(1-\rho^2)\)。这是你熟悉的 CFA 内容，本章 3.12.4 把它推到多元。
• 投影：把向量 \(y\) 投到 \(X\) 的列空间上得到拟合值 \(\hat y=Hy\)。投影的特点是"投两次和投一次一样"，即 \(H^2=H\)（幂等）。

符号提示：\(A[\{1..i\}]\) 或 \(A_{[i]}\) 表示取 \(A\) 的第 \(1\) 到 \(i\) 行、第 \(1\) 到 \(i\) 列；\(A[\{i+1\},\{1..i\}]\) 表示取第 \(i+1\) 行的前 \(i\) 个元素；\(O(n^2)\) 表示"运算量与 \(n^2\) 同一数量级"；\(\operatorname{nullity}\) 是零空间的维数（\(=n-\)秩）。

怎么读这一章。 核心必读：3.11.2 的"直观"部分、3.12.1–3.12.5，以及实战 1、2。第一遍可跳过：3.10 全部（Weyr 形，量化关系弱）、3.11.1（Littlewood 定理）、3.12.6（LPU 与三角等价，了解即可）、3.12.7 的 Lanczos 部分、实战 3 的前两段。建议顺序：3.12.1 → 3.12.3 → 3.12.4（与实战 1 对照着读）→ 3.12.5 → 实战 2，最后有兴趣再回头看 Weyr 形。

3.10 Weyr 标准形

3.10.1 动机：交换矩阵的结构

第 03a 章 3.6.4 节说：与单个 Jordan 块交换的矩阵是上三角 Toeplitz 矩阵；所以与非减次 Jordan 矩阵交换的矩阵都是上三角的。但一旦某个特征值有多个 Jordan 块（减次），情况就变复杂。原书的例子：

\(J\) 与 \(W\) 相似（\(W\) 就是后面要定义的 Weyr 形）。

• \(A\) 与 \(J\) 交换 ⇔ \(A=\begin{bmatrix}B&C\\D&E\end{bmatrix}\)，四块都是 \(2\times2\) 上三角 Toeplitz 矩阵——\(A\) 不一定是上三角的（\(D\) 可以非零）。
• \(A\) 与 \(W\) 交换 ⇔ \(A=\begin{bmatrix}B&C\\0&B\end{bmatrix}\)——分块上三角。

所以研究交换矩阵族时，Weyr 形比 Jordan 形更"整齐"。

推导拆解：第二条可以直接验证。记 \(N=W-\lambda I=\begin{bmatrix}0&I_2\\0&0\end{bmatrix}\)，\(A\) 与 \(W\) 交换等价于 \(A\) 与 \(N\) 交换（\(\lambda I\) 与一切交换）。分块相乘：\(AN=\begin{bmatrix}0&B\\0&D\end{bmatrix}\)，\(NA=\begin{bmatrix}D&E\\0&0\end{bmatrix}\)。令两者相等，逐块比较得 \(D=0\)、\(B=E\)，于是 \(A=\begin{bmatrix}B&C\\0&B\end{bmatrix}\)，\(B,C\) 任意。Jordan 形 \(J_2(\lambda)\oplus J_2(\lambda)\) 把两条链"并排"放，交换矩阵可以在两条链之间来回搬运，所以出现下三角块；Weyr 形把两条链"按层"放（先放所有链的第 1 层，再放第 2 层），交换矩阵只能把高层搬到低层，所以是分块上三角。

3.10.2 Weyr 块

回忆 Weyr 特征 \(w_k(A,\lambda)=\operatorname{rank}(A-\lambda I)^{k-1}-\operatorname{rank}(A-\lambda I)^k\)，即阶数 \(\ge k\) 的 Jordan 块个数，\(w_1\ge w_2\ge\cdots\ge w_q\ge1\)，\(q\) 是 \(\lambda\) 的指标。

定义（Weyr 块）　给定 \(\lambda\) 和非增正整数列 \(w=(w_1,\dots,w_q)\)，Weyr 块 \(W(w,\lambda)\) 是 \(q\times q\) 分块上三角块双对角矩阵

可以把它看成 Jordan 块的"分块版本"：对角块是纯量矩阵，超对角块是列满秩的 \(\begin{bmatrix}I\\0\end{bmatrix}\)。

秩的计算：由于 \(G_{w_{k-1},w_k}G_{w_k,w_{k+1}}=G_{w_{k-1},w_{k+1}}\)，每升一次幂，非零块整体上移一条块对角线：

于是 \(\operatorname{rank}(W-\lambda I)^{p-1}-\operatorname{rank}(W-\lambda I)^p=w_p\)：\(W(w,\lambda)\) 的 Weyr 特征恰好是 \(w\)。对角块个数 \(q\) 等于 \(\lambda\) 的指标。

推导拆解：用下面的 \(5\times5\) 例子（\(w=(2,2,1)\)）核对。\(N=W-\lambda I\) 只有三个 1，位于 \((1,3),(2,4),(3,5)\)，它把 \(e_3\to e_1\)、\(e_4\to e_2\)、\(e_5\to e_3\)。\(p=1\)：秩 3 \(=w_2+w_3=2+1\)；\(p=2\)：\(N^2\) 只剩 \(e_5\to e_1\)，秩 1 \(=w_3\)；\(p=3\)：\(N^3=0\)。相邻相减得 \(5-3=2\)、\(3-1=2\)、\(1-0=1\)，正好读回 \(w=(2,2,1)\)。"列满秩"的 \(G\) 块保证每一层的点都连到下一层的不同点上，不会"合并"导致秩额外下降。

Weyr 矩阵是特征值两两不同的 Weyr 块的直和。

例 (3.4.2.2)　第 03a 章例 (3.1.16a) 的 \(J\in M_{13}\)（Weyr 特征 \(6,5,2\)）对应

例　\(J=J_3(\lambda)\oplus J_2(\lambda)\)，\(w=(2,2,1)\)：

3.10.3 Weyr 标准形定理

定理 3.4.2.3　设 \(A\in M_n\) 的不同特征值为 \(\lambda_1,\dots,\lambda_d\)（任意给定次序），则存在非奇异 \(S\) 使

给定特征值次序后，这个 Weyr 矩阵由 \(A\) 唯一确定，称为 \(A\) 的 Weyr 标准形。\(A\) 实且特征值全实时 \(S\) 可取实。

证明一句话：\(W_A\) 与 \(A\) 对每个特征值的 Weyr 特征相同，而 Weyr 特征决定 Jordan 形（第 03a 章引理 3.1.18），所以两者相似于同一个 Jordan 矩阵。

Weyr 形与 Jordan 形的关系：

• 两者信息相同，呈现方式不同：Weyr 形直接显示 Weyr 特征（对角块大小），Jordan 形直接显示 Segre 特征（块大小）。
• \(W_A\) 与 \(J_A\) 置换相似（3.4.P8）。构造方法：画点图，第 \(k\) 行放 \(w_k\) 个点，逐行标号 \(1,2,\dots,n\)，再逐列读出标号得到置换。例：\(J_3(0)\oplus J_2(0)\)，\(w=(2,2,1)\)，点图逐行标号后逐列读出 \(\sigma=(1,3,5,2,4)\)，\(P=[e_1\ e_3\ e_5\ e_2\ e_4]\)，\(J=P^TWP\)（实战 3 验证）。
• 两者的非零元个数、非零非对角元个数都相同；Weyr 形同样在相似类中非零非对角元最少（3.4.P9）。
• 何时相同（3.4.P10–P11）：Weyr 形等于 Jordan 形 ⇔ 每个特征值要么只有一个 Jordan 块（非减次部分），要么所有块都是 \(1\times1\)（可对角化部分）。

3.10.4 与 Weyr 矩阵交换的矩阵

引理 3.4.2.4　设 \(F,F'\) 是同样分块（对角块 \(\lambda I_{n_1},\dots,\lambda I_{n_k}\)，\(n_1\ge\cdots\ge n_k\)）的块上三角矩阵，且 \(F'\) 的超对角块都列满秩。若 \(AF=F'A\)，则 \(A\) 共形分块上三角；若 \(A\) 还正规，则 \(A\) 分块对角。

证明思路：比较 \(AF=F'A\) 的块，从第一块列的最下方开始：\(F'_{k-1,k}A_{k1}=0\)，列满秩推出 \(A_{k1}=0\)；再向上、再向右，逐块推出下三角块全为零。

定理 3.4.2.5（Belitskii）　若 \(AB=BA\) 且 \(A=SW_AS^{-1}\)，则 \(S^{-1}BS\) 与 \(W_A\) 共形分块对角（不同特征值之间由 Sylvester 定理隔开），且每个对角块 \(B^{(\ell)}\) 与 \(W_A(\lambda_\ell)\) 的分块共形分块上三角。

原书进一步给出这些块之间的恒等式（3.4.2.12）：沿块对角线向上，\(B_{j-k-1,j-1}=\begin{bmatrix}B_{j-k,j}&\star\\0&\star\end{bmatrix}\)，即"上一层的块包含下一层的块作为左上角"。为了看清这种结构，原书引入标准分划（把 \(W_J\) 进一步细分，使每个对角块是纯量方阵、每个非对角块是单位阵或零矩阵的最粗分划），在 \(M_{13}\) 例子上逐块写出了与 \(W_J\) 交换的矩阵的形状。

定理 3.4.2.10(b)（O'Meara–Vinsonhaler）　设 \(\{A,A_1,A_2,\dots\}\) 是交换族，则存在同一个相似变换，把 \(A\) 化成 Weyr 标准形、同时把其余每个 \(A_i\) 化成上三角。

把"Weyr"换成"Jordan"，这个结论一般不成立。原书 3.4.P7 的反例：\(J=J_2(0)\oplus J_2(0)\) 与 \(A=\begin{bmatrix}0&I_2\\J_2(0)&0\end{bmatrix}\) 交换，但不存在相似变换使 \(J\) 保持 Jordan 形而 \(A\) 变成上三角。

中心化子的维数（3.4.P2–P4）　与 \(A\) 交换的矩阵全体 \(\mathcal C(A)\) 是一个代数，维数为

等号 ⇔ \(A\) 非减次。对 (3.1.16a) 的 \(J\)，\(\dim\mathcal C(J)=6^2+5^2+2^2=65\)。实战 3 用 Kronecker 积数值验证了这一点。

3.10.5 实 Jordan 形的其他推论与习题（选读）

• 3.4.P1：\(A\in M_n(\mathbf R)\) 满足 \(A^2=-I\)，则 \(n\) 为偶数，且实相似于 \(\begin{bmatrix}0&-I_{n/2}\\I_{n/2}&0\end{bmatrix}\)（"复结构"的标准形）。
• 3.4.P6：\(2\times2\) 实矩阵相似于 \(\begin{bmatrix}1&1\\-1&1\end{bmatrix}\) ⇔ 它形如 \(\begin{bmatrix}1+\alpha&(1+\alpha^2)/\beta\\-\beta&1-\alpha\end{bmatrix}\)，\(\beta\ne0\)。
• 推论 3.4.1.8：若 \(A=\begin{bmatrix}B&C\\0&0\end{bmatrix}\) 且 \(B\) 相似于实矩阵，则 \(A\) 也相似于实矩阵。

3.11 酉 Weyr 形与投影矩阵

3.11.1 Littlewood 定理

Jordan 形和 Weyr 形都需要一般的（可能很病态的）相似变换。如果只允许酉相似，能保留多少 Weyr 结构？

定理 3.4.3.1（Littlewood，要点）　任意 \(A\in M_n\) 酉相似于块上三角矩阵

其中每个特征值 \(\lambda_j\) 依次占据 \(q_j\)（它的指标）个对角块，这些对角块的大小恰为 \(\lambda_j\) 的 Weyr 特征 \(w_1(A,\lambda_j)\ge\cdots\ge w_{q_j}(A,\lambda_j)\)；同一特征值相邻两块之间的超对角块 \(F_{i,i+1}\) 是对角元为正实数的上三角矩阵。\(F\) 在"分块对角酉相似"意义下唯一：若还有满足同样条件的 \(F'\)，则 \(F'=UFU^*\)，\(U=U_1\oplus\cdots\oplus U_p\)。

证明：\(A=SW_AS^{-1}\)，对 \(S\) 做 QR 分解 \(S=QR\)，则 \(A=Q(RW_AR^{-1})Q^*\)；上三角相似 \(RW_AR^{-1}\) 保持对角块为纯量矩阵，并把超对角块 \(G\) 变成对角元为正的上三角块。唯一性用引理 3.4.2.4（酉矩阵是正规的，所以推出分块对角）。

这是 Schur 三角化的一个精化：Schur 形只保证对角线上是特征值，酉 Weyr 形还把 Jordan 结构"刻"进了分块形状里。

3.11.2 投影矩阵的酉标准形

推论 3.4.3.3　设 \(A^2=A\)（幂等，也叫投影），\(r=\operatorname{rank}A\)，奇异值中大于 1 的有 \(g\) 个：\(\sigma_1\ge\cdots\ge\sigma_g>1\)。则 \(A\) 酉相似于

推论：两个同阶投影酉相似 ⇔ 酉等价 ⇔ 奇异值相同。

证明：幂等矩阵的最小多项式整除 \(t(t-1)\)，可对角化，特征值 1 与 0 的指标都是 1，由 Littlewood 定理 \(A\) 酉相似于 \(\begin{bmatrix}I_r&F_{12}\\0&0\end{bmatrix}\)；对 \(F_{12}\) 做 SVD 再做置换相似，就拆成若干 \(\begin{bmatrix}1&s_i\\0&0\end{bmatrix}\) 块，其奇异值为 \(\sqrt{1+s_i^2}\) 和 0。

直观：

• 正交投影（\(A=A^*=A^2\)）的非零奇异值全是 1，没有 \(2\times2\) 块。OLS 的帽子矩阵 \(H=X(X^TX)^{-1}X^T\) 就是这种。
• 斜投影的每个 \(2\times2\) 块 \(\begin{bmatrix}1&s\\0&0\end{bmatrix}\) 是在一个平面内沿斜方向的投影，\(s=\sqrt{\sigma^2-1}\) 越大越"斜"，即投影方向与被投影子空间的夹角越小。
• 计量中的例子：恰好识别的工具变量估计对应斜投影 \(P=X(Z^TX)^{-1}Z^T\)（\(P^2=P\) 但 \(P\ne P^T\)）。单个回归元、单个工具时，\(\|P\|_2=1/|\cos\angle(x,z)|\)；弱工具意味着 \(x\) 与 \(z\) 几乎正交，奇异值巨大，投影极不稳定。实战 3 给出数值演示。

金融直觉：把投影想成"在墙上投影子"。正交投影是太阳在正上方，影子不会比物体长（奇异值 \(\le1\)）；斜投影是太阳很低，影子可以被拉得很长，太阳越低（投影方向与墙越接近平行）影子越长，这个"拉长倍数"就是大于 1 的奇异值 \(\sigma\)。 公式 \(\|P\|_2=1/|\cos\angle(x,z)|\) 的来源：单变量时 \(P=\dfrac{xz^T}{z^Tx}\) 是秩一矩阵，秩一矩阵 \(uv^T\) 的唯一非零奇异值是 \(\|u\|\|v\|\)，所以 \(\|P\|_2=\dfrac{\|x\|\|z\|}{|z^Tx|}=\dfrac1{|\cos\angle(x,z)|}\)。工具与内生变量的相关系数从 0.7 降到 0.1，噪声放大倍数从约 1.4 升到 10。这与 CFA 里"相关性低的对冲工具，对冲比率不稳定"是同一类现象。

3.4.P5（平方零矩阵）　\(A^2=0\)、秩 \(r\)、正奇异值 \(\sigma_1,\dots,\sigma_r\)，则 \(A\) 酉相似于 \(\bigoplus_i\begin{bmatrix}0&\sigma_i\\0&0\end{bmatrix}\oplus0_{n-2r}\)；两个平方零矩阵酉相似 ⇔ 奇异值相同（与第 02b 章 2.6.P23 一致）。

3.12 三角分解

3.12.1 为什么要三角分解

若 \(Ax=b\) 的系数矩阵是非奇异的上三角阵，回代（back substitution）即可：先由 \(a_{nn}x_n=b_n\) 求 \(x_n\)，再由第 \(n-1\) 个方程求 \(x_{n-1}\)，依次向上；下三角用前代（forward substitution）。每次只需 \(O(n^2)\) 次运算。若 \(A=LU\)（\(L\) 下三角、\(U\) 上三角），则先前代解 \(Ly=b\)，再回代解 \(Ux=y\)。

推导拆解：用一个 \(2\times2\) 例子把"消元 = LU"走一遍。\(A=\begin{bmatrix}4&2\\2&5\end{bmatrix}\)。高斯消元：第 2 行减去第 1 行的 \(\tfrac24=0.5\) 倍，得 \(U=\begin{bmatrix}4&2\\0&4\end{bmatrix}\)。把这个乘数 0.5 记在 \(L\) 的 \((2,1)\) 位置：\(L=\begin{bmatrix}1&0\\0.5&1\end{bmatrix}\)。验算 \(LU=\begin{bmatrix}4&2\\2&1+4\end{bmatrix}=A\)。所以 \(L\) 就是"消元时用过的乘数表"，\(U\) 是"消元的结果"。解 \(Ax=b=(8,13)^T\)：前代 \(Ly=b\) 得 \(y_1=8\)、\(y_2=13-0.5\cdot8=9\)；回代 \(Ux=y\) 得 \(x_2=9/4=2.25\)、\(x_1=(8-2\cdot2.25)/4=0.875\)。

更重要的是一次分解、多次求解：分解 \(A=LU\) 需要约 \(\frac23n^3\) 次浮点运算，之后每个新的右端 \(b\) 只需 \(O(n^2)\)。组合优化中用同一个协方差矩阵对多个期望收益向量求解、有限差分定价中每个时间步解同一个矩阵的方程组，都是这种情形。

3.12.2 LU 分解的存在性

定义 3.5.1　\(A=LU\)（\(L\) 下三角、\(U\) 上三角）称为 \(A\) 的 LU 分解。若 \(L\) 非奇异，总可以把 \(L\) 的对角元归一化为 1（单位下三角，unit lower triangular），把对角因子并入 \(U\)。

引理 3.5.2　若 \(A=LU\)，把三者都按左上角 \(k\times k\) 分块，则 \(A_{11}=L_{11}U_{11}\)。所以每个顺序主子矩阵也有 LU 分解，因子就是 \(L,U\) 的对应顺序主子矩阵。

定理 3.5.3　(a) \(A\) 有 \(L\) 非奇异的 LU 分解 ⇔ \(A\) 有行包含性质（row inclusion property）：对每个 \(i\)，第 \(i+1\) 行的前 \(i\) 个元素 \(A[\{i+1\},\{1..i\}]\) 是左上 \(i\times i\) 子矩阵 \(A[\{1..i\}]\) 各行的线性组合。 (b) \(A\) 有 \(U\) 非奇异的 LU 分解 ⇔ \(A\) 有对应的列包含性质。

必要性来自分块：\(A_{21}=L_{21}U_{11}=(L_{21}L_{11}^{-1})A_{11}\)。充分性用归纳，逐行构造 \(L\) 与 \(U\)。

白话解释："行包含性质"读起来拗口，意思其实是消元能进行下去。消元做到第 \(i+1\) 行时，要用前 \(i\) 行的组合把这一行的前 \(i\) 个元素消成 0。能做到的前提是：这一行的前 \(i\) 个元素，确实能由前 \(i\) 行的前 \(i\) 个元素线性组合出来。\(A_{[i]}\) 可逆时这一点自动成立（任何向量都能组合出来），所以非奇异情形下条件简化为"顺序主子式都非零"（推论 3.5.6）；只有 \(A_{[i]}\) 奇异时才需要真的检查。上面推导里的 \(L_{21}L_{11}^{-1}\) 就是那组组合系数。

推论 3.5.4　若 \(\operatorname{rank}A=k\) 且 \(A[\{1..j\}]\)（\(j=1,\dots,k\)）都非奇异，则 \(A\) 有 LU 分解，且任一因子可取单位三角；\(L,U\) 都非奇异 ⇔ \(k=n\)。

例 3.5.5　\(\begin{bmatrix}0&1\\1&0\end{bmatrix}\) 没有 LU 分解：若 \(A=LU\)，则 \(\ell_{11}u_{11}=a_{11}=0\)，\(L\) 或 \(U\) 奇异，与 \(A\) 非奇异矛盾。一般地，非奇异矩阵若有奇异的顺序主子矩阵，就没有 LU 分解。

两个提醒：

• 奇异矩阵的情况更微妙。原书给出一个 \(3\times3\) 矩阵，它既无行包含也无列包含性质，却有 LU 分解；而把它嵌入一个 \(4\times4\) 矩阵后就没有了。
• 即使要求 \(L\) 单位下三角，奇异矩阵的 LU 分解也可以不唯一：\(\begin{bmatrix}1&0\\a&1\end{bmatrix}\begin{bmatrix}0&1\\0&2-a\end{bmatrix}=\begin{bmatrix}0&1\\0&2\end{bmatrix}\) 对一切 \(a\) 成立。
• 3.5.P6：\((n,n)\) 元不影响 LU 分解是否存在。

3.12.3 LDU 分解与主元公式

推论 3.5.6　(a) 非奇异 \(A\) 有 LU 分解 ⇔ 所有顺序主子矩阵 \(A_{[i]}=A[\{1..i\}]\) 都非奇异。 (b) 此时 \(A=LDU\)，\(L\) 单位下三角，\(U\) 单位上三角，\(D=\operatorname{diag}(d_1,\dots,d_n)\)，且三者唯一，

即主元（pivot）= 相邻顺序主子式之比。证明：由引理 3.5.2，\(\det A_{[i]}=\det L_{[i]}\det D_{[i]}\det U_{[i]}=d_1\cdots d_i\)。

推导拆解：证明分三步。(1) 引理 3.5.2 说 \(A_{[i]}=L_{[i]}D_{[i]}U_{[i]}\)——左上角块只由三个因子的左上角块决定，因为三角矩阵右上（或左下）的零块让其余部分乘不进来。(2) 行列式可乘，且单位三角矩阵的行列式是对角元之积 \(=1\)，所以 \(\det A_{[i]}=\det D_{[i]}=d_1d_2\cdots d_i\)。(3) 相邻两式相除：\(d_i=\det A_{[i]}/\det A_{[i-1]}\)。 接上一个例子：\(\det A_{[1]}=4\)，\(\det A_{[2]}=16\)，所以 \(d_1=4\)、\(d_2=16/4=4\)，与消元得到的 \(U\) 的对角元 \(4,4\) 一致。这也说明唯一性：主元由 \(A\) 的子式完全决定，不依赖算法。

对称情形：若 \(A=A^T\) 且顺序主子式都非零，唯一性给出 \(U=L^T\)，即 \(A=LDL^T\)。若 \(A\) 是实对称正定的，所有 \(d_i>0\)，令 \(\tilde L=LD^{1/2}\) 得 Cholesky 分解 \(A=\tilde L\tilde L^T\)（第 02a 章由 QR 也推出过）。原书 3.5.P10 指出，对复对称矩阵（顺序主子矩阵非奇异）也有 \(A=LL^T\)，只是 \(L\) 可能是复的。

3.12.4 统计解读：主元就是逐次条件方差

设 \(\Sigma\) 是随机向量 \(x=(x_1,\dots,x_n)^T\) 的协方差矩阵（正定），\(\Sigma=LDL^T\)。

• \(d_i\) 是 \(x_i\) 对 \(x_1,\dots,x_{i-1}\) 做线性回归后的残差方差（条件方差）。理由：第 00 章 0.6 节，Schur 补就是条件协方差，而 \(\det\Sigma_{[i]}/\det\Sigma_{[i-1]}\) 正是 \(\Sigma_{[i]}\) 中 \((i,i)\) 元关于 \(\Sigma_{[i-1]}\) 的 Schur 补。
• 令 \(z=L^{-1}x\)，则 \(\operatorname{Cov}(z)=D\)：\(z\) 的分量互不相关，\(z_i\) 就是 \(x_i\) 的回归残差（"新息"）。所以 \(L^{-1}\) 第 \(i\) 行的相反数（去掉对角元）就是 \(x_i\) 对前 \(i-1\) 个变量的回归系数；\(L\) 的第 \(i\) 行是 \(x_i\) 在正交化新息 \(z_1,\dots,z_{i-1}\) 上的载荷。
• 这和第 02a 章"QR = 逐步正交化"是同一件事的两面：QR 对数据矩阵做 Gram–Schmidt，\(LDL^T\) 对协方差矩阵做同样的事。

推导拆解：用两资产情形把三条结论全部算出来。\(\Sigma=\begin{bmatrix}\sigma_1^2&\rho\sigma_1\sigma_2\\\rho\sigma_1\sigma_2&\sigma_2^2\end{bmatrix}\)。 第 1 步（消元）：第 2 行减去第 1 行的 \(\ell_{21}=\rho\sigma_1\sigma_2/\sigma_1^2=\rho\sigma_2/\sigma_1\) 倍。这个乘数正是 \(x_2\) 对 \(x_1\) 回归的斜率 \(\beta=\operatorname{Cov}/\operatorname{Var}(x_1)\)，也就是 CFA 里的最小方差对冲比率。 第 2 步（主元）：\(d_2=\sigma_2^2-\ell_{21}\cdot\rho\sigma_1\sigma_2=\sigma_2^2(1-\rho^2)\)，正是对冲后的残差方差。 第 3 步（新息）：\(L^{-1}=\begin{bmatrix}1&0\\-\ell_{21}&1\end{bmatrix}\)，\(z_2=x_2-\ell_{21}x_1\) 就是"多头 1 单位 \(x_2\)、空头 \(\ell_{21}\) 单位 \(x_1\)"的对冲组合收益，它与 \(x_1\) 不相关。 \(n\) 个资产时，同样的消元一列一列做下去，第 \(i\) 步就是"用前 \(i-1\) 个资产对冲第 \(i\) 个"。

量化用途：

• 风险归因的顺序分解：先放市场、再放行业、再放风格，\(d_i\) 给出每一步新增的"独立风险"；顺序不同，分解不同（与 QR 中列的顺序决定谁"吃掉"共同部分一样）。
• 蒙特卡洛：\(x=\tilde Lz\)（\(z\sim N(0,I)\)）生成相关正态向量。
• 对数似然：\(\log\det\Sigma=\sum\log d_i\)，\(x^T\Sigma^{-1}x=\|D^{-1/2}L^{-1}x\|^2\)，都不用显式求逆。
• 原书 3.5.P7 的一个漂亮例子：\(C_n=[1/\max\{i,j\}]\) 有 \(C_n=L_nL_n^T\)，\(L_n\) 的下三角元素 \(\ell_{ij}=1/i\)（\(i\ge j\)），所以 \(\det C_n=(1/n!)^2\)。

3.12.5 PLU 分解：任何矩阵都可以

不能 LU 分解的矩阵，先重排方程（行）就可以。

引理 3.5.7　若 \(A\in M_k\) 非奇异，存在置换矩阵 \(P\) 使 \(P^TA\) 的所有顺序主子式非零。（对 \(k\) 归纳：去掉最后一列，剩下 \(k-1\) 列无关，其中必有 \(k-1\) 行无关，把它们换到前面。）

定理 3.5.8（PLU 分解）　对每个 \(A\in M_n\)，存在置换矩阵 \(P\)、单位下三角 \(L\)、上三角 \(U\) 使 \(A=PLU\)。

解 \(Ax=b\) 时先解 \(Ly=P^Tb\)，再解 \(Ux=y\)。练习：每个 \(A\) 也可以写成 \(A=LUP\)（列置换）。

数值上，选主元不仅是为了存在性，更是为了稳定性（这是数值分析的常识，原书没有展开）。高斯消元若遇到很小的主元，乘数 \(\ell_{ik}=a_{ik}/a_{kk}\) 巨大，舍入误差被放大。部分选主元（partial pivoting）在每一步把当前列中绝对值最大的元素换到主元位置，保证 \(|\ell_{ik}|\le1\)。LAPACK 的 getrf（scipy.linalg.lu、numpy.linalg.solve 背后）都这样做，得到 \(PA=LU\) 型分解。实战 2 给出一个不选主元就完全算错的 \(2\times2\) 例子。对称正定矩阵做 Cholesky 不需要选主元，稳定且只需约 \(\frac13n^3\) 次运算。

白话解释：小主元的危害可以用"大数吃小数"理解。双精度浮点数只有约 16 位有效数字，\(10^{17}+1\) 存进去还是 \(10^{17}\)，那个 1 直接丢了。不选主元时，乘数 \(1/\epsilon\) 把第 2 行的原始信息乘到 \(10^{17}\) 量级，原来的"1"在加减中被吞掉，结果全错。选主元相当于先挑一个"块头大"的元素做除数，乘数不超过 1，就不会把误差放大。对称正定矩阵为什么不需要：它的主元 \(d_i\) 是条件方差，都是正数，而且 Cholesky 因子的每个元素都被对角元控制（\(\ell_{ij}^2\le a_{ii}\)），不会出现失控的大乘数。 （注：实战编号以文中实际为准——这个 \(2\times2\) 例子出现在实战 1 第 2 段，而不是实战 2。）

3.12.6 LPU 分解与三角等价（了解）

定理 3.5.11（LPU 分解）　对每个 \(A\in M_n\)，存在置换矩阵 \(P\)、单位下三角 \(L\)、上三角 \(U\) 使 \(A=LPU\)。若 \(A\) 非奇异，\(P\) 唯一。

构造：逐列消元，但每次用"该列中最上方的、尚未用过的非零行"作主元，只向下消元。\(n\) 步后得到一个"行置换后的上三角阵"。

唯一性的来源（3.5.10）：若 \(A=LBU\)，\(L,U\) 非奇异三角阵，则左上角 \(p\times q\) 子矩阵的秩相同：\(\operatorname{rank}A_{[p,q]}=\operatorname{rank}B_{[p,q]}\)。而置换矩阵 \(P\) 由所有 \(\operatorname{rank}P_{[p,q]}\) 唯一确定（3.5.P11）。

定义 3.5.12　若 \(A=LBU\)（\(L\) 非奇异下三角、\(U\) 非奇异上三角），称 \(A,B\) 三角等价。

定理 3.5.13　对非奇异 \(A,B\)，以下等价：(a) 二者三角等价于同一个（唯一的）置换矩阵；(b) 二者三角等价；(c) 对所有 \(p,q\)，\(\operatorname{rank}A_{[p,q]}=\operatorname{rank}B_{[p,q]}\)。所以置换矩阵是非奇异矩阵在三角等价下的标准形，左上角子矩阵的秩是完全不变量。

定理 3.5.14（LPDU 分解）　非奇异 \(A=LPDU\)，\(L\) 单位下三角，\(U\) 单位上三角，\(P\) 置换，\(D\) 非奇异对角；其中 \(P\) 和 \(D\) 唯一（\(L,U\) 不一定）。3.5.P3：每个非奇异矩阵单位三角等价于唯一的广义置换矩阵。

3.5.P12–P13 用 LPU 分解证明了第 00 章的互补零化度定律：\(A\) 非奇异、\(A^{-1}=[B_{ij}]\) 共形分块，则 \(\operatorname{nullity}A_{11}=\operatorname{nullity}B_{22}\)。

3.12.7 两个与算法相关的习题

3.5.P9（三对角矩阵）　\(A=\operatorname{tridiag}(-1,2,-1)\in M_n\)（二阶差分算子）：LU 分解中 \(L\) 的次对角元为 \(-\frac12,-\frac23,\dots,-\frac{n-1}n\)，\(U\) 的对角元为 \(2,\frac32,\dots,\frac{n+1}n\)、超对角元为 \(-1\)；\(\det A=n+1\)；特征值 \(\lambda_k=4\sin^2\frac{k\pi}{2(n+1)}\)，\(\lambda_1\to0\)、\(\lambda_n\to4\)，条件数随 \(n^2\) 增长。三对角矩阵的 LU 分解只需 \(O(n)\) 次运算（Thomas 算法），这是期权定价有限差分法每个时间步的核心计算。

3.5.P5（Lanczos 三对角化，了解）　若 Krylov 矩阵 \(X=[x\ Ax\ \cdots\ A^{n-1}x]\) 非奇异，则 \(X^{-1}AX\) 是 \(A\) 的特征多项式的友矩阵；用上三角矩阵修正后可得上 Hessenberg 形，再配合左 Krylov 矩阵的 LDU 分解得到三对角形；Hermite 情形给出三对角化算法。这是 Krylov 子空间方法（Lanczos、Arnoldi、共轭梯度）的代数源头，大规模稀疏协方差矩阵的特征计算用的就是这类方法。

量化实战

实战 1：协方差矩阵的 \(LDL^T\)：主元 = 条件方差，\(L^{-1}\) = 回归系数；选主元的必要性

场景：一个对冲组合包含股指、行业指数、个股、债券四类资产。按"股指 → 行业 → 个股 → 债券"的顺序分解风险，看每一步新增的独立方差。

import numpy as np
from scipy.linalg import ldl, lu

rng = np.random.default_rng(3)

# ---------- 1. LDL^T：主元 d_i = 顺序主子式之比 = 逐次条件方差 ----------
vols = np.array([0.20, 0.25, 0.30, 0.18])           # 年化波动：股指、行业、个股、债券
Corr = np.array([[1.0, 0.8, 0.6, -0.2],
                 [0.8, 1.0, 0.7, -0.1],
                 [0.6, 0.7, 1.0, 0.0],
                 [-0.2, -0.1, 0.0, 1.0]])
S = np.outer(vols, vols) * Corr
L, D, perm = ldl(S, lower=True)
print("置换（scipy 可能做对称选主元）:", perm, "  D 对角:", np.diag(D).round(6))
d = np.diag(D)
minors = [np.linalg.det(S[:k, :k]) for k in range(1, 5)]
ratios = [minors[0]] + [minors[k] / minors[k - 1] for k in range(1, 4)]
print("顺序主子式之比 det A_k / det A_{k-1}:", np.round(ratios, 6))
# 用 Cholesky 也能得到同样的 D：S = (L D^{1/2})(L D^{1/2})^T
C = np.linalg.cholesky(S)
print("Cholesky 对角元平方:", (np.diag(C) ** 2).round(6))

# 条件方差的回归解释：资产 k 对资产 1..k-1 回归的残差方差
T = 400_000
X = rng.multivariate_normal(np.zeros(4), S, size=T)
Linv = np.linalg.inv(L)
for k in range(1, 4):
    Z = X[:, :k]; y = X[:, k]
    b = np.linalg.lstsq(Z, y, rcond=None)[0]
    print("资产 %d | 前 %d 个：残差方差 %.6f vs d_%d = %.6f；回归系数 %s vs -L^{-1} 第 %d 行 %s"
          % (k + 1, k, np.var(y - Z @ b), k + 1, d[k], b.round(4), k + 1, (-Linv[k, :k]).round(4)))

# ---------- 2. 选主元的必要性 ----------
A0 = np.array([[0.0, 1.0], [1.0, 0.0]])
P, Lp, Up = lu(A0)
print("\n[[0,1],[1,0]] 没有 LU，但有 PLU：P =", P.ravel(), " L =", Lp.ravel(), " U =", Up.ravel())

def lu_nopivot(A):
    A = A.astype(float).copy(); n = len(A); L = np.eye(n)
    for k in range(n - 1):
        for i in range(k + 1, n):
            L[i, k] = A[i, k] / A[k, k]
            A[i, k:] -= L[i, k] * A[k, k:]
    return L, np.triu(A)

eps = 1e-17
Ae = np.array([[eps, 1.0], [1.0, 1.0]])
b = np.array([1.0, 2.0])
x_true = np.linalg.solve(Ae, b)
L1, U1 = lu_nopivot(Ae)
x_np = np.linalg.solve(U1, np.linalg.solve(L1, b))
P2, L2, U2 = lu(Ae)
x_pv = np.linalg.solve(U2, np.linalg.solve(L2, P2.T @ b))
print("ε=1e-17：不选主元 L =", L1.ravel(), "U =", U1.ravel())
print("  不选主元的解:", x_np, "  选主元的解:", x_pv, "  真解:", x_true.round(12))

关键输出：

置换（scipy 可能做对称选主元）: [0 1 2 3]   D 对角: [0.04     0.0225   0.0455   0.030335]
顺序主子式之比 det A_k / det A_{k-1}: [0.04     0.0225   0.0455   0.030335]
Cholesky 对角元平方: [0.04     0.0225   0.0455   0.030335]
资产 2 | 前 1 个：残差方差 0.022445 vs d_2 = 0.022500；回归系数 [1.0003] vs -L^{-1} 第 2 行 [1.]
资产 3 | 前 2 个：残差方差 0.045565 vs d_3 = 0.045500；回归系数 [0.1666 0.7327] vs -L^{-1} 第 3 行 [0.1667 0.7333]
资产 4 | 前 3 个：残差方差 0.030232 vs d_4 = 0.030335；回归系数 [-0.3156  0.0462  0.0994] vs -L^{-1} 第 4 行 [-0.3165  0.0475  0.0989]

[[0,1],[1,0]] 没有 LU，但有 PLU：P = [0. 1. 1. 0.]  L = [1. 0. 0. 1.]  U = [1. 0. 0. 1.]
ε=1e-17：不选主元 L = [1.e+00 0.e+00 1.e+17 1.e+00] U = [ 1.e-17  1.e+00  0.e+00 -1.e+17]
  不选主元的解: [0. 1.]   选主元的解: [1. 1.]   真解: [1. 1.]

读法：

1. 三种算法（scipy.linalg.ldl、顺序主子式之比、Cholesky 对角元平方）给出同一组主元，印证了 \(d_i=\det A_{[i]}/\det A_{[i-1]}\)。
2. 行业指数的方差 \(0.25^2=0.0625\)，对股指回归后剩 \(0.0225\)（\(=0.0625(1-0.8^2)\)），即 64% 的行业风险可由股指对冲；个股在股指、行业之后剩 0.0455（原方差 0.09 的一半）；债券与权益相关弱，剩余方差 0.0303 接近原值 0.0324。
3. 40 万个模拟样本的回归残差方差与 \(d_i\) 吻合，回归系数与 \(-L^{-1}\) 的行吻合（差异是抽样误差）。\(LDL^T\) 一次就给出全部"逐步对冲比率"和"对冲后残差风险"，不需要做任何回归。
4. \(\begin{bmatrix}0&1\\1&0\end{bmatrix}\) 没有 LU 分解，PLU 只是交换两行。\(\epsilon=10^{-17}\) 的矩阵不选主元时乘数达到 \(10^{17}\)，\(U\) 的 \((2,2)\) 元 \(1-10^{17}\) 在双精度下变成 \(-10^{17}\)，"1"被完全吞掉，解从 \((1,1)\) 变成了 \((0,1)\)；部分选主元先交换两行，结果正确。

实战 2：三对角 LU 与 Black–Scholes 有限差分定价

场景：用 Crank–Nicolson 格式解 Black–Scholes 偏微分方程。在对数价格 \(x=\ln S\) 与剩余期限 \(\tau\) 下，

空间离散后算子 \(\mathcal L\) 是三对角矩阵，每个时间步解 \((I-\theta\Delta\tau\mathcal L)V^{n+1}=(I+(1-\theta)\Delta\tau\mathcal L)V^n+\text{边界项}\)。先用原书 3.5.P9 的二阶差分矩阵验证三对角 LU 的结构。

import numpy as np, time
from scipy.linalg import solve_banded, lu_factor, lu_solve, lu
from scipy.stats import norm

# ---------- 1. 原书 3.5.P9：tridiag(-1, 2, -1) 的 LU、行列式与特征值 ----------
n = 6
A = 2 * np.eye(n) - np.eye(n, k=1) - np.eye(n, k=-1)
P, L, U = lu(A)
print("是否发生行交换:", not np.allclose(P, np.eye(n)))
print("L 的次对角元:", np.diag(L, -1).round(4), " 理论 -k/(k+1):", np.round([-k / (k + 1) for k in range(1, n)], 4))
print("U 的对角元  :", np.diag(U).round(4), " 理论 (k+1)/k:", np.round([(k + 1) / k for k in range(1, n + 1)], 4))
print("det = %.4f（理论 n+1 = %d）" % (np.linalg.det(A), n + 1))
lam = np.sort(np.linalg.eigvalsh(A))
print("特征值:", lam.round(4), " 理论 4sin^2(kπ/2(n+1)):", np.sort(4 * np.sin(np.arange(1, n + 1) * np.pi / (2 * (n + 1))) ** 2).round(4))
for m in [10, 100, 1000]:
    l = 4 * np.sin(np.array([1, m]) * np.pi / (2 * (m + 1))) ** 2
    print("  n=%4d  条件数 λ_max/λ_min = %.1f" % (m, l[1] / l[0]))

# ---------- 2. Black–Scholes 隐式差分：每个时间步解一个三对角方程 ----------
S0, K, r, sig, Tm = 100.0, 100.0, 0.03, 0.25, 1.0
M, N = 800, 400                                    # 空间网格、时间步
x = np.linspace(np.log(S0) - 5 * sig, np.log(S0) + 5 * sig, M + 1)
dx, dt = x[1] - x[0], Tm / N
V = np.maximum(np.exp(x) - K, 0.0)
a = 0.5 * sig ** 2 / dx ** 2; b = (r - 0.5 * sig ** 2) / (2 * dx)
lo, di, up = a - b, -2 * a - r, a + b               # 算子 L 的三条对角线
theta = 0.5                                        # Crank–Nicolson
ab = np.zeros((3, M - 1))                          # (I - θ dt L) 的带状存储
ab[0, 1:] = -theta * dt * up
ab[1, :] = 1 - theta * dt * di
ab[2, :-1] = -theta * dt * lo
t0 = time.perf_counter()
for j in range(1, N + 1):
    tau = j * dt
    Vi = V[1:-1]
    rhs = Vi + (1 - theta) * dt * (lo * V[:-2] + di * Vi + up * V[2:])
    right_new = np.exp(x[-1]) - K * np.exp(-r * tau)
    rhs[-1] += theta * dt * up * right_new          # 边界条件进入右端
    V[1:-1] = solve_banded((1, 1), ab, rhs)
    V[0], V[-1] = 0.0, right_new
t_band = time.perf_counter() - t0
price_fd = np.interp(np.log(S0), x, V)
d1 = (np.log(S0 / K) + (r + 0.5 * sig ** 2) * Tm) / (sig * np.sqrt(Tm)); d2 = d1 - sig * np.sqrt(Tm)
price_bs = S0 * norm.cdf(d1) - K * np.exp(-r * Tm) * norm.cdf(d2)
print("\nCrank–Nicolson 价格 = %.4f   Black–Scholes 公式 = %.4f   误差 = %.1e" % (price_fd, price_bs, price_fd - price_bs))

# 同一个三对角矩阵：带状求解 O(M) vs 稠密 LU 分解一次 + 反复回代 vs 每步稠密求解
Ad = np.diag(ab[1]) + np.diag(ab[0, 1:], 1) + np.diag(ab[2, :-1], -1)
rhs = np.ones(M - 1)
t0 = time.perf_counter(); [np.linalg.solve(Ad, rhs) for _ in range(N)]; t_dense = time.perf_counter() - t0
t0 = time.perf_counter(); lf = lu_factor(Ad); [lu_solve(lf, rhs) for _ in range(N)]; t_lu = time.perf_counter() - t0
print("%d 次求解耗时：带状 %.3fs（含全部定价循环）  稠密 LU 一次分解+回代 %.3fs  每次稠密 solve %.3fs"
      % (N, t_band, t_lu, t_dense))

关键输出（耗时与机器有关）：

是否发生行交换: False
L 的次对角元: [-0.5    -0.6667 -0.75   -0.8    -0.8333]  理论 -k/(k+1): [-0.5    -0.6667 -0.75   -0.8    -0.8333]
U 的对角元  : [2.     1.5    1.3333 1.25   1.2    1.1667]  理论 (k+1)/k: [2.     1.5    1.3333 1.25   1.2    1.1667]
det = 7.0000（理论 n+1 = 7）
特征值: [0.1981 0.753  1.555  2.445  3.247  3.8019]  理论 4sin^2(kπ/2(n+1)): [0.1981 0.753  1.555  2.445  3.247  3.8019]
  n=  10  条件数 λ_max/λ_min = 48.4
  n= 100  条件数 λ_max/λ_min = 4133.6
  n=1000  条件数 λ_max/λ_min = 406095.0

Crank–Nicolson 价格 = 11.3483   Black–Scholes 公式 = 11.3485   误差 = -1.8e-04
400 次求解耗时：带状 0.007s（含全部定价循环）  稠密 LU 一次分解+回代 0.074s  每次稠密 solve 1.026s

读法：

1. 二阶差分矩阵的 LU 因子、行列式、特征值都与原书 3.5.P9 的公式一致；它对角占优，部分选主元不会触发行交换。条件数大致按 \(n^2\) 增长（\(n\) 扩大 10 倍，条件数扩大约 100 倍），网格越细，方程越病态——不过三对角 LU 对这类对角占优矩阵是稳定的。
2. Crank–Nicolson 价格与 Black–Scholes 公式相差 \(1.8\times10^{-4}\)。同样的框架可以直接处理美式期权（每步取 \(\max\)）、局部波动率（系数随 \(x\) 变化）等没有闭式解的情形。
3. 速度：400 个时间步，带状求解（\(O(M)\)）全部定价只需几毫秒；稠密 LU 只分解一次、每步回代（\(O(M^2)\)）慢一个数量级；每步都重新做稠密求解（\(O(M^3)\)）再慢一个多数量级。"利用结构 + 分解一次多次回代"是数值实现的基本功。

实战 3：Weyr 形、中心化子与斜投影

import numpy as np
from scipy.linalg import block_diag, null_space

rng = np.random.default_rng(8)

def jordan_block(lam, k):
    return lam * np.eye(k) + np.diag(np.ones(k - 1), 1)

def weyr_block(w, lam=0.0):
    """W(w, λ)：对角块 λI_{w_i}，超对角块 G = [I; 0]"""
    n = sum(w); W = lam * np.eye(n); off = np.cumsum([0] + list(w))
    for i in range(len(w) - 1):
        W[off[i]:off[i] + w[i + 1], off[i + 1]:off[i + 2]] = np.eye(w[i + 1])
    return W

# ---------- 1. 原书 3.4.P8：J = J3(0)⊕J2(0) 与 Weyr 形置换相似 ----------
J = block_diag(jordan_block(0, 3), jordan_block(0, 2))
W = weyr_block([2, 2, 1])
print("Weyr 形 W_J(0) =\n", W.astype(int))
perm = [0, 2, 4, 1, 3]                       # σ = (1,3,5,2,4)
P = np.eye(5)[:, perm]
print("J = P^T W P ?", np.allclose(J, P.T @ W @ P))
print("rank W^p:", [int(np.linalg.matrix_rank(np.linalg.matrix_power(W, p))) for p in range(4)],
      " rank J^p:", [int(np.linalg.matrix_rank(np.linalg.matrix_power(J, p))) for p in range(4)])

# ---------- 2. 中心化子维数：dim{B: AB = BA} = Σ w_i^2（原书 3.4.P3–P4） ----------
def centralizer(A):
    n = A.shape[0]
    K = np.kron(np.eye(n), A) - np.kron(A.T, np.eye(n))   # vec(AB - BA) = K vec(B)
    return null_space(K)

blocks = [jordan_block(0, 3)] * 2 + [jordan_block(0, 2)] * 3 + [jordan_block(0, 1)]
J13 = block_diag(*blocks)
W13 = weyr_block([6, 5, 2])
print("\n(3.1.16a) 的中心化子维数: Jordan 形 %d, Weyr 形 %d, 理论 6²+5²+2² = %d"
      % (centralizer(J13).shape[1], centralizer(W13).shape[1], 36 + 25 + 4))
# 与 Weyr 形交换的随机矩阵：按 (6,5,2) 分块后是块上三角
Bw = (centralizer(W13) @ rng.normal(size=centralizer(W13).shape[1])).reshape(13, 13, order="F")
print("‖W B - B W‖ = %.1e；按 (6,5,2) 分块的下三角块范数: (2,1) %.1e  (3,1) %.1e  (3,2) %.1e"
      % (np.linalg.norm(W13 @ Bw - Bw @ W13), np.linalg.norm(Bw[6:11, :6]), np.linalg.norm(Bw[11:, :6]),
         np.linalg.norm(Bw[11:, 6:11])))
Bj = (centralizer(J13) @ rng.normal(size=centralizer(J13).shape[1])).reshape(13, 13, order="F")
print("与 Jordan 形交换的随机矩阵是上三角吗？", np.allclose(np.tril(Bj, -1), 0))
print("非减次矩阵（单个 J_13(0)）的中心化子维数:", centralizer(jordan_block(0, 13)).shape[1])

# ---------- 3. 幂等矩阵的酉标准形：正交投影 vs 斜投影（IV 估计） ----------
T = 200
z = rng.normal(size=(T, 1))
for name, strength in [("强工具", 1.0), ("弱工具", 0.1)]:
    xx = strength * z + rng.normal(size=(T, 1))
    H = xx @ np.linalg.solve(xx.T @ xx, xx.T)        # OLS 帽子矩阵：正交投影
    Pv = xx @ np.linalg.solve(z.T @ xx, z.T)         # 恰好识别 IV：斜投影 X (Z^T X)^{-1} Z^T
    sH = np.linalg.svd(H, compute_uv=False)[:2]; sP = np.linalg.svd(Pv, compute_uv=False)[:2]
    print("%s：H 幂等 %s 奇异值 %s；IV 投影幂等 %s 最大奇异值 %.3f = 1/|cos∠(x,z)| = %.3f"
          % (name, np.allclose(H @ H, H), sH.round(3), np.allclose(Pv @ Pv, Pv), sP[0],
             1 / abs((xx.T @ z).item() / np.linalg.norm(xx) / np.linalg.norm(z))))

关键输出：

Weyr 形 W_J(0) =
 [[0 0 1 0 0]
 [0 0 0 1 0]
 [0 0 0 0 1]
 [0 0 0 0 0]
 [0 0 0 0 0]]
J = P^T W P ? True
rank W^p: [5, 3, 1, 0]  rank J^p: [5, 3, 1, 0]

(3.1.16a) 的中心化子维数: Jordan 形 65, Weyr 形 65, 理论 6²+5²+2² = 65
‖W B - B W‖ = 7.5e-15；按 (6,5,2) 分块的下三角块范数: (2,1) 5.4e-15  (3,1) 1.4e-15  (3,2) 6.3e-16
与 Jordan 形交换的随机矩阵是上三角吗？ False
非减次矩阵（单个 J_13(0)）的中心化子维数: 13
强工具：H 幂等 True 奇异值 [1. 0.]；IV 投影幂等 True 最大奇异值 1.405 = 1/|cos∠(x,z)| = 1.405
弱工具：H 幂等 True 奇异值 [1. 0.]；IV 投影幂等 True 最大奇异值 8.230 = 1/|cos∠(x,z)| = 8.230

读法：

1. 用 Young 图读出的置换 \(\sigma=(1,3,5,2,4)\) 确实把 Weyr 形变成 Jordan 形；两者的秩序列相同。
2. 中心化子维数 65 与公式 \(\sum w_i^2\) 一致，相似的 \(J\) 与 \(W\) 当然给出相同的维数。与 Weyr 形交换的随机矩阵按 \((6,5,2)\) 分块后下三角块为零（Belitskii 定理），与 Jordan 形交换的矩阵则不是上三角的。非减次矩阵的中心化子维数恰为 \(n=13\)（等号情形），这时所有与之交换的矩阵都是它的多项式。
3. OLS 帽子矩阵是正交投影，非零奇异值为 1；恰好识别 IV 的投影矩阵 \(X(Z^TX)^{-1}Z^T\) 幂等但不对称，是斜投影，最大奇异值等于 \(1/|\cos\angle(x,z)|\)。工具变弱（\(x\) 与 \(z\) 几乎正交）时这个数从 1.4 跳到 8.2：估计量把噪声放大的倍数随之增大。这是"弱工具问题"的一个线性代数侧面，计量细节见第 05 册。

本章小结

Weyr 标准形用 Weyr 特征构造：每个特征值对应一个分块双对角的 Weyr 块，对角块是大小依次为 \(w_1\ge w_2\ge\cdots\) 的纯量矩阵，超对角块是列满秩的 \(\begin{bmatrix}I\\0\end{bmatrix}\)；它与 Jordan 形信息相同、置换相似，但与它交换的矩阵都是分块上三角的（Belitskii），因而交换族可以同时相似到"一个 Weyr 形 + 其余上三角"（O'Meara–Vinsonhaler），这是 Jordan 形做不到的；中心化子维数为 \(\sum w_i^2\ge n\)，等号恰为非减次。只允许酉相似时，Littlewood 定理给出酉 Weyr 形，由它推出投影矩阵的酉标准形：正交投影的非零奇异值全为 1，斜投影由大于 1 的奇异值刻画。三角分解部分：LU 分解存在 ⇔ 行包含性质，非奇异时 ⇔ 所有顺序主子式非零；此时 LDU 唯一，主元 \(d_i=\det A_{[i]}/\det A_{[i-1]}\)；对称正定时 \(\Sigma=LDL^T=\tilde L\tilde L^T\)，\(d_i\) 是逐次条件方差，\(L^{-1}\) 给出逐步回归系数。任何矩阵都有 PLU 与 LPU 分解；数值上必须选主元。非奇异矩阵在三角等价下的标准形是置换矩阵，左上角子矩阵的秩是完全不变量。三对角矩阵的 LU 只需 \(O(n)\)，是有限差分定价的核心。

练习

基础

1. 写出 \(J=J_4(\lambda)\oplus J_2(\lambda)\) 的 Weyr 特征与 Weyr 形，并用 Young 图求出把 Weyr 形变为 \(J\) 的置换。 答案要点：\(w=(2,2,1,1)\)；点图逐行标号 \(1,2/3,4/5/6\)，逐列读出 \((1,3,5,6,2,4)\)。
2. 求 \(A=\begin{bmatrix}4&2&2\\2&5&1\\2&1&6\end{bmatrix}\) 的 \(LDL^T\) 分解，并用顺序主子式验证主元。 答案要点：\(d_1=4\)；\(\det A_{[2]}=16\)，\(d_2=4\)；\(\det A=4(30-1)-2(12-2)+2(2-10)=116-20-16=80\)，\(d_3=80/16=5\)。
3. 说明 \(\begin{bmatrix}1&2\\2&4\end{bmatrix}\) 有 LU 分解但 \(U\) 奇异，\(\begin{bmatrix}0&1\\1&1\end{bmatrix}\) 没有 LU 分解；分别写出它们的 PLU 分解。
4. 设 \(\Sigma\) 是两只股票收益的协方差，\(\sigma_1=0.2\)、\(\sigma_2=0.3\)、\(\rho=0.6\)。求 \(LDL^T\)，解释 \(d_2\) 与最优对冲比率。 答案要点：\(\ell_{21}=\rho\sigma_2/\sigma_1=0.9\)（最优对冲比率），\(d_2=\sigma_2^2(1-\rho^2)=0.0576\)（对冲后残差方差）。
5. 证明幂等矩阵 \(A\) 是正交投影（\(A=A^*\)）当且仅当它的所有非零奇异值都等于 1。 提示：用推论 3.4.3.3 的酉标准形。

进阶

6. 证明：若 \(\Sigma=LDL^T\) 且 \(z=L^{-1}x\)，则 \(\operatorname{Cov}(z)=D\)，并由此证明 \(d_i\) 等于 \(x_i\) 对 \(x_1,\dots,x_{i-1}\) 的线性回归残差方差。 提示：\(x_i=z_i+\sum_{j<i}\ell_{ij}z_j\)，而 \(z_1,\dots,z_{i-1}\) 与 \(x_1,\dots,x_{i-1}\) 张成同一空间，\(z_i\) 与它们不相关。
7. 用 3.5.P7 的结论计算 \(C_4=[1/\max\{i,j\}]\) 的行列式，并用数值验证。 答案要点：\((1/4!)^2=1/576\)。
8. 证明：对角占优的三对角矩阵（\(|a_{ii}|>|a_{i,i-1}|+|a_{i,i+1}|\)）不选主元做 LU 时，所有主元非零且 \(|\ell_{i,i-1}|<1\)。 提示：对 \(i\) 归纳，证明 \(|u_{ii}|>|a_{i,i+1}|\)。
9. 设 \(A=\begin{bmatrix}1&s\\0&0\end{bmatrix}\)。求其奇异值，并说明它是沿哪个方向投影到哪条直线上；\(s\to\infty\) 时发生什么。 答案要点：奇异值 \(\sqrt{1+s^2}\) 与 0；投影到 \(e_1\) 轴，沿 \(\ker A=\operatorname{span}\{(s,-1)^T\}\) 方向；\(s\to\infty\) 时投影方向趋于与 \(e_1\) 轴平行，投影无界放大。
10. 编程：用 Cholesky 分解实现多元正态对数似然 \(-\frac12(\log\det\Sigma+x^T\Sigma^{-1}x+n\log2\pi)\)，不显式求逆；与 scipy.stats.multivariate_normal.logpdf 对比。

原书推荐习题

• 3.4.P1：\(A^2=-I\) 的实标准形（复结构）。
• 3.4.P4：中心化子维数 \(\sum w_i^2\)，Weyr 与 Segre 特征的联系。
• 3.4.P5：平方零矩阵的酉标准形。
• 3.4.P7：Jordan 形不能替代 Weyr 形的反例。
• 3.4.P8：用 Young 表构造 Weyr 形与 Jordan 形之间的置换相似。
• 3.5.P4：用高斯消元得到 LU。
• 3.5.P5：Lanczos 三对角化，与 Krylov 子空间方法相关。
• 3.5.P9：差分矩阵的 LU 与特征值（条件数直观）。
• 3.5.P10：对称矩阵的 \(LL^T\) 分解。
• 3.5.P11、P13：LPU 中置换矩阵的唯一性与互补零化度定律。

原书对照

第 3 章到此结束。原书第 4 章转向 Hermite 矩阵与对称矩阵：特征值的变分刻画（Courant–Fischer）、Weyl 不等式、交错定理、优超，这些是协方差矩阵估计误差分析和组合优化约束的理论基础。