主成分分析(PCA)与奇异值分解(SVD)的等价关系

你需要在数据降维时选择正确的工具。PCA和SVD是两个最常用的方法，很多人以为它们是不同的东西，但实际上它们可以互相转换。本文直接手把手教你理解并验证两者之间的等价关系，并用SVD实现PCA。

1. 理解主成分分析的核心逻辑

• 主成分分析的目标是将原始数据（可能有几十个特征）压缩成少数几个“综合指标”，同时尽量保留数据的波动信息。
• 核心操作：找到数据方差最大的方向（第一主成分），然后找到与它正交且方差次大的方向（第二主成分），以此类推。
• 关键步骤：计算数据的协方差矩阵（一个衡量各个特征之间如何一起变化的表格），然后对这个矩阵做特征分解。

2. 理解奇异值分解的核心逻辑

• 奇异值分解将一个数值矩阵 $A$ 分解成三个矩阵的乘积：$A = U \Sigma V^T$。
  • $U$ 和 $V$ 是正交矩阵（相当于旋转或反射），$\Sigma$ 是对角矩阵（对角线上是奇异值，代表各个方向上的缩放大小）。
• 关键点：$V$ 的列向量（右奇异向量）正好是 $A^T A$ 的特征向量，而 $\Sigma^2 / (n-1)$ 的对角线值等于协方差矩阵的特征值（当 $A$ 是中心化后的数据矩阵时）。

3. 建立等价关系（推导核心）

假设你有一个数据集，有 $n$ 个样本（行）和 $p$ 个特征（列）。准备数据矩阵 $X$（$n \times p$），已经对每一列做了中心化（减去该列的均值）。
现在：

1. 计算协方差矩阵：$C = \frac{1}{n-1} X^T X$（一个 $p \times p$ 的对称矩阵）。
2. 对 $X$ 做SVD：$X = U \Sigma V^T$。
   那么 $X^T X = V \Sigma^T U^T U \Sigma V^T = V \Sigma^2 V^T$（因为 $U$ 正交，$U^T U = I$）。
3. 对比特征分解：$C = \frac{1}{n-1} V \Sigma^2 V^T$。
   这说明 $V$ 的列向量就是协方差矩阵 $C$ 的特征向量（主成分方向）。而特征值 $\lambda_i = \frac{\sigma_i^2}{n-1}$，其中 $\sigma_i$ 是第 $i$ 个奇异值。

结论：对中心化后的数据矩阵做SVD，直接就能得到PCA的所有结果——不需要单独算协方差矩阵，也不需要做特征分解。

4. 实操步骤：用SVD计算PCA（以Python为例）

假设你有一个NumPy数组 X（每一行是一个样本，每一列是一个特征）。按照以下步骤实现PCA：

第一步：中心化数据

计算每一列的均值，减去均值。

import numpy as np

# 假设 X 是 (n, p) 的数组
mean = np.mean(X, axis=0)
X_centered = X - mean

第二步：执行SVD

使用NumPy的 linalg.svd 函数。

U, S, Vt = np.linalg.svd(X_centered, full_matrices=False)
# Vt 是 V 的转置，所以 V = Vt.T
V = Vt.T

• S 是一维数组，包含所有奇异值（按降序排列）。V 的列对应主成分方向（加载量）。

第三步：计算主成分得分

主成分得分是原始数据投影到主成分方向上的结果。计算 T = X_centered @ V（矩阵乘法）。

scores = X_centered @ V

• scores 的每一列是一个主成分，第一列是第一主成分得分。

第四步：获取解释方差比例

计算每个特征值：eigenvalues = S**2 / (n-1)，然后计算每个主成分的方差占比：explained_variance_ratio = eigenvalues / np.sum(eigenvalues)。

n = X_centered.shape[0]
eigenvalues = S**2 / (n - 1)
explained_variance_ratio = eigenvalues / np.sum(eigenvalues)

• explained_variance_ratio 告诉你前k个主成分能保留多少信息。

第五步：降维

保留前k个主成分：取 scores 的前k列，或者用 V[:, :k] 转换数据。

k = 2  # 保留2个主成分
X_reduced = X_centered @ V[:, :k]

5. 验证等价性（纯理论对照）

为了确认SVD结果与PCA完全一致，你可以做一次传统PCA（协方差矩阵特征分解）并比较：

• 传统PCA：C = (X_centered.T @ X_centered) / (n-1)，计算 eig_vals, eig_vecs = np.linalg.eigh(C)，注意eigh返回升序，需翻转。
• 对比：eig_vecs 的列（按特征值降序排列）应该与 V 的列只差一个符号（方向相反不影响PCA结果）。eig_vals 应该等于 S**2/(n-1)。

# 传统PCA（验证用）
C = (X_centered.T @ X_centered) / (n - 1)
eig_vals, eig_vecs = np.linalg.eigh(C)
# 降序排列
idx = np.argsort(eig_vals)[::-1]
eig_vecs = eig_vecs[:, idx]
eig_vals = eig_vals[idx]

# 验证特征向量方向（忽略符号）
print(np.abs(eig_vecs[:, :k]) - np.abs(V[:, :k]))
# 验证特征值
print(eig_vals[:k] - eigenvalues[:k])

如果两个矩阵相减的结果接近于 0，则等价关系成立。

6. 实际应用中的选择建议

• 当样本数远大于特征数（n >> p）时：直接对协方差矩阵做特征分解速度更快，因为只需要处理小矩阵 $p \times p$。此时PCA更经济。
• 当特征数远大于样本数（p >> n）时：SVD更稳定，因为它直接处理 $n \times p$ 的稀疏结构，而且不需要计算大型协方差矩阵。
• 当你已经用SVD做了其他分析时：可以直接从SVD结果中提取PCA结果，避免重复计算。

记住：无论你选择哪种方法，它们给出的主成分方向、得分、方差解释比例完全一致。两种方法本质上是同一个数学事实的不同表达方式。

7. 完整代码示例（一键执行）

import numpy as np

def pca_via_svd(X, n_components):
    # 中心化
    mean = np.mean(X, axis=0)
    X_centered = X - mean
    # SVD
    U, S, Vt = np.linalg.svd(X_centered, full_matrices=False)
    V = Vt.T
    # 得分
    scores = X_centered @ V
    # 解释方差比例
    n = X_centered.shape[0]
    eigenvalues = S**2 / (n - 1)
    explained_variance_ratio = eigenvalues / np.sum(eigenvalues)
    # 降维
    X_reduced = scores[:, :n_components]
    return X_reduced, explained_variance_ratio, V[:, :n_components]

# 使用示例
np.random.seed(42)
X_demo = np.random.randn(100, 5)  # 100个样本，5个特征
X_reduced, var_ratio, pcs = pca_via_svd(X_demo, 2)
print("前两个主成分解释方差比例:", var_ratio[:2].round(3))

运行上面的代码，你会得到与 sklearn.decomposition.PCA 完全一致的结果。至此，你已经掌握了PCA与SVD等价关系的核心以及实战用法。