PCA实战：用Python从零实现主成分分析（含特征值分解vs SVD对比）

PCA实战：用Python从零实现主成分分析（含特征值分解vs SVD对比）

主成分分析（PCA）是数据科学中最常用的降维技术之一，但很多初学者往往停留在调包使用的层面。本文将带你从零开始实现PCA算法，深入理解其数学本质，并对比特征值分解与奇异值分解（SVD）两种实现方式的性能差异和适用场景。

1. 环境准备与数据加载

在开始之前，我们需要准备Python环境和示例数据集。这里我们使用经典的Iris数据集作为演示，它包含150个样本，每个样本有4个特征（花萼长度、花萼宽度、花瓣长度、花瓣宽度）。

import numpy as np
from sklearn.datasets import load_iris
import matplotlib.pyplot as plt

# 加载Iris数据集
iris = load_iris()
X = iris.data
y = iris.target
feature_names = iris.feature_names

# 数据标准化
X_standardized = (X - np.mean(X, axis=0)) / np.std(X, axis=0)

为什么需要标准化？ PCA对数据的尺度非常敏感。如果特征的单位不同（比如一个特征以厘米为单位，另一个以米为单位），尺度较大的特征会主导主成分方向。标准化确保每个特征对结果的贡献是公平的。

2. PCA的数学基础

PCA的核心思想是通过线性变换将原始数据投影到新的坐标系中，使得：

1. 第一个主成分方向是数据方差最大的方向
2. 后续每个主成分都与前面的主成分正交（不相关）
3. 主成分按方差大小降序排列

2.1 协方差矩阵与特征值分解

PCA可以通过协方差矩阵的特征值分解来实现：

# 计算协方差矩阵
cov_matrix = np.cov(X_standardized.T)

# 特征值分解
eigen_values, eigen_vectors = np.linalg.eig(cov_matrix)

# 按特征值降序排列
sorted_index = np.argsort(eigen_values)[::-1]
sorted_eigenvalues = eigen_values[sorted_index]
sorted_eigenvectors = eigen_vectors[:, sorted_index]

print("特征值（方差贡献）:", sorted_eigenvalues)
print("解释方差比例:", sorted_eigenvalues / np.sum(sorted_eigenvalues))

2.2 主成分的几何解释

主成分可以理解为数据分布的最佳拟合直线（在更高维中是超平面）。第一个主成分方向是数据在该方向上投影方差最大的方向，第二个主成分是与第一个正交且剩余方差最大的方向，以此类推。

3. 两种PCA实现方式对比

3.1 特征值分解实现

基于特征值分解的PCA实现步骤如下：

def pca_eig(X, n_components=2):
    # 1. 标准化数据（已提前完成）
    # 2. 计算协方差矩阵
    cov_matrix = np.cov(X.T)
    
    # 3. 特征值分解
    eigen_values, eigen_vectors = np.linalg.eig(cov_matrix)
    
    # 4. 排序特征值和特征向量
    sorted_index = np.argsort(eigen_values)[::-1]
    sorted_eigenvectors = eigen_vectors[:, sorted_index]
    
    # 5. 选择前n个主成分
    components = sorted_eigenvectors[:, :n_components]
    
    # 6. 投影数据
    projected_data = X.dot(components)
    
    return projected_data, components

3.2 SVD实现

奇异值分解（SVD）提供了另一种计算PCA的方法：

def pca_svd(X, n_components=2):
    # 1. 标准化数据（已提前完成）
    # 2. 计算SVD
    U, S, Vt = np.linalg.svd(X)
    
    # 3. 选择前n个主成分
    components = Vt.T[:, :n_components]
    
    # 4. 投影数据
    projected_data = X.dot(components)
    
    return projected_data, components

3.3 两种方法的对比

实际建议：在大多数情况下，SVD是更好的选择，特别是当特征维度很高时。Scikit-learn的PCA实现默认就使用SVD方法。

4. 结果可视化与分析

让我们将Iris数据投影到前两个主成分并可视化：

# 使用特征值分解实现
projected_eig, components_eig = pca_eig(X_standardized)

# 使用SVD实现
projected_svd, components_svd = pca_svd(X_standardized)

# 可视化
plt.figure(figsize=(12, 5))

plt.subplot(1, 2, 1)
plt.scatter(projected_eig[:, 0], projected_eig[:, 1], c=y)
plt.title('PCA via Eigen Decomposition')
plt.xlabel('Principal Component 1')
plt.ylabel('Principal Component 2')

plt.subplot(1, 2, 2)
plt.scatter(projected_svd[:, 0], projected_svd[:, 1], c=y)
plt.title('PCA via SVD')
plt.xlabel('Principal Component 1')
plt.ylabel('Principal Component 2')

plt.tight_layout()
plt.show()

从可视化结果可以看出，两种方法得到的投影结果几乎相同（可能符号相反，但这不影响分析）。前两个主成分已经能够很好地分离三个类别的鸢尾花。

4.1 方差解释率分析

# 计算方差解释率
total_var = np.sum(sorted_eigenvalues)
explained_var_ratio = sorted_eigenvalues / total_var
cumulative_var_ratio = np.cumsum(explained_var_ratio)

plt.figure(figsize=(8, 5))
plt.bar(range(len(explained_var_ratio)), explained_var_ratio, alpha=0.5, 
        align='center', label='Individual explained variance')
plt.step(range(len(cumulative_var_ratio)), cumulative_var_ratio, 
         where='mid', label='Cumulative explained variance')
plt.ylabel('Explained variance ratio')
plt.xlabel('Principal components')
plt.legend(loc='best')
plt.title('Scree Plot')
plt.show()

方差解释率图（Scree Plot）帮助我们决定保留多少主成分。通常我们会选择累计解释方差达到一定阈值（如95%）的主成分数量。

5. 高级话题与实用技巧

5.1 核PCA简介

当数据不是线性可分时，标准PCA可能效果不佳。核PCA通过核技巧将数据映射到高维空间再进行PCA：

from sklearn.decomposition import KernelPCA

kpca = KernelPCA(n_components=2, kernel='rbf', gamma=15)
X_kpca = kpca.fit_transform(X_standardized)

plt.scatter(X_kpca[:, 0], X_kpca[:, 1], c=y)
plt.title('Kernel PCA (RBF Kernel)')
plt.show()

5.2 PCA在特征工程中的应用

PCA不仅用于可视化，在特征工程中也很有价值：

1. 降噪：去除方差小的成分可能去除噪声
2. 特征压缩：减少特征数量，提高模型效率
3. 去除相关性：主成分之间互不相关，适合某些模型

5.3 常见陷阱与注意事项

1. 分类问题中的PCA：在监督学习中，应该只在训练集上拟合PCA，然后转换测试集
2. 异常值影响：PCA对异常值敏感，可能需要先进行异常值处理
3. 解释性：主成分是原始特征的线性组合，可能难以解释
4. 信息损失：降维必然导致信息损失，需要权衡

6. 性能优化与大规模数据处理

当处理大规模数据时，可以考虑以下优化策略：

1. 增量PCA：适用于无法一次性加载到内存的大型数据集

   from sklearn.decomposition import IncrementalPCA
   ipca = IncrementalPCA(n_components=2, batch_size=10)
   X_ipca = ipca.fit_transform(X_standardized)
   

2. 随机化SVD：对于非常大的矩阵，可以使用随机算法近似计算SVD

   from sklearn.utils.extmath import randomized_svd
   U, S, Vt = randomized_svd(X_standardized, n_components=2)
   

3. 稀疏PCA：当数据稀疏时，可以使用稀疏PCA获得更稀疏的主成分

7. 实际案例：MNIST手写数字降维

让我们在一个更复杂的数据集上应用PCA——MNIST手写数字数据集：

from sklearn.datasets import fetch_openml
mnist = fetch_openml('mnist_784', version=1)
X_mnist = mnist.data / 255.0  # 归一化到[0,1]
y_mnist = mnist.target.astype(int)

# 使用SVD进行PCA
pca = PCA(n_components=2)
X_mnist_pca = pca.fit_transform(X_mnist)

plt.figure(figsize=(10, 8))
scatter = plt.scatter(X_mnist_pca[:, 0], X_mnist_pca[:, 1], c=y_mnist, 
                     cmap='tab10', alpha=0.5)
plt.colorbar(scatter)
plt.title('MNIST PCA Projection')
plt.xlabel('Principal Component 1')
plt.ylabel('Principal Component 2')
plt.show()

尽管只有两个主成分，我们仍然可以看到一些数字类别形成了相对清晰的簇，特别是数字0、1和7。这说明即使在高维数据中，PCA也能有效捕捉数据结构的主要模式。

7.1 重建图像与成分分析

我们可以通过选择不同数量的主成分来重建原始图像，观察信息保留情况：

def plot_digits_reconstruction(X, n_components):
    pca = PCA(n_components=n_components)
    X_pca = pca.fit_transform(X)
    X_reconstructed = pca.inverse_transform(X_pca)
    
    plt.figure(figsize=(8, 4))
    plt.subplot(1, 2, 1)
    plt.imshow(X[0].reshape(28, 28), cmap='gray')
    plt.title('Original')
    
    plt.subplot(1, 2, 2)
    plt.imshow(X_reconstructed[0].reshape(28, 28), cmap='gray')
    plt.title(f'Reconstructed (n={n_components})')
    plt.show()

# 使用不同数量的主成分重建
plot_digits_reconstruction(X_mnist[:100], 10)  # 使用前10个主成分
plot_digits_reconstruction(X_mnist[:100], 50)  # 使用前50个主成分

随着使用的主成分数量增加，重建图像的质量会逐渐提高。这个实验直观展示了PCA如何在保留大部分信息的同时大幅降低数据维度。

8. PCA与其他降维技术的比较

PCA是最基础的线性降维方法，还有其他多种降维技术各有特点：

选择哪种方法取决于具体需求。如果目标是探索性数据分析或可视化，t-SNE或UMAP可能更合适；如果目标是特征压缩或去噪，PCA通常是更好的选择。