用Python和NumPy理解奇异值分解：从Rayleigh商到Courant-Fischer定理的直观推导

在数据科学和机器学习领域，奇异值分解（SVD）是一个强大的工具，广泛应用于降维、推荐系统、自然语言处理等场景。然而，许多从业者仅停留在调用 sklearn.decomposition.TruncatedSVD 的层面，对其背后的数学原理知之甚少。本文将带你从Rayleigh商出发，通过Python和NumPy的实践，直观理解Courant-Fischer定理如何揭示SVD的本质。

1. Rayleigh商与特征值的极值性质

Rayleigh商是理解矩阵特征值和奇异值的重要桥梁。对于一个对称矩阵 M 和任意非零向量 x ，Rayleigh商定义为：

def rayleigh_quotient(M, x):
    return (x.T @ M @ x) / (x.T @ x)

这个简单的比值揭示了矩阵 M 在向量 x 方向上的"拉伸"程度。特别地，当 x 是 M 的特征向量时，Rayleigh商就等于对应的特征值。

我们可以用NumPy验证这一性质：

import numpy as np

# 生成一个对称矩阵
M = np.random.randn(5, 5)
M = M + M.T  # 确保对称

# 计算特征值和特征向量
eigvals, eigvecs = np.linalg.eigh(M)

# 验证Rayleigh商
for i in range(5):
    x = eigvecs[:, i]
    print(f"特征值: {eigvals[i]:.4f}, Rayleigh商: {rayleigh_quotient(M, x):.4f}")

输出将显示每个特征向量对应的Rayleigh商确实等于其特征值。

Rayleigh商的一个重要性质是极值性：对于对称矩阵，最大和最小特征值分别对应Rayleigh商的全局最大值和最小值。这一性质可以推广到中间特征值，这正是Courant-Fischer定理的核心内容。

2. Courant-Fischer定理的直观理解

Courant-Fischer定理以两种形式描述了特征值的极值特性：

2.1 子空间形式

对于n×n对称矩阵 M ，其特征值按降序排列λ₁ ≥ λ₂ ≥ ... ≥ λₙ，第k个特征值满足：

λₖ = max_{dim(S)=k} min_{x∈S} R(M,x) 
   = min_{dim(T)=n-k+1} max_{x∈T} R(M,x)

其中R(M,x)表示Rayleigh商，S和T是ℝⁿ的子空间。

2.2 正交补形式

另一种等价表述使用正交补：

λₖ = max_{w₁,...,w_{n-k}} min_{x⊥w₁,...,w_{n-k}} R(M,x)
   = min_{w₁,...,w_{k-1}} max_{x⊥w₁,...,w_{k-1}} R(M,x)

这两种形式虽然抽象，但揭示了特征值作为Rayleigh商在特定子空间约束下的极值性质。

3. 用NumPy验证Courant-Fischer定理

让我们通过代码验证Courant-Fischer定理的第一个等式。我们将构造不同的k维子空间S，计算min R(M,x)，然后寻找这些最小值中的最大值。

def verify_courant_fischer(M, k, trials=1000):
    n = M.shape[0]
    eigvals = np.linalg.eigvalsh(M)
    lambda_k = eigvals[-k]  # 因为eigvalsh返回升序
    
    max_min = -np.inf
    
    for _ in range(trials):
        # 随机生成k维子空间
        S = np.random.randn(n, k)
        S, _ = np.linalg.qr(S)  # 正交化
        
        # 计算子空间中的最小Rayleigh商
        min_rq = np.inf
        for __ in range(100):  # 采样子空间中的向量
            x = S @ np.random.randn(k)
            x = x / np.linalg.norm(x)
            rq = rayleigh_quotient(M, x)
            if rq < min_rq:
                min_rq = rq
        
        if min_rq > max_min:
            max_min = min_rq
    
    print(f"理论λ_{k}: {lambda_k:.6f}, 实验max-min: {max_min:.6f}")
    return abs(max_min - lambda_k) < 1e-4

运行这个函数，我们会发现实验得到的max-min值非常接近理论特征值λₖ，从而验证了Courant-Fischer定理。

4. 从Courant-Fischer到奇异值分解

奇异值分解(SVD)可以看作是对任意矩阵(不一定是方阵或对称)的特征值分解的推广。对于矩阵 A ，其奇异值σₖ与 AᵀA 的特征值λₖ满足σₖ=√λₖ。

Courant-Fischer定理可以推广到奇异值：

σₖ = max_{dim(S)=k} min_{x∈S, ||x||=1} ||Ax||
   = min_{dim(T)=n-k+1} max_{x∈T, ||x||=1} ||Ax||

这解释了为什么在PCA中，最大的奇异值对应的方向能保留数据最大的方差——因为它在所有k维子空间中最大化最小拉伸。

5. 应用实例：理解PCA中的方差最大化

主成分分析(PCA)常被描述为寻找数据最大方差方向的方法。通过Courant-Fischer定理，我们可以更深刻地理解这一点。

假设数据矩阵 X 已中心化，协方差矩阵为 C=XᵀX 。PCA的第一主成分方向是使样本方差最大的方向：

max_{||w||=1} wᵀCw

这正是Rayleigh商的最大值，对应 C 的最大特征值。类似地，第k主成分对应于在与前k-1个主成分正交的空间中寻找Rayleigh商的最大值，这正是Courant-Fischer定理所描述的。

我们可以用NumPy演示这一过程：

# 生成随机数据
X = np.random.randn(100, 10)  # 100个样本，10维
X = X - X.mean(axis=0)  # 中心化

# 计算协方差矩阵
C = X.T @ X / (X.shape[0]-1)

# 直接SVD
U, s, Vt = np.linalg.svd(X, full_matrices=False)

# 通过Courant-Fischer逐步寻找主成分
components = []
for k in range(1, 4):
    if k == 1:
        # 第一个主成分是Rayleigh商最大方向
        def objective(w):
            w = w / np.linalg.norm(w)
            return -(w.T @ C @ w)
        
        from scipy.optimize import minimize
        res = minimize(objective, np.random.randn(10))
        w1 = res.x / np.linalg.norm(res.x)
        components.append(w1)
    else:
        # 后续主成分需要在已找到的成分的正交补空间中寻找
        def constrained_objective(w, prev_components):
            w = w / np.linalg.norm(w)
            for pc in prev_components:
                w = w - (w @ pc) * pc
            w = w / np.linalg.norm(w)
            return -(w.T @ C @ w)
        
        res = minimize(constrained_objective, np.random.randn(10), 
                      args=(components,))
        wk = res.x
        for pc in components:
            wk = wk - (wk @ pc) * pc
        wk = wk / np.linalg.norm(wk)
        components.append(wk)

# 比较两种方法得到的主成分
print("SVD得到的右奇异向量(前3个):")
print(Vt[:3].T[:, :3])
print("\nCourant-Fischer逐步优化的主成分:")
print(np.array(components).T)

这个例子展示了Courant-Fischer定理如何指导我们逐步寻找数据中方差最大的方向，这正是PCA的核心思想。