哈喽，大家好~
今儿和大家详细的聊聊：核PCA（Kernel PCA）？
以及，我们为什么要用它？
简单的想象一下，你有一堆点，这些点在一个高维或者弯曲的低维流形上分布。我们常用的PCA（主成分分析）是线性的：它在原始空间里寻找一些“直线/平面”来捕捉数据最多的方差。
如果数据在一条弯曲的曲线上，比如圆环或 Swiss roll，线性PCA往往捕捉不到有意义的结构，你投影后看到的可能只是模糊的、不可分的点云。
核PCA的核心思想很直白但很聪明：先把数据“隐式地”映射到一个高维（甚至无限维）的特征空间，在那个空间里，原本非线性的结构可能变得线性可分；然后在那个高维空间里做PCA。
因为直接在高维空间里显式计算映射可能不现实（甚至不可行），我们用“核技巧（kernel trick）”实现：只通过数据之间的相似度（核函数）矩阵，就能完成在高维空间的PCA运算。
也就是，核PCA = 用核函数把数据隐式地升维到一个更容易做线性PCA的空间，然后在那个空间做PCA，最后把结果投影回低维去看。
核心逻辑
我们可以设数据点为 ，在特征映射 下对应 。
我们想在特征空间里做PCA，解决协方差矩阵的特征值问题：
求解 。在特征空间直接做非常难，但利用“核矩阵” 有捷径。
定义核矩阵 ，元素为 。可以证明我们可以先对 做中心化，然后求解：
得到的向量 给出特征空间的主方向
对新点 的投影为
其中要注意训练时与测试时核值的中心化（关键步骤，别忘了！）。
常用核函数：RBF（高斯）核 ，多项式核 等。RBF核非常常用，因为它能把邻近点的相似度放大，适合很多非线性流形。
优点与限制：
-
优点：能捕捉很多常见的非线性结构；实现上只需要核矩阵，不用显式映射。 -
限制：计算复杂度为 （核矩阵）和 （特征分解），在大数据上难以直接应用；另外核PCA的“逆映射”（pre-image）是难题，即如何把特征空间的投影精确地映射回原空间通常没有解析解，需要近似方法。
完整案例
我们用 PyTorch 完整实现核PCA，fit + transform + 对新样本映射~
其中代码注释非常的详细，大家可以详细的看看其中每一步的逻辑与含义~
import numpy as np
import torch
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
torch.manual_seed(42)
np.random.seed(42)
# 1) 生成 Swiss roll 数据（n x 3）
def make_swiss_roll(n_samples=1000, noise=0.0):
t = 1.5 * np.pi * (1 + 2 * np.random.rand(n_samples))
x = t * np.cos(t)
y = 21 * np.random.rand(n_samples)
z = t * np.sin(t)
X = np.vstack((x, y, z)).T
X += noise * np.random.normal(size=X.shape)
return X.astype(np.float32), t.astype(np.float32)
X_np, t = make_swiss_roll(n_samples=2800, noise=1)
X = torch.tensor(X_np) # (n, 3)
n = X.shape[0]
# 2) 定义核函数（RBF / Gaussian Kernel）
def rbf_kernel(X, Y=None, gamma=None):
# X: (n, d), Y: (m, d) or None -> if None use Y=X
if Y is None:
Y = X
X_norm = (X**2).sum(dim=1).view(-1, 1) # (n,1)
Y_norm = (Y**2).sum(dim=1).view(1, -1) # (1,m)
K = X_norm + Y_norm - 2.0 * X @ Y.t()
if gamma is None:
gamma = 1.0 / X.shape[1] # simple default
K = torch.exp(-gamma * K)
return K
# 3) 中心化核矩阵（训练集）
def center_kernel(K):
# K: (n,n)
n = K.shape[0]
one_n = torch.ones((n, n), dtype=K.dtype) / n
Kc = K - one_n @ K - K @ one_n + one_n @ K @ one_n
return Kc
# 4) KPCA fit：计算中心化核矩阵的特征分解，返回 alpha (n,k) 与 eigenvalues
def kpca_fit(K, n_components=2):
# K: centered kernel matrix
# 求解 K alpha = n * lambda * alpha
n = K.shape[0]
# PyTorch 的 eigh 需要对称矩阵的 numpy/tensor
vals, vecs = torch.linalg.eigh(K) # 升序返回
vals = vals.flip(0) # 降序
vecs = vecs.flip(1)
# 取前 k
vals_k = vals[:n_components].clone()
vecs_k = vecs[:, :n_components].clone()
# 归一化 alpha：要求 alpha_i^T K alpha_i = 1 -> vecs_k 已是K特征向量,需要归一化
# vecs_k 每列是 alpha (未缩放)
for i in range(n_components):
alpha = vecs_k[:, i]
norm = torch.sqrt((alpha.unsqueeze(0) @ K @ alpha.unsqueeze(1)).squeeze())
vecs_k[:, i] = alpha / (norm + 1e-12)
return vals_k, vecs_k # vals_k: (k,), vecs_k: (n,k)
# 5) 将训练点投影到主成分上（训练集上的变换）
def kpca_transform_train(Kc, alphas):
# Kc: centered kernel matrix (n,n)
# alphas: (n, k) 列向量是 alpha_i
# 投影 Z = Kc @ alphas 实际上等价于 (K @ alpha) 已包含 center 操作，Kc 更严谨
Z = Kc @ alphas
return Z # (n,k)
# 6) 将新样本投影到已训练的主成分上（注意中心化技巧）
def kpca_transform_new(X_train, X_new, alphas, K_train, gamma=None):
# X_train: (n,d), X_new: (m,d), alphas: (n,k), K_train: (n,n) original (uncentered) kernel on train
n = X_train.shape[0]
m = X_new.shape[0]
if gamma is None:
gamma = 1.0 / X_train.shape[1]
# 计算 k(x_new, X_train)
K_new = rbf_kernel(torch.tensor(X_new), X_train, gamma=gamma) # (m,n)
# 中心化： k_tilde(x, X) = k(x, X) - 1/n 1 K_train - 1/n K_x 1 + 1/n^2 1 K_train 1
one_n = torch.ones((n,), dtype=K_train.dtype) / n
K_row_mean = (K_train @ one_n).view(1, -1) # (1,n)
K_col_mean = K_row_mean.t() # (n,1)
K_total_mean = one_n @ K_row_mean.t() # scalar
# center each row of K_new:
K_new_centered = K_new - K_row_mean - K_new.mean(dim=1, keepdim=True) + K_total_mean
# 对应公式实际上是：K_new_centered = K_new - 1_n K_train - K_new 1_n + 1_n K_train 1_n （做广播）
# 投影值 Z_new = K_new_centered @ alphas -> (m,k)
Z_new = K_new_centered @ alphas
return Z_new.detach().numpy()
# 主流程
gamma = 1e-3 # RBF gamma（可调）
K = rbf_kernel(X, None, gamma=gamma) # 原始核矩阵（未中心化）
Kc = center_kernel(K) # 中心化核矩阵
# 特征分解
n_components = 2
eigvals, alphas = kpca_fit(Kc, n_components=n_components)
# 训练集投影
Z = kpca_transform_train(Kc, alphas) # (n, k), 训练点在主成分上的坐标
# 线性PCA做对比（直接在原始3D做PCA）
X_mean = X - X.mean(dim=0)
U, S, Vt = torch.linalg.svd(X_mean) # SVD: X_mean = U diag(S) Vt
Z_lin = (X_mean @ Vt.T)[:, :2].numpy() # 线性PCA的前2个主成分
# 可视化分析
# 1) 原始 Swiss roll（3D）
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
p = ax.scatter(X_np[:,0], X_np[:,1], X_np[:,2], c=t, cmap='rainbow', s=15, alpha=0.9)
ax.set_title('Original Swiss Roll (3D) ')
fig.colorbar(p, ax=ax, shrink=0.6, label='t parameter')
# 2) 核矩阵热力图（中心化后的 Kc）
plt.figure(figsize=(8,6))
plt.imshow(Kc.numpy(), cmap='viridis')
plt.title('Centered Kernel Matrix $K_c$')
plt.colorbar()
# 3) 特征值谱（降序）
plt.figure(figsize=(8,5))
vals_np = eigvals.numpy()
plt.plot(np.arange(1, len(vals_np)+1), vals_np, marker='o', color='crimson')
plt.title('Eigenvalues (top components) of Centered Kernel Matrix')
plt.xlabel('Component index')
plt.ylabel('Eigenvalue')
plt.grid(True, linestyle='--', alpha=0.6)
# 4) 投影结果对比：核PCA vs 线性PCA
plt.figure(figsize=(14,6))
plt.subplot(1,2,1)
plt.scatter(Z.detach().numpy()[:,0], Z.detach().numpy()[:,1], c=t, cmap='plasma', s=15)
plt.title('KPCA Projection (2D) with RBF kernel')
plt.xlabel('PC1 (kernel space)')
plt.ylabel('PC2 (kernel space)')
plt.colorbar(label='t parameter')
plt.subplot(1,2,2)
plt.scatter(Z_lin[:,0], Z_lin[:,1], c=t, cmap='plasma', s=15)
plt.title('Linear PCA Projection (2D)')
plt.xlabel('PC1 (linear)')
plt.ylabel('PC2 (linear)')
plt.colorbar(label='t parameter')
plt.show()
原始 Swiss roll：
展示数据真实的非线性结构，一个像卷轴一样的三维流形。
颜色按 渐变能看到点在卷轴上的顺序。
在原始空间用颜色可以显现点之间的拓扑顺序，但从 2D 投影观察时会丢失顺序信息（如果直接将 XYZ 投影到平面会混叠）。
核矩阵热力图：
显示训练样本在隐式特征空间中的两两相似度（做了中心化）。
沿对角线会看到高相似度（自相似），但重要的是在 off-diagonal 是否有明显的块状结构或渐变——这些表明在特征空间里数据按局部邻域或参数 t 排列成连续结构。RBF 的 gamma 大小会极大影响热力图的“颗粒度”：gamma 大 -> 相似度快速衰减 -> 更稀疏；gamma 小 -> 相似度更平滑。
特征值谱：
按降序显示 的特征值（与特征空间中主成分的“能量”相关）。大的特征值对应重要的主成分。
如果前 2 个或前少数几个特征值远大于其余，这说明数据在隐式特征空间里主要由这几个方向决定，说明降到低维很合理。Swiss roll 理想情况下前两个主成分会捕获“展开后的主流”信息。
KPCA vs Linear PCA：
把数据投影到 2 维平面，左是核PCA结果，右是线性PCA结果。点的颜色仍按 （原始参数）着色，便于看是否能把卷轴“展开”。
理想的核PCA图像会把颜色按顺序排列成一条接近单调曲线或流线（说明 KPCA 把 Swiss roll 展开了）；线性PCA 的颜色会混在一起，无法解开卷曲。这就清晰地说明了核PCA 能捕捉到非线性结构。
关于 gamma 的敏感性：
gamma（RBF 中的参数）类似于相似度尺度：较小 gamma -> 更宽的相似度尺度 -> 核矩阵更平滑；较大 gamma -> 相似度只在很近邻域内高。
建议实验：改变 gamma 值，重复上面的核矩阵热图与投影图，你会看到投影质量随 gamma 变化有较大不同。通常需要通过交叉验证或经验调参来选合适 gamma。
总结
总的来说，核PCA 的核心，隐式把数据升维 -> 在那个空间做线性PCA -> 利用核技巧避免显式计算映射。
关键数学量：核矩阵 、中心化 、解 ，投影为 （并且要对核做一致的中心化）。
核PCA 很强，但计算与内存代价高，参数（如 gamma）敏感，预影像求解复杂。实务中常结合近似核方法或随机特征映射。
最后
