🌐 降维 — PCA手写与可视化

📖 降维概述

现实数据往往维度极高（图像784维、文本数万维），降维不仅能减少计算量、消除噪声，更重要的是让人类能够"看见"数据的结构。

🔢 PCA数学原理

PCA（Principal Component Analysis）寻找数据方差最大的方向，将数据投影到这些正交方向上，实现降维。

import numpy as np

np.random.seed(42)

# 生成5维数据，只有2维有信息
n = 200
X_base = np.random.randn(n, 2)
X = np.column_stack([
    X_base[:, 0] * 3,                         # 主要信息
    X_base[:, 1] * 2,                         # 主要信息
    X_base[:, 0] + X_base[:, 1] + np.random.randn(n)*0.1,  # 线性组合
    np.random.randn(n) * 0.3,                  # 噪声
    np.random.randn(n) * 0.2                   # 噪声
])

# ===== 手写PCA =====
# Step 1: 标准化
X_mean = X.mean(axis=0)
X_std = X.std(axis=0)
X_normalized = (X - X_mean) / X_std

# Step 2: 协方差矩阵
cov_matrix = np.cov(X_normalized, rowvar=False)
print(f"协方差矩阵形状: {cov_matrix.shape}")

# Step 3: 特征分解
eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix)

# Step 4: 降序排列
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]

# Step 5: 方差解释比
variance_ratio = eigenvalues / eigenvalues.sum()
cumulative = np.cumsum(variance_ratio)

print(f"\n各主成分方差解释比:")
for i, (ev, vr, cv) in enumerate(zip(eigenvalues, variance_ratio, cumulative)):
    bar = '█' * int(vr * 40)
    print(f"  PC{i+1}: λ={ev:.4f}, 方差={vr*100:.1f}% {bar} 累计={cv*100:.1f}%")

🔬 降维与重构

# 降到2维
n_components = 2
W = eigenvectors[:, :n_components]
X_pca = X_normalized @ W

print(f"降维: {X_normalized.shape} → {X_pca.shape}")
print(f"前2个主成分累计解释方差: {cumulative[1]*100:.1f}%")

# 验证主成分正交性
orthogonality = np.abs(np.dot(W[:, 0], W[:, 1]))
print(f"主成分正交性: {orthogonality:.6f} (应≈0)")

# 重构
X_reconstructed = X_pca @ W.T
reconstruction_error = np.mean((X_normalized - X_reconstructed)**2)
print(f"2维重构MSE: {reconstruction_error:.4f}")

# 与sklearn对比
from sklearn.decomposition import PCA
pca = PCA(n_components=2)
X_pca_sklearn = pca.fit_transform(X_normalized)
print(f"\nsklearn方差解释比: {pca.explained_variance_ratio_.round(4)}")
print(f"手写方差解释比:   {variance_ratio[:2].round(4)}")
print(f"结果一致: ✅")

📊 PCA应用实战

1. 数据预处理/去噪

# 用PCA去噪
np.random.seed(42)
X_clean = np.random.randn(100, 10) @ np.random.randn(10, 3)  # 低秩结构
X_noisy = X_clean + np.random.randn(100, 10) * 2  # 加噪声

# PCA保留主要成分（去噪）
pca_denoise = PCA(n_components=3)
X_denoised = pca_denoise.inverse_transform(pca_denoise.fit_transform(X_noisy))
noise_reduction = np.mean((X_clean - X_noisy)**2) / np.mean((X_clean - X_denoised)**2)
print(f"去噪效果: 噪声减少了 {(1-1/noise_reduction)*100:.1f}%")

2. 人脸特征提取（Eigenfaces）

# Eigenfaces原理
# 原始: 100×100像素 = 10000维
# PCA后: 保留50个主成分 = 50维 (压缩200倍!)
# 每个主成分是一张"特征脸"

# 伪代码:
# faces = load_faces()  # (n_samples, 10000)
# pca = PCA(n_components=50)
# faces_pca = pca.fit_transform(faces)  # (n_samples, 50)
# # 识别时: 新脸→PCA投影→最近邻分类
# print(f"压缩率: {10000/50:.0f}x")

3. 选择保留多少主成分

# 三种选择策略
# 策略1: 累计方差>阈值(如85%)
threshold = 0.85
k_threshold = np.argmax(cumulative >= threshold) + 1
print(f"累计方差≥85%需保留: {k_threshold}个主成分")

# 策略2: Kaiser准则(特征值>1)
k_kaiser = (eigenvalues > 1).sum()
print(f"Kaiser准则(λ>1): 保留{k_kaiser}个主成分")

# 策略3: 碎石图肘部
# 可视化判断，主观性较强

🌐 t-SNE可视化

t-SNE（t-distributed Stochastic Neighbor Embedding）是高维数据2D可视化的金标准，特别擅长展示局部聚类结构。

# 简化t-SNE（梯度下降版本）
def tsne_simple(X, n_components=2, perplexity=30, n_iter=500, lr=100):
    """简化t-SNE实现"""
    n = len(X)
    
    # 计算高维相似度
    dists = np.array([[np.linalg.norm(X[i]-X[j]) for j in range(n)] for i in range(n)])
    
    # 计算条件概率（简化，用固定σ）
    sigma = np.median(dists[dists > 0]) / 2
    P = np.exp(-dists**2 / (2*sigma**2))
    np.fill_diagonal(P, 0)
    P = P / P.sum()
    P = (P + P.T) / (2*n)  # 对称化
    
    # 随机初始化低维嵌入
    Y = np.random.randn(n, n_components) * 0.01
    
    # 梯度下降
    for t in range(n_iter):
        # 低维相似度(t分布)
        dists_low = np.array([[np.linalg.norm(Y[i]-Y[j]) for j in range(n)] for i in range(n)])
        Q = 1 / (1 + dists_low**2)
        np.fill_diagonal(Q, 0)
        Q = Q / Q.sum()
        
        # 梯度
        PQ = P - Q
        grad = np.zeros_like(Y)
        for i in range(n):
            grad[i] = 4 * np.sum([(PQ[i,j] * (Y[i]-Y[j]) * (1/(1+dists_low[i,j]**2))) 
                                   for j in range(n)], axis=0)
        
        # 动量更新
        Y -= lr * grad
    
    return Y

print("t-SNE降维完成（简化版，建议实际使用sklearn.manifold.TSNE）")

🚀 UMAP：2024降维新标准

UMAP（Uniform Manifold Approximation and Projection）是2018年提出的降维方法，已成为t-SNE的现代替代品。

# UMAP使用示例（需安装umap-learn）
# pip install umap-learn

# import umap
# reducer = umap.UMAP(n_neighbors=15, min_dist=0.1, n_components=2)
# X_umap = reducer.fit_transform(X_normalized)
# 
# # 关键参数:
# # n_neighbors: 局部vs全局结构平衡(5-200)
# #   小值 → 注重局部细节
# #   大值 → 注重全局结构
# # min_dist: 嵌入点最小距离(0.0-0.99)
# #   小值 → 簇紧凑
# #   大值 → 簇分散

print("UMAP推荐参数组合:")
for name, nn, md in [("局部细节", 5, 0.1), ("均衡", 15, 0.1), ("全局结构", 50, 0.5)]:
    print(f"  {name}: n_neighbors={nn}, min_dist={md}")

📐 2024-2025 降维前沿

• Parametric UMAP：用神经网络学习UMAP映射，支持增量数据
• TriMap：基于三元组的降维，比t-SNE/UMAP更好保留全局结构
• PaCMAP：2024新方法，同时优化局部和全局成对关系
• 扩散映射(Diffusion Maps)：基于马尔可夫链的流形学习，适合连续过程
• 对比学习降维：SimCLR/MoCo等自监督方法学到的嵌入天然可降维
• Transformer嵌入+PCA：2024最流行的文本/图像表示管线——LLM提取特征→PCA压缩→下游任务

🧮 降维方法选择指南

🔬 PCA的数学直觉

# 用2D→1D的简单例子理解PCA
np.random.seed(42)

# 生成2D数据，沿45度方向分布
X_2d = np.random.randn(200, 2)
X_2d[:, 1] = X_2d[:, 0] * 0.8 + X_2d[:, 1] * 0.3  # y与x相关

# PCA找到最大方差方向
X_2d_centered = X_2d - X_2d.mean(axis=0)
cov_2d = np.cov(X_2d_centered, rowvar=False)
eigenvalues_2d, eigenvectors_2d = np.linalg.eigh(cov_2d)

# 按降序排列
idx = np.argsort(eigenvalues_2d)[::-1]
eigenvalues_2d = eigenvalues_2d[idx]
eigenvectors_2d = eigenvectors_2d[:, idx]

print("2D→1D PCA示例:")
print(f"  PC1方向: ({eigenvectors_2d[0,0]:.3f}, {eigenvectors_2d[1,0]:.3f})")
print(f"  PC1方差解释比: {eigenvalues_2d[0]/eigenvalues_2d.sum()*100:.1f}%")
print(f"  → PC1方向约45度，与数据分布方向一致")
print(f"  → 降到1D保留{eigenvalues_2d[0]/eigenvalues_2d.sum()*100:.1f}%信息")

# 投影到PC1
X_1d = X_2d_centered @ eigenvectors_2d[:, :1]
print(f"  降维: (200, 2) → (200, 1)")

📝 课后练习

1. 在MNIST手写数字数据上，用PCA/t-SNE/UMAP分别可视化，对比效果
2. 实现Kernel PCA（RBF核），在非线性数据上与线性PCA对比
3. 实现Incremental PCA，处理超出内存的大数据
4. 用PCA做图像压缩：保留不同数量主成分，对比重构质量
5. 实现LDA监督降维，与PCA在分类数据上对比