PCA 数值计算
数值计算了一些结果。
对于中心化数据矩阵B, B有d行n列，d代表特征个数，n代表样本个数，SVD分解为：
协方差矩阵为：
协方差矩阵的特征值λ_i与中心化数据矩阵B的奇异值s_i的关系：
解释方差比例：
import numpy as np def verify_variance_equals_eigenvalue(): """验证主成分方差等于特征值""" np.random.seed(42) # 生成示例数据 d, n = 5, 100 A = np.random.randn(d, n) * np.array([[2], [1], [0.5], [0.3], [0.1]]) # 中心化 mu = A.mean(axis=1, keepdims=True) B = A - mu # 方法1: 通过协方差矩阵的特征分解 S = B @ B.T / (n - 1) # 协方差矩阵 eigenvalues, eigenvectors = np.linalg.eigh(S) # 特征分解 # 按特征值降序排列 idx = eigenvalues.argsort()[::-1] eigenvalues = eigenvalues[idx] eigenvectors = eigenvectors[:, idx] print("方法1: 通过协方差矩阵特征分解") print("特征值:", eigenvalues.round(6)) # 方法2: 通过投影计算方差 variances_from_projection = [] for i in range(d): # 将数据投影到第i个特征向量方向 z = eigenvectors[:, i].T @ B variance = np.var(z, ddof=1) # 样本方差 variances_from_projection.append(variance) print("\n方法2: 通过投影计算方差") print("投影方差:", np.array(variances_from_projection).round(6)) # 方法3: 通过SVD U, S_svd, Vt = np.linalg.svd(B, full_matrices=False) eigenvalues_from_svd = (S_svd**2) / (n - 1) print("\n方法3: 通过SVD计算") print("从SVD得到的特征值:", eigenvalues_from_svd.round(6)) # 验证三种方法的一致性 print("\n验证一致性:") print("特征值与投影方差是否相等:", np.allclose(eigenvalues, variances_from_projection)) print("特征值与SVD特征值是否相等:", np.allclose(eigenvalues, eigenvalues_from_svd)) # 验证总方差守恒 total_variance_cov = np.trace(S) # 协方差矩阵的迹 total_variance_eigen = np.sum(eigenvalues) # 特征值之和 total_variance_data = np.sum(np.var(B, axis=1, ddof=1)) # 原始数据各特征方差之和 print("\n总方差验证:") print(f"协方差矩阵的迹: {total_variance_cov:.6f}") print(f"特征值之和: {total_variance_eigen:.6f}") print(f"原始数据各特征方差之和: {total_variance_data:.6f}") return eigenvalues, eigenvectors, B def geometric_interpretation(): """几何解释：主成分方差最大化""" np.random.seed(42) # 生成二维数据 n = 1000 theta = np.linspace(0, 2*np.pi, n) x = 3 * np.cos(theta) + np.random.randn(n) * 0.5 y = 1 * np.sin(theta) + np.random.randn(n) * 0.2 # 创建2×n的数据矩阵 A = np.vstack([x, y]) # 中心化 B = A - A.mean(axis=1, keepdims=True) # 计算协方差矩阵和特征分解 S = B @ B.T / (n - 1) eigenvalues, eigenvectors = np.linalg.eigh(S) idx = eigenvalues.argsort()[::-1] eigenvalues = eigenvalues[idx] eigenvectors = eigenvectors[:, idx] # 主成分方向 pc1, pc2 = eigenvectors[:, 0], eigenvectors[:, 1] print("二维数据示例:") print(f"特征值: {eigenvalues}") print(f"第一主成分方向: {pc1}") print(f"第二主成分方向: {pc2}") # 计算各方向投影方差 directions = np.linspace(0, np.pi, 180) # 测试180个方向 variances = [] for angle in directions: w = np.array([np.cos(angle), np.sin(angle)]) # 单位方向向量 z = w.T @ B variance = np.var(z, ddof=1) variances.append(variance) # 找到最大方差方向 max_idx = np.argmax(variances) max_angle = directions[max_idx] max_variance = variances[max_idx] print(f"\n通过搜索找到的最大方差方向: 角度={max_angle:.3f} rad") print(f"最大方差: {max_variance:.6f}") print(f"第一特征值: {eigenvalues[0]:.6f}") print(f"第一主成分方向的角度: {np.arctan2(pc1[1], pc1[0]):.3f} rad") # 验证：第一主成分方向的投影方差 z_pc1 = pc1.T @ B variance_pc1 = np.var(z_pc1, ddof=1) print(f"\n第一主成分方向的投影方差: {variance_pc1:.6f}") print(f"是否等于特征值: {np.isclose(variance_pc1, eigenvalues[0])}") return A, eigenvalues, eigenvectors, directions, variances if __name__ == "__main__": print("=" * 60) print("验证1: 代数推导验证") print("=" * 60) eigenvalues, eigenvectors, B = verify_variance_equals_eigenvalue() print("\n" + "=" * 60) print("验证2: 几何解释 - 方差最大化") print("=" * 60) A, eigvals, eigvecs, angles, var_list = geometric_interpretation() # 可视化几何解释 import matplotlib.pyplot as plt fig, axes = plt.subplots(1, 3, figsize=(15, 5)) # 1. 原始数据散点图 axes[0].scatter(A[0, :], A[1, :], alpha=0.5, s=10) axes[0].set_xlabel('X') axes[0].set_ylabel('Y') axes[0].set_title('原始数据') axes[0].grid(True, alpha=0.3) axes[0].set_aspect('equal') # 添加主成分方向 center = A.mean(axis=1) for i, (vec, val) in enumerate(zip(eigvecs.T, eigvals)): length = np.sqrt(val) * 3 axes[0].arrow(center[0], center[1], vec[0]*length, vec[1]*length, head_width=0.1, head_length=0.2, fc=f'C{i}', ec=f'C{i}', label=f'PC{i+1} (λ={val:.2f})') axes[0].legend() # 2. 各方向投影方差 axes[1].plot(angles, var_list, 'b-', label='投影方差') axes[1].axvline(x=angles[np.argmax(var_list)], color='r', linestyle='--', label=f'最大方差方向 ({angles[np.argmax(var_list)]:.2f} rad)') axes[1].axhline(y=eigvals[0], color='g', linestyle=':', label=f'第一特征值 ({eigvals[0]:.2f})') axes[1].set_xlabel('方向角度 (rad)') axes[1].set_ylabel('投影方差') axes[1].set_title('不同方向的投影方差') axes[1].legend() axes[1].grid(True, alpha=0.3) # 3. 特征值条形图 axes[2].bar(range(1, len(eigvals)+1), eigvals, alpha=0.7) axes[2].set_xlabel('主成分') axes[2].set_ylabel('特征值（方差）') axes[2].set_title('主成分方差（特征值）') axes[2].grid(True, alpha=0.3, axis='y') # 添加累计方差比例 cumulative = np.cumsum(eigvals) / np.sum(eigvals) for i, (val, cum) in enumerate(zip(eigvals, cumulative)): axes[2].text(i+1, val*1.05, f'{cum:.1%}', ha='center', fontsize=9) plt.tight_layout() plt.show()