线性代数（十四）：随机矩阵理论——混沌中的秩序

扔一百万次硬币，把结果排成一个 $1000 \times 1000$ 的对称矩阵，算出来的特征值居然会填满一个完美的半圆。一个理论上应该是单位阵的噪声协方差矩阵，特征值却会分布在一段区间上，而这段区间的宽度我甚至在看到任何数据之前就能预测出来。Wigner 矩阵的最大特征值的尾部分布无处不在——从晶体生长的高度涨落、随机排列中最长递增子序列的长度，到重核能级的能量分布，全都符合这个规律。随机矩阵理论（Random Matrix Theory，RMT）研究的就是这些规律为何会出现，以及如何加以利用。

本章你将学到
Wigner 矩阵和 Wishart 矩阵：驱动几乎所有现象的两大模型家族
半圆律和 Marchenko-Pastur 律：谱的“中心极限定理”
特征值排斥、Wigner 猜测、Tracy-Widom 边缘分布
自由概率论和 Stieltjes 变换的核心直觉
应用：MIMO 容量、协方差清洗、PCA 阈值选择、尖峰协方差与 BBP 相变
预备知识： 特征分解（第六章）、SVD（第九章）、基本概率（均值、方差、独立同分布、经典 CLT）。

第一个惊喜：随机矩阵为何“不随机”#

取一个 $1000 \times 1000$ 的矩阵，用独立的标准高斯变量填充所有元素，再将其对称化并除以 $\sqrt{n}$ 。计算其特征值并绘制直方图，你会发现无论使用什么随机种子，结果总是呈现出 $$[-2, 2]$$ 区间上的同一个半圆。换成均匀分布的 $\{-1, +1\}$ 元素？还是那个半圆。换成复埃尔米特高斯矩阵？半圆依然不变。

这并非魔法，而是经典中心极限定理在高维空间中的自然延伸。当你对上百万个弱相关的随机变量进行某种“平均”时（而一个特征值本质上就是矩阵元素的一种复杂平均），微观分布会被彻底抹平，最终极限仅由少数粗粒度统计量——均值、方差和对称性——决定。半圆律正是谱的中心极限定理。

在矩阵的世界里，变化的只是极限对象的形式：不再是实轴上的单个高斯分布，而是一整条特征值密度曲线，外加关于间距的精确描述和边缘处的涨落规律。RMT 正是研究这三者的数学工具。

两大模型家族#

Wigner 矩阵：对称噪声#

一个 Wigner 矩阵 $\mathbf{W} \in \mathbb{R}^{n\times n}$ 满足：

对角元 $w_{ii}$ 独立同分布，均值为 0，方差为 $\sigma_d^2$ ；
上三角元素 $w_{ij}$ （ $$i < j$$ ）独立同分布，均值为 0，方差为 $\sigma^2$ ；
下三角元素由对称性强制设定为 $w_{ji} = w_{ij}$ 。

当所有元素均为高斯分布时，就得到 高斯正交系综（GOE）。其名称中的“正交”源于分布对任意正交变换 $\mathbf{O}\mathbf{W}\mathbf{O}^\top$ 的不变性。对于复埃尔米特矩阵，对应的是 GUE（酉不变性）；对于四元数矩阵，则是 GSE（辛不变性）。GOE / GUE / GSE 这三者分别对应时间反演对称性存在、被磁场破坏、或因半整数自旋而破坏的情形——这也解释了为何这些矩阵在物理学中反复出现。

直观理解： 想象一个 Erdős–Rényi 社交网络，其中边权重 $w_{ij}$ 表示随机生成的“亲密度”分数。该亲密度矩阵的特征值揭示了其全局社群结构；而当亲密度完全随机时，这种结构恰恰是我们即将研究的普适规律。

Wishart 矩阵：样本协方差#

\mathbf{S} \;=\; \frac{1}{n}\,\mathbf{X}^\top \mathbf{X} \;\in\; \mathbb{R}^{p\times p}.

当 $n \gg p$ 时， $\mathbf{S}$ 是真实协方差的良好估计；但若 $$n$$ 和 $$p$$ 同时很大且比值 $\gamma = p/n$ 固定， $\mathbf{S}$ 就会严重偏离真实值——不过这种偏离是有结构、可预测的，由 Marchenko–Pastur 定律精确刻画。

直观理解： 假设你跟踪 $$p = 500$$ 只股票一年的每日收益率（约 $n \approx 252$ 个交易日）。你计算出的“协方差矩阵”包含约 125,000 个自由参数，却仅有约 126,000 个数据点。因此，即使真实协方差是单位阵，所得估计本质上仍是一个随机矩阵，其特征值会按照 RMT 规律扩散开来。

Wigner 半圆律#

定理陈述#

\hat{\mathbf{W}} \;=\; \frac{\mathbf{W}}{\sigma\sqrt{n}}.

f(x) \;=\; \frac{1}{2\pi}\sqrt{4 - x^2}, \qquad x \in [-2, 2],

而在 $$[-2, 2]$$ 之外， $$f(x) = 0$$ 。

为什么是半圆？三种视角#

1) 矩方法（严格的数学推导）。 计算 $m_k = \mathbb{E}[\frac{1}{n}\operatorname{tr}\hat{\mathbf{W}}^k]$ 。迹展开后对应 $$n$$ 个顶点上长度为 $$k$$ 的闭合路径，每条路径的权重由相应高斯矩乘积决定。由于独立性和零均值，只有两两配对且非交叉的路径能存活下来——这正是 Catalan 数 $C_{k/2}$ 所计数的结构。而 Catalan 数恰好是半圆分布的偶阶矩，因此极限分布必然是半圆。

\rho(\lambda_1, \dots, \lambda_n) \;\propto\; \prod_{i<j} |\lambda_i - \lambda_j|\;\exp\!\Big(-\tfrac{n}{4}\sum_i \lambda_i^2\Big),

它描述了 $$n$$ 个带电粒子在直线上受对数排斥力（Vandermonde 因子）和谐振子势约束的平衡态。排斥与约束之间的平衡结果，正是半圆分布。

3) 自由 CLT（代数视角）。 对称随机矩阵可视为大量“自由”的秩一扰动之和。在自由概率中，“独立变量求和”的对应概念是自由加法卷积，其对应的中心极限定理给出的极限分布是半圆而非高斯。第 7 节将展开讨论这一点。

为什么底层分布无关紧要？#

半圆律具有普适性：无论矩阵元素服从高斯、均匀、 $\pm 1$ ，还是任何均值为零且方差有限的分布，极限分布始终是同一个半圆。这与经典中心极限定理如出一辙——归一化后的和总收敛到高斯分布，与原始分布无关。

背后的机制也完全一致：在极限过程中，只有二阶矩得以保留，其余细节均被平均效应抹去。

代码验证#

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
import numpy as np
import matplotlib.pyplot as plt

n, repeats = 1500, 30
all_eigs = []
for _ in range(repeats):
    a = np.random.randn(n, n)
    w = (a + a.T) / np.sqrt(2 * n)        # 方差归一化：标准半圆
    all_eigs.append(np.linalg.eigvalsh(w))
all_eigs = np.concatenate(all_eigs)

x = np.linspace(-2, 2, 600)
plt.hist(all_eigs, bins=90, density=True, alpha=0.55, label="Empirical")
plt.plot(x, np.sqrt(np.maximum(4 - x**2, 0)) / (2 * np.pi),
         lw=2.5, label="Semicircle")
plt.legend(); plt.show()

实际上，当 $$n = 200$$ 时，模拟结果已与理论半圆几乎完美吻合。

Marchenko-Pastur 定律#

定义#

f(\lambda) \;=\; \frac{1}{2\pi\gamma\,\lambda}\sqrt{(\lambda_+ - \lambda)(\lambda - \lambda_-)},\qquad \lambda \in [\lambda_-, \lambda_+],

\lambda_\pm \;=\; (1 \pm \sqrt{\gamma})^2.

若 $\gamma > 1$ ，矩阵秩为 $$n < p$$ ，除区间 $[\lambda_-, \lambda_+]$ 内的特征值外，还有 $$p - n$$ 个精确为零的特征值。

这条密度的意义#

即使总体协方差是单位阵，有限样本也会导致特征值扩散：

$\gamma = 0.1$ ：边界 $$[0.47, 1.69]$$ ——轻微展宽；
$\gamma = 0.5$ ：边界 $$[0.09, 2.91]$$ ——最大样本特征值达到理论值的三倍；
$\gamma = 1.0$ ：边界 $$[0, 4]$$ ——谱触及零点，矩阵濒临奇异。

这一密度构成了 RMT 在统计学中所有定量应用的基础：任何严格大于 $\lambda_+$ 的样本特征值都显著偏离噪声；而落在 $[\lambda_-, \lambda_+]$ 内的特征值则与纯噪声无法区分。

不同纵横比下的 Marchenko-Pastur 密度：随着 γ 增大，谱范围变宽；当 γ=1 时，谱范围触及零

代码#

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
import numpy as np, matplotlib.pyplot as plt

n, p, repeats = 1000, 500, 30
gamma = p / n
all_eigs = []
for _ in range(repeats):
    x = np.random.randn(n, p)
    all_eigs.append(np.linalg.eigvalsh(x.T @ x / n))
all_eigs = np.concatenate(all_eigs)

lam_minus, lam_plus = (1 - np.sqrt(gamma))**2, (1 + np.sqrt(gamma))**2
xs = np.linspace(lam_minus + 1e-3, lam_plus - 1e-3, 500)
mp = np.sqrt((lam_plus - xs) * (xs - lam_minus)) / (2 * np.pi * gamma * xs)

plt.hist(all_eigs, bins=100, density=True, alpha=0.5, label="Empirical")
plt.plot(xs, mp, "r", lw=2, label="MP theory")
for e in (lam_minus, lam_plus):
    plt.axvline(e, ls="--", color="g")
plt.legend(); plt.show()

精细结构：排斥与边缘#

半圆律和 MP 律描述的是整体——宏观形状。RMT 还在两个更精细的尺度上给出了精确刻画。

间距分布：特征值互相排斥#

p(s) \;=\; \frac{\pi s}{2}\,\exp\!\Big(-\frac{\pi s^2}{4}\Big).

关键在于 $$p(0) = 0$$ ：特征值永不重合。它们主动排斥，且在 $$s=0$$ 处以二次方式趋零（GUE 则为三次——不同对称类对应不同排斥指数）。

相比之下，独立能级的间距服从指数分布 $p(s) = e^{-s}$ ，在 $$s=0$$ 处取得最大值。独立性所预期的聚集现象，恰恰是特征值坚决避免的。

边缘：Tracy-Widom 分布#

n^{2/3}\big(\lambda_{\max} - 2\big) \;\xrightarrow{d}\; \mathrm{TW}_1,

其中 $\mathrm{TW}_1$ 是 $\beta=1$ 的 Tracy-Widom 分布。该分布高度不对称：左尾超指数衰减（几乎不会看到 $\lambda_{\max}$ 明显小于 2），右尾按 $\exp(-\frac{2}{3} t^{3/2})$ 衰减（极端大值虽罕见但确实存在）。同一 Tracy-Widom 分布还出现在随机置换的最长递增子序列、晶体生长高度涨落，以及现代统计中大随机矩阵的最大奇异值问题中。

随机谱 vs 确定谱#

一张图胜过千言：随机 Wigner 谱呈现为平滑且强相关的整体；而确定性谱则是若干尖峰簇，每个簇内部仅有独立采样噪声。

应用#

无线通信：MIMO 容量随天线数量线性增长#

C \;=\; \sum_i \log_2\!\Big(1 + \frac{\mathrm{SNR}}{n_t}\,\lambda_i\Big),

其中 $\lambda_i$ 是 $\mathbf{H}\mathbf{H}^\dagger$ 的特征值， $\mathbf{H}$ 为随机信道矩阵。MP 律给出了 $\lambda_i$ 的极限密度，对其与 $\log_2(1 + \mathrm{SNR}\cdot\lambda)$ 的乘积积分即可得遍历容量。结论明确：容量与 $\min(n_t, n_r)$ 成线性关系。这正是现代手机、基站和 Wi-Fi 路由器普遍采用多天线设计的理论根基。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
import numpy as np

def mimo_capacity(n_r, n_t, snr_db, trials=500):
    snr = 10 ** (snr_db / 10)
    caps = []
    for _ in range(trials):
        H = (np.random.randn(n_r, n_t)
             + 1j * np.random.randn(n_r, n_t)) / np.sqrt(2)
        eigs = np.linalg.eigvalsh(H @ H.conj().T)
        caps.append(np.sum(np.log2(1 + snr / n_t * eigs)))
    return float(np.mean(caps))

print(f"4x4 @ 10 dB: {mimo_capacity(4, 4, 10):.2f} bits/s/Hz")
print(f"8x8 @ 10 dB: {mimo_capacity(8, 8, 10):.2f} bits/s/Hz")

金融：清洗样本协方差矩阵#

跟踪 $$p$$ 只股票在 $$n$$ 个交易日的数据，现实中纵横比 $\gamma = p/n$ 往往不小。MP 边界能精确识别哪些样本特征值属于噪声。标准流程如下：

对 $\mathbf{S}$ 做特征分解： $\mathbf{S} = \mathbf{U}\,\mathrm{diag}(\lambda_1, \dots, \lambda_p)\,\mathbf{U}^\top$ ；
估计噪声方差 $\sigma^2$ （例如取 MP 支撑区间内特征值的均值）；
计算边界 $\lambda_\pm = \sigma^2(1 \pm \sqrt{\gamma})^2$ ；
将 $[\lambda_-, \lambda_+]$ 内的所有特征值替换为其均值（或采用更精细的收缩方法）；
重构清洗后的协方差矩阵 $\tilde{\mathbf{S}} = \mathbf{U}\,\mathrm{diag}(\tilde\lambda_i)\,\mathbf{U}^\top$ 。

基于 $\tilde{\mathbf{S}}$ 构建的投资组合，样本外夏普比率通常提升 10%–30%，同时显著降低换手率。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
import numpy as np

def clean_covariance(returns):
    """把 MP 主体内的特征值替换为均值。"""
    n, p = returns.shape
    gamma = p / n
    S = np.cov(returns, rowvar=False)
    eigs, vecs = np.linalg.eigh(S)
    sigma2 = np.mean(eigs)
    lo, hi = sigma2 * (1 - np.sqrt(gamma))**2, sigma2 * (1 + np.sqrt(gamma))**2
    bulk = (eigs >= lo) & (eigs <= hi)
    eigs[bulk] = eigs[bulk].mean()
    return vecs @ np.diag(eigs) @ vecs.T

高维 PCA：保留多少主成分？#

在 PCA 流程中，当 $p \approx n$ 时，传统方法难以判断“多少主成分是信号”。RMT 给出了清晰答案：仅保留超过 $\lambda_+ = (1 + \sqrt{\gamma})^2$ 的特征值。低于 MP 边界的特征值在统计上无法与 i.i.d. 噪声区分。

1
2
3
4
5
6
7
8
import numpy as np

def pca_signal_count(X):
    n, p = X.shape
    gamma = p / n
    Xs = (X - X.mean(0)) / X.std(0)
    eigs = np.linalg.eigvalsh(Xs.T @ Xs / n)
    return int(np.sum(eigs > (1 + np.sqrt(gamma))**2))

尖峰协方差与 BBP 相变#

上述 PCA 准则在信号较弱时过于粗糙。更贴近现实的模型是尖峰协方差：总体协方差为 $\Sigma = \mathbf{I} + \sum_{k=1}^{r} (s_k - 1) \mathbf{v}_k\mathbf{v}_k^\top$ ，即单位阵叠加 $$r$$ 个“尖峰”。Baik-Ben Arous-Péché（BBP）相变指出：

若总体尖峰 $s_k > 1 + \sqrt{\gamma}$ ，对应样本特征值会脱离 MP 主体，落在 $s_k + \gamma s_k/(s_k - 1)$ ，且样本特征向量与 $\mathbf{v}_k$ 的夹角余弦为正；
若 $s_k \le 1 + \sqrt{\gamma}$ ，尖峰将完全不可见：它被埋入 MP 主体中，无论算法多么精巧，PCA 都无法恢复。

这是高维 PCA 的一条硬性不可能定理。它明确告诉你何时需要更多样本，何时问题在当前 $$n, p$$ 下统计上无解。

尖峰协方差：主体服从 MP 律，超临界尖峰从总体值偏移到外部，亚临界尖峰融进主体

神经网络初始化#

初始化权重时，若使 $\mathbf{W}^\top\mathbf{W}$ 的特征值集中在 1 附近，可逐层稳定激活值的方差。Xavier/Glorot 初始化规则 $\mathrm{Var}(w_{ij}) = 2 / (n_\text{in} + n_\text{out})$ 恰好使 $\mathbf{W}$ 的奇异值落在以 1 为中心的 MP 支撑范围内。正交初始化更进一步——它完全消除奇异值展宽，在深线性网络中实现更优的信号传播。

工具箱#

Stieltjes 变换#

m_\mu(z) \;=\; \int \frac{d\mu(\lambda)}{\lambda - z}, \qquad z \in \mathbb{C}^+.

f(\lambda) \;=\; -\frac{1}{\pi}\lim_{\eta \to 0^+}\operatorname{Im}\,m_\mu(\lambda + i\eta).

为何使用它？因为经验谱分布的 Stieltjes 变换恰等于 $\frac{1}{n}\operatorname{tr}(\mathbf{M} - z\mathbf{I})^{-1}$ ，而预解算子（resolvent）在代数上易于操作。现代 RMT 的大多数证明都在 Stieltjes 变换层面完成：写出 $$m(z)$$ 的自洽方程，求解，再反演。

对半圆律，方程为 $$m(z)^2 + z m(z) + 1 = 0$$ ，解得 $m(z) = (-z + \sqrt{z^2 - 4})/2$ 。Marchenko-Pastur 的方程同样简洁。

自由概率论一句话讲清楚#

在自由概率中，“非交换随机变量” $$a, b$$ 是自由独立的，若其交替混合矩满足特定的迹分解规则。Voiculescu 定理指出：大尺寸独立随机矩阵在渐近意义下是自由的。卷积的自由版本——自由加法卷积 $\boxplus$ ——将 $\mathbf{A}$ 和 $\mathbf{B}$ 的谱分布映射到 $\mathbf{A} + \mathbf{B}$ 的谱分布。其对应的中心极限定理给出半圆分布。实用启示是：对标量独立变量的所有操作（均值、方差、CLT、Berry-Esseen），在矩阵谱层面都有自由概率的对应物。

半圆律的证明骨架（矩方法）#

计算 $m_k = \mathbb{E}[\frac{1}{n}\operatorname{tr}\hat{\mathbf{W}}^k]$ 。展开迹后得到对 $\{1, \dots, n\}$ 上长度为 $$k$$ 的闭路径的求和，权重为 $\prod_e \mathbb{E}[w_e]$ 。独立性与零均值会消去所有含奇数次某条边的路径。当 $$k$$ 为偶数时，主导贡献来自两两配对且非交叉的路径——其数量恰为 Catalan 数 $C_{k/2}$ ，归一化后权重为 1。而 Catalan 数正是半圆分布的偶阶矩，故极限分布必为半圆。

练习题#

基础题#

写出 $3 \times 3$ GOE 矩阵的 $\mathbb{E}[\mathbf{W}]$ 和 $\mathrm{Cov}(w_{ij}, w_{kl})$ ，其中对角元方差为 2，非对角元方差为 1。
为何需要 $1/\sqrt{n}$ 归一化？若不归一化， $\lambda_{\max}$ 会如何变化？
设 $\gamma = 0.5$ ，计算 $\lambda_\pm$ 并草绘 MP 密度。其峰值位于何处？

计算题#

用代换 $x = 2\sin\theta$ 证明 $\int_{-2}^{2}\frac{1}{2\pi}\sqrt{4 - x^2}\,dx = 1$ 。
计算半圆分布的二阶和四阶矩，验证其是否等于 Catalan 数 $$C_1$$ 和 $$C_2$$ 。
设 $\mathbf{X}$ 元素独立同分布于 $\mathcal{N}(0, 1)$ ，证明 $\mathbb{E}[\mathbf{X}^\top\mathbf{X}] = n\mathbf{I}_p$ 。
推导 $2 \times 2$ 对称高斯矩阵两特征值的联合密度，并指出排斥因子。

编程题#

验证 Wigner 猜测：生成 200 个 $$n = 500$$ 的 GOE 矩阵，计算主体区域的最近邻间距，用局部均值归一化后，与 $$p(s)$$ 和 $e^{-s}$ 对比。
绘制 $\gamma \in \{0.1, 0.5, 1.0, 2.0\}$ 的 MP 直方图；当 $\gamma > 1$ 时，单独统计零特征值数量。
绘制 $2\times2$ 、 $4\times4$ 、 $8\times8$ 、 $16\times16$ MIMO 容量随 SNR 的变化曲线，验证高 SNR 下“天线翻倍，容量翻倍”的规律。
复现 BBP 相变：固定 $\gamma = 0.5$ ，让尖峰强度 $$s$$ 从 1 扫至 3，绘制最大样本特征值随 $$s$$ 的变化曲线，并标出临界点 $s = 1 + \sqrt{\gamma}$ 。

应用题#

一名量化研究员跟踪 100 只股票 200 天。计算 $\gamma$ ，求 $[\lambda_-, \lambda_+]$ ，并判断样本特征值 3.5 是信号还是噪声。
对 $8 \times 4$ 复高斯 MIMO 信道，写出容量表达式，在 30 dB 下估计容量，并预测升级至 $16 \times 8$ 的增益。
给定 1000 个样本和 500 个特征，计算 MP 阈值，并说明如何用于选择 PCA 主成分数量。

进阶题#

查阅 Tracy-Widom 分布及其在高维 MANOVA 中 Roy 最大根检验的作用。为何当 $$p/n$$ 不小时，经典卡方近似失效？
详述 BBP 相变现象。为何没有任何估计器能检测亚临界尖峰？
为何复杂原子核的能级间距服从 GOE 统计？时间反演对称性在 GOE 与 GUE 的选择中起何作用？

总结#

模型： Wigner 矩阵刻画对称噪声，Wishart 矩阵刻画样本协方差。其高斯版本 GOE、GUE、GSE 对应不同对称类。

主体规律： Wigner 谱服从半圆律，Wishart 谱服从 Marchenko-Pastur 分布。二者均具普适性：只要元素方差有限，具体分布无关紧要。

精细结构： 特征值间存在排斥效应（Wigner 猜测）；最大特征值涨落尺度为 $n^{-2/3}$ ，服从 Tracy-Widom 分布。

工具： Stieltjes 变换将谱问题转为代数方程；自由概率将经典 CLT 式推理推广至独立矩阵的谱分析。

为何重要： RMT 为 MIMO 容量、协方差清洗、PCA 阈值、尖峰模型 BBP 相变及神经网络初始化提供了精确且常无需参数的预测。核心启示始终如一：高维随机性背后隐藏着确定性的骨架，而我们可以利用它。

参考文献#

Bai, Z., & Silverstein, J. W. Spectral Analysis of Large Dimensional Random Matrices. Springer, 2010.
Anderson, G. W., Guionnet, A., & Zeitouni, O. An Introduction to Random Matrices. Cambridge University Press, 2010.
Mehta, M. L. Random Matrices. Academic Press, 2004.
Tao, T. Topics in Random Matrix Theory. AMS, 2012.
Tulino, A. M., & Verdú, S. Random Matrix Theory and Wireless Communications. Foundations and Trends, 2004.
Bouchaud, J.-P., & Potters, M. Financial Applications of Random Matrix Theory. arXiv:0910.1205 , 2009.
Couillet, R., & Debbah, M. Random Matrix Methods for Wireless Communications. Cambridge University Press, 2011.
Baik, J., Ben Arous, G., & Péché, S. “Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices.” Annals of Probability, 2005.