企业官网建设流程全解析

从零实现PCA：用Python代码拆解主成分分析的数学内核

主成分分析（PCA）是数据科学中最经典的降维技术之一，但大多数教程要么停留在数学推导的抽象层面，要么直接调用sklearn黑箱。本文将带你用NumPy从零实现PCA核心算法，通过可运行的Python代码揭示数学公式背后的计算本质，最后与sklearn的结果进行交叉验证。

1. 环境准备与数据生成

在开始之前，我们需要准备实验环境。建议使用Jupyter Notebook或Google Colab这类交互式环境，方便实时观察每一步的计算结果。

首先导入必要的库：

import numpy as np import matplotlib.pyplot as plt from sklearn.decomposition import PCA

为了直观理解PCA的效果，我们创建一个具有明显相关性的二维数据集：

np.random.seed(42) mean = [0, 0] cov = [[3, 2], [2, 3]] # 协方差矩阵 X = np.random.multivariate_normal(mean, cov, 100)

可视化原始数据：

plt.scatter(X[:, 0], X[:, 1], alpha=0.7) plt.axis('equal') plt.title('Original Data') plt.show()

2. PCA的数学核心：分步代码实现

2.1 数据标准化

PCA对数据的尺度敏感，因此第一步需要对数据进行中心化处理（减去均值）：

X_centered = X - np.mean(X, axis=0)

验证均值是否接近零：

print(f"Mean after centering: {np.round(np.mean(X_centered, axis=0), 6)}")

2.2 计算协方差矩阵

协方差矩阵是PCA的核心，它捕获了数据各维度之间的关系：

cov_matrix = np.cov(X_centered, rowvar=False) print("Covariance matrix:\n", cov_matrix)

2.3 特征值分解

PCA的魔力就藏在协方差矩阵的特征值和特征向量中：

eigenvalues, eigenvectors = np.linalg.eig(cov_matrix) print("Eigenvalues:", eigenvalues) print("Eigenvectors:\n", eigenvectors)

特征向量决定了主成分的方向，而特征值表示各主成分解释的方差量。我们可以按特征值大小降序排列：

# 获取排序索引 sorted_idx = np.argsort(eigenvalues)[::-1] eigenvalues = eigenvalues[sorted_idx] eigenvectors = eigenvectors[:, sorted_idx] print("Sorted eigenvalues:", eigenvalues) print("Sorted eigenvectors:\n", eigenvectors)

2.4 选择主成分

假设我们要降到1维，选择最大特征值对应的特征向量：

principal_component = eigenvectors[:, 0] print("First principal component:", principal_component)

2.5 数据投影

将数据投影到选定的主成分上：

X_pca = X_centered @ principal_component

可视化投影结果：

plt.figure(figsize=(10, 4)) # 原始数据 plt.subplot(1, 2, 1) plt.scatter(X_centered[:, 0], X_centered[:, 1], alpha=0.7) plt.axis('equal') plt.title('Original Data') # 投影后的数据 plt.subplot(1, 2, 2) plt.scatter(X_pca, np.zeros_like(X_pca), alpha=0.7) plt.title('Projected Data (1D)') plt.tight_layout() plt.show()

3. 与Scikit-learn的对比验证

为了验证我们的实现是否正确，下面使用sklearn的PCA进行对比：

sklearn_pca = PCA(n_components=1) X_sklearn = sklearn_pca.fit_transform(X) print("Sklearn components:\n", sklearn_pca.components_) print("Our components:\n", principal_component.reshape(1, -1))

比较两种方法的结果差异：

difference = np.abs(X_pca.reshape(-1, 1) - X_sklearn) print("Max difference:", np.max(difference))

注意：由于特征向量的方向不影响PCA的数学本质，两种实现可能得到方向相反但等价的向量，这属于正常现象。

4. PCA的进阶理解与应用技巧

4.1 方差解释率

每个主成分解释的方差比例是评估降维效果的重要指标：

explained_variance_ratio = eigenvalues / np.sum(eigenvalues) print("Explained variance ratio:", explained_variance_ratio)

通常我们会绘制碎石图（Scree Plot）来帮助选择主成分数量：

plt.plot(range(1, len(eigenvalues)+1), explained_variance_ratio, 'o-') plt.xlabel('Principal Component') plt.ylabel('Explained Variance Ratio') plt.title('Scree Plot') plt.show()

4.2 数据重构

了解如何从降维后的数据重构原始空间：

X_reconstructed = (X_pca.reshape(-1, 1) @ principal_component.reshape(1, -1)) + np.mean(X, axis=0)

计算重构误差：

reconstruction_error = np.mean(np.square(X - X_reconstructed)) print("Reconstruction MSE:", reconstruction_error)

4.3 高维数据可视化

对于高于3维的数据，PCA常被用于可视化：

from sklearn.datasets import load_iris iris = load_iris() X_iris = iris.data y_iris = iris.target # 我们的PCA实现 X_centered_iris = X_iris - np.mean(X_iris, axis=0) cov_iris = np.cov(X_centered_iris, rowvar=False) eigvals_iris, eigvecs_iris = np.linalg.eig(cov_iris) sorted_idx_iris = np.argsort(eigvals_iris)[::-1] eigvecs_iris = eigvecs_iris[:, sorted_idx_iris] X_pca_iris = X_centered_iris @ eigvecs_iris[:, :2] # 可视化 plt.scatter(X_pca_iris[:, 0], X_pca_iris[:, 1], c=y_iris) plt.xlabel('PC1') plt.ylabel('PC2') plt.title('Iris Dataset PCA Projection') plt.show()

5. PCA的实用注意事项

5.1 数据预处理

PCA对数据的尺度敏感，不同量纲的特征需要标准化：

from sklearn.preprocessing import StandardScaler scaler = StandardScaler() X_scaled = scaler.fit_transform(X)

5.2 特征值稳定性

对于小样本高维数据，协方差矩阵估计可能不稳定：

# 小样本高维数据示例 X_highdim = np.random.randn(10, 100) cov_highdim = np.cov(X_highdim, rowvar=False) print("Condition number:", np.linalg.cond(cov_highdim))

5.3 核PCA扩展

对于非线性数据，可以考虑核PCA：

from sklearn.decomposition import KernelPCA kpca = KernelPCA(n_components=2, kernel='rbf') X_kpca = kpca.fit_transform(X)

6. PCA在真实场景中的应用模式

6.1 数据压缩

评估压缩比和保留的信息量：

original_size = X.nbytes compressed_size = X_pca.nbytes + principal_component.nbytes compression_ratio = original_size / compressed_size print(f"Compression ratio: {compression_ratio:.1f}x")

6.2 噪声过滤

通过只保留主要成分来去除噪声：

# 添加噪声 X_noisy = X + np.random.normal(0, 0.5, X.shape) # 去噪 X_denoised = (X_noisy @ principal_component).reshape(-1, 1) @ principal_component.reshape(1, -1)

6.3 特征工程

PCA转换后的特征可以作为新特征：

from sklearn.ensemble import RandomForestClassifier from sklearn.model_selection import cross_val_score # 原始特征 scores_original = cross_val_score(RandomForestClassifier(), X, y_iris, cv=5) print("Original features accuracy:", np.mean(scores_original)) # PCA特征 scores_pca = cross_val_score(RandomForestClassifier(), X_pca_iris, y_iris, cv=5) print("PCA features accuracy:", np.mean(scores_pca))

7. PCA的数学本质深度解析

7.1 优化视角

PCA可以看作是在寻找使投影方差最大的方向：

# 验证第一主成分确实最大化投影方差 projection_variances = [] for angle in np.linspace(0, np.pi, 100): direction = np.array([np.cos(angle), np.sin(angle)]) projection = X_centered @ direction projection_variances.append(np.var(projection)) plt.plot(np.linspace(0, np.pi, 100), projection_variances) plt.axvline(np.arctan2(principal_component[1], principal_component[0]), color='r', linestyle='--') plt.xlabel('Angle (radians)') plt.ylabel('Projection Variance') plt.title('Variance vs Projection Direction') plt.show()

7.2 奇异值分解(SVD)关系

PCA也可以通过SVD实现：

U, s, Vt = np.linalg.svd(X_centered) print("From SVD:\n", Vt.T[:, 0]) print("From eigendecomposition:\n", principal_component)

7.3 统计解释

主成分实际上对应着数据的多元正态分布的等概率椭圆轴：

from scipy.stats import multivariate_normal x, y = np.mgrid[-4:4:.01, -4:4:.01] pos = np.dstack((x, y)) rv = multivariate_normal(mean, cov) plt.contour(x, y, rv.pdf(pos)) plt.scatter(X[:, 0], X[:, 1], alpha=0.3) plt.quiver(0, 0, principal_component[0], principal_component[1], angles='xy', scale_units='xy', scale=1, color='r') plt.axis('equal') plt.title('Principal Component as Normal Distribution Axis') plt.show()