1. 概述

主成分分析(Principal Component Analysis, PCA)是最经典的降维方法之一,通过线性变换将高维数据投影到低维空间,同时尽可能保留原始数据的方差信息。1

1.1 为什么需要PCA?

问题PCA的解决方案
维度灾难将高维数据映射到低维子空间
冗余特征识别并去除高度相关的特征方向
噪声干扰保留主要方差方向,过滤噪声
可视化将高维数据投影到2-3维空间
预处理作为其他算法的降维步骤

1.2 PCA与其他降维方法的对比

import numpy as np
from sklearn.decomposition import PCA, KernelPCA, IncrementalPCA
from sklearn.manifold import TSNE
import matplotlib.pyplot as plt
 
# 生成非线性流形数据(瑞士卷)
from sklearn.datasets import make_swiss_roll
X, color = make_swiss_roll(n_samples=1000, noise=0.1, random_state=42)
 
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
 
# PCA(线性投影)
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X)
axes[0, 0].scatter(X_pca[:, 0], X_pca[:, 1], c=color, cmap='viridis', s=1)
axes[0, 0].set_title(f'PCA (explained var: {pca.explained_variance_ratio_.sum():.2%})')
 
# Kernel PCA(RBF核)
kpca = KernelPCA(n_components=2, kernel='rbf', gamma=0.1)
X_kpca = kpca.fit_transform(X)
axes[0, 1].scatter(X_kpca[:, 0], X_kpca[:, 1], c=color, cmap='viridis', s=1)
axes[0, 1].set_title('Kernel PCA (RBF)')
 
# t-SNE(非线性流形学习)
tsne = TSNE(n_components=2, random_state=42, perplexity=30)
X_tsne = tsne.fit_transform(X)
axes[1, 0].scatter(X_tsne[:, 0], X_tsne[:, 1], c=color, cmap='viridis', s=1)
axes[1, 0].set_title('t-SNE')
 
# Incremental PCA
ipca = IncrementalPCA(n_components=2)
X_ipca = ipca.fit_transform(X)
axes[1, 1].scatter(X_ipca[:, 0], X_ipca[:, 1], c=color, cmap='viridis', s=1)
axes[1, 1].set_title('Incremental PCA')

2. PCA的数学推导

PCA可以通过两种等价的视角推导:最大方差视角最小重构误差视角2

2.1 最大方差视角

假设数据已中心化(均值为0),我们的目标是找到一组正交基向量(主轴),使得数据在这些方向上的投影方差最大。

单主成分

对于第一个主成分 (满足 ),投影后的方差为:

其中 是数据的协方差矩阵。

最大化方差约束在单位范数下,使用拉格朗日乘子法:

求导并设为零:

这正是特征值方程!因此 是协方差矩阵 的特征向量,对应的方差为特征值

多主成分

第二个主成分 不仅要最大化方差,还要与 正交()。这同样会导致 的特征向量,且

结论:PCA的主轴就是协方差矩阵的特征向量,按对应特征值从大到小排序。

2.2 最小重构误差视角

另一种等价的推导是最小化重构误差。设原始数据 ,我们将其投影到 维子空间:

其中 是投影系数。

重构误差:

最小化此目标函数同样得到 的列为协方差矩阵 的前 个特征向量。

2.3 SVD与PCA的关系

设数据矩阵 个样本, 维特征),中心化后的SVD为:

其中:

  • :左奇异向量
  • :奇异值矩阵
  • :右奇异向量

关键关系

因此,PCA的主成分(特征向量)就是数据矩阵的右奇异向量,奇异值的平方是对应的特征值。

import numpy as np
 
def pca_svd(X, n_components=None):
    """
    使用SVD进行PCA
    
    Args:
        X: 中心化后的数据矩阵 (n_samples, n_features)
        n_components: 保留的主成分数量
    
    Returns:
        components: 主成分方向 (n_components, n_features)
        explained_variance: 各主成分解释的方差
        explained_variance_ratio: 方差解释比例
    """
    # SVD分解
    U, s, Vt = np.linalg.svd(X, full_matrices=False)
    
    # 确定主成分数量
    if n_components is None:
        n_components = len(s)
    
    # 主成分方向(V的行)
    components = Vt[:n_components]
    
    # 解释的方差(奇异值的平方 / (n-1))
    explained_variance = s[:n_components]**2 / (X.shape[0] - 1)
    explained_variance_ratio = explained_variance / explained_variance.sum()
    
    return components, explained_variance, explained_variance_ratio
 
# 示例
np.random.seed(42)
X = np.random.randn(100, 5) @ np.random.randn(5, 3)  # 降维生成
X_centered = X - X.mean(axis=0)
 
components, var, var_ratio = pca_svd(X_centered, n_components=3)
print(f"主成分方向:\n{components}")
print(f"解释方差比例: {var_ratio}")

3. 概率PCA(Probabilistic PCA)

概率PCA是PCA的生成模型版本,提供了更灵活的框架和贝叶斯解释。3

3.1 生成模型

假设潜在变量 服从标准正态分布:

观测变量 由潜在变量线性生成,并加上各向同性噪声:

其中:

  • :负载矩阵(loading matrix)
  • :均值向量
  • :噪声

因此,观测模型为:

3.2 边缘分布

对潜在变量积分,得到观测数据的边缘分布:

其中 是数据的协方差矩阵。

关键洞察:当 时,,此时概率PCA退化为经典PCA。

3.3 后验分布

给定观测 ,潜在变量的后验分布为:

其中

3.4 EM算法求解

对于概率PCA,通常使用EM算法高效估计参数

E步:计算潜在变量的后验期望

M步:更新参数

import numpy as np
 
class ProbabilisticPCA:
    """
    概率PCA实现
    
    使用EM算法估计参数
    """
    def __init__(self, n_components, tol=1e-4, max_iter=100):
        self.n_components = n_components
        self.tol = tol
        self.max_iter = max_iter
    
    def fit(self, X):
        n, D = X.shape
        K = self.n_components
        
        # 中心化数据
        self.mean = np.mean(X, axis=0)
        X_centered = X - self.mean
        
        # 初始化W为随机正交矩阵
        W = np.random.randn(D, K)
        W, _ = np.linalg.qr(W)  # QR分解确保正交
        
        sigma2 = 1.0
        
        for iteration in range(self.max_iter):
            # E步
            M = W.T @ W + sigma2 * np.eye(K)
            M_inv = np.linalg.inv(M)
            
            # 后验均值和二阶矩
            Ez = M_inv @ W.T @ X_centered.T  # (K, n)
            EzzT = n * sigma2 * M_inv + Ez @ Ez.T  # (K, K)
            
            # M步
            W_new = (X_centered.T @ Ez.T) @ np.linalg.inv(EzzT)
            
            # 更新sigma2
            residual = X_centered - (W_new @ Ez).T
            sigma2_new = np.sum(residual**2) / (n * D) + np.trace(sigma2 * W_new @ M_inv @ W_new.T) * (D / n)
            
            # 检查收敛
            change = np.abs(sigma2_new - sigma2)
            if change < self.tol:
                print(f"EM converged at iteration {iteration}")
                break
            
            W = W_new
            sigma2 = sigma2_new
        
        self.W = W
        self.sigma2 = sigma2
        self.M = W.T @ W + sigma2 * np.eye(K)
        
        return self
    
    def transform(self, X):
        """投影到潜在空间"""
        X_centered = X - self.mean
        return (np.linalg.inv(self.M) @ self.W.T @ X_centered.T).T
    
    def inverse_transform(self, z):
        """从潜在空间重构"""
        return (self.W @ z.T).T + self.mean
 
# 示例
np.random.seed(42)
X = np.random.randn(500, 10) @ np.random.randn(10, 5) + np.random.randn(500, 5) * 0.1
 
ppca = ProbabilisticPCA(n_components=3)
ppca.fit(X)
 
z = ppca.transform(X[:10])
X_reconstructed = ppca.inverse_transform(z)
 
print(f"重构误差: {np.mean((X[:10] - X_reconstructed)**2):.4f}")

3.5 与经典PCA的关系

方面经典PCA概率PCA
参数估计直接SVD分解EM算法
噪声处理隐式忽略显式建模
新样本需要重新计算可直接投影
贝叶斯扩展不支持自然支持贝叶斯推断
计算复杂度

4. 核PCA(Kernel PCA)

核PCA将核技巧应用到PCA,实现非线性降维。4

4.1 动机

经典PCA只能学习线性子空间,对于非线性流形结构(如瑞士卷、螺旋)效果较差。

4.2 核技巧

核PCA通过核函数 将数据映射到高维特征空间 ,然后在特征空间中进行PCA:

4.3 数学推导

是特征空间中的数据中心化数据矩阵( 可能无穷大)。核PCA要求解的特征方程:

其中 是核矩阵

利用 ,得到:

其中 是特征向量。

4.4 常见核函数

核函数公式参数
线性核
多项式核
RBF/高斯核
拉普拉斯核
Sigmoid核

4.5 代码实现

import numpy as np
 
class KernelPCA:
    """
    核PCA实现
    
    支持多种核函数
    """
    def __init__(self, n_components, kernel='rbf', gamma=None, degree=3, coef0=1):
        self.n_components = n_components
        self.kernel = kernel
        self.gamma = gamma
        self.degree = degree
        self.coef0 = coef0
    
    def _compute_kernel(self, X):
        """计算核矩阵"""
        n = X.shape[0]
        K = np.zeros((n, n))
        
        if self.kernel == 'linear':
            K = X @ X.T
        
        elif self.kernel == 'rbf':
            if self.gamma is None:
                self.gamma = 1.0 / X.shape[1]
            # ||x - y||^2 = ||x||^2 + ||y||^2 - 2 x^T y
            X_sqnorms = np.sum(X**2, axis=1)
            K = np.exp(-self.gamma * (
                X_sqnorms[:, np.newaxis] + X_sqnorms[np.newaxis, :] - 2 * (X @ X.T)
            ))
        
        elif self.kernel == 'poly':
            K = (self.gamma * X @ X.T + self.coef0) ** self.degree
        
        elif self.kernel == 'sigmoid':
            K = np.tanh(self.gamma * X @ X.T + self.coef0)
        
        return K
    
    def _center_kernel_matrix(self, K):
        """中心化核矩阵"""
        n = K.shape[0]
        # K_centered = K - 1_n K - K 1_n + 1_n K 1_n
        ones_n = np.ones((n, n)) / n
        return K - ones_n @ K - K @ ones_n + ones_n @ K @ ones_n
    
    def fit(self, X):
        self.X = X
        
        # 计算并中心化核矩阵
        K = self._compute_kernel(X)
        self.K_centered = self._center_kernel_matrix(K)
        
        # 特征分解
        eigenvalues, eigenvectors = np.linalg.eigh(self.K_centered)
        
        # 按特征值从大到小排序
        idx = np.argsort(eigenvalues)[::-1]
        eigenvalues = eigenvalues[idx]
        eigenvectors = eigenvectors[:, idx]
        
        # 取前n_components个
        self.eigenvalues = eigenvalues[:self.n_components]
        self.eigenvectors = eigenvectors[:, :self.n_components]
        
        # 归一化特征向量
        self.alpha = self.eigenvectors / np.sqrt(np.abs(self.eigenvalues) + 1e-10)
        
        return self
    
    def transform(self, X_new):
        """投影新样本到核PCA空间"""
        n_train = self.X.shape[0]
        n_new = X_new.shape[0]
        
        # 计算测试集与训练集的核矩阵
        K_test = np.zeros((n_new, n_train))
        
        if self.kernel == 'linear':
            K_test = X_new @ self.X.T
        
        elif self.kernel == 'rbf':
            X_train_sqnorms = np.sum(self.X**2, axis=1)
            X_new_sqnorms = np.sum(X_new**2, axis=1)
            K_test = np.exp(-self.gamma * (
                X_new_sqnorms[:, np.newaxis] + X_train_sqnorms[np.newaxis, :] 
                - 2 * (X_new @ self.X.T)
            ))
        
        # 中心化
        ones_train = np.ones((n_train, n_train)) / n_train
        ones_new_train = np.ones((n_new, n_train)) / n_train
        
        K_test_centered = K_test - ones_new_train @ self.K_centered @ ones_train @ ones_train
        
        # 投影
        return K_test_centered @ self.alpha
 
# 示例
np.random.seed(42)
# 生成同心圆数据(线性不可分)
angles = np.random.uniform(0, 2 * np.pi, 500)
r1 = np.random.normal(1, 0.1, 250)
r2 = np.random.normal(2, 0.1, 250)
 
X1 = np.column_stack([r1 * np.cos(angles[:250]), r1 * np.sin(angles[:250])])
X2 = np.column_stack([r2 * np.cos(angles[250:]), r2 * np.sin(angles[250:])])
X = np.vstack([X1, X2])
y = np.array([0] * 250 + [1] * 250)
 
# 线性PCA
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X)
 
# 核PCA
kpca = KernelPCA(n_components=2, kernel='rbf', gamma=0.5)
X_kpca = kpca.fit_transform(X)

4.6 优缺点分析

方面优点缺点
非线性能捕捉复杂流形结构计算复杂度高
核选择灵活适应不同数据结构需要调参
可解释性核函数可解释特征向量难以解释
新样本可直接投影需要存储核矩阵

5. 增量PCA(Incremental PCA)

增量PCA用于处理大规模数据或在线学习场景。5

5.1 动机

经典PCA需要将整个数据集加载到内存中,对于大规模数据或流式数据不适用。

5.2 随机SVD算法

使用随机投影加速SVD分解:

import numpy as np
 
class IncrementalPCA:
    """
    增量PCA实现
    
    基于随机SVD的在线学习方法
    """
    def __init__(self, n_components, batch_size=100, n_iter=5):
        self.n_components = n_components
        self.batch_size = batch_size
        self.n_iter = n_iter
    
    def fit(self, X):
        n, D = X.shape
        
        # 初始化
        self.mean = np.zeros(D)
        self.V = np.zeros((D, self.n_components))
        self.S_sq = np.zeros(self.n_components)  # 奇异值平方
        
        n_samples = 0
        
        # 分批处理
        for start in range(0, n, self.batch_size):
            end = min(start + self.batch_size, n)
            X_batch = X[start:end]
            
            # 更新均值
            batch_mean = X_batch.mean(axis=0)
            delta = batch_mean - self.mean
            self.mean += delta * (len(X_batch) / (n_samples + len(X_batch)))
            n_samples += len(X_batch)
            
            # 中心化批次
            X_batch_centered = X_batch - batch_mean
            
            # QR分解获取当前基
            Q, R = np.linalg.qr(np.column_stack([self.V, X_batch_centered.T]))
            
            # 随机投影
            for _ in range(self.n_iter):
                Z = Q @ np.random.randn(Q.shape[1], self.n_components + D)
                Z, _ = np.linalg.qr(Z)
                Z = Q.T @ Z
                Z, _ = np.linalg.qr(Z)
                Q = Q @ Z if Z.shape[0] == Q.shape[1] else Q @ Z.T
            
            # 更新V和S
            B = Q.T @ X_batch_centered.T
            U_b, S_b, Vt_b = np.linalg.svd(B, full_matrices=False)
            
            k = min(self.n_components, len(S_b))
            self.V = Q @ U_b[:, :k]
            self.S_sq[:k] = S_b[:k]**2
            
            Q = Q @ Vt_b[:k].T if k < Vt_b.shape[0] else Q
        
        # 计算解释方差比例
        total_var = np.sum(self.S_sq)
        self.explained_variance_ratio_ = self.S_sq / total_var
        
        return self
    
    def transform(self, X):
        """投影到主成分空间"""
        X_centered = X - self.mean
        return X_centered @ self.V
    
    def fit_transform(self, X):
        self.fit(X)
        return self.transform(X)
 
# 示例:处理大规模数据
np.random.seed(42)
# 模拟100万个样本
X_large = np.random.randn(10000, 100) @ np.random.randn(100, 20)
 
ipca = IncrementalPCA(n_components=10, batch_size=500)
X_reduced = ipca.fit_transform(X_large)
 
print(f"原始维度: {X_large.shape}")
print(f"降维后: {X_reduced.shape}")
print(f"解释方差比例: {ipca.explained_variance_ratio_}")

5.3 场景选择

场景推荐方法
小数据集(<10K)经典PCA(SVD分解)
中等数据集(10K-1M)增量PCA或随机SVD
超大规模(>1M)分布式PCA或随机投影
在线学习增量PCA

6. PCA白化变换

白化(Whitening)是一种重要的预处理步骤,使数据具有单位协方差矩阵。

6.1 白化的定义

白化变换 满足:

即各维度不相关且方差为1。

6.2 PCA白化

基于PCA的白化:

其中 是特征向量矩阵, 是特征值对角矩阵。

def pca_whitening(X, n_components=None):
    """
    PCA白化实现
    
    白化后的数据具有:
    1. 零均值
    2. 单位方差
    3. 各维度不相关
    """
    # 中心化
    X_centered = X - X.mean(axis=0)
    
    # SVD分解
    U, s, Vt = np.linalg.svd(X_centered, full_matrices=False)
    
    if n_components is None:
        n_components = len(s)
    
    # 白化矩阵
    whitening_matrix = Vt[:n_components].T / np.sqrt(s[:n_components] + 1e-10)
    
    # 白化变换
    X_white = X_centered @ whitening_matrix
    
    return X_white, whitening_matrix, Vt[:n_components], s[:n_components]
 
# 示例
np.random.seed(42)
X = np.random.randn(500, 10)
 
X_white, W, V, s = pca_whitening(X, n_components=5)
 
print(f"白化后协方差矩阵:\n{np.cov(X_white.T):.4f}")
print(f"(对角元素接近1,其他元素接近0)")

6.3 ZCA白化

ZCA(Zero-phase Component Analysis)白化保持数据尽可能接近原始形式:

def zca_whitening(X, n_components=None, epsilon=1e-10):
    """
    ZCA白化实现
    
    特点:最小化与原始数据的重构误差
    """
    # 中心化
    X_centered = X - X.mean(axis=0)
    
    # 协方差矩阵的特征分解
    cov = np.cov(X_centered.T)
    eigenvalues, eigenvectors = np.linalg.eigh(cov)
    
    # 排序
    idx = np.argsort(eigenvalues)[::-1]
    eigenvalues = eigenvalues[idx]
    eigenvectors = eigenvectors[:, idx]
    
    if n_components is None:
        n_components = len(eigenvalues)
    
    # ZCA白化矩阵
    whitening_matrix = eigenvectors[:, :n_components] @ \
                       np.diag(1.0 / np.sqrt(eigenvalues[:n_components] + epsilon)) @ \
                       eigenvectors[:, :n_components].T
    
    X_white = X_centered @ whitening_matrix
    
    return X_white, whitening_matrix
 
# 可视化比较
import matplotlib.pyplot as plt
 
# 生成相关数据
np.random.seed(42)
X_correlated = np.random.randn(500, 2)
X_correlated[:, 1] = 0.8 * X_correlated[:, 0] + 0.2 * np.random.randn(500)
 
# PCA白化
X_pca_white, _ = pca_whitening(X_correlated)
 
# ZCA白化
X_zca_white, _ = zca_whitening(X_correlated)

6.4 在深度学习中的应用

应用作用
输入预处理标准化协方差,加速训练
白化神经元L2正则化的非线性版本
实例归一化类似ZCA的空间白化
对比学习协方差匹配作为正则化

7. 实践指南

7.1 何时使用PCA

适合的场景

  • 数据维度高于样本数
  • 特征之间存在强相关性
  • 作为下游算法的预处理步骤
  • 数据可视化

不适合的场景

  • 数据本身是离散的(如类别数据)
  • 需要保留可解释性
  • 数据方差主要分布在非线性流形上

7.2 选择主成分数量

def select_n_components(X, variance_threshold=0.95):
    """
    基于累积解释方差选择主成分数量
    """
    pca = PCA()
    pca.fit(X)
    
    cumsum = np.cumsum(pca.explained_variance_ratio_)
    n_components = np.searchsorted(cumsum, variance_threshold) + 1
    
    return n_components, cumsum
 
# 可视化解释方差
import matplotlib.pyplot as plt
 
plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
plt.plot(cumsum, 'bo-')
plt.axhline(y=0.95, color='r', linestyle='--')
plt.xlabel('Number of Components')
plt.ylabel('Cumulative Explained Variance')
plt.title('Scree Plot')
 
plt.subplot(1, 2, 2)
plt.bar(range(1, len(pca.explained_variance_ratio_) + 1), 
        pca.explained_variance_ratio_)
plt.xlabel('Principal Component')
plt.ylabel('Explained Variance Ratio')
plt.title('Individual Explained Variance')

7.3 常见陷阱

# 陷阱1:忘记中心化
# 错误做法
U, s, Vt = np.linalg.svd(X)  # X未中心化
 
# 正确做法
X_centered = X - X.mean(axis=0)
U, s, Vt = np.linalg.svd(X_centered)
 
# 陷阱2:对非高斯数据直接PCA
# 如果数据的分布不是椭球形,PCA可能不是最优的
# 考虑先进行非线性变换(如Box-Cox)
 
# 陷阱3:混淆训练集和测试集
# 正确做法:只在训练集上fit,transform测试集
pca = PCA(n_components=10)
pca.fit(X_train)  # 只在训练集上拟合
X_train_pca = pca.transform(X_train)
X_test_pca = pca.transform(X_test)  # 用相同的变换

8. 总结

PCA家族对比

方法线性/非线性概率解释计算复杂度适用场景
经典PCA线性通用降维
概率PCA线性生成模型
核PCA非线性非线性流形
增量PCA线性大规模数据
ZCA白化线性图像预处理

核心公式汇总

目标公式
协方差矩阵
PCA投影
重构
重构误差
PCA白化

参考资料


相关阅读

Footnotes

  1. Jolliffe, I. T. (2002). Principal Component Analysis (2nd ed.). Springer.

  2. Shlens, J. (2014). A Tutorial on Principal Component Analysis. arXiv:1404.1100.

  3. Tipping, M. E., & Bishop, C. M. (1999). Probabilistic Principal Component Analysis. Journal of the Royal Statistical Society: Series B, 61(3), 611-622.

  4. Schölkopf, B., Smola, A., & Müller, K. R. (1998). Nonlinear Component Analysis as a Kernel Eigenvalue Problem. Neural Computation, 10(5), 1299-1319.

  5. Weng, J., & Ahuja, N. (2000). Incremental Principal Component Analysis. Proceedings of the International Conference on Image Processing.