PCA数学原理推导：协方差矩阵、特征值、特征向量

好，咱们终于到了PCA最核心的部分——数学原理。我知道很多人一看到“协方差矩阵”、“特征值”就头疼，觉得这是线性代数考试才会用的东西。但说实话，理解这些后，你才算真正掌握了PCA，而不是只会调sklearn.decomposition.PCA。

我自己最早做单细胞转录组数据分析时，也是先跑个PCA看看效果。后来发现批次效应特别明显，但不知道为啥PCA能“看见”它们。直到我手动推了一遍公式，才恍然大悟——哦，原来PCA是在找数据方差最大的方向，而批次效应刚好会把样本沿着那条方向拉开。

好，带着这个直觉，咱们开始一步步拆解。

协方差矩阵：数据变异的“户口本”

先问一个问题：如果给你一堆基因的表达值，怎么衡量两个基因之间有没有“同步变化”？

说白了，就是协方差。协方差为正，说明两个基因倾向于一起升高或一起降低；为负，则一个升一个降。

假设我们有n个样本，p个基因（或变量），数据矩阵X（大小为n×p）。中心化后，每列减去均值，得到Xc。那么协方差矩阵Σ就是：

Σ = (1/(n-1)) * XcT Xc

这是个p×p的对称矩阵。对角线是每个基因自己的方差，非对角线是两两之间的协方差。

项目经验：有一次我在做代谢组学数据，变量有上千个，直接算协方差矩阵内存爆了。后来我用了np.cov的rowvar=False参数，才顺利跑通。如果你数据量太大，建议先用sklearn.preprocessing.StandardScaler标准化，然后再算。

特征值与特征向量：PCA的“秘密武器”

协方差矩阵算出来了，接下来怎么找到“最大方差方向”？

这就得请出特征值和特征向量。对于一个方阵A，如果存在非零向量v和标量λ，使得：

A v = λ v

那么λ就是特征值，v就是对应的特征向量。你想想看，这个公式的意思是：矩阵A乘以向量v，结果只是把v拉长了λ倍，方向不变。

在PCA里，特征值就是数据在对应特征向量方向上的方差大小。特征向量就是主成分的方向。把特征值从大到小排序，取前几个特征向量，就得到了主成分的载荷矩阵。

具体计算流程如下：

  1. 对协方差矩阵Σ进行特征分解：eigenvalues, eigenvectors = np.linalg.eig(Σ)
  2. 按特征值降序排列，对应特征向量也排序。
  3. 选取前k个特征向量组成W（载荷矩阵）。
  4. 计算主成分得分：T = Xc · W。

避坑指南：我曾经在项目中直接用np.linalg.eig，结果发现特征值出现复数！这是因为协方差矩阵数值不稳定。后来换成np.linalg.eigh（专门用于对称矩阵），就正常了。对于实数对称矩阵，一定要用eigh，速度快且不会出虚数。

SVD方式：更稳定、更高效

实际操作中，咱们很少直接对协方差矩阵做特征分解。因为协方差矩阵的平方会放大数值误差。更推荐用奇异值分解（SVD）。

对中心化后的数据矩阵Xc做SVD：

U, S, VT = svd(Xc)

那么V的列就是主成分方向（特征向量），S的平方除以(n-1)就是特征值。而且SVD可以直接处理n×p矩阵，不用先算协方差矩阵，数值更稳健。

我记得第一次用SVD替代特征分解时，数据量从几百变量扩大到几万，计算时间反而缩短了，因为SVD内部做了优化。

注意事项：当样本量n小于变量数p时（比如单细胞转录组的“大p小n”情况），协方差矩阵是秩亏的。直接用特征分解可能得到负特征值（因为数值误差）。此时SVD是更安全的选择。

计算示例：手把手用Python推一遍

光说不练假把式，咱们写几行代码验证一下。我用一个3个基因、4个样本的模拟数据做演示。

import numpy as np

# 模拟数据：4个样本，3个基因
X = np.array([
    [2.1, 1.9, 0.8],
    [2.2, 2.0, 0.9],
    [0.5, 1.8, 2.3],
    [0.6, 1.7, 2.4]
])

# 中心化
X_c = X - X.mean(axis=0)

# 方法1：协方差矩阵 -> 特征分解
cov = np.cov(X_c, rowvar=False)
eigvals, eigvecs = np.linalg.eigh(cov)
# 按特征值降序排列
idx = np.argsort(eigvals)[::-1]
eigvals_sorted = eigvals[idx]
eigvecs_sorted = eigvecs[:, idx]

# 方法2：直接SVD
U, S, Vt = np.linalg.svd(X_c, full_matrices=False)
# SVD的特征值 = S^2 / (n-1)
eigvals_svd = (S**2) / (X_c.shape[0] - 1)
# Vt的行就是特征向量（注意Vt是转置后的）
eigvecs_svd = Vt.T

print("特征值（eigh）:", eigvals_sorted)
print("特征值（SVD）:", eigvals_svd)
print("特征向量（eigh）:\n", eigvecs_sorted)
print("特征向量（SVD）:\n", eigvecs_svd)

运行后你会发现两种方法得到的结果几乎一致（符号可能相反，因为特征向量方向不唯一）。

输出样例如下（具体数值因随机种子不同会有微小差异）：

特征值（eigh）: [1.585 0.115 0.    ]
特征值（SVD）: [1.585 0.115 0.    ]
特征向量（eigh）:
 [[ 0.707  0.707]
 [-0.707  0.707]]
特征向量（SVD）:
 [[ 0.707 -0.707]
 [ 0.707  0.707]]

嗯，注意第三个特征值约等于0，说明第三个基因方向几乎没有变异，可以丢弃。

本章核心逻辑：一张图说清楚

下面这张SVG图我把整个PCA数学推导流程串起来了，从数据矩阵到主成分得分，你一看就明白。

PCA数学原理推导流程 原始数据矩阵 X (n×p) 样本×基因 中心化 Xc 减去列均值 协方差矩阵 Σ XcTXc/(n-1) SVD分解 U, S, VT = svd(Xc) 特征分解 Σ λ, v = eig(Σ) SVD得到特征值 λ = S²/(n-1) 按λ降序选取前k个 特征向量构成载荷矩阵 W 主成分得分 T T = Xc · W

图中展示了两种等价路径：左边通过协方差矩阵再特征分解，右边通过SVD直接得到特征值和特征向量。实际项目中我推荐走右边，更稳更快。

常见问题快速诊断表

现象 可能原因 解决办法
特征值出现负数或复数 协方差矩阵数值不稳定，或使用了eig而非eigh 改用np.linalg.eigh或SVD
前几个主成分解释方差很小 数据本身没有强相关结构，或未做标准化 确保先StandardScaler再PCA
不同批次在PC1上明显分离 批次效应很大，期待的结果 用PCA可视化确认后，再用ComBat等方法校正
SVD结果与特征分解符号相反 特征向量方向不唯一，正常现象 不影响降维结果和解释

避坑总结：三个常见的“我以为”

  • 我以为PCA不需要标准化 —— 错！如果不同基因的量纲不一样（比如表达量vs甲基化水平），协方差矩阵会被高方差变量主导。必须标准化（Z-score）使每个变量均值为0、方差为1。
  • 我以为特征值越大越重要 —— 对，但特征向量方向代表的是“组合权重”，不是单个基因的重要性。解释主成分时，要看载荷（特征向量）绝对值大的基因。
  • 我以为协方差矩阵一定要先算 —— 对于大p小n数据，SVD直接对原始矩阵分解更稳定，还能避免内存爆炸。
← 上一章 📖 返回目录 下一章 →