第 08b 章 本原矩阵、随机矩阵与 Birkhoff 定理

对应原书：Horn & Johnson《Matrix Analysis》第 2 版，第 8 章的 8.5 节 Primitive matrices、8.6 节 A general limit theorem、8.7 节 Stochastic and doubly stochastic matrices（书 p.540–553，PDF p.560–573）。Perron–Frobenius 定理与不可约性见第 08a 章。

上一章的 Perron–Frobenius 定理给出了不可约非负矩阵的正 Perron 向量，但没有保证 \((A/\rho)^m\) 收敛——周期性会让幂永远振荡。这一章补上最后一块：本原性（不可约 + 非周期）恰好是收敛的条件，并且可以纯粹从零/非零模式判定。随后把理论落到最常用的一类非负矩阵上：随机矩阵（Markov 链转移矩阵）与双随机矩阵。对量化读者，这一章回答：评级迁移矩阵多期推算会收敛到什么、多快；regime 模型的平稳分布与持续期；PageRank 为什么要加阻尼；以及分配类组合优化问题为什么有整数最优解（Birkhoff 定理）。

学习目标

读完本章，你应当能够：

1. 理解本原矩阵的定义与三种等价刻画：唯一最大模特征值、某个幂为正、各节点闭路长度的最大公约数为 1。
2. 掌握本原指数的上界（Wielandt 界 \(n^2-2n+2\)、最短圈界、Holladay–Varga 界），会用反复平方的布尔运算高效判定本原性。
3. 理解不可约非负矩阵的 Cesàro 极限定理与 \(O(1/N)\) 收敛速度，以及它与 Markov 链遍历定理的关系。
4. 掌握随机矩阵与双随机矩阵的基本性质，理解 Birkhoff 定理的构造性证明，以及"凸函数在双随机矩阵集合上的最大值在置换矩阵处取得"。
5. 在量化中应用：评级迁移与 regime 链的多期推算、平稳分布与混合速度、PageRank 式中心性、分配问题的整数性与 von Neumann 迹不等式。

读前导读

这一章在解决什么问题

这一章分成互相独立的两半。

前一半（8.6–8.7）回答：Markov 链什么时候"忘掉"初始状态？ 上一章说不可约的非负矩阵有正的 Perron 向量，但幂不一定收敛：如果状态只能"牛 → 熊 → 牛 → 熊"地交替跳，分布会永远来回摆。本章给出精确条件：除了互通（不可约），还要"非周期"，两者合起来叫本原。本原的链从任何初始状态出发，最终都收敛到同一个平稳分布，速度由第二大特征值的模 \(|\lambda_2|\) 决定。你做 regime 模型（例如隐马尔可夫的牛/熊/震荡）、评级迁移的多年推算、PageRank 式的网络中心性时，用的都是这套结论。CFA 信用部分的"累计违约概率"和"违约前期望年数"，本章给出它们的矩阵公式。

后一半（8.8）是关于双随机矩阵的组合优化结论。 双随机矩阵是行和、列和都为 1 的非负矩阵，可以看作"软分配"：第 \(i\) 笔订单按比例分给各个通道。Birkhoff 定理说每个软分配都是若干个"硬分配"（一一对应，即置换矩阵）的加权平均。推论很实用：线性目标的分配问题，放松成连续变量求解，最优解自动是整数解。你在 CFA 里见过的线性规划"最优解在可行域的角点上"，这里的角点恰好是置换矩阵。

需要先想起来的数学

1. Markov 链与转移矩阵。 \(P_{ij}\) 是"今天在状态 \(i\)、明天在状态 \(j\)"的条件概率，每行和为 1。分布用行向量写：\(\pi_{t+1}^T=\pi_t^TP\)。两期转移是 \(P^2\)，因为"今天 \(i\) 后天 \(j\)"要对中间状态 \(k\) 求和 \(\sum_kP_{ik}P_{kj}\)（全概率公式）。见 第 00 册第 07 章 概率中的分析工具。

2. 几何级数。 \(\sum_{k\ge0}q^k=1/(1-q)\)，\(|q|<1\)。矩阵版 \(\sum_kQ^k=(I-Q)^{-1}\) 用于吸收链。例：每年违约概率恒为 5%，存活到第 \(k\) 年的概率 \(0.95^k\)，期望存活年数 \(\sum_k0.95^k=1/0.05=20\)。见 第 00 册第 04 章 级数与收敛。

3. 最大公约数 gcd。 \(\gcd(4,6)=2\)，\(\gcd(3,4)=1\)。本章用它判断周期：回到起点的步数如果全是 2 的倍数，链就有周期 2。

4. 凸组合与极点。 \(t_1P_1+\cdots+t_NP_N\)（\(t_i\ge0\)、\(\sum t_i=1\)）叫凸组合，即加权平均。"极点"是凸集中不能写成其他点平均的点，如三角形的三个顶点。线性函数在凸多面体上的最大、最小值一定在某个顶点取到。见本册 附录 A 与 第 00 册第 05 章 多元微积分与优化。

5. 大 O 记号。 "误差为 \(O(1/N)\)"指误差不超过某个常数乘以 \(1/N\)：\(N\) 增大 10 倍，误差缩小约 10 倍。与之对比，"几何速度"\(|\lambda_2|^N\) 每多一步就乘一个小于 1 的因子，快得多。见 第 00 册第 07 章。

本章符号速查。 \(\gamma(A)\) 是本原指数。\(e\) 是全 1 向量，\(\pi\) 是平稳分布（不是圆周率）。\(\sigma\) 在 Birkhoff 部分表示一个置换（把 \(1,\dots,n\) 重新排列的一一对应），\(a_{1\sigma(1)}\cdots a_{n\sigma(n)}\) 是沿这个置换取出的 \(n\) 个元素之积，每行每列各取一个；在 von Neumann 部分 \(\sigma_i(A)\) 又是奇异值，注意区分。\(\lceil x\rceil\) 是向上取整。\(\operatorname{Re}\) 是复数实部，\(\operatorname{tr}\) 是迹（对角元之和）。\(a^\downarrow\) 是把向量分量从大到小排好。

怎么读这一章

核心必读：8.6.1（本原的定义与三种刻画）、8.7 的直观段落、8.8.1 全部（随机矩阵、Markov 链字典、吸收链）、8.8.2 的 Birkhoff 定理陈述与 Corollary 8.7.4。实战 1 和实战 3 是本章对量化最直接的价值。

第一次可以只看结论：8.6.2 本原指数的各种上界（知道"判本原只看零模式，并且有现成上界"即可）、8.6.3–8.6.4、8.7 的证明、8.8.3 von Neumann 迹定理的证明（结论要记住，第 07b 章的最佳逼近都靠它）、8.8.4 的习题罗列。

建议顺序：8.6.1 → 8.8.1 → 实战 1 → 实战 3 → 8.7 → 8.8.2 → 实战 4 → 其余。

8.6 本原矩阵

8.6.1 定义与刻画

Definition 8.5.0。 非负矩阵 \(A\) 称为本原的（primitive），若它不可约，且只有一个最大模特征值（即 \(\rho(A)\) 本身）。

第 08a 章证明 Perron 定理极限 \((A/\rho)^m\to xy^T\) 时只用到两点：\(\rho\) 单重且有正的左右特征向量、其余特征值的模严格小于 \(\rho\)。本原矩阵恰好满足这两点，所以：

Theorem 8.5.1。 \(A\) 非负本原，\(x,y\) 为右、左 Perron 向量（\(y^Tx=1\)），则

至此 Perron 定理的全部结论从正矩阵推广到了本原矩阵。

Theorem 8.5.2（Frobenius）。 非负 \(A\) 本原 ⇔ 存在 \(m\ge1\) 使 \(A^m>0\)。

证明：若 \(A^m>0\)，则图中任意两点间有长 \(m\) 的路，强连通，不可约；再对 \(A^m\) 用 Perron 定理，\(A\) 只有一个最大模特征值。反过来由 8.5.1，\((A/\rho)^m\to xy^T>0\)，充分大的 \(m\) 必使 \(A^m>0\)。练习：不可约且 \(A^m>0\) ⇒ 对所有 \(p>m\) 都有 \(A^p>0\)（不可约矩阵每列有非零元）。

白话解释："某个幂全正"的意思是：存在一个固定步数 \(m\)，从任何状态出发、恰好走 \(m\) 步，都有正概率到达任何状态。不可约只要求"总能走到"，不要求"同一个步数都能走到"。 用 regime 链举例。牛、熊两个状态，如果每期必然切换（\(P=\begin{bmatrix}0&1\\1&0\end{bmatrix}\)），从牛出发奇数步一定在熊，偶数步一定在牛，没有哪个步数能同时到达两者，所以不本原，分布永远在 \((1,0)\) 和 \((0,1)\) 之间跳。只要牛市有一点"留在原地"的概率（对角元为正），一步之后就可能在牛也可能在熊，两步之后四个方向都有正概率，\(P^2>0\)，本原。

Theorem 8.5.3（图论判据）。 \(A\) 不可约非负，\(L_i\) 为从节点 \(i\) 出发回到 \(i\) 的所有闭路长度，\(g_i=\gcd(L_i)\)。则 \(A\) 本原 ⇔ 所有 \(g_i=1\)。

证明：本原时 \(A^k>0\) 对所有大 \(k\) 成立，相邻长度的闭路都存在，\(g_i=1\)。若不本原、有 \(k>1\) 个最大模特征值，由第 08a 章 Corollary 8.4.7，\(k\nmid m\) 时 \(A^m\) 对角元全为 0，即不存在长度不被 \(k\) 整除的闭路，\(g_i\ge k\)。（原书此处引用编号写作 8.4.8，按内容应为 8.4.7。）

推导拆解："gcd 为 1 ⇒ 大步数的闭路都存在"这一步，证明里一笔带过，用一个例子补上。设节点 \(i\) 有长 3 和长 4 的两条闭路。绕前者 \(a\) 圈、后者 \(b\) 圈，得到长 \(3a+4b\) 的闭路。\(3a+4b\) 能凑出 \(3,4,6,7,8,9,10,\dots\)，从 6 起每个整数都能凑出（6=3+3，7=3+4，8=4+4，之后每个数加 3 即可）。这是"两个互素整数能凑出所有足够大的整数"的一般事实。所以足够大的 \(m\) 都有长 \(m\) 的闭路，\((A^m)_{ii}>0\)；再借助不可约把"回到 \(i\)"扩展到"到达任何 \(j\)"。 反过来，如果闭路长度只有 3 和 6（gcd 为 3），凑出的永远是 3 的倍数，长 4、5、7 的闭路不存在，对应的 \(A^m\) 对角元为 0，永远不会全正。实战 2 第 4 部分正是这两个例子。

Romanovsky 定理。 不可约非负矩阵的所有 \(g_i\) 相等，且恰等于最大模特征值的个数 \(k\)——周期（period）。Markov 链语言：不可约链要么非周期（本原），要么周期为 \(k\)；周期为 \(k\) 的链只能在 \(k\) 的倍数步回到出发点。

8.5.P17：本原性只依赖零元的位置，与非零元的大小无关。8.5.P18：对称不可约非负矩阵本原 ⇔ \(-\rho(A)\) 不是特征值；对无向图邻接矩阵，本原 ⇔ 连通且不是二部图。

8.6.2 本原指数

使 \(A^k>0\) 的最小 \(k\) 叫本原指数（index of primitivity），记 \(\gamma(A)\)。

Lemma 8.5.4。 不可约且对角元全为正 ⇒ \(A^{n-1}>0\)，所以本原，\(\gamma(A)\le n-1\)。

证明：\(A\ge\alpha(I+B/\alpha)\)，\(B\) 为去掉对角的部分（不可约），\((I+B/\alpha)^{n-1}>0\)（第 08a 章 Lemma 8.4.1）。

Lemma 8.5.5。 本原矩阵的所有幂 \(A^m\) 仍本原。（不可约矩阵的幂可以可约：\(\begin{bmatrix}0&1\\1&0\end{bmatrix}^2=I\)。）

Theorem 8.5.6。 本原 ⇒ 存在 \(k\le(n-1)n^n\) 使 \(A^k>0\)（粗糙的界，证明依次让每个节点有自环）。

Theorem 8.5.7（最短圈界）。 设 \(\Gamma(A)\) 的最短圈长为 \(s\)，则 \(\gamma(A)\le n+s(n-2)\)。

证明思路：写 \(n+s(n-2)=(n-s)+s(n-1)\)，考虑 \(A^{n-s}(A^s)^{n-1}\)。最短圈上的 \(s\) 个节点在 \(\Gamma(A^s)\) 中有自环；圈外每个节点都能用恰好 \(n-s\) 步走到圈上（不够就绕圈补足）；有自环的节点在 \(A^s\) 的图中用恰好 \(n-1\) 步可达任何节点（在自环上"等待"凑步数）。拼起来每个 \((i,j)\) 元都为正。

Corollary 8.5.8（Wielandt 定理）。 非负 \(A\in M_n\) 本原 ⇔ \(A^{n^2-2n+2}>0\)。

证明：本原且 \(n>1\) 时最短圈长 \(s\le n-1\)（若唯一的圈长都是 \(n\) 的倍数，由 8.5.3 不本原），代入 8.5.7：\(\gamma\le n+(n-1)(n-2)=n^2-2n+2=(n-1)^2+1\)。Wielandt 矩阵（8.5.P4：\(a_{12}=a_{23}=\cdots=a_{n-1,n}=a_{n1}=a_{n2}=1\)，其余为 0）的图只有长 \(n\) 与 \(n-1\) 两种圈，它的本原指数恰好是 \(n^2-2n+2\)，界是紧的（实战 2 验证）。

Theorem 8.5.9（Holladay–Varga）。 不可约且有 \(d\ge1\) 个正对角元 ⇒ \(\gamma(A)\le2n-d-1\)。

原书例：\(A=\begin{bmatrix}0&1\\1&1\end{bmatrix}\) 本原，特征值 \((1\pm\sqrt5)/2\)；最短圈界给 \(\gamma\le2+1\cdot0=2\)，Holladay–Varga 给 \(\gamma\le2\cdot2-1-1=2\)，实际 \(A^2=\begin{bmatrix}1&1\\1&2\end{bmatrix}>0\)，\(\gamma=2\)。

8.6.3 实用判定：反复平方

判定本原性只需要零/非零模式（布尔矩阵乘法），并利用"一旦为正就一直为正"：

• 判不可约：\(n=10\) 时不必算 \((I+A)^9\)（8 次乘法），算 \((I+A)^2,(I+A)^4,(I+A)^8,(I+A)^{16}\) 即可（4 次），因为 \((I+A)^{16}\) 的零模式包含于 \((I+A)^9\) 的零模式。
• 判本原：算 \(A^2,A^4,\dots,A^{128}\)（7 次乘法），\(128\ge82=n^2-2n+2\)，检查 \(A^{128}\) 是否全正。直接按 Wielandt 界要 81 次乘法。

一般地，判本原需要 \(\lceil\log_2(n^2-2n+2)\rceil\) 次布尔矩阵乘法。

8.6.4 习题中的其他结论

• 8.5.P2：本原时 \(\lim_m(a_{ij}^{(m)})^{1/m}=\rho(A)\) 对每个 \((i,j)\) 成立。
• 8.5.P3：本原矩阵之积不一定本原。
• 8.5.P8–P12：幂等的不可约非负矩阵是秩一正矩阵；\(\lim(A/\rho)^m\) 存在不需要本原（如 \(I\)）；但不可约且极限存在 ⇒ 本原。
• 8.5.P13：\(A\) 不可约非负 ⇒ \(A+\epsilon I\) 本原（加正对角元）。所以不可约矩阵是本原矩阵的极限。Markov 链中的"惰性化"（lazy chain：以 1/2 概率原地不动）就是用这个技巧消除周期性。
• 8.5.P14：组合对称（\(a_{ij}>0\iff a_{ji}>0\)）的本原矩阵满足 \(A^{2n-2}>0\)。
• 8.5.P15：\(n\) 为素数时，不可约非负非奇异矩阵要么本原，要么相似于 \(t^n-\rho^n\) 的友矩阵（全部特征值最大模）。
• 8.5.P16（幂法）：\(x^{(0)}>0\)，\(y^{(m+1)}=Ax^{(m)}\)，\(x^{(m+1)}=y^{(m+1)}/\sum_iy_i^{(m+1)}\)。\(A\) 本原时 \(x^{(m)}\to\) Perron 向量，\(\sum_iy_i^{(m+1)}\to\rho(A)\)，速度由 \(|\lambda_2|/\rho\) 决定；不本原时（如 \(\begin{bmatrix}0&1\\1&0\end{bmatrix}\)）可能振荡。这就是 PageRank 计算的数学基础（实战 3）。

8.7 一般极限定理：Cesàro 平均

不可约但不本原时 \((A/\rho)^m\) 振荡，但时间平均收敛。两个引理（练习）：

• \(\theta\in(0,2\pi)\) 时 \(\frac1N\sum_{m=1}^Ne^{im\theta}=\frac{e^{i\theta}-e^{i(N+1)\theta}}{N(1-e^{i\theta})}\to0\)；
• \(\rho(B)<1\) 时 \(\frac1N\sum_{m=1}^NB^m=\frac1N(B-B^{N+1})(I-B)^{-1}\to0\)。

Theorem 8.6.1。 \(A\) 不可约非负，\(n\ge2\)，\(x,y\) 为右、左 Perron 向量（\(y^Tx=1\)），则

且存在常数 \(C\) 使误差的 \(\infty\)-范数不超过 \(C/N\)。

证明：若有 \(k\) 个最大模特征值，它们是 \(\rho e^{2\pi ij/k}\)，都单重（第 08a 章 8.4.6），所以 \(A/\rho=S([1]\oplus[e^{i\theta}]\oplus\cdots\oplus[e^{i(k-1)\theta}]\oplus B)S^{-1}\)，\(\rho(B)<1\)。平均以后，单位根对应的项是 \(O(1/N)\)，\(B\) 对应的项也是 \(O(1/N)\)，只剩 \(S([1]\oplus0)S^{-1}=xy^T\)。8.6.P1 用 \(\begin{bmatrix}0&1\\1&0\end{bmatrix}\) 说明 \(O(1/N)\) 不能改进（奇数 \(N\) 时误差恰为 \(\frac1{2N}\) 量级）。8.6.P2：不可约时每个 \((i,j)\) 元在无穷多个 \(m\) 上为正，但对周期矩阵也在无穷多个 \(m\) 上为零。

推导拆解：两个引理为什么成立？第一个是等比数列求和：\(\sum_{m=1}^Ne^{im\theta}\) 是首项 \(e^{i\theta}\)、公比 \(e^{i\theta}\) 的等比数列和，等于 \(\frac{e^{i\theta}(1-e^{iN\theta})}{1-e^{i\theta}}\)。分子的模不超过 2（两个模为 1 的复数之差），分母是与 \(N\) 无关的非零常数（\(\theta\ne0\)），所以和是有界的，除以 \(N\) 后就是 \(O(1/N)\)。第二个是矩阵版的同一公式，\(\rho(B)<1\) 保证 \(B^{N+1}\to0\)，和同样有界。 直观：\(e^{im\theta}\) 是在单位圆上转圈的点，转来转去互相抵消，平均值趋于 0。唯一不转的是特征值 1（\(\theta=0\)），所以平均之后只剩它对应的部分 \(xy^T\)。周期链的振荡就这样被平均掉了。

直观：周期 Markov 链的 \(m\) 步转移概率不收敛，但"在各状态停留的时间比例"收敛到平稳分布——这是遍历定理的矩阵版本。量化中用一条长模拟路径的时间平均估计平稳量（例如 regime 模型的长期平均波动率），对周期链也是合法的，只是误差按 \(1/N\) 而不是几何速度下降。

8.8 随机矩阵与双随机矩阵

8.8.1 随机矩阵

定义。 非负矩阵 \(A\) 若 \(Ae=e\)（每行和为 1），称为**（行）随机矩阵（row stochastic matrix）；转置后的 \(e^TA=e^T\) 为列随机矩阵**（第 08a 章的人口迁移模型）。每一行是 \(n\) 个状态上的一个概率分布，Markov 链的一步转移矩阵就是行随机矩阵：\(P_{ij}=\Pr(X_{t+1}=j\mid X_t=i)\)，分布按行向量演化 \(\pi_{t+1}^T=\pi_t^TP\)。

基本性质：

• \(Ae=e\) 且 \(e>0\)，由第 08a 章 8.3.4，\(\rho(A)=1\)。
• 随机矩阵集合是紧凸集，对乘法封闭（8.7.P1，半群）。
• 幂有界（元素都在 \([0,1]\)），所以每个模为 1 的特征值都半单（8.7.P2）。
• 8.7.P3（对角相似化归）：若非负 \(A\) 有正特征向量 \(x\)，则 \(\rho(A)^{-1}D^{-1}AD\)（\(D=\operatorname{diag}(x)\)）是随机矩阵。许多关于"有正特征向量的非负矩阵"的问题因此都可以化为随机矩阵的问题。

Markov 链的 Perron–Frobenius 字典。

有吸收态的链。 信用评级迁移矩阵通常把"违约"设为吸收态：\(P=\begin{bmatrix}Q&r\\0&1\end{bmatrix}\)，\(Q\) 是非违约评级之间的子矩阵（次随机，行和 \(<1\)）。\(P\) 可约，\(\pi=(0,\dots,0,1)\) 是平凡的平稳分布；有意义的量来自瞬时块 \(Q\)：\(m\) 年累计违约概率是 \(P^m\) 的最后一列，"至今未违约"的概率约按 \(\rho(Q)^m\) 衰减，违约前的期望年数是 \((I-Q)^{-1}e\)（基本矩阵 \((I-Q)^{-1}=\sum Q^k\ge0\) 正是第 08a 章的 M 矩阵）。实战 1 演示。

推导拆解：为什么违约前期望年数是 \((I-Q)^{-1}e\)？ 第一步，\(Q^k\) 的 \((i,j)\) 元是"从评级 \(i\) 出发，第 \(k\) 年处于未违约评级 \(j\)"的概率（违约后出不来，所以只要还在 \(Q\) 的状态里就说明没违约）。第二步，\((Q^ke)_i=\sum_j(Q^k)_{ij}\) 是第 \(k\) 年仍未违约的概率，即存活概率 \(S_i(k)\)。第三步，期望存活年数 \(=\sum_{k\ge0}S_i(k)\)（每年存活就"贡献"一年，期望是逐年存活概率之和）。第四步，\(\sum_kQ^ke=(I-Q)^{-1}e\)，用了 Neumann 级数（\(\rho(Q)<1\) 保证收敛）。 金融直觉：标量版本你很熟：年违约率恒为 \(h\) 时，存活概率 \((1-h)^k\)，期望存活 \(\sum_k(1-h)^k=1/h\)，这与年金现值公式是同一个几何级数。矩阵版本允许评级在存活期间上下迁移，\(\rho(Q)\) 相当于"综合后的长期年存活率"。实战 1 中 \(\rho(Q)=0.982\)，对应长期每年约 1.8% 的违约强度。\(P^m\) 最后一列的累计违约概率曲线，则是信用利差期限结构和 IFRS 9 整个存续期预期损失的输入。

8.8.2 双随机矩阵

定义。 \(A\) 与 \(A^T\) 都随机，即非负且 \(Ae=e\)、\(e^TA=e^T\)，称双随机矩阵（doubly stochastic）。例子：置换矩阵；正交随机矩阵 \([|u_{ij}|^2]\)（\(U\) 实正交或酉，第 04a 章 Schur 优超定理与第 06 章 Hoffman–Wielandt 定理中出现过）。

练习：\(n\ge2\)，若双随机矩阵某个 \(a_{ii}=1\)，则第 \(i\) 行、第 \(i\) 列其余元素为 0，\(A\) 置换相似于 \([1]\oplus B\)，\(B\) 双随机。

Lemma 8.7.1。 \(A\) 双随机且 \(A\ne I\)，则存在非恒等置换 \(\sigma\) 使 \(a_{1\sigma(1)}\cdots a_{n\sigma(n)}>0\)。

证明：反设所有非恒等置换的乘积都为 0，则行列式展开中只剩对角项，\(p_A(t)=\prod(t-a_{ii})\)，对角元就是特征值；1 是特征值，所以某 \(a_{ii}=1\)，剥掉这一行列后对剩下的双随机矩阵重复，最终 \(A=I\)，矛盾。

Theorem 8.7.2（Birkhoff 定理）。 \(A\) 双随机 ⇔ 存在置换矩阵 \(P_1,\dots,P_N\) 与正数 \(t_1,\dots,t_N\)，\(\sum t_i=1\)，使

且可取 \(N\le n^2-n+1\)（8.7.P15 改进为 \(n^2-2n+2\)）。

构造性证明：若 \(A\) 不是置换矩阵，由 8.7.1（或对 \(A\) 本身用同样论证）找一个"正对角线"\(\sigma\)，令 \(\alpha_1=\min_ia_{i\sigma(i)}\)，\(P_1\) 为 \(\sigma\) 对应的置换矩阵，

\(A_1\) 仍双随机，且至少多一个零元。重复至多 \(n^2-n\) 次就剩下一个置换矩阵。每一步"剥离一个置换矩阵的正倍数并制造新零元"；找 \(\sigma\) 就是在二分图上找一个完美匹配（实战 4 用 linear_sum_assignment 实现）。8.7.P15 的改进来自维数：双随机矩阵满足 \(2n-1\) 个独立线性约束，集合是 \(\mathbf R^{(n-1)^2}\) 中的凸多面体，由 Carathéodory 定理（原书附录 B，本册附录 A.2）任一点是至多 \((n-1)^2+1\) 个顶点的凸组合。

推导拆解：手算一个 \(3\times3\) 例子，不做归一化，直接逐步剥离。 \(A=\begin{bmatrix}.5&.5&0\\.5&.25&.25\\0&.25&.75\end{bmatrix}\)，行和、列和都是 1。 第 1 步：对角线 \((1,1),(2,2),(3,3)\) 上的元素 \(.5,.25,.75\) 都为正，最小是 \(.25\)。减去 \(.25I\)，剩 \(\begin{bmatrix}.25&.5&0\\.5&0&.25\\0&.25&.5\end{bmatrix}\)，行和、列和都变成 \(.75\)，并且多了一个零元。 第 2 步：取置换 \(1\to2,2\to1,3\to3\)，元素 \(.5,.5,.5\)，减去 \(.5P_{(12)}\)，剩 \(\begin{bmatrix}.25&0&0\\0&0&.25\\0&.25&0\end{bmatrix}\)。 第 3 步：剩下的恰好是 \(.25P_{(23)}\)（置换 \(1\to1,2\to3,3\to2\)）。 结果：\(A=.25I+.5P_{(12)}+.25P_{(23)}\)，权重和为 1。原文的 \(A_1=(A-\alpha_1P_1)/(1-\alpha_1)\) 只是把每一步的剩余部分重新缩放成双随机，做法相同。 金融直觉：把 \(A\) 读成"3 笔订单分给 3 个通道"的软分配方案，分解就告诉你：25% 的时间按原样分配，50% 的时间交换前两笔订单的通道，25% 的时间交换后两笔。按这些概率随机选择硬分配，期望效果与软分配完全一样，而每次执行都是一笔订单只走一个通道。

历史：Kőnig 在 1916 年已对非负有理元素的情形得到等价结论，所以也称 Birkhoff–Kőnig 定理。8.7.P7–P8：置换矩阵恰好是双随机矩阵集合的极点。8.7.P11：分解不唯一。

Corollary 8.7.4。 凸（凹）实值函数在双随机矩阵集合上的最大值（最小值）在某个置换矩阵处取得。

证明：最大值点 \(A=\sum t_iP_i\)，\(f(A)\le\sum t_if(P_i)\le\max_if(P_i)\)。

这是原书附录 B（本册附录 A.2）"凸函数在紧凸集上的最大值在极点取得"的特例，也是第 06 章 Hoffman–Wielandt 定理证明的关键。线性函数既凸又凹，所以线性目标的分配问题 \(\min_{X\text{ 双随机}}\sum c_{ij}x_{ij}\) 的 LP 松弛总有一个置换矩阵最优解——分配问题天然具有整数性（实战 4）。

白话解释：证明只有一行，关键是"凸函数的平均值不小于平均处的函数值"（凸函数定义，即 Jensen 不等式）。\(f(A)=f(\sum t_iP_i)\le\sum t_if(P_i)\)，而加权平均不超过其中的最大者。所以在任何软分配处取到的值，都被某个硬分配超过或追平。 实务意义：分配问题本来是整数规划（每笔订单只能选一个通道，\(n!\) 种方案），一般很难。但因为可行域的角点全是置换矩阵，放松成连续变量的线性规划就能直接得到整数解，用匈牙利算法 \(O(n^3)\) 即可求解。注意这只对"每行每列和为 1"这种约束结构成立；加上额外约束（如通道容量不同）就可能失去整数性。

8.8.3 双次随机矩阵与 von Neumann 迹定理

非负且行和、列和都 \(\le1\) 的矩阵叫双次随机（doubly substochastic）。

Lemma 8.7.5。 双次随机矩阵被某个双随机矩阵逐元素控制：\(A\le S\)。（逐个增大"行和与列和都不足 1"的位置上的元素。）

练习：\(U,V\) 酉，则 \([|u_{ij}v_{ji}|]\) 双次随机（Cauchy–Schwarz）。

Theorem 8.7.6（von Neumann 迹定理）。 \(A,B\in M_n\)：

证明：\(A=V_1\Sigma_AW_1^*\)，\(B=V_2\Sigma_BW_2^*\)，令 \(U=W_1^*V_2\)、\(V=W_2^*V_1\)（酉），

\([|u_{ij}v_{ji}|]\) 双次随机，被某双随机矩阵控制；右端关于双随机矩阵是线性函数，由 8.7.4 最大值在置换矩阵 \(\pi\) 处取得，得 \(\le\sum_i\sigma_i(A)\sigma_{\pi(i)}(B)\le\sum_i\sigma_i(A)\sigma_i(B)\)（重排不等式）。这就是第 07b 章所有最佳逼近定理（Eckart–Young、Procrustes）的"发动机"。

白话解释：最后一步"\(\sum_i\sigma_i(A)\sigma_{\pi(i)}(B)\le\sum_i\sigma_i(A)\sigma_i(B)\)"是重排不等式：两列数配对相乘求和，大配大、小配小时最大。金融例子：手上有三笔资金 100、50、10，三个策略的预期收益率 8%、5%、2%，要一对一分配，总收益最大的方案是 100 配 8%、50 配 5%、10 配 2%。 整个定理说的是：\(\operatorname{tr}(AB)\) 可以看作 \(A\) 与 \(B\) 的"对齐程度"，旋转 \(A\)、\(B\) 的方向能改变它，但最多只能做到把两者的奇异向量对齐、奇异值大配大。实战 4 第 3 部分的数值正是这个上界。

8.8.4 习题中的其他结论

• 8.7.P5：双随机矩阵满足 \(\sigma_1(A)=\rho(A)=1\)（谱范数为 1）。证明：\(\|A\|_2\le\sum t_i\|P_i\|_2=1\)。8.7.P6：任何由置换不变向量范数诱导的矩阵范数在双随机矩阵上取值 1。
• 8.7.P9：双随机矩阵不能恰有 \(n+1\) 个正元。
• 8.7.P10：\(2\times2\) 双随机矩阵都形如 \(\begin{bmatrix}a&1-a\\1-a&a\end{bmatrix}\)。
• 8.7.P12：\(A\) 双随机、对称、半正定 ⇒ \(A^{1/2}e=e\)；\(A^{1/2}\) 一般不一定非负，但 \(n=2\) 时非负。
• 8.7.P13：酉不变矩阵范数满足 \(\|A\|\le\|I\|\)（\(A\) 双随机），用 Ky Fan 占优定理（第 07b 章）。
• 8.7.P14：可约的双随机矩阵置换相似于 \(A_1\oplus A_2\)（两块都双随机）——可约即完全可分解。

量化实战

实战 1：评级迁移与市场 regime——吸收、平稳分布与混合速度

场景一：一年期评级迁移矩阵（AAA、A、BBB、BB、B、违约 D），D 为吸收态。求多年累计违约概率、违约前期望年数。场景二：三状态市场 regime 链（牛、震荡、熊），求平稳分布与收敛速度。场景三：一条周期为 2 的链，演示只有 Cesàro 平均收敛。

import numpy as np
from scipy.sparse.csgraph import connected_components

np.set_printoptions(suppress=True)
# ---------- 1. 一年期评级迁移矩阵（行随机，虚构但形状典型）：AAA, A, BBB, BB, B, D ----------
R = ["AAA", "A", "BBB", "BB", "B", "D"]
P = np.array([[0.92, 0.07, 0.008, 0.001, 0.0005, 0.0005],
              [0.02, 0.90, 0.065, 0.010, 0.003, 0.002],
              [0.002, 0.05, 0.87, 0.06, 0.012, 0.006],
              [0.001, 0.005, 0.07, 0.82, 0.08, 0.024],
              [0.000, 0.002, 0.01, 0.07, 0.83, 0.088],
              [0.0, 0.0, 0.0, 0.0, 0.0, 1.0]])
assert np.allclose(P.sum(1), 1)
nc, lab = connected_components(P > 0, directed=True, connection='strong')
print("强连通分量:", nc, "标签", lab, "→ 可约：{AAA..B} 是瞬时类，D 是吸收态")
for m in [1, 5, 10, 30]:
    print("%2d 年累计违约概率:" % m, dict(zip(R[:5], np.linalg.matrix_power(P, m)[:5, 5].round(4).tolist())))
Q = P[:5, :5]                                       # 瞬时块
print("瞬时块谱半径 ρ(Q) = %.4f → 'P(仍未违约)' 约按 ρ(Q)^m 衰减；平均违约前年数 (I-Q)^{-1}1 =" % max(abs(np.linalg.eigvals(Q))),
      dict(zip(R[:5], np.linalg.solve(np.eye(5) - Q, np.ones(5)).round(1).tolist())))

# ---------- 2. 不可约 + 本原的状态链：市场 regime（牛/震荡/熊） ----------
S = np.array([[0.90, 0.08, 0.02],
              [0.10, 0.80, 0.10],
              [0.05, 0.15, 0.80]])
w, V = np.linalg.eig(S.T)
pi = np.real(V[:, np.argmax(w.real)]); pi /= pi.sum()      # 左 Perron 向量 = 平稳分布
lam = np.sort(abs(w))[::-1]
print("\nregime 链：平稳分布 π =", pi.round(4), " 次特征值模 |λ2| = %.4f" % lam[1])
for m in [5, 20, 50]:
    err = np.abs(np.linalg.matrix_power(S, m) - np.outer(np.ones(3), pi)).max()
    print("  m=%2d: max|S^m - 1π^T| = %.2e   |λ2|^m = %.2e" % (m, err, lam[1] ** m))
print("  半衰期 ≈ ln 2 / (-ln|λ2|) = %.1f 期" % (np.log(2) / -np.log(lam[1])))

# ---------- 3. 周期链：只有 Cesàro 平均收敛（原书 8.6.1） ----------
C = np.array([[0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1], [0.5, 0, 0.5, 0]])   # 周期 2
print("\n周期链特征值模:", np.round(np.sort(abs(np.linalg.eigvals(C)))[::-1], 4))
wc, Vc = np.linalg.eig(C.T); pc = np.real(Vc[:, np.argmin(abs(wc - 1))]); pc /= pc.sum()
for N in [10, 100, 1000]:
    avg = sum(np.linalg.matrix_power(C, m) for m in range(1, N + 1)) / N
    print("  N=%4d: C^N 与极限差 %.3f；Cesàro 平均与 1π^T 差 %.4f（≈ C/N）"
          % (N, np.abs(np.linalg.matrix_power(C, N) - np.outer(np.ones(4), pc)).max(),
             np.abs(avg - np.outer(np.ones(4), pc)).max()))

关键输出：

强连通分量: 2 标签 [1 1 1 1 1 0] → 可约：{AAA..B} 是瞬时类，D 是吸收态
 1 年累计违约概率: {'AAA': 0.0005, 'A': 0.002, 'BBB': 0.006, 'BB': 0.024, 'B': 0.088}
 5 年累计违约概率: {'AAA': 0.0051, 'A': 0.0174, 'BBB': 0.0463, 'BB': 0.1387, 'B': 0.3302}
10 年累计违约概率: {'AAA': 0.0192, 'A': 0.0531, 'BBB': 0.1185, 'BB': 0.2725, 'B': 0.4979}
30 年累计违约概率: {'AAA': 0.1738, 'A': 0.2795, 'BBB': 0.4028, 'BB': 0.5777, 'B': 0.7458}
瞬时块谱半径 ρ(Q) = 0.9821 → 'P(仍未违约)' 约按 ρ(Q)^m 衰减；平均违约前年数 (I-Q)^{-1}1 = {'AAA': 78.1, 'A': 67.6, 'BBB': 56.9, 'BB': 42.2, 'B': 27.4}

regime 链：平稳分布 π = [0.4464 0.3393 0.2143]  次特征值模 |λ2| = 0.8306
  m= 5: max|S^m - 1π^T| = 2.18e-01   |λ2|^m = 3.95e-01
  m=20: max|S^m - 1π^T| = 1.47e-02   |λ2|^m = 2.44e-02
  m=50: max|S^m - 1π^T| = 5.65e-05   |λ2|^m = 9.34e-05
  半衰期 ≈ ln 2 / (-ln|λ2|) = 3.7 期

周期链特征值模: [1.     1.     0.7071 0.7071]
  N=  10: C^N 与极限差 0.354；Cesàro 平均与 1π^T 差 0.0458（≈ C/N）
  N= 100: C^N 与极限差 0.333；Cesàro 平均与 1π^T 差 0.0044（≈ C/N）
  N=1000: C^N 与极限差 0.333；Cesàro 平均与 1π^T 差 0.0004（≈ C/N）

读法：

1. 评级矩阵有两个强连通分量：五个评级互通但都通向 D，D 只通向自己。\(P\) 可约，长期看所有债券都会违约（30 年时 B 级的累计违约率 75%），所以"平稳分布"没有信息量。有用的是瞬时块：\(\rho(Q)=0.982\) 决定了"存活概率"的长期衰减速度，基本矩阵 \((I-Q)^{-1}\) 给出违约前的期望年数（AAA 约 78 年、B 约 27 年）。多期违约概率曲线是信用利差期限结构、CDS 定价和 IFRS 9 预期损失的输入。
2. regime 链本原（对角元全正，Lemma 8.5.4），平稳分布 \(\pi=(0.446,0.339,0.214)\) 就是左 Perron 向量：长期看约 45% 的时间处于牛市。\(S^m\) 收敛到 \(e\pi^T\) 的速度与 \(|\lambda_2|^m=0.83^m\) 一致，"记忆"的半衰期约 3.7 期。\(|\lambda_2|\) 越接近 1，regime 越持续，用短样本估计平稳分布越不可靠。
3. 周期 2 的链有两个模为 1 的特征值，\(C^N\) 与极限的差始终是 0.33，永不收敛；Cesàro 平均的误差按 \(1/N\) 下降，与 Theorem 8.6.1 的 \(O(1/N)\) 一致。

实战 2：本原性与本原指数的判定

import numpy as np

def bool_mult(X, Y):
    return (X.astype(int) @ Y.astype(int)) > 0

def primitivity_index(A, max_power=None):
    """最小的 k 使 A^k > 0（只看零模式）；不本原返回 None。"""
    n = A.shape[0]; B = A > 0; P = B.copy()
    for k in range(1, (max_power or n * n - 2 * n + 2) + 1):
        if P.all():
            return k
        P = bool_mult(P, B)
    return None

def is_primitive_by_squaring(A):
    """反复平方：A^{2^j} 的零模式，2^j ≥ n^2-2n+2 时检查是否全正（Wielandt 界）。"""
    n = A.shape[0]; P = A > 0; mults = 0; p = 1
    while p < n * n - 2 * n + 2:
        P = bool_mult(P, P); p *= 2; mults += 1
    return P.all(), mults

# ---------- 1. Wielandt 矩阵：本原指数恰好达到上界 n^2-2n+2（原书 8.5.P4） ----------
for n in [3, 5, 8, 10]:
    W = np.zeros((n, n))
    for i in range(n - 1):
        W[i, i + 1] = 1
    W[n - 1, 0] = W[n - 1, 1] = 1
    print("n=%2d Wielandt 矩阵：γ = %3d，界 n²-2n+2 = %3d；反复平方判定本原=%s，用 %d 次布尔乘法"
          % (n, primitivity_index(W), n * n - 2 * n + 2, *is_primitive_by_squaring(W)))

# ---------- 2. 原书例 [[0,1],[1,1]]：γ = 2 ----------
F = np.array([[0, 1], [1, 1]])
print("\n[[0,1],[1,1]]：γ =", primitivity_index(F), "；A^2 =", np.linalg.matrix_power(F, 2).tolist())

# ---------- 3. Holladay–Varga：d 个正对角元 ⇒ γ ≤ 2n-d-1（原书 8.5.9） ----------
rng = np.random.default_rng(1)
n = 10
ring = np.roll(np.eye(n), 1, axis=1)                     # 单向环：不可约但周期 n
for d in [0, 1, 3, 10]:
    A = ring.copy(); A[np.arange(d), np.arange(d)] = 1
    print("环 + %2d 个自环：γ = %s，Holladay–Varga 界 2n-d-1 = %s"
          % (d, primitivity_index(A, 200), (2 * n - d - 1) if d > 0 else "—（不本原）"))

# ---------- 4. 图论判据：各节点闭路长度的 gcd（原书 8.5.3） ----------
C = np.zeros((6, 6))
for i in range(5): C[i, i + 1] = 1
C[5, 0] = 1; C[2, 0] = 1                                  # 闭路长度 6 与 3 → gcd 3：不本原
print("\n闭路长度 {3,6} 的图：本原指数 =", primitivity_index(C, 200),
      "；最大模特征值个数 =", int(np.sum(np.isclose(abs(np.linalg.eigvals(C)), max(abs(np.linalg.eigvals(C)))))))
C[3, 0] = 1                                               # 再加一条长度 4 的闭路 → gcd(3,4,6)=1
print("再加长度 4 的闭路：本原指数 =", primitivity_index(C, 200))

关键输出：

n= 3 Wielandt 矩阵：γ =   5，界 n²-2n+2 =   5；反复平方判定本原=True，用 3 次布尔乘法
n= 5 Wielandt 矩阵：γ =  17，界 n²-2n+2 =  17；反复平方判定本原=True，用 5 次布尔乘法
n= 8 Wielandt 矩阵：γ =  50，界 n²-2n+2 =  50；反复平方判定本原=True，用 6 次布尔乘法
n=10 Wielandt 矩阵：γ =  82，界 n²-2n+2 =  82；反复平方判定本原=True，用 7 次布尔乘法

[[0,1],[1,1]]：γ = 2 ；A^2 = [[1, 1], [1, 2]]
环 +  0 个自环：γ = None，Holladay–Varga 界 2n-d-1 = —（不本原）
环 +  1 个自环：γ = 18，Holladay–Varga 界 2n-d-1 = 18
环 +  3 个自环：γ = 16，Holladay–Varga 界 2n-d-1 = 16
环 + 10 个自环：γ = 9，Holladay–Varga 界 2n-d-1 = 9

闭路长度 {3,6} 的图：本原指数 = None ；最大模特征值个数 = 3
再加长度 4 的闭路：本原指数 = 13

读法：

1. Wielandt 矩阵的本原指数恰好等于 \(n^2-2n+2\)，说明 Wielandt 界是紧的；\(n=10\) 时反复平方只需 7 次布尔乘法，与原书的计算一致。
2. 单向环（周期 \(n\)）不本原；只要加一个自环就变成本原，且环加自环的例子恰好达到 Holladay–Varga 界 \(2n-d-1\)。Markov 链中"惰性化"（给每个状态加一点停留概率）消除周期性，就是 \(d=n\) 的情形，\(\gamma\le n-1\)。
3. 闭路长度只有 3 与 6 时 gcd 为 3，不本原，恰有 3 个最大模特征值（Romanovsky 定理）；加一条长 4 的闭路后 gcd 变成 1，本原。

实战 3：PageRank 式中心性——阻尼、本原性与收敛速度

场景：在资金流、持股或"注意力"有向网络上计算节点的重要性。随机游走矩阵 \(P\)（行随机）往往可约：有悬挂节点（没有出边）、有只进不出的"陷阱"子图。PageRank 的做法是混入均匀跳转：\(G=\alpha P+(1-\alpha)\frac1nee^T\)，\(G\) 的元素全为正，因而本原，左 Perron 向量唯一且为正，幂法收敛，次特征值满足 \(|\lambda_2(G)|\le\alpha\)。

import numpy as np
from scipy.sparse.csgraph import connected_components

rng = np.random.default_rng(3)
# ---------- 有向"资金流/持股"网络：Adj[i, j] = 1 表示 i 把权重（资金、注意力）指向 j ----------
n = 40
Adj = (rng.uniform(size=(n, n)) < 0.08).astype(float)
np.fill_diagonal(Adj, 0)
Adj[rng.choice(36, 15, replace=False), 0] = 1          # 节点 0：被很多节点指向的"枢纽"
Adj[37, :] = 0                                        # 节点 37：没有出边（悬挂节点）
Adj[38, :] = 0; Adj[39, :] = 0; Adj[38, 39] = Adj[39, 38] = 1   # 38↔39：只进不出的"陷阱"
Adj[5, 38] = 1
out = Adj.sum(1)
nc, lab = connected_components(Adj > 0, directed=True, connection='strong')
print("强连通分量 %d 个 → 随机游走矩阵可约" % nc)

# 行随机转移矩阵；悬挂节点均匀跳转
P = np.where(out[:, None] > 0, Adj / np.maximum(out, 1)[:, None], 1.0 / n)

def pagerank(P, alpha, tol=1e-10, maxit=100_000):
    n = P.shape[0]
    G = alpha * P + (1 - alpha) / n                    # Google 矩阵：元素全正 → 本原（原书 8.5.2）
    x = np.full(n, 1.0 / n)
    for k in range(maxit):
        x_new = x @ G                                  # x_{k+1}^T = x_k^T G（左 Perron 向量的幂法）
        if np.abs(x_new - x).sum() < tol:
            return x_new, k + 1, np.sort(abs(np.linalg.eigvals(G)))[::-1][1]
        x = x_new

for alpha in [0.5, 0.85, 0.95, 0.99]:
    x, it, l2 = pagerank(P, alpha)
    print("α=%.2f：幂法 %5d 次收敛，|λ2(G)| = %.4f ≤ α，陷阱节点 38、39 的总权重 %.3f，枢纽 0 的权重 %.3f"
          % (alpha, it, l2, x[38] + x[39], x[0]))

# ---------- 不加阻尼（α=1）：可约 → 质量被陷阱吸走 ----------
x = np.full(n, 1.0 / n)
for _ in range(5000):
    x = x @ P
print("α=1：5000 步后陷阱节点 38、39 的总权重 %.3f（其余节点的'中心性'全被吸走）" % (x[38] + x[39]))

关键输出：

强连通分量 6 个 → 随机游走矩阵可约
α=0.50：幂法    28 次收敛，|λ2(G)| = 0.5000 ≤ α，陷阱节点 38、39 的总权重 0.067，枢纽 0 的权重 0.095
α=0.85：幂法   116 次收敛，|λ2(G)| = 0.8500 ≤ α，陷阱节点 38、39 的总权重 0.140，枢纽 0 的权重 0.136
α=0.95：幂法   364 次收敛，|λ2(G)| = 0.9500 ≤ α，陷阱节点 38、39 的总权重 0.294，枢纽 0 的权重 0.121
α=0.99：幂法  1857 次收敛，|λ2(G)| = 0.9900 ≤ α，陷阱节点 38、39 的总权重 0.658，枢纽 0 的权重 0.060
α=1：5000 步后陷阱节点 38、39 的总权重 1.000（其余节点的'中心性'全被吸走）

读法：

1. 网络有 6 个强连通分量，\(P\) 可约。不加阻尼时，随机游走最终被"只进不出"的 38↔39 吸收，所有其他节点的中心性都变成 0——这是 Frobenius 标准形中闭类（常返类）吸收质量的直接后果，结果毫无意义。
2. 加阻尼后 \(G>0\)，由 Perron 定理有唯一正的平稳分布。由于 \(P\) 有多个闭类（特征值 1 多重），\(|\lambda_2(G)|\) 恰好等于 \(\alpha\)，幂法的迭代次数大致按 \(\log(\text{tol})/\log\alpha\) 增长：\(\alpha=0.85\) 时 116 次，\(\alpha=0.99\) 时 1857 次。
3. \(\alpha\) 是"跟随网络结构"与"均匀先验"之间的权衡：\(\alpha\) 太大，结果被陷阱结构主导（38、39 占 66%），枢纽节点 0 的重要性反而被稀释；\(\alpha=0.85\) 是常见折中。在量化中，这类中心性用于识别资金流网络中的核心标的、供应链网络中的关键企业，以及构造"网络动量/溢出"因子。

实战 4：Birkhoff 分解、分配问题的整数性与 von Neumann 迹不等式

场景：(1) 把一个双随机矩阵按 Birkhoff 定理的构造性证明分解成置换矩阵的凸组合。(2) 5 笔大单分配给 5 个执行通道（经纪商/算法），\(c_{ij}\) 为预估冲击成本，求总成本最小的一一分配：LP 松弛的最优解自动是置换矩阵。(3) 数值验证 von Neumann 迹不等式。

import numpy as np
from scipy.optimize import linear_sum_assignment, linprog

rng = np.random.default_rng(6)

def sinkhorn(M, iters=2000):
    """交替归一化行和列：正矩阵 → 双随机矩阵（Sinkhorn–Knopp）。"""
    A = M.copy()
    for _ in range(iters):
        A /= A.sum(1, keepdims=True); A /= A.sum(0, keepdims=True)
    return A

def birkhoff(A, tol=1e-12):
    """原书 8.7.2 的构造性证明：反复剥离一个置换矩阵。"""
    A = A.copy(); terms = []
    while A.max() > tol:
        # 找一个正"对角线"：在正元素上的完美匹配（代价取 -log 时就是最大乘积匹配）
        cost = np.where(A > tol, -np.log(np.maximum(A, tol)), 1e9)
        r, c = linear_sum_assignment(cost)
        t = A[r, c].min()
        Pm = np.zeros_like(A); Pm[r, c] = 1
        terms.append((t, Pm)); A -= t * Pm
        A[np.abs(A) < tol] = 0
    return terms

n = 5
D = sinkhorn(rng.uniform(0.1, 1, (n, n)) * (rng.uniform(size=(n, n)) < 0.7) + 1e-3 * np.eye(n))
print("双随机? 行和", D.sum(1).round(6), " 列和", D.sum(0).round(6))
terms = birkhoff(D)
recon = sum(t * Pm for t, Pm in terms)
print("Birkhoff 分解：%d 个置换矩阵（上界 n²-2n+2 = %d），权重和 %.6f，重构误差 %.1e"
      % (len(terms), n * n - 2 * n + 2, sum(t for t, _ in terms), np.abs(recon - D).max()))
for t, Pm in terms[:3]:
    print("   权重 %.4f  置换 %s" % (t, Pm.argmax(1).tolist()))

# ---------- 分配问题：LP 松弛的最优解自动是置换矩阵（原书 8.7.4） ----------
# 例：5 笔大单分给 5 个执行通道，c_ij = 第 i 笔单走通道 j 的预估冲击成本（bp）
Cost = rng.uniform(2, 20, (n, n)).round(1)
A_eq = np.vstack([np.kron(np.eye(n), np.ones(n)), np.kron(np.ones(n), np.eye(n))])   # 行和、列和 = 1
res = linprog(Cost.ravel(), A_eq=A_eq, b_eq=np.ones(2 * n), bounds=(0, 1), method="highs")
Xlp = res.x.reshape(n, n)
r, c = linear_sum_assignment(Cost)
print("\nLP 松弛最优值 %.1f；匈牙利算法最优值 %.1f；LP 解是 0-1 置换矩阵? %s"
      % (res.fun, Cost[r, c].sum(), np.allclose(Xlp, Xlp.round())))

# ---------- von Neumann 迹不等式（原书 8.7.6）：Re tr(AB) ≤ Σ σ_i(A) σ_i(B) ----------
A = rng.standard_normal((n, n)); B = rng.standard_normal((n, n))
sA, sB = np.linalg.svd(A, compute_uv=False), np.linalg.svd(B, compute_uv=False)
best = max(np.trace(A @ Q @ B @ Q2) for Q, Q2 in
           [(np.linalg.qr(rng.standard_normal((n, n)))[0], np.linalg.qr(rng.standard_normal((n, n)))[0]) for _ in range(20000)])
U1, _, V1t = np.linalg.svd(A); U2, _, V2t = np.linalg.svd(B)
Qopt, Q2opt = V1t.T @ U2.T, V2t.T @ U1.T                 # 对齐奇异向量
print("\nΣ σ_i(A)σ_i(B) = %.4f；随机正交旋转下 tr(A Q B Q') 的最大值 %.4f；对齐奇异向量后 %.4f"
      % ((sA * sB).sum(), best, np.trace(A @ Qopt @ B @ Q2opt)))

关键输出：

双随机? 行和 [1. 1. 1. 1. 1.]  列和 [1. 1. 1. 1. 1.]
Birkhoff 分解：9 个置换矩阵（上界 n²-2n+2 = 17），权重和 1.000000，重构误差 1.1e-16
   权重 0.2820  置换 [3, 1, 4, 0, 2]
   权重 0.2004  置换 [3, 0, 4, 2, 1]
   权重 0.2029  置换 [3, 1, 2, 0, 4]

LP 松弛最优值 34.2；匈牙利算法最优值 34.2；LP 解是 0-1 置换矩阵? True

Σ σ_i(A)σ_i(B) = 26.3813；随机正交旋转下 tr(A Q B Q') 的最大值 19.3325；对齐奇异向量后 26.3813

读法：

1. Sinkhorn 迭代（交替归一化行、列）把正矩阵变成双随机矩阵——它在最优传输、"把打分矩阵变成软分配"中常用。Birkhoff 分解用了 9 个置换矩阵，少于理论上界 17，重构误差在机器精度。每一步找"正对角线"就是一次二分图完美匹配。分解的一个直接用途：把"软"的分配方案（例如每笔订单按比例分到各通道）实现为若干个"硬"分配方案的随机混合，期望效果完全相同。
2. 分配问题的 LP 松弛（只要求 \(X\) 双随机）最优值与匈牙利算法完全相同，LP 解本身就是 0-1 置换矩阵：线性目标在双随机多面体上的最优解落在极点（Corollary 8.7.4）。配对交易中的资产配对、交易与对手方的一一匹配、因子与因子的匹配（跨期对齐的离散版本）都属于这一类。
3. 两万次随机正交旋转得到的 \(\operatorname{tr}(AQBQ')\) 最大只有 19.3，而把两个矩阵的奇异向量对齐后恰好达到上界 \(\sum\sigma_i(A)\sigma_i(B)=26.38\)。这正是第 07b 章 Procrustes 问题的解法原理。

本章小结

本原矩阵是不可约且只有一个最大模特征值的非负矩阵，等价于某个幂为正，又等价于每个节点的闭路长度最大公约数为 1；本原性只依赖零模式。本原时 \((A/\rho)^m\to xy^T>0\)，Perron 定理全部成立。本原指数有一系列上界：对角元全正时 \(\le n-1\)；最短圈长为 \(s\) 时 \(\le n+s(n-2)\)；一般地 \(\le n^2-2n+2\)（Wielandt，Wielandt 矩阵说明它紧）；有 \(d\) 个正对角元时 \(\le2n-d-1\)（Holladay–Varga）。实际判定用反复平方的布尔运算。不可约而不本原（周期 \(k>1\)）时幂振荡，但 Cesàro 平均以 \(O(1/N)\) 的速度收敛到 \(xy^T\)。随机矩阵的谱半径为 1，单位圆上的特征值半单；Markov 链的平稳分布是左 Perron 向量，本原（遍历）链 \(P^m\to e\pi^T\)，速度由 \(|\lambda_2|\) 决定；有吸收态的链要看瞬时块和基本矩阵 \((I-Q)^{-1}\)。双随机矩阵恰好是置换矩阵的凸包（Birkhoff 定理，构造性证明每步剥离一个置换矩阵），所以凸函数在其上的最大值、线性函数的最优值都在置换矩阵处取得——分配问题的整数性和 von Neumann 迹不等式都由此而来。

练习

基础

1. 判断下列 0-1 矩阵是否不可约、是否本原，并求本原指数：\(\begin{bmatrix}0&1\\1&0\end{bmatrix}\)，\(\begin{bmatrix}1&1\\1&0\end{bmatrix}\)，\(\begin{bmatrix}0&1&0\\0&0&1\\1&1&0\end{bmatrix}\)。 答案要点：第一个不可约不本原（周期 2）；第二个本原，\(\gamma=2\)；第三个是 \(n=3\) 的 Wielandt 矩阵，本原，\(\gamma=5\)。
2. 对 regime 转移矩阵 \(\begin{bmatrix}0.95&0.05\\0.10&0.90\end{bmatrix}\)，求平稳分布、次特征值和"牛市"状态的期望持续期。 答案要点：\(\pi=(2/3,1/3)\)；\(\lambda_2=0.85\)；牛市持续期 \(1/0.05=20\) 期。
3. 证明：行随机矩阵的每个特征值满足 \(|\lambda|\le1\)，且 \(e\) 是特征值 1 的右特征向量。
4. 说明为什么"惰性化" \(\frac12(I+P)\) 总是本原（\(P\) 不可约），并比较它与 \(P\) 的平稳分布。 答案要点：对角元为正 + 不可约 ⇒ 本原（Lemma 8.5.4）；平稳分布相同。
5. 用 Birkhoff 定理证明双随机矩阵的谱范数等于 1。 提示：8.7.P5，\(\|A\|_2\le\sum t_i\|P_i\|_2=1\)，且 \(Ae=e\)。
6. 证明：\(2\times2\) 双随机矩阵都形如 \(\begin{bmatrix}a&1-a\\1-a&a\end{bmatrix}\)，并写出它的 Birkhoff 分解。 答案要点：\(aI+(1-a)\begin{bmatrix}0&1\\1&0\end{bmatrix}\)。

进阶

7. 设评级链 \(P=\begin{bmatrix}Q&r\\0&1\end{bmatrix}\)，\(\rho(Q)<1\)。证明 \(P^m\) 的最后一列收敛到 \(e\)（最终必然违约），并证明从状态 \(i\) 出发到违约的期望步数是 \(((I-Q)^{-1}e)_i\)。 提示：\(P^m=\begin{bmatrix}Q^m&(I+Q+\cdots+Q^{m-1})r\\0&1\end{bmatrix}\)；\((I-Q)^{-1}=\sum Q^k\) 的 \((i,j)\) 元是访问 \(j\) 的期望次数。
8. 证明 Google 矩阵 \(G=\alpha P+(1-\alpha)\frac1nee^T\) 的除 1 以外的特征值都是 \(\alpha\) 乘以 \(P\) 的某个特征值，从而 \(|\lambda_2(G)|\le\alpha\)。 提示：取非奇异 \(S=[e\ \ S_1]\)。因 \(Pe=e\)，\(S^{-1}PS=\begin{bmatrix}1&\star\\0&P_1\end{bmatrix}\)，\(P_1\) 的特征值是 \(\lambda_2,\dots,\lambda_n\)；又 \(S^{-1}e=e_1\)，所以 \(S^{-1}(\frac1nee^T)S=\begin{bmatrix}1&\star\\0&0\end{bmatrix}\)。于是 \(S^{-1}GS=\begin{bmatrix}1&\star\\0&\alpha P_1\end{bmatrix}\)。
9. 证明 Wielandt 矩阵（\(n\ge3\)）的图只有长 \(n\) 和 \(n-1\) 两种圈，并据此说明它本原。
10. 用 Corollary 8.7.4 证明：对任意实向量 \(a,b\)，\(\max_{X\text{ 双随机}}a^TXb=\sum_ia_i^\downarrow b_i^\downarrow\)（重排不等式），并用它解释"按信号强弱排序配对"的最优性。
11. 实现 Sinkhorn 迭代，证明它保持矩阵的"交叉比" \(a_{ij}a_{kl}/(a_{il}a_{kj})\) 不变，并讨论当正矩阵改为有零元的非负矩阵时它何时会失败。 提示：每一步都是左右乘正对角阵；需要存在"支撑"在零模式上的双随机矩阵（完全不可分解性）。

原书推荐习题

• 8.5.P4：Wielandt 矩阵，界的紧性；8.5.P13：\(A+\epsilon I\) 本原；8.5.P16：幂法收敛性；8.5.P18：对称矩阵本原 ⇔ \(-\rho\) 不是特征值（无向图非二部）。
• 8.6.P1：\(O(1/N)\) 界的最优性。
• 8.7.P2、P3：随机矩阵特征值半单、对角相似化归为随机矩阵；8.7.P5：双随机矩阵谱范数为 1；8.7.P15：Birkhoff 项数 \(n^2-2n+2\)（Carathéodory 定理）。

原书对照

原书正文到此结束。附录 A–F 的背景知识（复数、凸集与凸函数、代数基本定理、特征值的连续性、紧性与 Weierstrass 定理、典范对）见本册附录（A_附录_凸性与连续性等背景.md）。