特征值分解与 SVD 的推导和应用场景？

核心概念

特征值分解（Eigenvalue Decomposition, EVD）是将一个方阵分解为由其特征值和特征向量表示的三个矩阵的乘积。其核心思想是找到一组特殊的向量（特征向量），当矩阵对这些向量进行线性变换时，仅仅是对其进行缩放，而方向保持不变，缩放的比例就是特征值。EVD 揭示了矩阵变换的内在“主轴”方向和拉伸强度，但它仅适用于方阵。

奇异值分解（Singular Value Decomposition, SVD）是一种更通用的矩阵分解方法，适用于任意形状的 $m \times n$ 矩阵。它将一个矩阵分解为一个酉矩阵、一个对角矩阵和一个酉矩阵的转置的乘积。SVD 揭示了矩阵所代表的线性变换的几何本质：它将输入空间中的一组标准正交基，映射为输出空间中的另一组标准正交基，并伴随着沿这些基方向的缩放，缩放因子即为奇异值。

原理与推导

特征值分解 (EVD)

1. 定义与公式

对于一个 $n \times n$ 的方阵 $A$，如果存在一个标量 $\lambda$ 和一个非零向量 $\mathbf{v}$，使得： $A\mathbf{v} = \lambda\mathbf{v}$ 那么，$\lambda$ 称为矩阵 $A$ 的一个特征值（Eigenvalue），$\mathbf{v}$ 称为对应于 $\lambda$ 的特征向量（Eigenvector）。

2. 推导

• 特征方程: 上述定义可以改写为 $(A - \lambda I)\mathbf{v} = \mathbf{0}$，其中 $I$ 是单位矩阵。由于 $\mathbf{v}$ 是非零向量，这意味着矩阵 $(A - \lambda I)$ 是奇异的（不可逆），其行列式必须为零： $\det(A - \lambda I) = 0$ 这个方程被称为特征方程，它是一个关于 $\lambda$ 的 $n$ 次多项式，求解这个方程可以得到 $n$ 个特征值（可能包含重根或复数）。

• 矩阵分解形式: 如果矩阵 $A$ 有 $n$ 个线性无关的特征向量 $\{\mathbf{v}_1, \dots, \mathbf{v}_n\}$，对应的特征值为 $\{\lambda_1, \dots, \lambda_n\}$。我们可以将这些特征向量作为列向量构成一个矩阵 $P = [\mathbf{v}_1, \dots, \mathbf{v}_n]$。那么： $AP = A[\mathbf{v}_1, \dots, \mathbf{v}_n] = [A\mathbf{v}_1, \dots, A\mathbf{v}_n] = [\lambda_1\mathbf{v}_1, \dots, \lambda_n\mathbf{v}_n]$ 同时，我们定义一个对角矩阵 $\Lambda$，其对角线元素为特征值： $\Lambda = \text{diag}(\lambda_1, \dots, \lambda_n)$ 那么，$[\lambda_1\mathbf{v}_1, \dots, \lambda_n\mathbf{v}_n] = P\Lambda$。 因此，我们得到 $AP = P\Lambda$。如果 $P$ 可逆（即特征向量线性无关），则可得特征值分解： $A = P\Lambda P^{-1}$

• 对称矩阵的特殊情况: 在机器学习中，我们经常处理对称矩阵（如协方差矩阵）。对于实对称矩阵 $A$，其特征值均为实数，且不同特征值对应的特征向量是正交的。我们可以将这些正交向量归一化为标准正交向量，此时 $P$ 变为一个正交矩阵 $Q$（满足 $Q^T Q = I$，即 $Q^{-1} = Q^T$）。分解形式变为： $A = Q\Lambda Q^T$ 这被称为谱分解（Spectral Decomposition）。

3. 几何解释

特征值分解可以看作是为矩阵 $A$ 所代表的线性变换找到了一个“最自然”的坐标系。在这个坐标系下（由特征向量张成的空间），变换 $A$ 的作用非常简单，仅仅是在各个坐标轴（特征向量方向）上进行独立的缩放，缩放因子就是对应的特征值。

4. 复杂度

求解一个 $n \times n$ 稠密矩阵的特征值分解通常使用如 QR 算法等迭代方法，其时间复杂度约为 $O(n^3)$。

奇异值分解 (SVD)

1. 定义与公式

对于任意一个 $m \times n$ 的矩阵 $A$，其奇异值分解形式为： $A = U\Sigma V^T$ 其中：

• $U$ 是一个 $m \times m$ 的酉矩阵（或正交矩阵，如果 $A$ 是实矩阵），其列向量称为左奇异向量。
• $\Sigma$ 是一个 $m \times n$ 的对角矩阵，其对角线上的元素 $\sigma_i$ 称为奇异值（Singular Values），且非负，通常按降序排列 $\sigma_1 \ge \sigma_2 \ge \dots \ge 0$。
• $V$ 是一个 $n \times n$ 的酉矩阵（或正交矩阵），其列向量称为右奇异向量。$V^T$ 是 $V$ 的共轭转置（或转置）。

2. 推导 (与 EVD 的联系)

SVD 的推导非常巧妙，它通过构造两个对称半正定矩阵 $A^TA$ 和 $AA^T$ 来与 EVD 建立联系。

• 构造 $A^TA$: 这是一个 $n \times n$ 的对称矩阵。 $A^TA = (U\Sigma V^T)^T (U\Sigma V^T) = (V\Sigma^T U^T)(U\Sigma V^T)$ 由于 $U$ 是酉矩阵，$U^TU = I$。所以： $A^TA = V(\Sigma^T\Sigma)V^T$ 这里的 $\Sigma^T\Sigma$ 是一个 $n \times n$ 的对角矩阵，其对角元素为 $\sigma_i^2$。这正是对 $A^TA$ 的特征值分解！

  • 结论1: $V$ 的列向量（右奇异向量）是 $A^TA$ 的特征向量。
  • 结论2: $\Sigma^T\Sigma$ 的对角元素（奇异值的平方 $\sigma_i^2$）是 $A^TA$ 的特征值。

• 构造 $AA^T$: 这是一个 $m \times m$ 的对称矩阵。 $AA^T = (U\Sigma V^T)(U\Sigma V^T)^T = (U\Sigma V^T)(V\Sigma^T U^T)$ 由于 $V$ 是酉矩阵，$V^TV = I$。所以： $AA^T = U(\Sigma\Sigma^T)U^T$ 这里的 $\Sigma\Sigma^T$ 是一个 $m \times m$ 的对角矩阵，其对角元素也是 $\sigma_i^2$。

  • 结论3: $U$ 的列向量（左奇异向量）是 $AA^T$ 的特征向量。

• 总结推导步骤:

  1. 计算 $n \times n$ 的对称矩阵 $A^TA$。
  2. 对 $A^TA$ 进行特征值分解，得到特征值 $\lambda_i$ 和特征向量 $\mathbf{v}_i$。
  3. 奇异值 $\sigma_i = \sqrt{\lambda_i}$。
  4. 右奇异向量矩阵 $V$ 由 $\mathbf{v}_i$ 构成。
  5. 左奇异向量 $\mathbf{u}_i$ 可以通过关系 $A\mathbf{v}_i = \sigma_i\mathbf{u}_i$ 计算得出，即 $\mathbf{u}_i = \frac{1}{\sigma_i}A\mathbf{v}_i$（对于 $\sigma_i > 0$）。

3. 几何解释

SVD 将任意线性变换 $A$ 分解为三个基本操作：

1. $V^T$ (旋转): 首先，在输入空间 $\mathbb{R}^n$ 中进行一次坐标系旋转（或反射），将标准基向量对齐到右奇异向量 $V$ 的方向上。
2. $\Sigma$ (缩放): 接着，沿着新的坐标轴（$V$ 的方向）进行缩放，每个轴的缩放比例由对应的奇异值 $\sigma_i$ 决定。如果 $m \neq n$ 或矩阵是奇异的，某些方向可能会被压缩至零或维度发生变化。
3. $U$ (旋转): 最后，在输出空间 $\mathbb{R}^m$ 中进行另一次坐标系旋转，将缩放后的向量对齐到左奇异向量 $U$ 的方向上，得到最终的变换结果。

4. 复杂度

对于一个 $m \times n$ 的稠密矩阵，SVD 的计算复杂度通常是 $O(\min(m^2n, mn^2))$。

代码实现

code_drill: 手写幂迭代求最大特征值

幂迭代（Power Iteration）是一种简单而有效的迭代算法，用于计算矩阵绝对值最大的特征值及其对应的特征向量。

 1import numpy as np
 2
 3def power_iteration(A, num_simulations: int = 100):
 4    """
 5    使用幂迭代法计算矩阵A的最大特征值和对应的特征向量。
 6
 7    参数:
 8    A (np.ndarray): 输入的方阵 (n, n)。
 9    num_simulations (int): 迭代次数。
10
11    返回:
12    tuple: (最大特征值, 对应的特征向量)
13    """
14    # 为什么这样做: 随机初始化一个向量b_k。理论上，只要这个向量不与第二大特征向量
15    # 构成的子空间正交，算法就能收敛。随机选择几乎总能保证这一点。
16    n = A.shape[0]
17    b_k = np.random.rand(n)
18
19    for _ in range(num_simulations):
20        # 为什么这样做: 这是幂迭代的核心步骤。每次将向量b_k左乘矩阵A。
21        # 如果b_k可以表示为特征向量的线性组合，每次相乘会使得最大特征值
22        # 对应的分量在幅度上增长得最快，逐渐主导整个向量。
23        b_k1 = np.dot(A, b_k)
24
25        # 为什么这样做: 归一化。如果不归一化，向量的模会趋向于无穷大（如果|λ_max|>1）
26        # 或零（如果|λ_max|<1），导致数值溢出或下溢。归一化后，向量b_k会
27        # 收敛到最大特征值对应的特征向量。
28        b_k1_norm = np.linalg.norm(b_k1)
29        b_k = b_k1 / b_k1_norm
30
31    # 为什么这样做: 使用瑞利商(Rayleigh quotient)计算特征值。
32    # 当b_k已经收敛到特征向量v时，Av ≈ λv。
33    # 那么 v.T * A * v ≈ v.T * λ * v = λ * (v.T * v)。
34    # 因为v是单位向量，v.T * v = 1，所以 λ ≈ v.T * A * v。
35    eigenvalue = np.dot(b_k.T, np.dot(A, b_k))
36
37    return eigenvalue, b_k
38
39# --- 示例 ---
40# 创建一个对称矩阵，确保特征值为实数且有明确的最大值
41A = np.array([[4, 1, 1],
42              [1, 3, -1],
43              [1, -1, 2]])
44
45# 使用手写的幂迭代
46max_eigenvalue_manual, max_eigenvector_manual = power_iteration(A)
47print("--- 手写幂迭代结果 ---")
48print(f"最大特征值: {max_eigenvalue_manual:.6f}")
49print(f"对应特征向量: {max_eigenvector_manual}")
50
51# 使用NumPy内置函数进行验证
52eigenvalues, eigenvectors = np.linalg.eig(A)
53max_idx = np.argmax(np.abs(eigenvalues))
54max_eigenvalue_np = eigenvalues[max_idx]
55max_eigenvector_np = eigenvectors[:, max_idx]
56
57# 确保符号一致以便比较
58if np.sign(max_eigenvector_manual[0]) != np.sign(max_eigenvector_np[0]):
59    max_eigenvector_np = -max_eigenvector_np
60
61print("\n--- NumPy linalg.eig 验证结果 ---")
62print(f"最大特征值: {max_eigenvalue_np:.6f}")
63print(f"对应特征向量: {max_eigenvector_np}")

code_drill: 用 SVD 压缩一张图像并可视化不同秩的重建效果

SVD 的一个重要应用是低秩近似。通过保留最大的 k 个奇异值及其对应的奇异向量，我们可以得到原矩阵的最佳 k 秩近似，从而实现数据压缩。

 1import numpy as np
 2import matplotlib.pyplot as plt
 3from skimage import data
 4from skimage.color import rgb2gray
 5from skimage.transform import resize
 6
 7def svd_image_compression(image, ranks):
 8    """
 9    使用SVD对图像进行压缩和重建，并可视化结果。
10
11    参数:
12    image (np.ndarray): 输入的灰度图像矩阵。
13    ranks (list of int): 用于重建的不同秩(k)的列表。
14    """
15    # 为什么这样做: 对图像矩阵进行SVD分解。
16    # U包含图像的行空间信息，V^T包含列空间信息，s包含能量（重要性）信息。
17    # 奇异值s已经按降序排列，最大的奇异值对应图像中最重要的模式。
18    U, s, Vt = np.linalg.svd(image, full_matrices=False)
19    
20    # full_matrices=False 进行经济SVD，U为(m,k), s为(k,), Vt为(k,n)，其中k=min(m,n)
21    # 这在计算上更高效，且不影响重建结果。
22
23    fig, axes = plt.subplots(1, len(ranks) + 1, figsize=(15, 5))
24    
25    # 显示原始图像
26    axes[0].imshow(image, cmap='gray')
27    axes[0].set_title("原始图像\nRank={}".format(np.linalg.matrix_rank(image)))
28    axes[0].axis('off')
29
30    for i, k in enumerate(ranks):
31        # 为什么这样做: 这就是低秩近似的核心。我们只保留前k个最大的奇异值
32        # 及其对应的左右奇异向量。这相当于抓住了图像中最主要的k个“成分”。
33        U_k = U[:, :k]
34        s_k = s[:k]
35        Vt_k = Vt[:k, :]
36
37        # 为什么这样做: 通过截断后的U, s, V^T重建图像。
38        # np.diag(s_k)将一维奇异值数组转换为对角矩阵。
39        # 重建的矩阵是原矩阵在Frobenius范数意义下的最佳k秩近似。
40        reconstructed_image = U_k @ np.diag(s_k) @ Vt_k
41        
42        # 计算压缩率
43        original_size = image.shape[0] * image.shape[1]
44        # 存储U_k, s_k, Vt_k所需的元素数量
45        compressed_size = U_k.size + s_k.size + Vt_k.size
46        compression_ratio = compressed_size / original_size * 100
47
48        axes[i+1].imshow(reconstructed_image, cmap='gray')
49        axes[i+1].set_title(f"重建 (k={k})\n压缩率: {compression_ratio:.2f}%")
50        axes[i+1].axis('off')
51        
52    plt.tight_layout()
53    plt.show()
54
55# --- 示例 ---
56# 加载一张示例图片并转换为灰度图
57# 为什么这样做: SVD直接作用于2D矩阵，所以需要先将RGB三通道图像转为单通道灰度图。
58# 同时缩小尺寸以加速计算。
59plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
60plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
61
62# 使用scikit-image自带的图像
63original_image = data.camera()
64# 调整尺寸以加快SVD计算
65image_gray = resize(original_image, (256, 256), anti_aliasing=True)
66
67# 定义不同的秩进行比较
68ranks_to_test = [5, 20, 50, 100]
69
70svd_image_compression(image_gray, ranks_to_test)

工程实践

• 使用场景:

  • EVD:
    • PCA (主成分分析): 核心就是对数据协方差矩阵进行特征值分解。特征向量是主成分方向，特征值是沿该方向的方差。
    • 谱聚类 (Spectral Clustering): 对图拉普拉斯矩阵进行特征值分解，利用其特征向量进行聚类。
    • 马尔可夫链: 分析状态转移矩阵的稳态分布，与最大特征值为1对应的特征向量有关。
  • SVD:
    • 降维与数据压缩: 如图像压缩、特征降维。截断SVD（Truncated SVD）是实现这一目标的主要手段。
    • 推荐系统: 矩阵分解模型（如FunkSVD）通过SVD的思想来发现用户和物品的隐向量，填补评分矩阵中的缺失值。
    • NLP中的LSA: 潜在语义分析（Latent Semantic Analysis）对词-文档矩阵进行SVD，发现词与文档之间的潜在主题。
    • 数值稳定性: 用于计算伪逆（Moore-Penrose Pseudoinverse），$A^+ = V\Sigma^+U^T$，其中 $\Sigma^+$ 是将 $\Sigma$ 的非零元素取倒数后转置得到。这对于求解病态或非方阵的线性方程组 $A\mathbf{x}=\mathbf{b}$ 非常鲁棒。

• 超参数选择:

  • 在使用截断SVD或PCA时，最重要的超参数是保留的秩 k。
  • 经验法则:
    1. 能量保持: 绘制奇异值（或其平方）的累积和曲线（Scree Plot）。选择一个 k，使得保留的奇异值能量（平方和）占总能量的90%~99%。$\text{Energy}(k) = \frac{\sum_{i=1}^k \sigma_i^2}{\sum_{i=1}^r \sigma_i^2}$。
    2. 肘部法则 (Elbow Method): 在奇异值大小的折线图上，寻找曲线斜率变化最剧烈的“肘部”，该点之后的奇异值迅速减小，可被视为噪声。
    3. 任务驱动: 根据下游任务的性能（如分类准确率、重建误差）来交叉验证选择最佳的 k。

• 性能 / 显存 / 吞吐 的权衡:

  • 完整SVD/EVD: 计算成本高 ($O(n^3)$)，不适用于大矩阵。
  • 截断SVD: 对于只需要前 k 个奇异值/向量的场景，应使用专门的迭代算法，如 scipy.sparse.linalg.svds，其复杂度大致为 $O(mnk)$，远快于完整SVD。
  • 随机SVD (Randomized SVD): 对于超大规模矩阵，随机SVD是一种非常高效的近似算法。它通过随机投影将大矩阵降维到一个小矩阵，对小矩阵进行SVD，然后映射回原空间。速度极快，且精度通常很高。sklearn.utils.extmath.randomized_svd 是一个很好的实现。

• 常见坑和调试技巧:

  • 中心化: 在将SVD用于PCA时，必须先对数据进行中心化（减去均值），否则第一个主成分会和数据的均值高度相关，无法正确反映方差最大的方向。
  • $A^TA$ 的数值问题: 直接计算 $A^TA$ 来求SVD可能会导致数值不稳定，因为矩阵的条件数会被平方 ($\text{cond}(A^TA) = \text{cond}(A)^2$)，加剧病态问题。专业的SVD库（如LAPACK，被NumPy/SciPy调用）使用更稳定的算法，如Golub-Kahan双对角化，直接在 $A$ 上操作。

常见误区与边界情况

• EVD vs. SVD 的混淆:

  • 适用范围: EVD仅限方阵；SVD适用于任意矩阵。
  • 基的正交性: EVD的特征向量不总是正交的（仅当矩阵为正规矩阵，如对称/反对称/酉矩阵时才正交）；SVD的左、右奇异向量各自构成标准正交基。
  • 值的性质: EVD的特征值可以是复数；SVD的奇异值总是实数且非负。
  • 关系: 对于实对称矩阵 $A=Q\Lambda Q^T$，其SVD为 $A=U\Sigma V^T$。此时，$U$ 和 $V$ 与 $Q$ 非常相似（可能差一个符号），奇异值 $\sigma_i = |\lambda_i|$。

• 幂迭代的边界情况:

  • 收敛失败: 如果绝对值最大的特征值不唯一（例如，$\lambda_1=5, \lambda_2=-5$），幂迭代可能不会收敛到一个单一向量，而是在两个向量之间振荡。
  • 初始向量: 如果初始随机向量恰好与最大特征值对应的特征向量正交，算法会收敛到第二大特征值。在浮点数运算中，这种情况发生的概率极低。

• SVD的唯一性:

  • 奇异值 $\sigma_i$ 是唯一确定的。
  • 如果所有奇异值都是唯一的正数，那么左右奇异向量也唯一确定（最多差一个正负号）。
  • 如果存在重复的奇异值，对应的奇异向量在其张成的子空间内可以任意旋转，不再唯一。

• 常见面试追问:

  • 问: PCA和SVD有什么关系？
  • 答: 对数据矩阵 $X$（已中心化）进行PCA，等价于对 $X$ 进行SVD。$X$ 的右奇异向量 $V$ 就是主成分方向。具体来说，协方差矩阵 $C = \frac{1}{n-1}X^TX$，对 $C$ 做EVD得到 $C = P\Lambda P^{-1}$。同时，$X=U\Sigma V^T$，则 $X^TX = V\Sigma^T\Sigma V^T$。对比可知，PCA的主成分（$P$的列）就是SVD的右奇异向量（$V$的列）。
  • 问: 为什么在推荐系统中，我们不直接用标准的SVD，而是用矩阵分解算法（如ALS, FunkSVD）？
  • 答: 推荐系统中的评分矩阵通常是极其稀疏的。标准SVD无法处理含有大量缺失值的矩阵。而矩阵分解算法（如ALS）是专门为稀疏矩阵设计的迭代优化算法，它只在已知的评分上计算损失并更新隐向量，从而有效地为所有用户和物品找到隐向量表示。
  • 问: 奇异值的大小代表什么物理意义？
  • 答: 奇异值衡量了矩阵在对应奇异向量方向上的“拉伸”或“重要性”程度。在图像压缩中，它代表了某个模式对图像整体贡献的大小。在PCA中，它与主成分方向上的数据方差成正比。在信息论角度，大奇异值对应了信号，小奇异值往往对应了噪声。