第 08a 章 非负矩阵与 Perron–Frobenius 定理
对应原书:Horn & Johnson《Matrix Analysis》第 2 版,第 8 章 Positive and Nonnegative Matrices 的 8.0–8.4 节(书 p.517–540,PDF p.537–560)。本原矩阵、随机矩阵与 Birkhoff 定理见第 08b 章。
第 07 章(07a–07d)的"正"指二次型为正;这一章的"正"指每个元素为正。两者是完全不同的概念:\(\begin{bmatrix}1&2\\2&1\end{bmatrix}\) 元素全正但不是正定矩阵,\(\begin{bmatrix}2&-1\\-1&2\end{bmatrix}\) 正定却有负元素。元素非负的矩阵在量化里随处可见:Markov 链的转移概率矩阵、信用评级迁移矩阵、产业间投入产出系数矩阵、银行间敞口矩阵、资金流与持股网络的邻接矩阵。它们共同的数学结构由 Perron–Frobenius 理论刻画:谱半径本身就是一个特征值,对应一个非负(往往为正)的特征向量,它决定了长期分布、增长率、中心性和系统性风险的放大倍数。
学习目标
读完本章,你应当能够:
- 区分逐元素序 \(A\ge B\) 与第 07d 章的 Loewner 序 \(A\succeq B\),熟练使用 \(|Ax|\le|A||x|\) 与"谱半径关于逐元素序单调"。
- 用行和、列和以及任意正向量加权(Collatz–Wielandt 界)估计非负矩阵的谱半径,理解"正特征向量必属于谱半径"。
- 陈述并理解 Perron 定理(正矩阵)的六条结论及其证明思路,知道收敛速度由次特征值决定。
- 知道对一般非负矩阵哪些结论保留、哪些失效;掌握不可约性、\((I+A)^{n-1}>0\) 判据与 Perron–Frobenius 定理,以及最大模特征值的循环(单位根)结构。
- 在量化中应用:Leontief 投入产出模型与 M 矩阵、风险传染网络的谱半径与关键连边识别、特征向量中心性。
读前导读
这一章在解决什么问题
一句话:一个所有元素都非负的矩阵反复作用很多次之后,结果长什么样?
你在 CFA 信用风险部分见过评级迁移矩阵:一年后 AAA 有 90% 留在 AAA、8% 降到 AA……把它连乘 10 次就是 10 年迁移矩阵。问题是:连乘很多次之后会稳定吗?稳定到什么分布?多快稳定?这类问题在量化里反复出现:Markov 链的长期分布、银行间损失一轮轮传染会不会放大、产业链上一个需求冲击最终带动多少总产出、网络里哪个节点最重要(PageRank)。
这些矩阵有个共同点:元素都是比例、概率、敞口,不可能为负。Perron–Frobenius 定理告诉你,正是"非负"这个看似平凡的条件,带来了很强的结构:最大的那个增长率 \(\rho(A)\) 一定是个正实数(不会是复数或负数),它对应的方向(Perron 向量)各分量都为正,可以直接解释成"长期分布"或"重要性得分"。
如果你学过凯恩斯乘数 \(1/(1-c)\)(\(c\) 是边际消费倾向),本章的 Leontief 逆 \((I-A)^{-1}=I+A+A^2+\cdots\) 就是它的矩阵版:一轮支出带来下一轮支出,乘数存在的条件从"\(c<1\)"变成"\(\rho(A)<1\)"。
需要先想起来的数学
1. 特征值与特征向量。 \(Ax=\lambda x\):矩阵作用在 \(x\) 上只是把它放大 \(\lambda\) 倍,方向不变。例:\(\begin{bmatrix}0.9&0.3\\0.1&0.7\end{bmatrix}\begin{bmatrix}3\\1\end{bmatrix}=\begin{bmatrix}3\\1\end{bmatrix}\),所以 \((3,1)\) 是特征值 1 的特征向量。反复作用 \(m\) 次,沿特征方向的分量被放大 \(\lambda^m\) 倍,所以模最大的特征值主导长期行为。见 第 00 册第 06 章 线性代数速成。
2. 谱半径。 \(\rho(A)=\max|\lambda_i|\),所有特征值绝对值(复数时取模)的最大者。它是矩阵"每作用一次大约放大多少倍"的长期平均。\(\rho<1\) 则 \(A^m\to0\),\(\rho>1\) 则一般会爆炸。
3. 几何级数。 \(1+a+a^2+\cdots=1/(1-a)\) 当且仅当 \(|a|<1\)。年金现值公式就是它。矩阵版 \(\sum_kA^k=(I-A)^{-1}\) 成立的条件是 \(\rho(A)<1\),叫 Neumann 级数。见 第 00 册第 04 章 级数与收敛。
4. 极限与"子列收敛"。 \(A^m\) 随 \(m\to\infty\) 趋于某个矩阵,是逐元素的极限。8.4.1 的证明里用到"紧性":在一个有界闭集(例如所有分量非负、和为 1 的向量)里任取一串点,总能挑出一个收敛的子序列。见 第 00 册第 01 章 函数极限与连续。
5. 有向图。 把矩阵看成网络:\(a_{ij}>0\) 就画一条从 \(i\) 到 \(j\) 的箭头。"强连通"指从任何节点出发沿箭头都能走到任何其他节点。\((A^k)_{ij}>0\) 当且仅当存在一条恰好 \(k\) 步的路径 \(i\to j\)。
本章符号速查。 \(x>0\) 指向量每个分量都 \(>0\),\(x\ge0\) 指每个分量 \(\ge0\)(不是"向量大于零"这种模糊说法)。\(e\) 是全 1 向量,\(Ae\) 就是行和向量。\(J_n\) 是全 1 矩阵。\(|x|\)、\(|A|\) 是逐元素取绝对值。\(\sigma(A)\) 是全部特征值的集合(谱),不是标准差。\(\operatorname{adj}\) 是伴随矩阵,满足 \(A\operatorname{adj}A=(\det A)I\)。\(E_{ij}\) 是只有 \((i,j)\) 位置为 1 的矩阵。"代数单重"指特征值作为特征多项式的根只出现一次;"几何单重"指它的特征向量方向唯一;"半单"指 Jordan 块都是 \(1\times1\),没有 \(\begin{bmatrix}1&1\\0&1\end{bmatrix}\) 那样会让幂线性增长的结构。
怎么读这一章
核心必读:8.1 两城市例子(把它读透,全章的结论在它身上都看得见)、8.2.3–8.2.4(谱半径单调和行和界,实务中最常用)、8.3.2 Perron 定理的六条陈述、8.5.1–8.5.2(不可约与 Perron–Frobenius 定理)、8.5.3 灵敏度公式。
第一次可以只看结论:8.3.1 的证明链条(Lemma 8.2.1 建议读,其余跳过)、8.3.3 和 8.5.5 的习题罗列、8.4.1 中 8.3.5 半单性、8.5.4 Wielandt 定理的证明。8.4.3 的 M 矩阵结合实战 2 读,效果最好。
建议顺序:8.1 → 8.2 → 8.3.2 → 实战 1 → 8.5.1–8.5.3 → 实战 3 → 回头读 8.4 和实战 2 → 8.5.4。
8.1 引言:人口迁移模型
原书用一个人口迁移模型引出全章。\(n\) 个城市,每天早上城市 \(j\) 当前人口的固定比例 \(a_{ij}\) 迁往城市 \(i\)(\(a_{jj}\) 是留下的比例)。第 \(m\) 天的人口向量满足
这样的 \(A\) 叫列随机矩阵。(量化中 Markov 链习惯用行随机矩阵 \(P\)、行向量分布 \(\pi_{t+1}=\pi_tP\),两者互为转置,第 08b 章统一处理。)长期人口分布取决于 \(A^m\) 的渐近行为。
两城市详解。 设 \(a_{21}=\alpha\)(城市 1 迁往城市 2 的比例),\(a_{12}=\beta\):
特征值为 \(1\) 与 \(1-\alpha-\beta\);\(0\le\alpha,\beta\le1\) 时 \(|1-\alpha-\beta|\le1=\rho(A)\),谱半径本身就是特征值。对应 1 的特征向量 \(x=(\beta,\alpha)^T\) 非负。若 \(\alpha,\beta\) 不同时为 0 也不同时为 1,则 \(|1-\alpha-\beta|<1\),
均衡分布与初始分布无关,正比于特征向量 \(x\);\(A^m\) 的极限是秩一矩阵,每列都正比于 \(x\)。两种例外:\(\alpha=\beta=0\) 时 \(A=I\)(两城互不往来,极限依赖初值);\(\alpha=\beta=1\) 时 \(A=\begin{bmatrix}0&1\\1&0\end{bmatrix}\),两城每天整体互换,\(A^m\) 不收敛,但 Cesàro 平均 \(\frac1m\sum_{k=1}^mA^k\to\begin{bmatrix}.5&.5\\.5&.5\end{bmatrix}\) 收敛。
推导拆解:两个特征值怎么来的?第一步,每列和为 1 意味着 \(e^TA=e^T\)(\(e=(1,1)^T\)),所以 1 是 \(A^T\) 的特征值,也就是 \(A\) 的特征值。第二步,特征值之和等于迹:\(1+\lambda_2=(1-\alpha)+(1-\beta)\),得 \(\lambda_2=1-\alpha-\beta\)。第三步,验证 \(x=(\beta,\alpha)^T\):\(A x\) 的第一个分量是 \((1-\alpha)\beta+\beta\alpha=\beta\),第二个是 \(\alpha\beta+(1-\beta)\alpha=\alpha\),确实不变。 第四步,为什么极限与初值无关?把初值分解成"沿 \(x\) 的部分 + 沿第二特征向量的部分",每天前者不变、后者乘以 \(1-\alpha-\beta\)。\(|1-\alpha-\beta|<1\) 时后者指数衰减,只剩前者。总人口 \(p_1+p_2\) 守恒,所以极限由总人口和 \(x\) 的方向完全决定。 金融直觉:把两城换成"正常/违约前预警"两个信用状态,\(\alpha\) 是每期恶化概率,\(\beta\) 是修复概率。长期有 \(\alpha/(\alpha+\beta)\) 的债务人处于预警状态,向这个比例收敛的速度是 \((1-\alpha-\beta)^m\)。两个状态之间来回越频繁(\(\alpha+\beta\) 越大),越快忘掉初始状态。
这个例子预示了本章要对一般 \(n\) 证明的五条结论:
- 谱半径 \(\rho(A)\) 本身是特征值(不只是某个特征值的模);
- 对应特征向量可取为非负,在"不可约"条件下为正;
- 若 \(A\) 的元素全正,\(\rho(A)\) 是单特征值,且严格大于其他所有特征值的模;
- 若 \(A>0\),\(\lim(A/\rho(A))^m\) 存在,是每列正比于 \(x\) 的秩一矩阵;
- 即使有零元素,Cesàro 平均 \(\frac1m\sum_{k=1}^m(A/\rho(A))^k\) 仍收敛(第 08b 章)。
8.0.P1 提醒:\(\begin{bmatrix}1&1\\0&1\end{bmatrix}\) 的谱半径为 1,幂却无界——"谱半径为 1"本身不保证幂有界,非负矩阵的结构才是关键。8.0.P2 给出一族 \(A_\epsilon\),\(\lim A_\epsilon^m=xy^T\) 在 \(\epsilon\to0\) 时发散:极限矩阵对矩阵元素可以非常敏感。
8.2 逐元素序与基本不等式
8.2.1 记号
对实矩阵,\(A\ge0\) 表示所有 \(a_{ij}\ge0\)(非负矩阵,nonnegative matrix),\(A>0\) 表示所有 \(a_{ij}>0\)(正矩阵,positive matrix);\(A\ge B\) 表示 \(A-B\ge0\)。\(|A|=[|a_{ij}|]\) 是逐元素绝对值。注意与第 07d 章 Loewner 序 \(\succeq\) 的区别:这里比较的是元素,不是二次型。
基本事实(8.1.1–8.1.17):\(|A+B|\le|A|+|B|\);\(|AB|\le|A||B|\),\(|A^m|\le|A|^m\);\(0\le A\le B\)、\(0\le C\le D\) ⇒ \(AC\le BD\);\(A>0\)、\(x\ge0\) 且 \(x\ne0\) ⇒ \(Ax>0\);\(A\ge0\)、\(x>0\)、\(Ax=0\) ⇒ \(A=0\);\(\|A\|_F=\||A|\|_F\)。
8.2.2 三角不等式及其等号
Proposition 8.1.8。 (a) \(|Ax|\le|A||x|\);(b) 若 \(A\ge0\) 有一行全正且 \(|Ax|=A|x|\),则存在 \(\theta\) 使 \(e^{-i\theta}x=|x|\);(c) 若 \(x>0\) 且 \(Ax=|A|x\),则 \(A=|A|\)。
(b) 的依据是复数三角不等式的等号条件(原书附录 A,本册附录 A.1):\(|\sum z_k|=\sum|z_k|\) 当且仅当所有 \(z_k\) 在同一条射线上。这条"等号 ⇒ 相位一致"的论证是本章证明的基本技巧。
8.2.3 谱半径的单调性
Theorem 8.1.18。 \(B\ge0\)、\(|A|\le B\) ⇒ \(\rho(A)\le\rho(|A|)\le\rho(B)\)。
证明:\(|A^m|\le|A|^m\le B^m\),于是 \(\|A^m\|_F\le\||A|^m\|_F\le\|B^m\|_F\),开 \(m\) 次方令 \(m\to\infty\),用 Gelfand 公式 \(\rho(A)=\lim\|A^m\|^{1/m}\)(第 05b 章)。
推导拆解:第一步,\(|A^m|\le|A|^m\) 是把 \(|AB|\le|A||B|\) 用 \(m-1\) 次(每个元素是乘积之和,和的绝对值不超过绝对值之和)。\(|A|^m\le B^m\) 是因为非负矩阵相乘保持逐元素序。第二步,Frobenius 范数只看元素绝对值的平方和,逐元素小则范数小。第三步,开 \(m\) 次方再取极限。 金融直觉:Gelfand 公式可以用"年化收益"来理解。\(\|A^m\|\) 是 \(m\) 期之后的总放大倍数,\(\|A^m\|^{1/m}\) 是每期的几何平均放大倍数,正如 \(m\) 年累计收益开 \(m\) 次方得到年化收益。期数越长,初期的波动被平均掉,剩下的就是长期增长率 \(\rho(A)\)。这样看,定理就是说"每条边都更大的网络,长期增长率不会更小"。
推论:(8.1.19) \(0\le A\le B\Rightarrow\rho(A)\le\rho(B)\);(8.1.20) 非负矩阵任一主子矩阵的谱半径不超过整体,\(\max_ia_{ii}\le\rho(A)\)。非负性不可去:\(\begin{bmatrix}1&1\\-1&-1\end{bmatrix}\) 幂零,\(\rho=0<1=a_{11}\)。
量化读法:传染矩阵中任何一条敞口增加,系统的放大倍数(谱半径)都不会减少;只看一部分机构(主子矩阵)得到的放大倍数是整体的下界。
8.2.4 行和、列和夹逼
Lemma 8.1.21。 \(A\ge0\) ⇒ \(\rho(A)\le\) 最大行和,\(\rho(A)\le\) 最大列和;行和全相等时 \(\rho(A)\) 就等于行和(\(e\) 是特征向量)。
Theorem 8.1.22。 \(A\ge0\):
下界的证明(出人意料地简单):设最小行和 \(\alpha>0\),把每行等比例缩小成行和恰为 \(\alpha\) 的矩阵 \(B\),\(0\le B\le A\),\(\rho(B)=\alpha\),由单调性 \(\alpha\le\rho(A)\)。
用对角相似 \(S=\operatorname{diag}(x)\)(\(x>0\))不改变谱,把同样的论证用于 \(S^{-1}AS=[a_{ij}x_j/x_i]\):
Theorem 8.1.26(加权界)。 \(A\ge0\),任意 \(x>0\):
白话解释:\((Ax)_i/x_i\) 是"节点 \(i\) 在一步之后变成原来的几倍"。如果你随便给一个正的初始分布 \(x\),各节点的一步增长倍数有高有低,真正的长期增长率 \(\rho(A)\) 一定夹在最低和最高之间。\(x\) 越接近 Perron 向量,各节点的倍数越接近,区间越窄;\(x\) 恰好是 Perron 向量时,所有节点倍数相同,都等于 \(\rho(A)\)。 为什么能用 \(S^{-1}AS\)?对角相似 \([a_{ij}x_j/x_i]\) 相当于把每个节点换一种计量单位(例如节点 \(i\) 用"千元"、节点 \(j\) 用"万元"计)。换单位改变了行和,但不改变长期增长率(特征值在相似变换下不变)。所以可以挑一个最有利的单位制来夹逼。 (Collatz–Wielandt 型比较)。** \(x>0\),\(\alpha x\le Ax\le\beta x\) ⇒ \(\alpha\le\rho(A)\le\beta\);严格不等号给出严格结论。
Corollary 8.1.30。 非负矩阵的正特征向量必然对应 \(\rho(A)\):\(x>0\)、\(Ax=\lambda x\) ⇒ \(\lambda=\rho(A)\)。
Corollary 8.1.31。 若 \(A\ge0\) 有正特征向量,则
8.1.P7(幂迭代的一步)。 取 \(x=e\) 得行和界;取 \(x=Ae\)(行和向量)得到更紧的界;不断迭代 \(x\leftarrow Ax\),上下界收敛到 \(\rho(A)\)——这就是幂迭代,并且每一步都给出可验证的上下界(实战 1)。
Corollary 8.1.33(幂的增长)。 若 \(A\ge0\) 有正特征向量 \(x\),则 \(A^m\) 的每个行和夹在 \(\frac{\min x_k}{\max x_k}\rho(A)^m\) 与 \(\frac{\max x_k}{\min x_k}\rho(A)^m\) 之间;\((A/\rho(A))^m\) 的元素一致有界。
8.1.P11:\(\rho(A)\ge(a_{1\sigma(1)}\cdots a_{n\sigma(n)})^{1/n}\),任一"环"上元素的几何平均是谱半径的下界。
8.3 正矩阵:Perron 定理
Perron (1907) 对元素全正的矩阵建立了最干净的理论。
8.3.1 证明的主线
Lemma 8.2.1。 \(A>0\),\(Ax=\lambda x\),\(|\lambda|=\rho(A)\),则 \(|x|>0\) 且 \(A|x|=\rho(A)|x|\)。
证明:令 \(z=A|x|>0\),\(z\ge|Ax|=\rho(A)|x|\),\(y=z-\rho(A)|x|\ge0\)。若 \(y\ne0\),则 \(Ay>0\),即 \(Az>\rho(A)z\),由 8.1.29 得 \(\rho(A)>\rho(A)\),矛盾。
推导拆解:逐步看每一行用了什么。 ① \(z=A|x|>0\):\(A\) 每个元素为正,\(|x|\) 非负且不全为零,所以每个分量都是正数之和(8.2.1 节"\(A>0\)、\(x\ge0\) 非零 ⇒ \(Ax>0\)")。 ② \(z\ge|Ax|\):三角不等式 \(|Ax|\le|A||x|=A|x|\)(\(A\) 本身非负,\(|A|=A\))。而 \(|Ax|=|\lambda x|=\rho(A)|x|\)。 ③ 于是 \(y=z-\rho(A)|x|\ge0\)。假设 \(y\) 不全为零,左乘 \(A\) 得 \(Ay>0\),即 \(Az-\rho(A)A|x|=Az-\rho(A)z>0\)。 ④ \(Az>\rho(A)z\) 且 \(z>0\),说明每个节点的一步增长倍数都严格大于 \(\rho(A)\),由 Collatz–Wielandt 下界得 \(\rho(A)>\rho(A)\),不可能。所以 \(y=0\),即 \(A|x|=\rho(A)|x|\),且 \(|x|=z/\rho(A)>0\)。 要点:把任意最大模特征向量取绝对值,它就自动变成谱半径的正特征向量。这是"非负性带来的刚性"的第一次体现。
Theorem 8.2.2。 \(A>0\) ⇒ 存在正向量 \(x,y\) 使 \(Ax=\rho(A)x\)、\(y^TA=\rho(A)y^T\)。由此 Collatz–Wielandt 公式对正矩阵成立。
Lemma 8.2.3 → Theorem 8.2.4(严格占优)。 若 \(|\lambda|=\rho(A)\),由 8.1.8(b) 特征向量可以旋转成正向量,再由 8.1.30 得 \(\lambda=\rho(A)\)。所以正矩阵的谱圆上只有 \(\rho(A)\) 一个特征值。
Theorem 8.2.5(几何单重)。 若 \(p,q>0\) 都是 \(\rho(A)\) 的特征向量,令 \(\beta=\min q_i/p_i\),\(r=q-\beta p\ge0\) 有零分量;若 \(r\ne0\) 则 \(Ar=\rho(A)r>0\) 矛盾,故 \(q=\beta p\)。
Theorem 8.2.7(代数单重与极限)。 由 \(y^Tx>0\) 得 \(\rho(A)\) 代数单重(第 01 章:左右特征向量不正交 ⇒ 单重)。于是 \(A=S([\rho(A)]\oplus B)S^{-1}\),\(\rho(B)<\rho(A)\),
8.3.2 Perron 定理
定义。 满足 \(Ax=\rho(A)x\)、\(\sum x_i=1\) 的正向量 \(x\) 叫(右)Perron 向量,\(\rho(A)\) 叫 Perron 根;\(A^T\) 的对应特征向量 \(y\) 按 \(y^Tx=1\) 归一化后叫左 Perron 向量。
Theorem 8.2.8(Perron 定理)。 \(A\in M_n\),\(A>0\),则
- (a) \(\rho(A)>0\);
- (b) \(\rho(A)\) 是代数单重特征值;
- (c) 存在唯一的 \(x>0\),\(Ax=\rho(A)x\),\(\sum x_i=1\);
- (d) 存在唯一的 \(y>0\),\(y^TA=\rho(A)y^T\),\(y^Tx=1\);
- (e) 其余特征值都满足 \(|\lambda|<\rho(A)\);
- (f) \((A/\rho(A))^m\to xy^T\)。
白话解释:(f) 是最有用的一条。它说对任何初始向量 \(p_0\),\(A^mp_0\approx\rho(A)^m\,x\,(y^Tp_0)\)。三个因子各有分工:\(\rho(A)^m\) 是总量的增长,\(x\) 是长期的形状(各节点占比,与初值无关),\(y^Tp_0\) 是初值的加权规模,\(y_i\) 告诉你"初始时放在节点 \(i\) 的一单位,长期贡献多少"。 金融直觉:在评级迁移、资金流动这类模型里,\(x\) 回答"长期大家分布在哪里",\(y\) 回答"从哪里出发影响最大"。实战 3 把 \(x\) 叫脆弱性、\(y\) 叫系统重要性,就是这个分工。对称矩阵时 \(x\) 与 \(y\) 方向相同,非对称网络里两者可以差别很大。
收敛速度(8.2.10–8.2.11)。 \(\|(A/\rho(A))^m-xy^T\|_\infty\le Cr^m\),\(r\) 可取 \((|\lambda_2|/\rho(A),1)\) 中任意数,\(|\lambda_2|\) 是除 \(\rho(A)\) 外最大的特征值模,叫次特征值(secondary eigenvalue)。一个只用元素就能算的上界(Ostrowski 型):
元素越"均匀",收敛越快。它通常很保守(实战 1 中实际比值 0.12,界 0.91),但不需要算任何特征值。Markov 链的混合速度、幂迭代(PageRank)的收敛速度都由次特征值比决定。
8.3.3 相关结论
- Theorem 8.2.9(Fan):\(B\ge0\) 且 \(b_{ij}\ge|a_{ij}|\)(\(i\ne j\)),则 \(A\) 的特征值都在 \(\bigcup_i\{z:|z-a_{ii}|\le\rho(B)-b_{ii}\}\) 中——用 \(B\) 的 Perron 向量作权重的 Geršgorin 定理(第 06 章),比普通 Geršgorin 圆盘更紧。
- 8.2.P5:\(A>B>0\) ⇒ \(\rho(A)>\rho(B)\)(严格单调)。8.2.P15:只有 \(0\le A\le B\)、\(A\ne B\) 时不一定严格(\(\begin{bmatrix}0&1\\0&0\end{bmatrix}\) 与 \(\begin{bmatrix}0&2\\0&0\end{bmatrix}\)),但 \(B>0\) 时严格。
- 8.2.P6:\(\rho(A)=\sum_{i,j}a_{ij}x_j\)(\(x\) 为归一化 Perron 向量)。
- 8.2.P7:\(n\ge2\) 时正矩阵的逆不可能非负;非负矩阵的逆非负 ⇔ 它是置换矩阵乘正对角阵。
- 8.2.P8:对任意正的左右特征向量,\((A/\rho)^m\to(y^Tx)^{-1}xy^T\)。
- 8.2.P9:若最小行和或最大行和等于 \(\rho(A)\),则所有行和相等。
- 8.2.P13:\(\rho(A)=\lim_m(\operatorname{tr}A^m)^{1/m}\)。
- 8.2.P14:\(\operatorname{adj}(\rho(A)I-A)>0\),它的每列是右 Perron 向量的正倍数、每行是左 Perron 向量的正倍数——已知 \(\rho(A)\) 时不必解方程就能得到 Perron 向量。8.2.P11 给出 \(\operatorname{adj}(\rho I-A)=\gamma xy^T\),\(\gamma>0\)。
8.4 一般非负矩阵
8.4.1 能保留什么
对元素可以为零的非负矩阵,Perron 定理只有一部分能通过取极限保留下来:
Theorem 8.3.1。 \(A\ge0\) ⇒ \(\rho(A)\) 是 \(A\) 的特征值,且存在非负非零向量 \(x\) 使 \(Ax=\rho(A)x\)。
证明:\(A(\epsilon)=A+\epsilon J_n>0\)(\(J_n\) 为全 1 阵),取其 Perron 向量 \(x(\epsilon)\),\(\|x(\epsilon)\|_1=1\)。单位单纯形紧,取收敛子列 \(x(\epsilon_k)\to x\ge0\);\(\rho(A(\epsilon))\) 随 \(\epsilon\) 单调递减且 \(\ge\rho(A)\),极限 \(\rho\) 满足 \(Ax=\rho x\),于是 \(\rho\le\rho(A)\),故 \(\rho=\rho(A)\)。
Theorem 8.3.2。 \(A\ge0\),\(x\ge0\) 非零,\(Ax\ge\alpha x\) ⇒ \(\rho(A)\ge\alpha\)。
Corollary 8.3.3(Collatz–Wielandt 的 max-min 半边)。
但另一半(min-max)对一般非负矩阵不成立:\(A=\begin{bmatrix}1&0\\0&2\end{bmatrix}\)、\(x=(1,0)^T\) 满足 \(Ax\le1\cdot x\),而 \(\rho(A)=2\)。这个矩阵没有正的左或右特征向量。
Theorem 8.3.4。 若 \(A\ge0\) 有正的(左或右)特征向量,对应特征值 \(\lambda\ge0\),则 \(\lambda=\rho(A)\)。特别地,行随机矩阵(\(Ae=e\))与列随机矩阵的谱半径都是 1,幂有界。
Theorem 8.3.5。 若 \(A\ge0\) 有正的左特征向量:(a) \(Ax\ge\rho(A)x\) ⇒ \(Ax=\rho(A)x\);(b) \(A\ne0\) 时 \(\rho(A)>0\),且所有模为 \(\rho(A)\) 的特征值都是半单的(Jordan 块为 \(1\times1\))。
术语提醒。 一般非负矩阵的 Perron 根仍有定义,但它的特征向量即使归一化也未必唯一(如 \(A=I\),任何非负向量都是),所以一般非负矩阵没有良定义的"Perron 向量"。
8.4.2 Frobenius 标准形与状态分类(8.3.P8)
任何非负矩阵要么不可约,要么置换相似于分块上三角形
每个对角块不可约(可以是 \(1\times1\) 零块),且 \(\sigma(A)=\sigma(A_1)\cup\cdots\cup\sigma(A_k)\)。对角块就是有向图 \(\Gamma(A)\) 的强连通分量。对 Markov 链,这正是状态分类:常返类(闭的强连通分量)与瞬时状态;违约这样的吸收态是 \(1\times1\) 块 \([1]\)。所以评级迁移矩阵(含违约吸收态)是可约的,Perron–Frobenius 的"正特征向量、单重"结论要分块应用(第 08b 章)。
8.4.3 习题中的重要概念
- 8.3.P9 本质非负矩阵(essentially nonnegative,Metzler matrix):非对角元非负的实矩阵 \(A\)。因为 \(\lambda I+A\ge0\)(\(\lambda\) 足够大),由 8.3.1 它有实特征值 \(r(A)\),满足 \(r(A)\ge\operatorname{Re}\lambda_i\) 对所有特征值,称主导特征值。连续时间 Markov 链的生成元(转移速率矩阵 \(Q\))、线性常微分方程 \(\dot x=Ax\) 的稳定性分析都用到它:\(r(A)<0\) ⇔ 系统渐近稳定。
- 8.3.P10:\(\rho(A)\le\lambda_{\max}(\frac12(A+A^T))\)。
- 8.3.P11:幻方矩阵(\(1,\dots,n^2\) 排成行、列、对角线和相等)的谱半径是 \(\frac12n(n^2+1)\)。
- 8.3.P15 M 矩阵(M-matrix):非对角元 \(\le0\) 且所有实特征值为正的矩阵 \(A\) 满足 \(A^{-1}\ge0\)。证明:\(\mu=\max a_{ii}\),\(B=\mu I-A\ge0\),\(\mu-\rho(B)\) 是 \(A\) 的实特征值,故 \(\mu>\rho(B)\),于是 Neumann 级数 \(A^{-1}=\mu^{-1}\sum_{k\ge0}(B/\mu)^k\ge0\)。Leontief 模型的 \(I-A\)(\(\rho(A)<1\))就是 M 矩阵(实战 2)。
金融直觉:M 矩阵的证明核心是 Neumann 级数,它就是年金/永续年金公式的矩阵版。标量时 \(\frac1{1-a}=1+a+a^2+\cdots\),\(|a|<1\) 才收敛,每一项都非负,和自然非负。矩阵时每一项 \(A^k\ge0\),和 \((I-A)^{-1}\ge0\) 也自然非负,收敛条件换成 \(\rho(A)<1\)。经济读法:\(A^k d\) 是第 \(k\) 轮的间接需求,总产出是各轮之和;每轮的规模按 \(\rho(A)\) 的速度衰减,就像永续年金每期按贴现因子衰减。
- 8.3.P16 单调矩阵(monotone matrix):\(A\) 非奇异且 \(A^{-1}\ge0\) ⇔ (\(Ax\ge Ay\Rightarrow x\ge y\))。M 矩阵都是单调矩阵。经济含义:最终需求增加,各部门产出不会减少。
原书说明:8.3.P12(b) 中的"positive"。 这道题讨论 \(r>\rho(A)\) 时的 \((rI-A)^{-1}\) 与 \(\operatorname{adj}(\rho(A)I-A)\)。由 Neumann 级数 \((rI-A)^{-1}=r^{-1}\sum_k(A/r)^k\ge0\)。原书该小问写作 \((rI-A)^{-1}\) 为 "positive",对一般非负矩阵这只在 \(A\) 不可约时成立;\(A\) 可约时只能保证非负。反例:\(A=I_2\) 时 \((1.5I-A)^{-1}=2I\);\(A=\begin{bmatrix}1&1\\0&1\end{bmatrix}\) 时 \((1.5I-A)^{-1}=\begin{bmatrix}2&4\\0&2\end{bmatrix}\),都有零元(实战 3 末尾数值验证)。本教材按"非负"处理:令 \(r\downarrow\rho(A)\),\(\det(rI-A)(rI-A)^{-1}=\operatorname{adj}(rI-A)\to\operatorname{adj}(\rho(A)I-A)\),得到 \(\operatorname{adj}(\rho(A)I-A)\ge0\),这正是 (c) 要的结论,不受影响。不可约时 8.4.P23 给出更强的 \(\operatorname{adj}(\rho I-A)=c\,xy^T>0\)。
8.5 不可约非负矩阵:Perron–Frobenius 定理
8.5.1 不可约性
非负矩阵 \(A\) 称为不可约(irreducible),若它不能经置换相似化为 \(\begin{bmatrix}B&C\\0&D\end{bmatrix}\)(\(B,D\) 为方阵)的形式;等价地,有向图 \(\Gamma(A)\)(\(a_{ij}>0\) 时有边 \(i\to j\))强连通——从任一节点出发都能到达任一节点(第 06 章 6.2 节)。启发式原则:对"无零元素"矩阵成立的结论,常常可以推广到不可约矩阵。
Lemma 8.4.1。 \(A\ge0\) 不可约 ⇔ \((I+A)^{n-1}>0\)。
直观:\((I+A)^{n-1}=\sum_k\binom{n-1}{k}A^k\) 的 \((i,j)\) 元为正 ⇔ 存在长度 \(\le n-1\) 的路 \(i\to j\);强连通图中任意两点间的最短路长度不超过 \(n-1\)。实际计算只关心零/非零模式(布尔运算),并可反复平方以减少乘法次数(第 08b 章)。
推导拆解:为什么 \((A^2)_{ij}>0\) 等于"有两步路"?\((A^2)_{ij}=\sum_ka_{ik}a_{kj}\),每项非负,所以和为正当且仅当某个 \(k\) 让 \(a_{ik}>0\) 且 \(a_{kj}>0\),即存在 \(i\to k\to j\)。\(A^k\) 同理。\((I+A)^{n-1}\) 展开后包含 \(I,A,\dots,A^{n-1}\) 的正系数组合,所以它的 \((i,j)\) 元为正等于"存在 0 到 \(n-1\) 步的路"。不经过重复节点的路最多 \(n-1\) 步,所以够用。 白话解释:可约的银行网络,就是存在一群银行"只借出不借入"或"只借入不借出",冲击从外面传不进去(或传不出来)。例如三家银行,1 欠 2、2 欠 3,但 3 不欠任何人:从 3 出发走不到 1,网络可约。评级迁移矩阵因为违约是吸收态(进去出不来),也是可约的。
Lemma 8.4.3。 若 \(A\ge0\) 且某个 \(A^m>0\),则 \(\rho(A)\) 是唯一的最大模特征值,为正且代数单重(对 \(A^m\) 用 Perron 定理)。这类矩阵叫本原矩阵(第 08b 章)。
8.5.2 定理与证明
Theorem 8.4.4(Perron–Frobenius 定理)。 \(A\in M_n\) 不可约非负,\(n\ge2\),则
- (a) \(\rho(A)>0\);
- (b) \(\rho(A)\) 是代数单重特征值;
- (c) 存在唯一的 \(x>0\),\(Ax=\rho(A)x\),\(\sum x_i=1\);
- (d) 存在唯一的 \(y>0\),\(y^TA=\rho(A)y^T\),\(y^Tx=1\)。
证明:(a) 不可约矩阵没有零行,由行和下界 \(\rho(A)>0\)(8.1.25)。(b) \(\rho(I+A)=\rho(A)+1\);若 \(\rho(A)\) 多重,则 \((1+\rho(A))^{n-1}\) 是正矩阵 \((I+A)^{n-1}\) 的多重特征值,与 Perron 定理矛盾。(c) 由 8.3.1 取非负特征向量 \(x\),则 \((I+A)^{n-1}x=(1+\rho(A))^{n-1}x\),左边是正矩阵乘非负非零向量,为正,故 \(x>0\)。(d) 对 \(A^T\)。
与 Perron 定理相比,缺少两条:"\(\rho(A)\) 是唯一最大模特征值"和"\((A/\rho(A))^m\) 收敛"。反例:\(\begin{bmatrix}0&1\\1&0\end{bmatrix}\) 不可约,特征值 \(\pm1\),幂振荡。这两条要靠"本原性"才能恢复(第 08b 章)。
由 (c)(d),8.1.30–8.1.33、8.3.4–8.3.5 的结论都适用于不可约非负矩阵,特别是 Collatz–Wielandt 公式的两半都成立。
8.4 节几个配套结论:
- 8.4.P15、P21:不可约非负矩阵的任何非负特征向量都是 Perron 向量的正倍数;有两个线性无关的非负特征向量的非负矩阵必可约。
- 8.4.P14:\(A\) 不可约、\(B\ge0\) 非零 ⇒ \(\rho(A+B)>\rho(A)\)(严格单调)。
- 8.4.P3:不可约只是存在正特征向量的充分条件:可约的 \(\begin{bmatrix}1&0\\1&0\end{bmatrix}\) 有正特征向量 \((1,1)^T\),可约的 \(\begin{bmatrix}1&1\\0&0\end{bmatrix}\) 没有。
- 8.4.P5:不可约性不是相似不变量,\(AB\) 不可约时 \(BA\) 可以可约。
8.5.3 谱半径的灵敏度(8.4.P13)
设 \(A\) 不可约非负,\(x\)、\(y\) 为右、左 Perron 向量,\(y^Tx=1\)。\(\rho(A)\) 是单重特征值,所以对元素可微,且
推导拆解:设 \(A\) 沿某个方向变化,\(A(t)x(t)=\rho(t)x(t)\)。两边对 \(t\) 求导(乘积法则):\(A'x+Ax'=\rho'x+\rho x'\)。左乘 \(y^T\):\(y^TA'x+y^TAx'=\rho'y^Tx+\rho\,y^Tx'\)。因为 \(y^TA=\rho y^T\),左边第二项 \(=\rho\,y^Tx'\),与右边第二项抵消。剩下 \(\rho'=y^TA'x/y^Tx=y^TA'x\)(归一化 \(y^Tx=1\))。只改 \(a_{ij}\) 时 \(A'=E_{ij}\),\(y^TE_{ij}x=y_ix_j\)。左特征向量的妙处就在于把未知的 \(x'\) 消掉了。 正号来自 Perron–Frobenius:\(x,y\) 都是正向量,乘积必正。
这是单特征值一阶扰动公式 \(d\lambda=y^T(dA)x/y^Tx\) 在 \(dA=E_{ij}\) 时的特例。勘误说明:精读笔记把它记作 \(x_iy_j\);按"\(x\) 为右、\(y\) 为左 Perron 向量"的约定,正确的是 \(y_ix_j\)(两者只有在 \(A\) 对称时才一致),实战 3 用有限差分验证了 \(y_ix_j\)。8.4.P14 的一种证法就是把这个导数沿 \(A+tB\) 积分。
量化读法:在传染/放大模型中,"哪条连边对系统放大倍数影响最大"由 \(y_ix_j\) 回答——连边起点的左 Perron 分量与终点的右 Perron 分量之积。
8.5.4 最大模特征值的循环结构
Theorem 8.4.5(Wielandt)。 \(A\ge0\) 不可约,\(A\ge|B|\),\(\lambda=e^{i\varphi}\rho(B)\) 是 \(B\) 的最大模特征值。若 \(\rho(A)=\rho(B)\),则存在对角酉阵 \(D\) 使 \(B=e^{i\varphi}DAD^{-1}\)。
取 \(B=A\),就得到最大模特征值的结构:
Theorem 8.4.6。 \(A\ge0\) 不可约,恰有 \(k\) 个最大模特征值,则
- (a) \(A\) 相似于 \(e^{2\pi ip/k}A\),\(p=0,1,\dots,k-1\);
- (b) 整个谱(连同 Jordan 结构)在绕原点旋转 \(2\pi/k\) 下不变;
- (c) 最大模特征值恰为 \(\rho(A)e^{2\pi ip/k}\),\(p=0,\dots,k-1\),每个都代数单重。
证明要点:最大模特征值的辐角集合 \(\mathcal S\) 在模 \(2\pi\) 加法下封闭(\(A\) 相似于 \(e^{i\varphi}A\) 与 \(e^{i\psi}A\) ⇒ 相似于 \(e^{i(\varphi+\psi)}A\));有限的、对加法封闭的辐角集合只能是 \(k\) 次单位根的辐角 \(\{2\pi p/k\}\)。单重性来自 8.4.4(b) 与旋转对称。
推论:
- 任意圆周 \(|z|=r>0\) 上的特征值个数是 \(k\) 的倍数,所以 \(k\) 整除非零特征值的个数。例:不可约非负 \(3\times3\) 矩阵不可能有特征值 \(1,i,-i\)(\(k\) 只能是 1 或 3,都不符;或看 \(\operatorname{tr}A^2=1-1-1<0\),与非负矛盾)。
- Corollary 8.4.7:\(k>1\) 时 \(A\) 的对角元全为 0;更一般地,\(k\nmid m\) 时 \(A^m\) 的对角元全为 0(因为 \(\operatorname{tr}A^m=e^{2\pi im/k}\operatorname{tr}A^m\))。所以有一个正对角元的不可约非负矩阵,\(\rho(A)\) 必是唯一最大模特征值;但这不是必要条件:\(J_3-I\) 对角全零,特征值 \(2,-1,-1\)。
- 循环块形(8.4.8):\(k>1\) 时存在置换 \(P\) 使
图的节点分成 \(k\) 类,边只从第 \(j\) 类指向第 \(j+1\) 类,循环往复。这样的矩阵叫指数为 \(k\) 的循环矩阵(cyclic of index \(k\))(8.4.P9),特征多项式形如 \(t^r(t^k-\rho^k)(t^k-\mu_2^k)\cdots\)(8.4.P10)。Markov 链中这就是周期为 \(k\) 的链。
白话解释:循环块形就是"节点分成 \(k\) 组,每一步都必须从一组跳到下一组"。例如 \(k=2\):资金只在"银行 → 基金"和"基金 → 银行"之间流动,同类之间不直接往来。那么从银行出发,奇数步后一定在基金、偶数步后一定在银行,分布永远在两组之间来回摆,\(A^m\) 不会收敛。\(-\rho\) 这个特征值就是这种来回摆动的数学表现。只要任何一个节点有"自环"(对角元为正,能原地停留一步),节奏就被打乱,周期消失,这就是 Corollary 8.4.7 的含义。
8.5.5 其他习题结论
- 8.4.P12(Cauchy 根界):多项式 \(p(t)=t^n+a_{n-1}t^{n-1}+\cdots+a_0\) 的根的模不超过 \(\tilde p(t)=t^n-|a_{n-1}|t^{n-1}-\cdots-|a_0|\) 的唯一正根(对友矩阵用 Perron–Frobenius)。可用于快速判断 AR(\(p\)) 特征多项式的根是否在某个圆内。
- 8.4.P17–P18(非负矩阵的最佳秩一逼近):\(A\ge0\),若 \(AA^T\) 或 \(A^TA\) 不可约,则 Frobenius 范数下的最佳秩一逼近唯一且非负:\(\sqrt{\rho(AA^T)}\,vw^T\),\(v,w\) 是正的首奇异向量。这解释了为什么对非负数据(成交量矩阵、持仓矩阵)做 SVD 时首个奇异向量可以取成全正——它代表"整体规模"因子。\(I_n\)(\(n>1\))的最佳秩一逼近不唯一。
- 8.4.P22:\(\mathbf R^n\) 中两两夹角大于 \(\pi/2\) 的向量至多 \(n+1\) 个。
- 8.4.P25:对任意非负 \(A\),\(\limsup_m(\operatorname{tr}A^m)^{1/m}=\rho(A)\)(极限可以不存在,如 \(\begin{bmatrix}0&1\\1&0\end{bmatrix}\))。
- 8.4.P26:\(A\) 不可约 ⇔ \(\rho(A)I-A\) 没有零主子式 ⇔ \(\operatorname{adj}(\rho(A)I-A)=cxy^T\) 为正。
量化实战
实战 1:迁移极限、谱半径的可验证界与 Perron 收敛
import numpy as np
rng = np.random.default_rng(0)
rho = lambda M: max(abs(np.linalg.eigvals(M)))
# ---------- 1. 原书 8.0 两城市迁移:A = [[1-α, β], [α, 1-β]](列随机) ----------
a, b = 0.1, 0.3
A = np.array([[1 - a, b], [a, 1 - b]])
print("特征值:", np.round(np.sort(np.linalg.eigvals(A).real), 4), " (1 与 1-α-β = %.1f)" % (1 - a - b))
print("A^50 =\n", np.round(np.linalg.matrix_power(A, 50), 4))
print("理论极限 [[β,β],[α,α]]/(α+β) =\n", np.array([[b, b], [a, a]]) / (a + b))
P = np.array([[0., 1.], [1., 0.]]) # α=β=1:两城每天整体互换
cesaro = sum(np.linalg.matrix_power(P, k) for k in range(1, 1001)) / 1000
print("α=β=1 时 A^m 在 I 与 P 之间振荡;Cesàro 平均 =\n", cesaro)
# ---------- 2. 谱半径的行和夹逼与 Collatz–Wielandt 界(8.1.22、8.1.26、8.1.P7) ----------
n = 6
M = rng.uniform(0, 1, (n, n)) * (rng.uniform(size=(n, n)) < 0.6)
r = M.sum(1)
print("\nρ(M) = %.4f;行和界 [%.4f, %.4f];列和界 [%.4f, %.4f]"
% (rho(M), r.min(), r.max(), M.sum(0).min(), M.sum(0).max()))
x = M @ np.ones(n) # 取 x = Ae(幂迭代一步)
q = (M @ x) / x
print("x = Ae 的 Collatz–Wielandt 界 [%.4f, %.4f]" % (q.min(), q.max()))
for k in range(30): # 继续幂迭代,界不断收紧
x = M @ x; x /= x.sum()
q = (M @ x) / x
print("幂迭代 30 步后的界 [%.4f, %.4f]" % (q.min(), q.max()))
# ---------- 3. Perron 定理:正矩阵的幂与收敛速度(8.2.7、8.2.10–8.2.11) ----------
Apos = rng.uniform(0.2, 1.0, (n, n))
lam = np.linalg.eigvals(Apos); idx = np.argsort(-abs(lam))
r0, r2 = lam[idx[0]].real, abs(lam[idx[1]])
w, V = np.linalg.eig(Apos); x = np.abs(V[:, np.argmax(w.real)].real); x /= x.sum()
w2, U = np.linalg.eig(Apos.T); y = np.abs(U[:, np.argmax(w2.real)].real); y /= y @ x
kap = Apos.min() / Apos.max()
print("\nρ = %.4f,次特征值模比 |λ2|/ρ = %.4f ≤ (1-κ²)/(1+κ²) = %.4f(κ = min/max = %.3f)"
% (r0, r2 / r0, (1 - kap**2) / (1 + kap**2), kap))
for m in [1, 3, 6, 10]:
err = np.abs(np.linalg.matrix_power(Apos / r0, m) - np.outer(x, y)).max()
print("m = %2d ||(A/ρ)^m - x y^T||_max = %.2e (|λ2|/ρ)^m = %.2e" % (m, err, (r2 / r0) ** m))
关键输出:
特征值: [0.6 1. ] (1 与 1-α-β = 0.6)
A^50 =
[[0.75 0.75]
[0.25 0.25]]
理论极限 [[β,β],[α,α]]/(α+β) =
[[0.75 0.75]
[0.25 0.25]]
α=β=1 时 A^m 在 I 与 P 之间振荡;Cesàro 平均 =
[[0.5 0.5]
[0.5 0.5]]
ρ(M) = 2.4081;行和界 [1.5450, 3.6626];列和界 [1.5359, 3.7045]
x = Ae 的 Collatz–Wielandt 界 [2.0579, 2.6911]
幂迭代 30 步后的界 [2.4081, 2.4081]
ρ = 4.3759,次特征值模比 |λ2|/ρ = 0.1188 ≤ (1-κ²)/(1+κ²) = 0.9135(κ = min/max = 0.213)
m = 1 ||(A/ρ)^m - x y^T||_max = 1.10e-01 (|λ2|/ρ)^m = 1.19e-01
m = 3 ||(A/ρ)^m - x y^T||_max = 1.53e-03 (|λ2|/ρ)^m = 1.67e-03
m = 6 ||(A/ρ)^m - x y^T||_max = 1.73e-06 (|λ2|/ρ)^m = 2.80e-06
m = 10 ||(A/ρ)^m - x y^T||_max = 4.94e-10 (|λ2|/ρ)^m = 5.58e-10
读法:
- 两城市例子中 \(A^{50}\) 已经与理论极限一致到 4 位小数(收敛速度 \(0.6^m\));\(\alpha=\beta=1\) 时只有 Cesàro 平均收敛。
- 行和、列和给出的区间很宽;取 \(x=Ae\) 后区间收窄到 \([2.06,2.69]\),继续幂迭代 30 步后上下界重合于 \(\rho=2.4081\)。在大型稀疏网络上,这种带上下界的幂迭代比调用稠密特征值算法便宜得多,而且每一步都能报告误差范围。
- 正矩阵的 \((A/\rho)^m\) 以 \((|\lambda_2|/\rho)^m\) 的速度收敛到秩一矩阵 \(xy^T\)。Ostrowski 型界 0.91 远比实际的 0.12 保守,但它只用最大、最小元素就能算出。
实战 2:Leontief 投入产出模型——M 矩阵与冲击传导
场景:行业轮动和宏观因子研究中,常需要估计"某行业需求冲击如何沿产业链传导"。Leontief 模型:\(a_{ij}\) 为生产 1 元第 \(j\) 部门产品需要投入的第 \(i\) 部门产品(直接消耗系数),总产出 \(x\) 满足 \(x=Ax+d\)(\(d\) 为最终需求),所以 \(x=(I-A)^{-1}d\)。只有当 \((I-A)^{-1}\ge0\) 时,任何非负需求都对应非负产出——这正是 M 矩阵理论(8.3.P15–P16)。
import numpy as np
rho = lambda M: max(abs(np.linalg.eigvals(M)))
sectors = ["上游原材料", "能源", "中游制造", "建筑", "消费", "金融服务"]
# 直接消耗系数 a_ij:生产 1 元第 j 部门产品需要投入的第 i 部门产品(虚构数据)
A = np.array([[0.10, 0.15, 0.30, 0.25, 0.05, 0.00],
[0.12, 0.10, 0.12, 0.08, 0.04, 0.02],
[0.05, 0.08, 0.20, 0.30, 0.15, 0.03],
[0.01, 0.02, 0.01, 0.02, 0.01, 0.03],
[0.01, 0.01, 0.02, 0.02, 0.10, 0.05],
[0.04, 0.05, 0.06, 0.05, 0.08, 0.15]])
print("列和(每 1 元产出的中间投入):", A.sum(0).round(2), " ρ(A) = %.4f < 1" % rho(A))
# ---------- 1. Leontief 逆 (I-A)^{-1} = Σ A^k ≥ 0:M 矩阵(原书 8.3.P15) ----------
I = np.eye(6)
Linv = np.linalg.inv(I - A)
print("(I-A)^{-1} 最小元素 = %.4f(非负;A 不可约时全正)" % Linv.min())
S, Ak = I.copy(), I.copy()
for k in range(1, 60):
Ak = Ak @ A; S += Ak
print("Neumann 级数 60 项与逆矩阵的差: %.2e" % np.abs(S - Linv).max())
print("产出乘数(列和:最终需求 +1 元带动的总产出):")
for s_, m in zip(sectors, Linv.sum(0)):
print(" %-6s %.3f" % (s_, m))
# ---------- 2. 冲击传导:建筑业最终需求下降 100 ----------
dd = np.zeros(6); dd[3] = -100
dx = Linv @ dd
print("\n建筑需求 -100 → 各部门总产出变化:", dict(zip(sectors, dx.round(1).tolist())))
print("直接效应 %.1f,第 1 轮间接 %.1f,第 2 轮 %.1f ……"
% (dd.sum(), (A @ dd).sum(), (A @ A @ dd).sum()))
# ---------- 3. 可行性 ⇔ ρ(A) < 1:所有投入系数等比例放大 ----------
for c in [2.0, 2.5]:
A_bad = c * A
print("\n投入系数放大 %.1f 倍: ρ = %.4f;(I-A)^{-1} 最小元素 = %.3f" % (c, rho(A_bad), np.linalg.inv(I - A_bad).min()))
关键输出:
列和(每 1 元产出的中间投入): [0.33 0.41 0.71 0.72 0.43 0.28] ρ(A) = 0.4453 < 1
(I-A)^{-1} 最小元素 = 0.0193(非负;A 不可约时全正)
Neumann 级数 60 项与逆矩阵的差: 4.44e-16
产出乘数(列和:最终需求 +1 元带动的总产出):
上游原材料 1.582
能源 1.734
中游制造 2.290
建筑 2.380
消费 1.817
金融服务 1.489
建筑需求 -100 → 各部门总产出变化: {'上游原材料': -48.0, '能源': -22.2, '中游制造': -45.6, '建筑': -103.9, '消费': -4.8, '金融服务': -13.4}
直接效应 -100.0,第 1 轮间接 -72.0,第 2 轮 -36.5 ……
投入系数放大 2.0 倍: ρ = 0.8905;(I-A)^{-1} 最小元素 = 0.199
投入系数放大 2.5 倍: ρ = 1.1131;(I-A)^{-1} 最小元素 = -5.794
读法:
- 所有列和都小于 1,由列和上界 \(\rho(A)<1\) 立即可知模型可行,无需算特征值。\(\rho(A)=0.45\)。
- Leontief 逆是 Neumann 级数 \(\sum A^k\):第 \(k\) 项是"第 \(k\) 轮间接需求"。建筑需求下降 100,直接效应 −100,第一轮间接 −72(向上游原材料、中游制造采购减少),第二轮 −36.5……总计约 −238,即建筑业的产出乘数 2.38。逆矩阵所有元素为正(\(A\) 不可约),所以任何一个部门的需求冲击都会传到所有部门。
- 产出乘数排序(建筑、中游制造最高)可用于判断"稳增长政策"的受益链条与传导强度;乘积 \(x_iy_j\) 型灵敏度(8.5.3 节)则告诉我们哪个投入系数的变化对整体乘数影响最大。
- 系数放大到 \(\rho(A)>1\) 时,\((I-A)^{-1}\) 出现负元素:需求增加会让某些部门的"产出"减少,模型失去经济意义(对应经济无法自我维持)。\(\rho(A)<1\iff I-A\) 是非奇异 M 矩阵 \(\iff(I-A)^{-1}\ge0\)。
实战 3:风险传染网络——谱半径、关键连边与中心性
场景:8 家银行,\(w_{ij}\) 表示银行 \(i\) 对银行 \(j\) 的敞口占银行 \(i\) 资本的比例。线性化的损失级联(DebtRank 类模型)\(\ell^{(t+1)}=W\ell^{(t)}\):\(\rho(W)<1\) 时冲击逐轮衰减,总放大倍数约为 \(1/(1-\rho)\);\(\rho\) 越接近 1 越危险。监管关心:哪些机构最重要、削减哪条敞口最能降低系统放大倍数。
import numpy as np
from scipy.sparse.csgraph import connected_components
rng = np.random.default_rng(42)
n = 8
banks = [f"B{i}" for i in range(n)]
# ---------- 1. 银行间敞口传染矩阵 W:w_ij = 银行 i 对银行 j 的敞口 / 银行 i 的资本 ----------
W = rng.uniform(0.0, 0.25, (n, n)) * (rng.uniform(size=(n, n)) < 0.35)
np.fill_diagonal(W, 0)
W[0, 1] = 0.6; W[1, 2] = 0.5; W[2, 0] = 0.55 # 三家大行之间的互相持有(一个环)
# 不可约检验:(I+W)^{n-1} > 0 ⇔ 有向图强连通(原书 8.4.1)
irr = (np.linalg.matrix_power(np.eye(n) + (W > 0), n - 1) > 0).all()
ncomp, lab = connected_components(W > 0, directed=True, connection='strong')
print("不可约? %s;强连通分量个数 %d,标签 %s" % (irr, ncomp, lab))
lam, V = np.linalg.eig(W); k = np.argmax(lam.real)
r = lam[k].real; x = np.abs(V[:, k].real); x /= x.sum()
lam2, U = np.linalg.eig(W.T); y = np.abs(U[:, np.argmax(lam2.real)].real); y /= y @ x
print("ρ(W) = %.4f(DebtRank 类模型中 ρ<1 表示损失级联会衰减;越接近 1 放大越强)" % r)
print("放大倍数 1/(1-ρ) = %.2f" % (1 / (1 - r)))
# ---------- 2. 谱半径对每条边的灵敏度:∂ρ/∂w_ij = y_i x_j(原书 8.4.P13;y 左、x 右 Perron 向量,y^T x = 1) ----------
Sens = np.outer(y, x) * (W > 0)
i, j = np.unravel_index(np.argmax(Sens), Sens.shape)
h = 1e-6; W2 = W.copy(); W2[i, j] += h
fd = (max(abs(np.linalg.eigvals(W2))) - r) / h
print("\n最关键的边 w[%s,%s]:解析导数 y_i x_j = %.4f,有限差分 %.4f" % (banks[i], banks[j], y[i] * x[j], fd))
top = np.argsort(-Sens, axis=None)[:4]
print("灵敏度前 4 的边:", [(banks[t // n], banks[t % n], round(float(Sens.flat[t]), 4)) for t in top])
W3 = W.copy(); W3[i, j] *= 0.5
print("把该边敞口砍半:ρ 从 %.4f 降到 %.4f" % (r, max(abs(np.linalg.eigvals(W3)))))
# ---------- 3. 特征向量中心性:右 Perron 向量 x("脆弱性")与左 Perron 向量 y("系统重要性") ----------
print("\n右 Perron 向量 x(受传染程度):", dict(zip(banks, x.round(3).tolist())))
print("左 Perron 向量 y(系统重要性,展示时归一为和 1):", dict(zip(banks, (y / y.sum()).round(3).tolist())))
# ---------- 4. 原书 8.3.P12(b) 的说明:可约非负矩阵 (rI-A)^{-1} 只能保证非负 ----------
for name, A in [("I_2", np.eye(2)), ("[[1,1],[0,1]]", np.array([[1., 1.], [0., 1.]]))]:
rr = 1.5
print("A = %-14s r = %.1f > ρ(A)=1: (rI-A)^{-1} =" % (name, rr), np.linalg.inv(rr * np.eye(2) - A).round(3).tolist())
关键输出:
不可约? True;强连通分量个数 1,标签 [0 0 0 0 0 0 0 0]
ρ(W) = 0.7132(DebtRank 类模型中 ρ<1 表示损失级联会衰减;越接近 1 放大越强)
放大倍数 1/(1-ρ) = 3.49
最关键的边 w[B0,B1]:解析导数 y_i x_j = 0.2714,有限差分 0.2714
灵敏度前 4 的边: [('B0', 'B1', 0.2714), ('B1', 'B2', 0.2616), ('B2', 'B0', 0.2561), ('B3', 'B0', 0.2177)]
把该边敞口砍半:ρ 从 0.7132 降到 0.6166
右 Perron 向量 x(受传染程度): {'B0': 0.228, 'B1': 0.211, 'B2': 0.225, 'B3': 0.117, 'B4': 0.07, 'B5': 0.037, 'B6': 0.07, 'B7': 0.041}
左 Perron 向量 y(系统重要性,展示时归一为和 1): {'B0': 0.215, 'B1': 0.195, 'B2': 0.188, 'B3': 0.16, 'B4': 0.108, 'B5': 0.023, 'B6': 0.108, 'B7': 0.002}
A = I_2 r = 1.5 > ρ(A)=1: (rI-A)^{-1} = [[2.0, 0.0], [0.0, 2.0]]
A = [[1,1],[0,1]] r = 1.5 > ρ(A)=1: (rI-A)^{-1} = [[2.0, 4.0], [0.0, 2.0]]
读法:
- 网络强连通,Perron–Frobenius 定理适用:\(\rho(W)=0.71\) 是单重特征值,左右 Perron 向量都为正。
- 解析灵敏度 \(y_ix_j\) 与有限差分完全一致。最关键的三条边恰好是三家大行之间的互持环 \(B0\to B1\to B2\to B0\):环会让损失反复回流。只把 \(w_{B0,B1}\) 砍半,放大倍数 \(1/(1-\rho)\) 就从 3.49 降到 2.61。这给出了"针对性降低关联"的量化依据,比"按敞口大小排序"更准确,因为它考虑了网络位置。
- 右 Perron 向量 \(x\) 衡量在系统性级联中各机构受损的相对程度(脆弱性),左 Perron 向量 \(y\) 衡量初始冲击发生在某机构时对整个系统的影响(系统重要性)。两者在非对称网络中可以差别很大:B7 几乎不影响别人(\(y\) 只有 0.002),但自身仍会被波及。同样的思路用于资产相关网络、供应链网络的特征向量中心性。
- 最后两行是 8.3.P12(b) 说明中的反例:可约非负矩阵的 \((rI-A)^{-1}\) 有零元,只是非负而非正。
本章小结
非负矩阵用逐元素序比较,\(|Ax|\le|A||x|\) 和 Gelfand 公式给出谱半径关于逐元素序单调:\(|A|\le B\Rightarrow\rho(A)\le\rho(B)\)。谱半径夹在最小与最大行和(列和)之间,用任意正向量加权得到 Collatz–Wielandt 界 \(\min_i(Ax)_i/x_i\le\rho(A)\le\max_i(Ax)_i/x_i\),幂迭代让它不断收紧;正特征向量一定属于谱半径。Perron 定理说正矩阵的谱半径为正、代数单重、严格大于其他特征值的模,有唯一的正的左右 Perron 向量,且 \((A/\rho)^m\to xy^T\),速度由次特征值比决定。一般非负矩阵只保留"\(\rho(A)\) 是特征值、有非负特征向量"和 max-min 刻画;有正左特征向量时最大模特征值半单。Frobenius 标准形把任意非负矩阵分解成不可约块,对应有向图的强连通分量和 Markov 链的状态分类。不可约(\((I+A)^{n-1}>0\),图强连通)时,Perron–Frobenius 定理恢复了单重性与正的左右 Perron 向量,但最大模特征值可以有 \(k\) 个,均匀分布在半径 \(\rho\) 的圆上,整个谱在旋转 \(2\pi/k\) 下不变,矩阵可置换成循环块形。\(\partial\rho/\partial a_{ij}=y_ix_j\) 给出谱半径对每个元素的灵敏度。M 矩阵 \(I-A\)(\(\rho(A)<1\))的逆非负,是 Leontief 模型可行性的条件;Metzler 矩阵有实的主导特征值。
| 概念/公式 | 表达式 | 用途 |
|---|---|---|
| 逐元素序 | \(A\ge B\iff a_{ij}\ge b_{ij}\) | 与 Loewner 序区分 |
| 谱半径单调 | \(\vert A\vert \le B\Rightarrow\rho(A)\le\rho(B)\) | 连边增加,放大不减 |
| 行和界 | \(\min_i\sum_ja_{ij}\le\rho\le\max_i\sum_ja_{ij}\) | 快速判断 \(\rho<1\) |
| Collatz–Wielandt | \(\min_i\frac{(Ax)_i}{x_i}\le\rho\le\max_i\frac{(Ax)_i}{x_i}\) | 带误差界的幂迭代 |
| 正特征向量 | \(x>0\),\(Ax=\lambda x\Rightarrow\lambda=\rho(A)\) | 识别 Perron 根 |
| Perron 定理 | \(A>0\):\(\rho\) 单重、严格占优,\(x,y>0\),\((A/\rho)^m\to xy^T\) | 长期分布 |
| 次特征值界 | \(\vert \lambda_2\vert /\rho\le(1-\kappa^2)/(1+\kappa^2)\) | 收敛速度 |
| 一般非负 | \(\rho(A)\) 是特征值,有 \(x\ge0\) | 极限论证 |
| Frobenius 标准形 | 分块上三角,对角块不可约 | 状态分类、吸收态 |
| 不可约判据 | \((I+A)^{n-1}>0\) | 图强连通 |
| Perron–Frobenius | 不可约:\(\rho>0\) 单重,\(x,y>0\) | 中心性、平稳分布 |
| 循环结构 | 最大模特征值 \(=\rho e^{2\pi ip/k}\) | 周期链 |
| 灵敏度 | \(\partial\rho/\partial a_{ij}=y_ix_j\) | 关键连边 |
| M 矩阵 | \(\rho(A)<1\Rightarrow(I-A)^{-1}=\sum A^k\ge0\) | Leontief、产出乘数 |
| Metzler 矩阵 | 主导特征值 \(r(A)\ge\operatorname{Re}\lambda_i\) | 连续时间链、ODE 稳定 |
练习
基础
- 对两城市迁移矩阵 \(A=\begin{bmatrix}0.9&0.3\\0.1&0.7\end{bmatrix}\),求特征值、Perron 向量和 \(\lim A^m\),并估计 \(A^m\) 与极限的误差何时小于 \(10^{-6}\)。 答案要点:特征值 1、0.6;\(x=(0.75,0.25)^T\);误差约 \(0.6^m\),\(m\ge28\)。
- 不算特征值,判断 \(\begin{bmatrix}0.2&0.5&0.1\\0.3&0.1&0.4\\0.1&0.2&0.3\end{bmatrix}\) 的谱半径是否小于 1。 提示:行和 0.8、0.8、0.6,最大行和 \(<1\)。
- 举例说明非负矩阵的非负特征向量不一定对应谱半径(8.3.P4),但正特征向量一定对应。 答案要点:\(\operatorname{diag}(1,2)\) 的 \(e_1\) 是非负特征向量,对应 1 而非 \(\rho=2\)。
- 判断 \(\begin{bmatrix}0&1&0\\0&0&1\\1&0&0\end{bmatrix}\) 是否不可约,求它的最大模特征值个数 \(k\),并验证 \(k\) 整除非零特征值个数。 答案要点:循环置换矩阵,不可约;特征值为 3 个三次单位根,\(k=3\)。
- 设 Leontief 系数矩阵每列和都 \(<1\)。证明 \((I-A)^{-1}\ge0\),并说明列和的经济含义。 提示:列和界 + Neumann 级数;列和是每 1 元产出的中间投入成本,小于 1 意味着有增加值。
- 证明:若 \(A\ge0\) 不可约且有一个正对角元,则 \(\rho(A)\) 是唯一的最大模特征值。 提示:Corollary 8.4.7。
进阶
- 证明谱半径的灵敏度公式 \(\partial\rho/\partial a_{ij}=y_ix_j\)(\(A\) 不可约,\(y^Tx=1\)),并说明为什么它总是正的。 提示:对 \(A(t)x(t)=\rho(t)x(t)\) 求导,左乘 \(y^T\)。
- 用 \(A(\epsilon)=A+\epsilon J\) 与紧性证明 Theorem 8.3.1,并指出为什么这个方法不能保住"单重"和"正特征向量"。 提示:极限过程中重根可以合并,正分量可以趋于零。考虑 \(A=I\)。
- 设 \(Q\) 是连续时间 Markov 链的生成元(非对角元 \(\ge0\),行和为 0)。证明 \(Q\) 是 Metzler 矩阵、主导特征值为 0,且 \(e^{tQ}\ge0\)。 提示:\(Q+\lambda I\ge0\),\(e^{tQ}=e^{-\lambda t}e^{t(Q+\lambda I)}\)。
- (非负最佳秩一逼近)设 \(A\ge0\) 且 \(AA^T\) 不可约。证明 \(A\) 的首个左、右奇异向量可以取成正向量,并解释在"股票 × 日期"的成交额矩阵中这个秩一分量的含义。 提示:对 \(AA^T\)、\(A^TA\) 用 Perron–Frobenius;它是"股票规模 × 市场整体活跃度"的乘积结构。
- 证明 8.3.P12 的修正版本:\(A\ge0\)、\(r>\rho(A)\) 时 \((rI-A)^{-1}\ge0\);\(A\) 不可约时 \((rI-A)^{-1}>0\)。给出可约时有零元的例子。 提示:Neumann 级数;不可约时 \(\sum_kA^k/r^k\ge c(I+A)^{n-1}>0\)。
原书推荐习题
- 8.0.P2:显式计算 \(A_\epsilon^m\) 的极限与 \(xy^T\);8.1.P7:用 \(x=Ae\) 改进行和界(幂迭代一步)。
- 8.2.P5、P8、P14:严格单调、未归一化特征向量的极限、由伴随矩阵求 Perron 向量;8.2.P11:单重性的特征多项式证明。
- 8.3.P8:Frobenius 标准形;8.3.P9:Metzler 矩阵;8.3.P15–P16:M 矩阵与单调矩阵。
- 8.4.P13:Perron 根对元素的导数;8.4.P17:非负矩阵最佳秩一逼近;8.4.P25:\(\limsup(\operatorname{tr}A^m)^{1/m}=\rho(A)\)。
原书对照
| 本章小节 | 原书小节 | 书页 | PDF 页 |
|---|---|---|---|
| 8.1 引言:人口迁移 | 8.0 Introduction(含习题 8.0) | 517–519 | 537–539 |
| 8.2 逐元素序与基本不等式 | 8.1 Inequalities and generalities(含习题 8.1.P1–P11) | 519–524 | 539–544 |
| 8.3 Perron 定理 | 8.2 Positive matrices(含习题 8.2.P1–P16) | 524–529 | 544–549 |
| 8.4 一般非负矩阵(含 8.3.P12(b) 说明) | 8.3 Nonnegative matrices(含习题) | 529–533 | 549–553 |
| 8.5 Perron–Frobenius 与循环结构 | 8.4 Irreducible nonnegative matrices(含习题 8.4.P1–P26) | 533–540 | 553–560 |
下一章(第 08b 章)在不可约之上再加"本原"(非周期)条件,恢复 \((A/\rho)^m\) 的收敛,并系统讨论 Markov 链转移矩阵、PageRank 与双随机矩阵。