【统计学习方法死磕系列Ⅱ】主成分分析

一、引言：从高维到低维

在机器学习中，我们经常面对高维数据——一张 64x64 的人脸图片有 4096 维，一篇有 10000 个词汇的文档有 10000 维。高维数据不仅带来计算的”维度灾难”，还隐藏一个更深层的问题：数据的内在维度通常远低于观测维度。

想象你拍摄一个旋转的硬币——虽然每张照片是 1024x1024 像素（约 100 万维），但硬币的变化实际上只由两个因素决定：旋转角度和硬币的面值。也就是说，数据虽然存在于 100 万维空间，却几乎落在某个 2 维流形上。

主成分分析（Principal Component Analysis, PCA）的目标就是发现这种低维结构：找到一组新的坐标轴（主成分），使得数据在这些轴上的投影方差最大化，从而实现降维。

本文将介绍 PCA 算法的步骤、两种经典推导（最大方差与最小重构误差）、与 SVD 的关系，以及核 PCA 和特征脸等应用。

二、PCA 算法步骤（快速预览）

假设有 $m$ 条 $n$ 维数据。

将原始数据按列组成 $n$ 行 $m$ 列矩阵 $X$
将 $X$ 的每一行（代表一个属性字段）进行零均值化，即减去这一行的均值
求出协方差矩阵 $C = \frac{1}{m-1} X X^T$
求出协方差矩阵的特征值及对应的特征向量
将特征向量按对应特征值大小从上到下按行排列成矩阵，取前 $k$ 行组成矩阵 $P$
$Y = PX$ 即为降维到 $k$ 维后的数据

这六个步骤构成了 PCA 的标准流程。接下来我们将从数学原理出发，深入理解每一步的动机和推导。

三、PCA 的两种等价推导

PCA 有两种经典的推导角度，殊途同归：

3.1 推导一：最大方差视角

目标：找到一个投影方向 $w$（$|w| = 1$），使得数据在该方向上的投影方差最大化。

给定中心化后的数据 $X = [x_1, x_2, \dots, x_m] \in \mathbb{R}^{n \times m}$，其中 $\sum_i x_i = 0$。

数据在方向 $w$ 上的投影（标量）为 $w^T x_i$。投影方差为：

$$ \frac{1}{m-1} \sum_{i=1}^{m} (w^T x_i)^2 = \frac{1}{m-1} \sum_i w^T x_i x_i^T w = w^T C w $$

其中 $C = \frac{1}{m-1} \sum_i x_i x_i^T = \frac{1}{m-1} X X^T$ 是协方差矩阵。

因此，问题变为约束优化：

$$ \max_{w} w^T C w, \quad \text{s.t. } \|w\| = 1 $$

拉格朗日乘子法：

$$ \mathcal{L}(w, \lambda) = w^T C w - \lambda (w^T w - 1) $$

对 $w$ 求导并置零：

$$ \frac{\partial \mathcal{L}}{\partial w} = 2Cw - 2\lambda w = 0 \implies Cw = \lambda w $$

即 $w$ 必须是协方差矩阵 $C$ 的特征向量，$\lambda$ 是特征值。

代入目标函数：

$$ w^T C w = w^T (\lambda w) = \lambda w^T w = \lambda $$

因此，要最大化方差，就应选择最大特征值对应的特征向量作为第一主成分。

推广到 k 个主成分：第二主成分在与第一主成分正交的约束下最大化方差，依此类推。由 Rayleigh 商定理可知，前 $k$ 个主成分就是协方差矩阵最大的 $k$ 个特征值对应的特征向量。

3.2 推导二：最小重构误差视角

目标：找到一组正交基，使得用这组基的线性组合重构数据时，误差最小。

设 $U_k = [u_1, u_2, \dots, u_k] \in \mathbb{R}^{n \times k}$ 是 $k$ 个正交基向量（$U_k^T U_k = I$）。

每个数据点 $x_i$ 在低维空间的表示为 $z_i = U_k^T x_i$，重构为 $\hat{x}_i = U_k z_i = U_k U_k^T x_i$。

重构误差为：

$$ J = \frac{1}{m} \sum_{i=1}^{m} \|x_i - U_k U_k^T x_i\|^2 $$

展开并化简：

$$ J = \frac{1}{m} \sum_i \|x_i\|^2 - \frac{1}{m} \sum_i \|U_k^T x_i\|^2 $$

第一项是常数（数据的总能量），因此最小化 $J$ 等价于最大化投影坐标的模长之和：

$$ \max_{U_k} \sum_i \|U_k^T x_i\|^2 = \max_{U_k} \sum_i \text{tr}(U_k^T x_i x_i^T U_k) $$

也就是最大化 $\text{tr}(U_k^T C U_k)$。这与最大方差视角完全一致，解也是 $C$ 的前 $k$ 个最大特征值对应的特征向量。

3.3 两种视角的统一

视角	优化目标	数学形式	解
最大方差	最大化投影方差	$\max_w w^T C w$	$C$ 的最大特征值对应的特征向量
最小误差	最小化重构误差	$\min_{U_k} \sum \|x_i - U_k U_k^T x_i\|^2$	$C$ 的前 k 个最大特征值对应的特征向量

两种视角在数学上等价，因为它们本质上共同受限于同一个协方差矩阵 $C$ 的谱结构。这也符合直觉：最大化了投影方差，自然就最小化了重构误差（因为总方差 = 投影方差 + 重构误差）。

四、协方差矩阵与特征分解

4.1 协方差矩阵

协方差矩阵 $C$ 是 PCA 的核心。它是一个 $n \times n$ 的对称半正定矩阵：

$$ C = \frac{1}{m-1} \sum_{i=1}^{m} (x_i - \bar{x})(x_i - \bar{x})^T $$

$C$ 的对角线元素 $C_{jj}$ 是特征 $j$ 的方差，非对角线元素 $C_{ij}$ 是特征 $i$ 和 $j$ 的协方差。

4.2 特征分解

由于 $C$ 是对称半正定矩阵，由谱定理可知它可以对角化：

$$ C = U \Lambda U^T = \sum_{i=1}^{n} \lambda_i u_i u_i^T $$

其中：

$\Lambda = \text{diag}(\lambda_1, \lambda_2, \dots, \lambda_n)$，$\lambda_1 \geq \lambda_2 \geq \dots \geq \lambda_n \geq 0$
$U = [u_1, u_2, \dots, u_n]$ 是正交矩阵，$u_i$ 是 $C$ 的特征向量
$\lambda_i$ 是第 $i$ 个主成分的方差

4.3 方差解释率

第 $i$ 个主成分的方差贡献率：

$$ r_i = \frac{\lambda_i}{\sum_{j=1}^{n} \lambda_j} $$

前 $k$ 个主成分的累计方差贡献率：

$$ R_k = \frac{\sum_{i=1}^{k} \lambda_i}{\sum_{j=1}^{n} \lambda_j} $$

通常选择 $k$ 使 $R_k$ 达到某个阈值（如 95%），这保证了降维后的数据保留了原始数据 95% 的能量。

五、坡碎石图（Scree Plot）

坡碎石图是选择主成分数量 $k$ 的可视化工具。它将特征值从大到小绘制在图上，寻找”肘部”。

import numpy as np
import matplotlib.pyplot as plt

def scree_plot(eigenvalues, threshold=0.95):
    """Plot scree plot with cumulative variance"""
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
    
    # Left: Eigenvalue magnitude
    ax1.plot(range(1, len(eigenvalues) + 1), eigenvalues, 'bo-')
    ax1.set_xlabel('Principal Component')
    ax1.set_ylabel('Eigenvalue')
    ax1.set_title('Scree Plot')
    
    # Right: Cumulative variance
    cumsum = np.cumsum(eigenvalues) / np.sum(eigenvalues)
    ax2.plot(range(1, len(eigenvalues) + 1), cumsum, 'ro-')
    ax2.axhline(y=threshold, color='g', linestyle='--', 
                label=f'{threshold*100}% threshold')
    ax2.set_xlabel('Number of Components')
    ax2.set_ylabel('Cumulative Explained Variance')
    ax2.set_title('Cumulative Variance Explained')
    ax2.legend()
    
    plt.tight_layout()
    plt.show()

六、PCA 与 SVD 的关系

在实际计算中，通常不直接对协方差矩阵做特征分解，而是使用 SVD（奇异值分解），尤其在样本维度很高的情况下。

6.1 通过 SVD 计算 PCA

设中心化后的数据矩阵 $X \in \mathbb{R}^{n \times m}$（$n$ 为特征数，$m$ 为样本数）。

SVD 分解：

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

其中 $U \in \mathbb{R}^{n \times n}$ 是左奇异向量矩阵，$\Sigma \in \mathbb{R}^{n \times m}$ 是对角奇异值矩阵，$V \in \mathbb{R}^{m \times m}$ 是右奇异向量矩阵。

则协方差矩阵为：

$$ C = \frac{1}{m-1} X X^T = \frac{1}{m-1} U \Sigma V^T V \Sigma^T U^T = \frac{1}{m-1} U \Sigma^2 U^T $$

对比 $C = U \Lambda U^T$，可知：

左奇异向量 $U$ 就是主成分方向
奇异值的平方 $\sigma_i^2$ 正比于特征值：$\lambda_i = \frac{\sigma_i^2}{m-1}$

降维后的数据为：

$$ Z = U_k^T X = \Sigma_k V_k^T $$

6.2 SVD vs EVD 的数值比较

方面	协方差矩阵 EVD	SVD
矩阵大小	$n \times n$	$n \times m$
当 $n \gg m$ 时	计算 $n \times n$ 矩阵非常昂贵	只需计算 $U$ 的前 $k$ 列
数值稳定性	协方差矩阵可能病态	SVD 更稳定
适用性	小到中等维度	任意维度

实际建议：永远用 SVD 来实现 PCA。sklearn.decomposition.PCA 用的就是 SVD。

七、PCA 的完整 Python 实现

7.1 从零实现（使用 EVD）

import numpy as np

class PCA:
    def __init__(self, n_components):
        self.n_components = n_components
        self.components_ = None  # shape: (n_components, n_features)
        self.mean_ = None
        self.explained_variance_ = None
        self.explained_variance_ratio_ = None
    
    def fit(self, X):
        n_samples, n_features = X.shape
        
        # 1. Center the data
        self.mean_ = np.mean(X, axis=0)
        X_centered = X - self.mean_
        
        # 2. Compute covariance matrix
        C = (X_centered.T @ X_centered) / (n_samples - 1)
        
        # 3. Eigendecomposition
        eigenvalues, eigenvectors = np.linalg.eigh(C)
        
        # 4. Sort by eigenvalue descending
        idx = np.argsort(eigenvalues)[::-1]
        eigenvalues = eigenvalues[idx]
        eigenvectors = eigenvectors[:, idx]
        
        # 5. Keep top k components
        self.components_ = eigenvectors[:, :self.n_components].T
        self.explained_variance_ = eigenvalues[:self.n_components]
        total_var = np.sum(eigenvalues)
        self.explained_variance_ratio_ = eigenvalues[:self.n_components] / total_var
        
        return self
    
    def transform(self, X):
        X_centered = X - self.mean_
        return X_centered @ self.components_.T
    
    def fit_transform(self, X):
        self.fit(X)
        return self.transform(X)
    
    def inverse_transform(self, Z):
        return Z @ self.components_ + self.mean_

7.2 使用 sklearn

from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline

# Note: StandardScaler is recommended before PCA
# so features with large scales don't dominate
pipeline = Pipeline([
    ('scaler', StandardScaler()),
    ('pca', PCA(n_components=0.95))  # keep 95% variance
])

Z = pipeline.fit_transform(X)
print(f"Original dim: {X.shape[1]}, Reduced dim: {Z.shape[1]}")

7.3 为什么 PCA 前要标准化

PCA 对变量的尺度非常敏感。如果不标准化，方差大的变量（即使不包含有价值信息）会主导主成分。例如，将”年收入（0-1000000）”和”年龄（0-100）”一起做 PCA 而不标准化，第一个主成分将几乎完全由年收入主导。

标准化（Z-score：$x’ = \frac{x - \mu}{\sigma}$）使所有特征具有相同尺度，确保 PCA 的公平性。

7.4 何时不标准化

特征已经在同一量纲和相同尺度上（如都是像素值 $[0, 255]$）
特征的相对方差本身包含有意义的信号

八、核 PCA（Kernel PCA）

8.1 动机

标准 PCA 只能发现线性结构的主成分。当数据在原始空间中高度非线性时（如瑞士卷数据集），标准 PCA 失效。

核 PCA 通过核技巧将数据隐式映射到高维特征空间，在特征空间中做线性 PCA。

8.2 数学原理

设映射 $\phi: \mathbb{R}^n \to \mathcal{H}$ 将数据映射到特征空间。我们不直接计算 $\phi(x_i)$，而是使用核函数 $K(x_i, x_j) = \langle \phi(x_i), \phi(x_j) \rangle$。

核 PCA 的核心是：在特征空间中，PCA 等价于对中心化的核矩阵 $\tilde{K}$ 做特征分解：

$$ \tilde{K} = K - \mathbf{1}_m K - K \mathbf{1}_m + \mathbf{1}_m K \mathbf{1}_m $$

其中 $K_{ij} = K(x_i, x_j)$，$\mathbf{1}_m$ 是 $m \times m$ 的全 1 矩阵除以 $m$。

然后对 $\tilde{K}$ 特征分解：

$$ \tilde{K} v = \lambda v $$

归一化后（$|v| = 1/\sqrt{\lambda}$），新数据 $x$ 在第 $k$ 个主成分上的投影为：

$$ z_k = \sum_{i=1}^{m} v_{ki} K(x_i, x) $$

8.3 Python 实现

from sklearn.decomposition import KernelPCA

kpca = KernelPCA(n_components=2, kernel='rbf', gamma=0.1)
Z = kpca.fit_transform(X)

# Common kernels: 'linear', 'poly', 'rbf', 'sigmoid', 'cosine'

九、应用：特征脸（Eigenfaces）

Eigenfaces 是 PCA 在人脸识别中的经典应用（Turk & Pentland, 1991）。

9.1 思想

将每张 $h \times w$ 的人脸图像拉成 $d = hw$ 维的向量。用 PCA 找到人脸空间的主成分——即特征脸（eigenfaces）。每张人脸可以表示为若干张特征脸的线性组合。

9.2 实现

from sklearn.decomposition import PCA
from sklearn.datasets import fetch_lfw_people

# Load face data
lfw = fetch_lfw_people(min_faces_per_person=70, resize=0.4)
X = lfw.data  # (n_samples, h*w)

# PCA
n_components = 150
pca = PCA(n_components=n_components, whiten=True)
Z = pca.fit_transform(X)

# Eigenfaces = pca.components_
# Each component is of shape (h*w,), can be reshaped to (h, w)

# Reconstruction
X_reconstructed = pca.inverse_transform(Z)

# Recognition: 1-NN in PCA space
from sklearn.neighbors import KNeighborsClassifier
knn = KNeighborsClassifier(n_neighbors=1)
knn.fit(Z_train, y_train)
predictions = knn.predict(Z_test)

十、PCA 的局限性与变体

10.1 主要局限

局限	说明	替代方案
只捕获线性结构	非线性数据效果差	核 PCA, t-SNE, UMAP
对尺度敏感	特征尺度不一则偏斜	标准化
主成分难解释	每个 PC 是原始特征的线性组合	稀疏 PCA
假设高斯分布	数据非高斯时效果不一定最优	ICA
无监督	不利用标签信息	LDA（监督降维）

10.2 相关变体

方法	核心思想	与 PCA 的区别
概率 PCA（PPCA）	用概率模型框架解释 PCA	有生成模型、可处理缺失值
稀疏 PCA	加 L1 正则化使主成分稀疏	主成分更具可解释性
鲁棒 PCA（RPCA）	将数据分解为低秩 + 稀疏	对离群点鲁棒
ICA（独立分量分析）	最大化统计独立性	非高斯、比 PCA 更高阶
LDA（线性判别分析）	最大化类间/类内方差比	有监督降维

十一、面试高频问题

Q1：PCA 和 LDA 的核心区别是什么？何时用 PCA，何时用 LDA？

PCA 是无监督的降维方法，寻找方差最大的方向，不考虑标签。LDA 是有监督的降维方法，寻找使类间分离最大、类内分散最小的方向。

方面	PCA	LDA
监督	无监督	有监督
目标	最大化方差	最大化 Fisher 比率
最大降维维数	$\min(n, m)$	$c - 1$（$c$ 为类别数）
适用场景	数据探索、降噪、特征提取	分类任务的特征降维

Q2：为什么要对数据做中心化（减去均值）？

中心化是 PCA 的数学前提。PCA 寻找的是过原点的方向，协方差矩阵的定义隐含了中心化（$C = \mathbb{E}[(x-\mu)(x-\mu)^T]$）。如果不中心化，第一个主成分将主要指向数据中心而非方差最大方向。

事实上，PCA 的推导中通过拉格朗日乘子法得到 $Cw = \lambda w$，其中 $C$ 本身就以均值为中心定义，因此中心化是推导的自然结果，而非额外的技巧。

Q3：为什么用 SVD 而非 EVD 来计算 PCA？

效率：当 $n \gg m$（特征数远大于样本数）时，协方差矩阵 $C$ 是 $n \times n$ 的，EVD 计算量极高。而 SVD 可以只求 $X$ 的前 $m$ 个奇异向量，更高效。
数值稳定性：构造协方差矩阵 $X X^T$ 时可能丢失精度（条件数平方恶化）。SVD 直接在 $X$ 上操作，数值更稳定。
内存：当 $n$ 很大（如图像数据百万维），$n \times n$ 的协方差矩阵根本无法存储。SVD 可以通过随机化方法近似。

Q4：PCA 选择了 k 个主成分，怎样评估这个 k 是否合适？

累计方差贡献率：选 $k$ 使 $R_k \geq 95%$（或其他阈值）
坡碎石图：找特征值递减曲线的”肘部”
下游任务性能：用降维后的特征做分类/回归，选使任务性能最高的 $k$
重构误差：画 $\frac{|X - \hat{X}|^2}{|X|^2}$ vs $k$，选误差下降的拐点
交叉验证：在概率 PCA 框架下，可以用对数似然做模型选择

Q5：PCA 假设数据服从高斯分布吗？如果数据不服从高斯分布会怎样？

标准 PCA 不假设数据服从高斯分布——它仅是寻找方差最大的方向，对任何分布的数据都”能工作”。但在以下意义上，高斯假设是隐含的：

PCA 只使用了二阶统计量（协方差/方差），而高斯分布完全由均值和方差决定，所以 PCA 对于高斯数据是最优的降维方法
概率 PCA 的框架显式假设了高斯先验
如果数据高度非高斯（如存在簇结构），PCA 可能产生无意义的主成分——例如两个球形簇的 PCA 第一主成分可能会穿过两个簇之间的空白区域，而这在实际意义上是无用的

对于非高斯数据，ICA 或 t-SNE/UMAP 可能是更好的选择。

Q6：PCA 和自编码器（Autoencoder）有什么联系和区别？

这是一个经常被问到的问题，因为两者的目标非常相似——都是学习数据的低维表示。

联系：当自编码器使用线性激活函数且隐藏层无 bias 时，其学习到的权重矩阵恰好张成了与 PCA 主成分相同的子空间。换句话说，线性自编码器等价于 PCA（Baldi & Hornik, 1989）。具体而言，如果编码器为 $h = Wx$，解码器为 $\hat{x} = W’h$，用 MSE 损失训练，$W$ 的列向量将收敛到数据协方差矩阵的前 $k$ 个主成分方向。

区别：

方面	PCA	Autoencoder
映射	线性（$z = U_k^T x$）	可非线性（通过激活函数）
可解释性	高（主成分有明确的方差解释）	低（隐空间维度是黑箱）
优化	闭式解（特征分解）	梯度下降（非凸，可能局部最优）
参数数量	固定（$n \times k$）	可变（隐藏层可任意设计）
对非线性数据	差（只能找线性结构）	好（深度网络可学习非线性流形）

实践建议：如果数据是近似线性结构或需要可解释性，用 PCA。如果数据高度非线性且重点在于重建精度（如去噪自编码器用于图像去噪），用自编码器。有时可以先做 PCA 作为基线，再用自编码器看能否进一步降低重构误差。

十二、总结

主成分分析是数据降维的”瑞士军刀”——简单、优雅、高效。它通过协方差矩阵的特征分解（或等价地，数据矩阵的 SVD）来发现数据方差最大的正交方向，从而在保留最多信息的前提下降低数据维度。

PCA 的核心洞见是：在大量观测变量的表象之下，数据的内在结构往往只需要少数几个维度就能刻画。PCA 不仅是一种降维工具，更是一种数据探索的方法——通过特征脸（eigenfaces）、主成分载荷图等可视化手段，帮助我们理解数据的底层结构。

在实际工作中，建议在建模前先对数据跑一遍 PCA，用前两个主成分做散点图观察数据的聚集模式、离群点、线性可分性等。这往往能为后续的模型选择和特征工程提供宝贵的洞察。

参考文献：

Pearson, K. (1901). On lines and planes of closest fit to systems of points in space. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 2(11), 559-572.
Hotelling, H. (1933). Analysis of a complex of statistical variables into principal components. Journal of Educational Psychology, 24(6), 417.
Turk, M., & Pentland, A. (1991). Eigenfaces for recognition. Journal of Cognitive Neuroscience, 3(1), 71-86.
Scholkopf, B., Smola, A., & Muller, K. R. (1998). Nonlinear component analysis as a kernel eigenvalue problem. Neural Computation, 10(5), 1299-1319.
Tipping, M. E., & Bishop, C. M. (1999). Probabilistic principal component analysis. Journal of the Royal Statistical Society: Series B, 61(3), 611-622.