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()