奇异值分解(SVD):线性代数的瑞士军刀

奇异值分解(SVD):线性代数的瑞士军刀

本节摘要:SVD 是线性代数的瑞士军刀——每个矩阵都有一个,每个数据科学家都需要一个。特征分解只对方阵且要求满秩特征向量,而 SVD 对任意形状、任意秩的矩阵都成立,无任何附加条件。它把矩阵拆成三个因子,揭示矩阵对空间做的几何动作:A = U·Σ·Vᵀ——Vᵀ 旋转输入、Σ 沿轴缩放、U 旋转到输出。本节从几何直觉出发,讲清左奇异向量(输出方向)、奇异值(缩放因子)、右奇异向量(输入方向),以及 A·vᵢ = σᵢ·uᵢ 的坐标级图像;给出外积形式与 Eckart-Young 定理(截断 SVD 是最优低秩近似);用从零的幂迭代实现 SVD;演示 SVD 在图像压缩(800×600 图存成 k×1401 个值)、推荐系统(Netflix 大奖的潜在因子)、NLP 潜在语义分析、去噪、伪逆求解最小二乘上的真实威力;最后点明 PCA 本质就是 SVD——sklearn 的 PCA 内部就用 SVD 而非协方差特征分解,因为后者平方条件数、损失精度。

对应原课程:Phase 01 · Lesson 11 · singular-value-decomposition(原英文 phases/01-math-foundations/11-singular-value-decomposition/docs/en.md)。前置:第 1、2、3 节(线性代数)。

学习目标

阅读完本节,你应当能够:

1. 用幂迭代从零实现 SVD,解释 U、Σ、Vᵀ 的几何含义。
2. 应用截断 SVD 做图像压缩,权衡压缩率与重建误差。
3. 用 SVD 算 Moore-Penrose 伪逆,求解超定最小二乘系统。
4. 把 SVD 接上 PCA、推荐系统(潜在因子)、NLP 潜在语义分析。

一、问题与直觉

你有一个 1000×2000 的矩阵(用户-电影评分、文档-词频、图像像素),要压缩它、去噪、找隐藏结构或解最小二乘。特征分解只对方阵且要求满秩特征向量;SVD 对任意矩阵、任意形状、任意秩都成立,无任何条件。它把矩阵拆成三个因子,揭示矩阵对空间做了什么。它是线性代数里最通用、最有用的分解。

1.1 SVD 的几何含义

每个矩阵,无论形状,都依次做三件事:旋转、缩放、旋转。SVD 让这分解显式:

A = U · Σ · Vᵀ
   m×n     m×m   m×n   n×n
  (任意)  (旋转) (缩放) (旋转)
​

• Vᵀ 在输入空间(n 维)旋转;
• Σ 沿每个轴缩放(拉伸或压缩);
• U 把结果旋转到输出空间(m 维)。

想象:SVD 接过一个矩阵,告诉你「它把一个输入球面先按 Vᵀ 旋转、再用 Σ 拉成椭球、最后用 U 旋转到输出」。奇异值就是椭球各轴的长度。

1.2 完整分解

对 m×n 矩阵 A:

A = U·Σ·Vᵀ,其中:
  U     m×m,正交(UᵀU = I)
  Σ     m×n,对角(奇异值在对角)
  V     n×n,正交(VᵀV = I)
奇异值 σ₁ ≥ σ₂ ≥ … ≥ σᵣ > 0,r = rank(A)
​

U 的列是左奇异向量,V 的列是右奇异向量,Σ 对角元是奇异值(非负、降序)。

1.3 三组件的几何含义

• 右奇异向量(V 的列):输入空间 ℝⁿ 的标准正交基,是被矩阵映到输出空间正交方向的「自然输入坐标」。
• 奇异值(Σ 对角):缩放因子,σᵢ 告诉你矩阵沿第 i 个右奇异向量拉伸多少;σ 为 0 表示该方向被压扁。
• 左奇异向量(U 的列):输出空间 ℝᵐ 的标准正交基;第 i 个左奇异向量是第 i 个右奇异向量(缩放后)的落点方向。

关系:A·vᵢ = σᵢ·uᵢ——矩阵把第 i 个右奇异向量按 σᵢ 缩放,映到第 i 个左奇异向量。这给了「矩阵做什么」的逐坐标图像。

1.4 外积形式

SVD 可写成秩 1 矩阵之和:

A = σ₁·u₁·v₁ᵀ + σ₂·u₂·v₂ᵀ + … + σᵣ·uᵣ·vᵣᵀ
​

每一项 σᵢ·uᵢ·vᵢᵀ 是秩 1 矩阵(外积)。第一项捕获最重要的单一模式,第二项次之,依此类推。截断求和即任意秩的最优近似(Eckart-Young 定理)。

1.5 与特征分解的关系

A 的奇异值/向量直接来自 AᵀA 与 AAᵀ 的特征值/向量:

• 右奇异向量 V 是 AᵀA 的特征向量,σᵢ² 是 AᵀA 的特征值;
• 左奇异向量 U 是 AAᵀ 的特征向量,其特征值也是 σᵢ²。

三个推论:

1. 奇异值恒实非负(半正定矩阵特征值的平方根)。
2. 你可以经 AᵀA 特征分解算 SVD,但这平方了条件数、损失精度,专用 SVD 算法避免之。
3. A 方阵且对称半正定时,SVD 与特征分解是同一件事。

1.6 截断 SVD:低秩近似

Eckart-Young-Mirsky 定理:对 A 的最佳秩 k 近似(无论 Frobenius 还是谱范数),就是保留前 k 个奇异值及对应向量:

A_k = U_k · Σ_k · V_kᵀ
近似误差(谱范数)   = σ_{k+1}
近似误差(Frobenius) = √(σ_{k+1}² + … + σᵣ²)
​

这不仅是「一个好的」近似,而是可证最优的秩 k 近似——没有别的秩 k 矩阵比它更接近 A。若奇异值衰减快,小 k 就能捕获大部分矩阵;衰减慢则无低秩结构。

1.7 图像压缩

灰度图是像素强度矩阵。800×600 图有 48 万值。SVD 用远少于此近似它:

原始:800×600 = 480000 值
秩 k 近似:U_k(800×k) + Σ_k(k) + V_k(600×k) = k·(800+600+1) = k·1401 值
  k=10:  14010 值  (2.9%)
  k=50:  70050 值  (14.6%)
  k=100: 140100 值 (29.2%)
​

关键洞见:自然图像奇异值衰减快——前几个捕获大结构(形状、梯度),后面的捕获细节与噪声。秩 50 截断常能给出肉眼几乎无差、却省 85% 存储的图。

1.8 推荐系统

Netflix 大奖让 SVD 出名。用户-电影评分矩阵大部分项缺失,但该矩阵有低秩——用户口味并不完全独立,有少数潜在因子(动作 vs 文艺、老 vs 新、烧脑 vs 直观)解释大部分偏好。SVD 把(填充后的)评分矩阵分解为 U(用户在潜在因子空间的画像)、Σ(各因子重要性)、Vᵀ(电影画像)。预测评分 = 用户画像与电影画像的点积(按奇异值加权),低秩近似填上缺失项。实战用 Simon Funk 的增量 SVD 或 ALS(交替最小二乘)直接处理缺失,但核心都是 SVD 潜在因子分解。

1.9 NLP 潜在语义分析(LSA)

LSA(或 LSI)把 SVD 用于「词-文档」矩阵。秩 k=2 截断后:每个文档变成 2D「概念空间」里的点,每个词也变成同空间里的点;同主题文档聚拢,同义词聚拢(因它们出现在相似文档里,SVD 把它们归到同一潜在维度)。LSA 是最早从原始文本捕获语义相似度的成功方法之一,现代词嵌入(Word2Vec、GloVe)可视为其思想后裔。

1.10 去噪

带噪数据的信号集中在前几个奇异值,噪声分散在所有奇异值上。截断去掉噪声底。当奇异值有明显「间隙」时,间隙之上是信号(保留)、之下是噪声(丢弃)。这用于信号处理、科学测量、数据清洗——任何被加性噪声污染的矩阵,截断 SVD 都是分离信噪的原则性方法。

1.11 SVD 求伪逆

Moore-Penrose 伪逆 A⁺ 把矩阵求逆推广到非方阵与奇异矩阵,SVD 让它变得平凡:

若 A = U·Σ·Vᵀ,则 A⁺ = V·Σ⁺·Uᵀ
其中 Σ⁺:转置 Σ,把每个非零对角元 σᵢ 换成 1/σᵢ,零保持零。
​

伪逆解最小二乘:若 Ax = b 无精确解(超定系统),则 x = A⁺b 是最小二乘解(最小化 ‖Ax−b‖)。结果与正规方程 (AᵀA)⁻¹Aᵀb 相同,但数值更稳。

1.12 数值稳定性优势

经 AᵀA 算特征分解会平方奇异值(AᵀA 特征值是 σᵢ²),从而平方条件数、放大误差:

A 奇异值 [1000, 1, 0.001],条件数 10⁶
AᵀA 特征值 [10⁶, 1, 10⁻⁶],条件数 10¹²  → 多丢 6 位精度
​

现代 SVD 算法(Golub-Kahan 双对角化)直接在 A 上算,从不形成 AᵀA。这就是为什么总该用 np.linalg.svd(A) 而非 np.linalg.eig(A.T @ A)。

1.13 与 PCA 的关系

PCA 本质就是 SVD 在中心化数据上的应用——这不是类比,是字面意义上的同一计算。

中心化数据矩阵 X(n_samples × n_features):
  协方差 C = (1/(n−1))·XᵀX
  X = U·Σ·Vᵀ(SVD)
  XᵀX = V·Σ²·Vᵀ
  C   = (1/(n−1))·V·Σ²·Vᵀ
故主成分恰是右奇异向量 V,每成分解释方差 = σᵢ²/(n−1)。
sklearn 的 PCA 用 SVD 而非特征分解实现——更快更稳。
​

第 10 节学的降维,全是 SVD 在底下。PCA 是 SVD 在 ML 里最常见的应用。

二、从零实现

完整源码见 phases/01-math-foundations/11-singular-value-decomposition/code/(Python svd.py,Julia svd.jl)。

2.1 用幂迭代从零算 SVD

思路:对 AᵀA 用幂迭代找最大奇异值与向量,再「缩减」矩阵,重复找下一个。

def power_iteration(M, num_iters=100):
    v = np.random.randn(M.shape[1]); v /= np.linalg.norm(v)
    for _ in range(num_iters):
        v = M @ v; v /= np.linalg.norm(v)
    eigenvalue = v @ M @ v
    return eigenvalue, v

def svd_from_scratch(A, k=None):
    m, n = A.shape
    k = k or min(m, n)
    sigmas, us, vs = [], [], []
    A_residual = A.copy().astype(float)
    for _ in range(k):
        eigenvalue, v = power_iteration(A_residual.T @ A_residual, num_iters=200)
        if eigenvalue < 1e-10: break          # 剩余方向已被压扁
        sigma = np.sqrt(eigenvalue)
        u = A_residual @ v / sigma
        sigmas.append(sigma); us.append(u); vs.append(v)
        A_residual -= sigma * np.outer(u, v)  # 缩减,找下一个
    return np.column_stack(us), np.array(sigmas), np.column_stack(vs)
​

设计要点:幂迭代对最大特征向量线性收敛;每轮「缩减排掉已找到的秩 1 成分」,等价于找剩余矩阵的最大奇异方向。教育性强但远慢于 LAPACK 的 Golub-Kahan。

2.2 图像压缩

def compress_image_svd(image, k):
    U, S, Vt = np.linalg.svd(image, full_matrices=False)
    return U[:, :k] @ np.diag(S[:k]) @ Vt[:k, :]

for k in [1, 5, 10, 20, 50]:
    compressed = compress_image_svd(image, k)
    error = np.linalg.norm(image - compressed) / np.linalg.norm(image)
    ratio = k * (rows + cols + 1) / (rows * cols)
    print(f"k={k}  error={error:.4f}  storage={ratio:.1%}")
​

2.3 伪逆解最小二乘

U, S, Vt = np.linalg.svd(A, full_matrices=False)
A_pinv = Vt.T @ np.diag(1.0 / S) @ U.T
x_svd = A_pinv @ b                          # = 最小二乘解
​

三、框架对比

NumPy/LAPACK 的 SVD 是工业级实现:

U, S, Vt = np.linalg.svd(A, full_matrices=False)   # 直接在 A 上,不形成 AᵀA
​

full_matrices=False(瘦 SVD)在 m≫n 或 n≫m 时省内存。sklearn 的 PCA 内部正是 svd(X_centered, full_matrices=False)。

💡 永远用 np.linalg.svd(A) 而非 eig(A.T @ A)——后者平方条件数,高条件数矩阵会丢多位精度。

四、可复用产物

• outputs/skill-svd.md:一份关于「何时、如何在实际项目里用 SVD」的技能说明。
• Julia 版 code/svd.jl 演示同样概念,用 Julia 原生 svd() 与 LinearAlgebra。

源码见 phases/01-math-foundations/11-singular-value-decomposition/code/。

五、练习

1. (Easy) 不用幂迭代,改用 AᵀA 特征分解从零实现完整 SVD(V 与奇异值来自特征分解,U = AVΣ⁻¹),与幂迭代版、NumPy 比较精度。
2. (Medium) 加载真实灰度图,在秩 1、5、10、25、50、100 压缩,算每秩压缩率与相对误差,找视觉可接受的秩。
3. (Medium) 搭微型推荐系统:10×8 用户-电影矩阵,缺失项用行均值填,SVD 后秩 3 重建,预测缺失评分,验证合理性。
4. (Hard) 造 100×50 文档-词矩阵含 3 个合成主题(每主题 5 词)加噪声,SVD 后验证前 3 个奇异值远大于其余;把文档投到 3D 潜在空间,验证同主题文档聚拢。
5. (Hard) 造秩 3 的 50×40 干净矩阵,加不同水平高斯噪声(σ=0.1、0.5、1.0、2.0);对每噪声水平扫 k=1~40,对干净矩阵测重建误差,观察最优 k 如何随噪声变化。

本节要点回顾

1. SVD 对任意矩阵成立——任意形状、任意秩、无附加条件,是最通用的矩阵分解。
2. 几何三步:A = U·Σ·Vᵀ = 旋转(Vᵀ)+ 缩放(Σ)+ 旋转(U);奇异值是椭球轴长。
3. A·vᵢ = σᵢ·uᵢ:右奇异向量是输入方向,奇异值是缩放因子,左奇异向量是输出落点。
4. 外积形式把 A 拆成秩 1 之和,每项加一层结构,首项最重要。
5. Eckart-Young 定理:截断 SVD 是可证最优的秩 k 近似。
6. 图像压缩:自然图奇异值衰减快,小秩即可省 85% 存储而肉眼无差。
7. 推荐系统与 LSA:SVD 揭示潜在因子(口味、主题),填缺失评分、捕语义相似度。
8. 去噪:信号集中在前奇异值、噪声分散在所有,SVD 截断分离信噪。
9. 伪逆 A⁺ = V·Σ⁺·Uᵀ 解最小二乘,比正规方程更稳。
10. PCA 本质是 SVD:sklearn 的 PCA 用 SVD 实现,因为 AᵀA 特征分解平方条件数、损精度。

下一节,我们把矩阵推广到任意维——张量运算:深度学习里 [batch, channel, height, width] 这样的多维数组如何 reshape、广播、einsum,以及它们在 PyTorch/JAX 里的高效实现。