随机矩阵理论中Marchenko-Pastur分布与样本协方差矩阵谱

拿到一个高维数据矩阵，想判断它的特征值里哪些是真实信号、哪些只是噪声，最直接的工具就是Marchenko-Pastur分布（简称MP分布）。下面按可执行顺序拆成七个阶段，每一步都能用纸笔或几行代码完成。

阶段一：把原始数据变成样本协方差矩阵

1. 整理 数据矩阵 $X$，让每一行是一个样本、每一列是一个特征，矩阵尺寸记为 $n \times p$，其中 $n$ 是样本数，$p$ 是特征数。

2. 中心化 每一列：用该列所有元素减去该列均值，使每一列的均值为零。这一步不能省，否则最大的特征值会被均值方向污染。

3. 构造 样本协方差矩阵：

$$ S = \frac{1}{n} X^{T} X $$

这里分母用 $n$ 而不是 $n-1$，这是随机矩阵理论的标准约定，和统计软件默认的 $n-1$ 只差一个常数倍，不影响谱的形状判断。

4. 计算 矩阵 $S$ 的全部特征值，得到一个长度为 $p$ 的实数列表。由于 $S$ 是对称半正定矩阵，特征值全部非负。

5. 排列 特征值从小到大，记为 $\lambda_1 \le \lambda_2 \le \cdots \le \lambda_p$。这个列表就是后面要反复对照的“经验谱”。

阶段二：明确“谱”到底在描述什么

谱就是特征值的集合。把每个特征值看成一根竖线，统计落在每个小区间里的特征值个数，再除以总数 $p$，就得到经验谱分布。它回答的问题是：随机抽一个特征值，它落在某个区间里的概率有多大。

样本协方差矩阵的特征值会随数据一起随机波动。随机矩阵理论关心的不是某一次抽样的具体值，而是当 $n$ 和 $p$ 同时变得很大、并且比值趋于固定时，这堆特征值的整体分布会稳定到一个确定的形状。这个极限形状就是MP分布。

阶段三：套用Marchenko-Pastur公式

计算 长宽比：

$$ c = \frac{p}{n} $$

这个比值是MP分布唯一的形状参数。它衡量“特征数相对样本数有多多”。

写出 MP分布的密度函数。假设数据每个元素独立同分布、均值为零、方差为 $\sigma^2$，则当 $n, p \to \infty$ 且 $p/n \to c$ 时，经验谱分布收敛到：

$$ f_{MP}(x) = \frac{1}{2\pi c \sigma^2 x}\sqrt{(b-x)(x-a)}, \quad a \le x \le b $$

其中上下边界为：

$$ a = \sigma^2 \left(1 - \sqrt{c}\right)^2, \qquad b = \sigma^2 \left(1 + \sqrt{c}\right)^2 $$

注意 当 $c > 1$ 时，也就是特征数比样本数还多，矩阵 $S$ 不满秩，会在 $x=0$ 处额外堆积质量 $1 - 1/c$ 的特征值。这时密度公式只在 $[a, b]$ 上有效，零点单独处理。

记住 三个最常用的结论：

1. 谱的支撑区间是 $[a, b]$，区间外理论上不该有特征值。
2. 谱的均值等于 $\sigma^2$，与 $c$ 无关。
3. 谱的方差等于 $\sigma^4 c$，比值越大，谱越“胖”。

阶段四：用三个数快速核对

不必画图也能判断数据是否符合MP律。计算 下面三个量，与理论值逐一对照。

1. 核对 最大值。经验谱的最大特征值应接近 $b = \sigma^2(1+\sqrt{c})^2$。如果明显超出，说明存在强信号方向或数据不独立。

2. 核对 最小值。$c < 1$ 时经验谱的最小特征值应接近 $a = \sigma^2(1-\sqrt{c})^2$；$c \ge 1$ 时应接近 0。

3. 核对 平均值。把所有特征值取平均，结果应接近 $\sigma^2$。若偏离很大，先检查是否漏了中心化，或者数据方差本身不是常数。

估计 $\sigma^2$ 的实用做法：如果确信数据以噪声为主，直接用特征值均值作为 $\sigma^2$ 的估计，再代入边界公式。

阶段五：用代码做一次完整验证

运行 下面这段脚本，直接对比经验谱与理论边界：

import numpy as np

rng = np.random.default_rng(42)
n, p = 2000, 500
sigma2 = 1.0
c = p / n

X = rng.normal(0, np.sqrt(sigma2), size=(n, p))
X = X - X.mean(axis=0)
S = X.T @ X / n

eigvals = np.linalg.eigvalsh(S)
a = sigma2 * (1 - np.sqrt(c)) ** 2
b = sigma2 * (1 + np.sqrt(c)) ** 2

print(f"c = {c:.4f}")
print(f"理论边界: [{a:.4f}, {b:.4f}]")
print(f"经验最小: {eigvals.min():.4f}")
print(f"经验最大: {eigvals.max():.4f}")
print(f"经验均值: {eigvals.mean():.4f}")
print(f"理论均值: {sigma2:.4f}")

观察 输出中的经验最小、最大值应与理论边界相差很小；经验均值应接近理论均值。若偏差随 $n$ 增大而缩小，说明符合MP律。

调整 把 p 改成 n 的 1.2 倍再跑一次，会看到 $c > 1$，零点附近出现大量接近零的特征值，验证不满秩时的额外质量。

阶段六：把MP分布用于信号筛选

这是MP分布最实用的落地场景：从噪声里挑出真实方向。

1. 设 零假设为“数据全是噪声”，此时特征值应全部落在 $[a, b]$ 内。

2. 标记 所有大于 $b$ 的特征值。它们超出噪声谱的上边界，对应真实信号方向，保留。

3. 标记 所有小于 $a$ 的特征值（$c<1$ 时）。它们同样偏离噪声谱，通常提示数据存在相关性结构或方差非均匀，需要单独排查。

4. 丢弃 落在 $[a, b]$ 内的特征值对应的方向。在零假设下它们与随机噪声不可区分，用来做降维或建模只会引入虚假结构。

5. 重算 剔除信号方向后，用剩余特征值重新估计 $\sigma^2$，再更新边界 $a, b$，迭代一两轮可让阈值更稳。

这个方法在PCA降维、金融协方差矩阵估计、无线信号检测里都直接可用，替代了凭经验画“碎石图”的粗糙做法。

阶段七：避开四个常见误用

1. 忘记中心化。均值不为零会让第一个特征值远大于 $b$，被误判为强信号。先减列均值再做一切。

2. 样本量过小。$n$ 只有几十时，经验谱本身波动极大，边界附近判断不可靠。经验上 $n$ 至少几百，边界才有参考价值。

3. 数据重尾。MP分布假设元素独立同分布且四阶矩有限。金融收益率这类重尾数据，谱的上边界会被少数极端值抬高，直接套公式会高估 $b$。此时要么对数据做截尾，要么改用更稳健的谱估计。

4. 把 $c$ 当成固定值。$c = p/n$ 随样本量变化，阈值 $b$ 也跟着变。每次换数据集都要重新算 $c$ 和边界，不能沿用上一次的数值。

完成 上述七步后，手里就有了一个不依赖主观判断的谱诊断流程：从矩阵构造、参数计算、边界核对，到信号筛选和误用排查，全部可由数值计算闭环验证。