SVD分解

SVD分解

小象实战讲义 · 人工智能数学基础

在机器学习和数据分析中，我们常常需要处理高维、大规模的矩阵数据。如何有效地压缩数据、去除噪声或提取核心特征，是实际应用中的关键问题。奇异值分解（Singular Value Decomposition, SVD）作为一种强大的矩阵分解工具，为解决这些问题提供了统一的数学框架。本节我们将深入理解SVD的定义、性质及其在图像压缩、降噪和主成分分析（PCA）中的核心应用。

💡 核心导读

• SVD的定义与存在性：理解任意矩阵均可分解为三个特定矩阵的乘积形式，并掌握其矩阵维数的对应关系。
• 奇异值与特征值的关系：掌握奇异值来源于矩阵 $A^TA$ 或 $AA^T$ 的特征值平方根，并理解其几何与代数意义。
• 分解的几何解释与低秩近似：从加权和的角度理解SVD，掌握利用前 $r$ 个奇异值进行矩阵低秩近似的原理与方法。
• 核心应用实践：学习SVD在图像压缩、图像降噪中的具体实现，并理解其与主成分分析（PCA）的内在等价关系。

SVD的定义与基本性质

对于任意一个 $m \times n$ 的实矩阵 $A$，都存在一个奇异值分解（SVD），使得 $A$ 可以表示为三个矩阵的乘积：

$$A = U \Sigma V^T$$

其中：

• $U$ 是一个 $m \times m$ 的正交矩阵，即 $U^T U = U U^T = I_m$。其列向量 $\mathbf{u}_1, \dots, \mathbf{u}_m$ 称为 $A$ 的左奇异向量。
• $V$ 是一个 $n \times n$ 的正交矩阵，即 $V^T V = V V^T = I_n$。其列向量 $\mathbf{v}_1, \dots, \mathbf{v}_n$ 称为 $A$ 的右奇异向量。
• $\Sigma$ 是一个 $m \times n$ 的对角矩阵（更准确地说，是矩形对角矩阵），其形式取决于 $m$ 和 $n$ 的大小关系：
  • 若 $m \ge n$，则 $\Sigma$ 形如 $\begin{bmatrix} D \ 0 \end{bmatrix}$，其中 $D = \text{diag}(\sigma_1, \sigma_2, \dots, \sigma_n)$。
  • 若 $m \le n$，则 $\Sigma$ 形如 $\begin{bmatrix} D & 0 \end{bmatrix}$，其中 $D = \text{diag}(\sigma_1, \sigma_2, \dots, \sigma_m)$。

矩阵 $\Sigma$ 主对角线（即行标与列标相等的元素）上的非零元素 $\sigma_1, \sigma_2, \dots, \sigma_p$（其中 $p = \min(m, n)$）称为矩阵 $A$ 的奇异值。我们约定奇异值按降序排列：

$$\sigma_1 \ge \sigma_2 \ge \dots \ge \sigma_p \ge 0$$

奇异值与矩阵 $A^TA$ 和 $AA^T$ 的特征值有直接关系。$A^TA$ 是一个 $n \times n$ 的半正定矩阵，其特征值 $\lambda_i$ 均为非负实数。奇异值 $\sigma_i$ 正是这些特征值的平方根：

$$\sigma_i = \sqrt{\lambda_i(A^TA)}, \quad i = 1, \dots, n$$

同理，$AA^T$ 的特征值也包含这些 $\lambda_i$（另有 $|m-n|$ 个零特征值）。因此，奇异值可以看作是矩阵 $A$ 的“拉伸”强度的度量。

SVD的推导：与特征值分解的联系

SVD的存在性及其性质可以通过特征值分解来理解和证明。假设 $A = U \Sigma V^T$ 成立，我们来考察 $A^TA$ 和 $AA^T$。

1. 考察 $A^TA$： $$\begin{aligned} A^TA &= (U \Sigma V^T)^T (U \Sigma V^T) \ &= V \Sigma^T U^T U \Sigma V^T \ &= V \Sigma^T \Sigma V^T \quad (\text{因为 } U^T U = I) \ &= V \Lambda V^T \end{aligned}$$ 其中 $\Lambda = \Sigma^T \Sigma$ 是一个 $n \times n$ 的对角矩阵。由于 $\Sigma$ 是广义对角矩阵，$\Sigma^T \Sigma$ 的计算结果是一个标准的对角方阵，其对角线元素正是 $\sigma_i^2$。 上式 $A^TA = V \Lambda V^T$ 正是实对称矩阵 $A^TA$ 的正交对角化。因此：

• $V$ 的列向量 $\mathbf{v}_i$ 是 $A^TA$ 的单位特征向量，对应于特征值 $\lambda_i = \sigma_i^2$。
• $V$ 就是使得 $A^TA$ 对角化的正交矩阵。

2. 考察 $AA^T$： $$\begin{aligned} AA^T &= (U \Sigma V^T) (U \Sigma V^T)^T \ &= U \Sigma V^T V \Sigma^T U^T \ &= U \Sigma \Sigma^T U^T \quad (\text{因为 } V^T V = I) \ &= U \tilde{\Lambda} U^T \end{aligned}$$ 其中 $\tilde{\Lambda} = \Sigma \Sigma^T$ 是一个 $m \times m$ 的对角矩阵，其前 $p$ 个对角线元素也是 $\sigma_i^2$，其余为0。 同理，上式表明 $U$ 的列向量 $\mathbf{u}_i$ 是 $AA^T$ 的单位特征向量，对应于非零特征值 $\sigma_i^2$。

这个推导过程也揭示了计算SVD的一种理论方法：分别计算 $A^TA$ 和 $AA^T$ 的特征值与特征向量，然后取特征值的平方根得到奇异值，并将对应的特征向量排列成 $V$ 和 $U$。需要注意的是，由于特征向量在符号选择（$\pm$）和特征子空间内正交基选取上的自由度，SVD分解本身不是唯一的。但在奇异值降序排列的约定下，奇异值本身是唯一的。

SVD的几何解释与低秩近似

将矩阵 $U$ 和 $V$ 按列分块：$U = [\mathbf{u}_1, \dots, \mathbf{u}_m]$， $V = [\mathbf{v}_1, \dots, \mathbf{v}_n]$。为简化讨论，假设 $m \ge n$，则 $\Sigma = \begin{bmatrix} D \ 0 \end{bmatrix}$，其中 $D = \text{diag}(\sigma_1, \dots, \sigma_n)$。

将SVD展开为矩阵乘法的和形式： $$\begin{aligned} A &= U \Sigma V^T = [\mathbf{u}_1, \dots, \mathbf{u}_m] \begin{bmatrix} \sigma_1 & & \ & \ddots & \ & & \sigma_n \ & 0 & \end{bmatrix} \begin{bmatrix} \mathbf{v}_1^T \ \vdots \ \mathbf{v}_n^T \end{bmatrix} \ &= \sigma_1 \mathbf{u}_1 \mathbf{v}_1^T + \sigma_2 \mathbf{u}_2 \mathbf{v}_2^T + \dots + \sigma_n \mathbf{u}_n \mathbf{v}_n^T \end{aligned}$$

这里，每一项 $\mathbf{u}_i \mathbf{v}_i^T$ 是一个 $m \times n$ 的矩阵，其秩为1。因此，SVD将原矩阵 $A$ 分解为 $n$ 个秩为1的矩阵的加权和，权重就是奇异值 $\sigma_i$。

由于奇异值按降序排列，$\sigma_1$ 对应的分量 $\sigma_1 \mathbf{u}_1 \mathbf{v}_1^T$ 对 $A$ 的“贡献”最大，$\sigma_2$ 次之，依此类推。在许多实际矩阵（如图像数据、用户-物品评分矩阵）中，奇异值衰减得非常快，前10%甚至1%的奇异值之和就可能占到全部奇异值之和的99%以上。

这启发我们进行低秩近似：只保留前 $r$（$r \ll \min(m, n)$）个最大的奇异值及其对应的左右奇异向量，用它们的加权和来近似原矩阵 $A$。

$$A \approx A_r = \sum_{i=1}^{r} \sigma_i \mathbf{u}_i \mathbf{v}_i^T = U_r \Sigma_r V_r^T$$

其中：

• $U_r$ 是 $U$ 的前 $r$ 列，维度为 $m \times r$。
• $\Sigma_r$ 是 $\Sigma$ 的左上角 $r \times r$ 子矩阵，即 $\text{diag}(\sigma_1, \dots, \sigma_r)$。
• $V_r$ 是 $V$ 的前 $r$ 列，维度为 $n \times r$，因此 $V_r^T$ 的维度为 $r \times n$。

这种近似在Frobenius范数意义下是最优的，即 $A_r$ 是所有秩不超过 $r$ 的矩阵中，与 $A$ 的Frobenius距离最小的那个。

低秩近似带来了两大核心应用：

1. 数据压缩：存储原始矩阵 $A$ 需要 $m \times n$ 个元素。存储近似矩阵 $A_r$ 只需要存储 $U_r$、$\Sigma_r$ 和 $V_r$，总计 $mr + r + nr = r(m + n + 1)$ 个元素。当 $r$ 很小时，存储量大幅降低。
2. 噪声抑制：信号或数据中的主要信息往往蕴含在较大的奇异值对应的分量中，而噪声或细微扰动则体现在较小的奇异值分量里。通过舍弃这些小奇异值分量，可以达到滤除噪声、保留主体结构的目的。

应用实践：图像压缩与降噪

图像可以看作一个像素值矩阵。对图像矩阵进行SVD，并利用其低秩近似特性，可以实现图像的有损压缩和降噪。

图像压缩

压缩过程：

1. 将灰度图像读入为一个 $m \times n$ 的矩阵 $A$（每个元素代表一个像素的亮度）。
2. 对 $A$ 进行SVD分解：$A = U \Sigma V^T$。
3. 选择截断秩 $r$，计算近似矩阵 $A_r = U_{:, 1:r} \cdot \Sigma_{1:r, 1:r} \cdot V_{:, 1:r}^T$。
4. 将 $A_r$ 的值重新缩放到合法的像素值范围（如0-255），并保存为图像。

压缩比约为：$\frac{r(m + n + 1)}{m \times n}$。$r$ 越小，压缩比越高，但图像质量损失也越大。

图像降噪

降噪过程与压缩类似，其核心思想是：噪声通常对应于较小的奇异值。通过选择一个合适的截断秩 $r$，舍弃那些代表噪声的小奇异值分量，即可在尽量保留图像主体信息的同时抑制噪声。

选择 $r$ 是关键：$r$ 太小会损失有用信息（如图像边缘变得模糊），$r$ 太大则降噪效果不明显。通常需要根据具体图像和噪声特性进行权衡。

下面的代码示例演示了如何利用 numpy 对一张示例图像进行SVD压缩与重构。

import numpy as np
import matplotlib.pyplot as plt

# 1. 生成示例图像（用数学函数合成一张灰度测试图并叠加噪声，无需联网下载）
np.random.seed(42)
m, n = 200, 200
gx = np.linspace(0, 4 * np.pi, n)
gy = np.linspace(0, 4 * np.pi, m)
X, Y = np.meshgrid(gx, gy)
A = 128 + 80 * np.sin(X) * np.cos(Y) + 40 * np.sin(2 * X + Y)
A = A + 10 * np.random.randn(m, n)  # 叠加少量随机噪声，模拟真实图像
A = np.clip(A, 0, 255).astype(np.float64)
print(f"原始图像尺寸: {m} x {n}， 总像素数: {m * n}")

# 2. 对图像矩阵进行SVD分解
U, S, Vt = np.linalg.svd(A, full_matrices=False)
# 注意：np.linalg.svd返回的S是奇异值的一维向量，Vt = V^T
p = len(S)  # 奇异值个数，p = min(m, n)

# 3. 定义不同秩r进行压缩和重构
ranks = [1, 5, 20, 50]
fig, axes = plt.subplots(1, len(ranks) + 1, figsize=(15, 4))

# 显示原图
axes[0].imshow(A, cmap='gray')
axes[0].set_title(f'原图\n{m}x{n}')
axes[0].axis('off')

for idx, r in enumerate(ranks, start=1):
    # 使用前r个奇异值重构图像
    # A_r = U[:, :r] @ np.diag(S[:r]) @ Vt[:r, :]
    # 更高效的乘法方式：
    A_reconstructed = (U[:, :r] * S[:r]) @ Vt[:r, :]
    
    # 确保像素值在合理范围内
    A_reconstructed = np.clip(A_reconstructed, 0, 255)
    
    # 计算存储量
    original_size = m * n
    compressed_size = r * (m + n + 1)  # U_r: m*r, S_r: r, V_r^T: r*n
    compression_ratio = compressed_size / original_size
    
    # 显示重构图像
    axes[idx].imshow(A_reconstructed, cmap='gray')
    axes[idx].set_title(f'r={r}\n压缩比: {compression_ratio:.2%}')
    axes[idx].axis('off')

plt.tight_layout()
plt.show()

# 4. 打印奇异值衰减情况（前20个）
print("前20个奇异值:")
print(S[:20])
print(f"\n前10个奇异值之和 / 总和: {np.sum(S[:10]) / np.sum(S):.2%}")
print(f"前20个奇异值之和 / 总和: {np.sum(S[:20]) / np.sum(S):.2%}")

SVD与主成分分析（PCA）的关系

主成分分析（PCA）是一种经典的数据降维方法，其目标是为原始数据寻找一组新的正交基（主成分），使得数据在这些新方向上的方差最大。SVD是实现PCA最稳定、最常用的数值算法。

设定：假设我们有 $n$ 个样本，每个样本有 $p$ 个特征，构成数据矩阵 $X \in \mathbb{R}^{n \times p}$。通常我们假设 $X$ 已经中心化，即每个特征的均值为0。

PCA的目标：寻找一个投影矩阵 $W \in \mathbb{R}^{p \times k}$（$k < p$），其列向量是单位正交的（即 $W^T W = I_k$），使得投影后的数据 $T = X W \in \mathbb{R}^{n \times k}$ 的方差尽可能大。$T$ 的列称为主成分得分。

SVD的桥梁作用：对中心化后的数据矩阵 $X$ 进行“精简”SVD（即不计算完整的 $U$ 和 $V$）： $$X = U \Sigma V^T, \quad \text{其中 } U \in \mathbb{R}^{n \times p}, \Sigma \in \mathbb{R}^{p \times p}, V \in \mathbb{R}^{p \times p}$$ （这里假设 $n \ge p$，且 $\Sigma$ 是方阵；若 $n < p$，形式类似，核心结论不变）。

1. 协方差矩阵：样本协方差矩阵为 $\frac{1}{n-1} X^T X$。根据SVD，有： $$X^T X = V \Sigma^T U^T U \Sigma V^T = V \Sigma^2 V^T$$ 这正是 $X^T X$ 的特征值分解。$V$ 的列向量就是 $X^T X$ 的特征向量，即我们要求的主成分方向。因此，$W = V$（取前 $k$ 列）。

2. 主成分得分：将 $W = V$ 代入 $T = X W$，得到： $$T = X V = U \Sigma V^T V = U \Sigma$$ 因此，主成分得分矩阵 $T$ 就是左奇异向量矩阵 $U$ 与奇异值矩阵 $\Sigma$ 的乘积。

3. 方差：第 $i$ 个主成分的方差为 $\frac{1}{n-1} \sigma_i^2$。奇异值 $\sigma_i$ 的降序排列保证了第一主成分方差最大，第二主成分次之，依此类推。

结论：对数据矩阵 $X$ 进行SVD，可以直接得到PCA的所有结果：

• 右奇异矩阵 $V$ 的列就是主成分方向。
• $U \Sigma$ 的列就是主成分得分。
• $\sigma_i^2 / (n-1)$ 就是各主成分的方差。

因此，在实际计算中，我们通常直接对数据矩阵调用SVD算法来执行PCA，这比先计算协方差矩阵再求特征值分解更数值稳定、更高效。

📝 动手练一练

1. 手动验证SVD：给定矩阵 $A = \begin{bmatrix} 1 & 2 \ 0 & 1 \ 2 & 0 \end{bmatrix}$。 a) 计算 $A^TA$ 和 $AA^T$。 b) 分别计算这两个矩阵的特征值和单位特征向量。 c) 根据特征值平方根得到奇异值，并根据特征向量构造 $V$ 和 $U$，验证 $A = U \Sigma V^T$ 是否成立（允许符号差异）。

2. 低秩近似与误差：使用上一题中的矩阵 $A$，计算其秩为1的近似矩阵 $A_1$。计算原始矩阵 $A$ 与近似矩阵 $A_1$之间的Frobenius范数误差 $||A - A_1||_F$。然后，尝试解释为什么这个误差等于 $\sigma_2$（假设有2个奇异值 $\sigma_1, \sigma_2$）？

参考答案：

1. a) $A^TA = \begin{bmatrix} 5 & 2 \ 2 & 5 \end{bmatrix}$, $AA^T = \begin{bmatrix} 5 & 2 & 2 \ 2 & 1 & 0 \ 2 & 0 & 4 \end{bmatrix}$。 b) $A^TA$ 的特征值：$\lambda_1=7, \lambda_2=3$。对应单位特征向量可分别为 $\mathbf{v}_1 = \frac{1}{\sqrt{2}}[1, 1]^T$, $\mathbf{v}_2 = \frac{1}{\sqrt{2}}[1, -1]^T$。故 $V = \frac{1}{\sqrt{2}}\begin{bmatrix} 1 & 1 \ 1 & -1 \end{bmatrix}$，奇异值 $\sigma_1=\sqrt{7}, \sigma_2=\sqrt{3}$。 c) 根据 $A \mathbf{v}_i = \sigma_i \mathbf{u}_i$ 可计算 $\mathbf{u}_1, \mathbf{u}_2$，并通过正交化补全 $U$ 的第三列。最终可验证分解成立。
2. 秩1近似 $A_1 = \sigma_1 \mathbf{u}_1 \mathbf{v}_1^T$。误差 $||A - A_1||$F^2 = \sum{i=2}^{p} \sigma_i^2。对于本题 $p=2$，所以 $||A - A_1||_F = \sigma_2 = \sqrt{3}$。这是因为SVD的低秩近似在Frobenius范数意义下是最优的，舍弃的分量对应的奇异值平方和就是近似误差的平方。

本章小结

本节我们系统学习了奇异值分解（SVD）这一线性代数中的核心工具。

• 核心定义：任何矩阵 $A_{m \times n}$ 都可分解为 $A = U \Sigma V^T$，其中 $U, V$ 为正交矩阵，$\Sigma$ 为广义对角矩阵，其正对角线元素（奇异值）降序排列。
• 数学本质：奇异值是 $A^TA$ 或 $AA^T$ 特征值的平方根；左右奇异向量分别是这两个矩阵的特征向量。SVD将矩阵分解为一系列秩1矩阵的加权和。
• 关键应用：利用奇异值衰减快的特性进行低秩近似，是数据压缩和噪声抑制的数学基础。SVD也是实现主成分分析（PCA）的稳定算法。

行动清单

1. 动手验证：运行讲义中的图像SVD示例代码，调整截断秩 $r$，直观感受压缩比与图像质量的关系。
2. 概念串联：画一张关系图，将矩阵 $A$、$A^TA$/$AA^T$、特征值/特征向量、奇异值/奇异向量、低秩近似、PCA主成分/得分这些概念连接起来，厘清其内在逻辑。
3. 拓展思考：除了图像，SVD还广泛应用于推荐系统（协同过滤）、自然语言处理（潜在语义分析LSA）等领域。尝试选择一个你感兴趣的领域，搜索“SVD在[某领域]的应用”，了解其如何解决实际问题。

— 小象教研组