数值计算了一些结果。对于中心化数据矩阵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(axis1, keepdimsTrue) 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, ddof1) # 样本方差 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_matricesFalse) 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, axis1, ddof1)) # 原始数据各特征方差之和 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(axis1, keepdimsTrue) # 计算协方差矩阵和特征分解 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, ddof1) 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, ddof1) 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, :], alpha0.5, s10) axes[0].set_xlabel(X) axes[0].set_ylabel(Y) axes[0].set_title(原始数据) axes[0].grid(True, alpha0.3) axes[0].set_aspect(equal) # 添加主成分方向 center A.mean(axis1) 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_width0.1, head_length0.2, fcfC{i}, ecfC{i}, labelfPC{i1} (λ{val:.2f})) axes[0].legend() # 2. 各方向投影方差 axes[1].plot(angles, var_list, b-, label投影方差) axes[1].axvline(xangles[np.argmax(var_list)], colorr, linestyle--, labelf最大方差方向 ({angles[np.argmax(var_list)]:.2f} rad)) axes[1].axhline(yeigvals[0], colorg, linestyle:, labelf第一特征值 ({eigvals[0]:.2f})) axes[1].set_xlabel(方向角度 (rad)) axes[1].set_ylabel(投影方差) axes[1].set_title(不同方向的投影方差) axes[1].legend() axes[1].grid(True, alpha0.3) # 3. 特征值条形图 axes[2].bar(range(1, len(eigvals)1), eigvals, alpha0.7) axes[2].set_xlabel(主成分) axes[2].set_ylabel(特征值方差) axes[2].set_title(主成分方差特征值) axes[2].grid(True, alpha0.3, axisy) # 添加累计方差比例 cumulative np.cumsum(eigvals) / np.sum(eigvals) for i, (val, cum) in enumerate(zip(eigvals, cumulative)): axes[2].text(i1, val*1.05, f{cum:.1%}, hacenter, fontsize9) plt.tight_layout() plt.show()
PCA 数值计算
发布时间：2026/5/28 23:30:19
数值计算了一些结果。对于中心化数据矩阵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(axis1, keepdimsTrue) 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, ddof1) # 样本方差 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_matricesFalse) 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, axis1, ddof1)) # 原始数据各特征方差之和 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(axis1, keepdimsTrue) # 计算协方差矩阵和特征分解 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, ddof1) 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, ddof1) 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, :], alpha0.5, s10) axes[0].set_xlabel(X) axes[0].set_ylabel(Y) axes[0].set_title(原始数据) axes[0].grid(True, alpha0.3) axes[0].set_aspect(equal) # 添加主成分方向 center A.mean(axis1) 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_width0.1, head_length0.2, fcfC{i}, ecfC{i}, labelfPC{i1} (λ{val:.2f})) axes[0].legend() # 2. 各方向投影方差 axes[1].plot(angles, var_list, b-, label投影方差) axes[1].axvline(xangles[np.argmax(var_list)], colorr, linestyle--, labelf最大方差方向 ({angles[np.argmax(var_list)]:.2f} rad)) axes[1].axhline(yeigvals[0], colorg, linestyle:, labelf第一特征值 ({eigvals[0]:.2f})) axes[1].set_xlabel(方向角度 (rad)) axes[1].set_ylabel(投影方差) axes[1].set_title(不同方向的投影方差) axes[1].legend() axes[1].grid(True, alpha0.3) # 3. 特征值条形图 axes[2].bar(range(1, len(eigvals)1), eigvals, alpha0.7) axes[2].set_xlabel(主成分) axes[2].set_ylabel(特征值方差) axes[2].set_title(主成分方差特征值) axes[2].grid(True, alpha0.3, axisy) # 添加累计方差比例 cumulative np.cumsum(eigvals) / np.sum(eigvals) for i, (val, cum) in enumerate(zip(eigvals, cumulative)): axes[2].text(i1, val*1.05, f{cum:.1%}, hacenter, fontsize9) plt.tight_layout() plt.show()
