三维点云处理（四）：Kernel PCA 核主成分分析
在上一章中，我们推导了标准 PCA 的数学原理。然而，现实世界中的点云数据往往呈现出非线性分布——例如螺旋状点云、环形轨迹或任意弯曲的曲面。对于这类数据，线性 PCA 的降维效果非常有限。
Kernel PCA（核主成分分析）通过**核技巧（Kernel Trick）**将数据隐式地映射到高维特征空间，使得在原空间中非线性可分的数据在映射后的高维空间中变为线性可分，再执行标准 PCA。
一、从线性到非线性：PCA 的局限性
1.1 线性 PCA 的失效场景
在二维平面中，考虑三组同心圆分布的点云数据（如图所示）。对于这种向心对称的非线性分布，数据在各个方向上的总体方差几乎一致，不存在一条直线能够将外圈、中圈和内圈线性分隔开。
线性 PCA 只能找到数据总体方差最大的直线方向。当数据沿曲线、椭圆或同心圆环分布时，在原始维度的任何直线投影都会丢失其内在的非线性拓扑结构。
1.2 核方法的核心思想
核方法的基本思路是：通过非线性映射 将数据从原始空间映射到高维（甚至无限维）特征空间 ，然后在 中执行线性 PCA。
根据 Cover 定理（Cover's Theorem），复杂的高维非线性模式在更高的维特征空间中几乎总是线性可分的。
二、数学推导：从协方差到核矩阵
2.1 特征空间中的 PCA 优化问题
设映射函数为 ，且有 个数据点 。假设映射后的数据已中心化：
特征空间中的协方差矩阵为：
PCA 的目标是求解 的特征向量 ：
2.2 特征向量存在于数据张成的空间中
由于 是 个 的外积之和，其特征向量 必然位于 所张成的子空间中。因此存在一组系数 使得：
2.3 核矩阵的引入
将 的展开式代入特征方程 ：
对两边左乘 （对任意 ），并引入核函数 ：
用矩阵形式表示为：
其中核矩阵 定义为：
2.4 简化与求解
消去 （假设 可逆），得到简化形式：
这意味着 是核矩阵 的特征向量， 是对应的特征值。
数学推导的核心结论总结
我们不需要显式知道高维映射函数 具体长什么样。只要我们能定义一种衡量两个点域之间相似度的核函数 （例如高斯核），我们就能直接构建核矩阵 ，对其进行特征值分解，从而完成在高维空间的 PCA。这就是著名的核技巧 (Kernel Trick)。
三、常用核函数与维基百科经典实例
核函数 隐含地定义了一个特征映射 。选择不同的核函数相当于选择不同的特征空间。
| 核函数 | 数学形式 | 参数 | 适用场景 |
|---|---|---|---|
| 线性核 | 无 | 退化为标准 PCA | |
| 多项式核 | , | 已知多项式阶数关系 | |
| 高斯核 (RBF) | （带宽） | 通用场景，最常用 | |
| Sigmoid 核 | , | 神经网络视角 |
3.1 经典同心圆案例：多项式核（Polynomial Kernel）
考虑维基百科中三组同心圆点云数据，使用阶数 、常数项 的多项式核函数：
在特征空间求解 Kernel PCA 后的投影结果如下图所示：
物理与数学直觉：
- 展开二次多项式核 ，这等价于将二维平面数据点 显式映射到了 6 维特征空间：
- 在这个 6 维特征空间中，同心圆的半径平方项 变为一个线性维度！
- 当在这个高维特征空间中求解 PCA 并将数据投影到第一主成分（PC1）上时，三个在二维平面中交叉嵌套的同心圆被完全拉开到不同的数值坐标区间上。如上图所示，仅凭第一主成分坐标就可以完美地线性分割这三组同心圆。
3.2 高斯核（Gaussian / RBF Kernel）
高斯核（又称径向基核函数 RBF Kernel）定义为：
在相同的同心圆点云数据上应用高斯核 Kernel PCA 的投影结果如下图所示：
物理与数学直觉：
- 高斯核本质上是对数据点之间局部空间相似度的度量：当两点重合（）时，相似度达到最大值 ；随着两点间距离增大约束衰减，相似度呈指数级下降趋近于 。
- 映射对应的特征空间是一个无限维希尔伯特空间（Hilbert Space）。在此无限维特征空间中，不同半径的独立同心圈簇被映射到了相互正交的基向量方向上。
- 如上图所示，高斯核 KPCA 投影后，外圈、中圈、内圈的三组点云分别紧密地凝聚在低维平面的不同离散区域中，展现出极强的非线性流形解耦能力。
四、中心化处理
4.1 特征空间中的中心化
前面的推导假设 ，但这在现实中并不成立。特征空间的中心化公式为：
4.2 中心化后的核矩阵
中心化后的核矩阵 可以通过原始核矩阵 直接计算：
其中 是 的矩阵，每个元素均为 。
用代码表示即为：
N = K.shape[0]
one_n = np.ones((N, N)) / N
K_centered = K - one_n @ K - K @ one_n + one_n @ K @ one_n五、新数据点的投影
5.1 投影公式
对于一个新数据点 ，其到第 个主成分 的投影为：
其中 是核矩阵 的第 个特征向量（需要归一化为 ）。
5.2 重构（Pre-image Problem）
与标准 PCA 不同，Kernel PCA 的逆映射（Pre-image）问题通常没有闭式解——在特征空间中找到一个点后，很难找到原始空间中对应的点。这需要使用数值优化方法（如梯度下降）来近似求解：
六、Python 实现
import numpy as np
from scipy.linalg import eigh
def rbf_kernel(X, Y=None, sigma=1.0):
"""
计算高斯核 (RBF Kernel) 矩阵。
k(x, y) = exp(-||x - y||² / (2σ²))
:param X: N x d 的输入矩阵
:param Y: M x d 的输入矩阵 (None 表示 Y = X)
:param sigma: 高斯核带宽参数
:return: N x M 的核矩阵
"""
if Y is None:
Y = X
# 计算成对距离平方: ||x - y||² = ||x||² + ||y||² - 2 x·y
X_norm_sq = np.sum(X ** 2, axis=1).reshape(-1, 1)
Y_norm_sq = np.sum(Y ** 2, axis=1).reshape(1, -1)
sq_dists = X_norm_sq + Y_norm_sq - 2 * np.dot(X, Y.T)
return np.exp(-sq_dists / (2 * sigma ** 2))
def kernel_pca(X, n_components=2, kernel='rbf', sigma=1.0, degree=3, coef0=1):
"""
执行 Kernel PCA 降维。
:param X: N x d 的输入数据矩阵，每行为一个数据点
:param n_components: 保留的主成分数量
:param kernel: 核函数类型 ('rbf', 'poly', 'linear')
:param sigma: RBF 核的带宽参数
:param degree: 多项式核的阶数
:param coef0: 多项式核的常数项
:return: (alphas, lambdas, projected)
alphas — 特征向量矩阵 (N x n_components)
lambdas — 特征值向量 (n_components,)
projected — 投影后的数据 (N x n_components)
"""
N = X.shape[0]
# 1. 计算核矩阵
if kernel == 'rbf':
K = rbf_kernel(X, sigma=sigma)
elif kernel == 'poly':
K = (np.dot(X, X.T) + coef0) ** degree
elif kernel == 'linear':
K = np.dot(X, X.T)
else:
raise ValueError(f"Unsupported kernel: {kernel}")
# 2. 中心化核矩阵
one_n = np.ones((N, N)) / N
K_centered = K - one_n @ K - K @ one_n + one_n @ K @ one_n
# 3. 特征值分解（取最大的 n_components 个特征对）
eigenvalues, eigenvectors = eigh(K_centered)
# eigh 返回升序，我们取最后 n_components 个（即最大的）
idx = np.argsort(eigenvalues)[::-1][:n_components]
lambdas = eigenvalues[idx]
alphas = eigenvectors[:, idx]
# 4. 归一化特征向量: ||α_k||² = 1 / λ_k
for k in range(n_components):
alphas[:, k] /= np.sqrt(lambdas[k])
# 5. 计算投影坐标
projected = K_centered @ alphas
return alphas, lambdas, projected
# ────── 使用示例 ──────
if __name__ == "__main__":
import matplotlib.pyplot as plt
# 生成同心圆数据（线性不可分）
np.random.seed(42)
n_samples = 200
# 外圈
theta_outer = np.random.uniform(0, 2 * np.pi, n_samples)
r_outer = 2.0 + 0.1 * np.random.randn(n_samples)
outer = np.column_stack([r_outer * np.cos(theta_outer),
r_outer * np.sin(theta_outer)])
# 内圈
theta_inner = np.random.uniform(0, 2 * np.pi, n_samples)
r_inner = 1.0 + 0.1 * np.random.randn(n_samples)
inner = np.column_stack([r_inner * np.cos(theta_inner),
r_inner * np.sin(theta_inner)])
X = np.vstack([outer, inner]) # 400 x 2
# Kernel PCA 降维
alphas, lambdas, projected = kernel_pca(X, n_components=2, kernel='rbf', sigma=0.5)
# 可视化对比
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
axes[0].scatter(X[:n_samples, 0], X[:n_samples, 1], c='blue', s=15, alpha=0.6, label='Outer')
axes[0].scatter(X[n_samples:, 0], X[n_samples:, 1], c='red', s=15, alpha=0.6, label='Inner')
axes[0].set_title("Original Space (concentric circles)")
axes[0].set_xlabel("x"); axes[0].set_ylabel("y")
axes[0].legend(); axes[0].axis('equal')
axes[1].scatter(projected[:n_samples, 0], projected[:n_samples, 1],
c='blue', s=15, alpha=0.6, label='Outer')
axes[1].scatter(projected[n_samples:, 0], projected[n_samples:, 1],
c='red', s=15, alpha=0.6, label='Inner')
axes[1].set_title(f"Kernel PCA (RBF, σ=0.5)\n"
f"λ₁={lambdas[0]:.1f}, λ₂={lambdas[1]:.1f}")
axes[1].set_xlabel("PC1"); axes[1].set_ylabel("PC2")
axes[1].legend()
plt.tight_layout()
plt.show()6.2 工程实战：基于 Sklearn 处理 3D Swiss Roll (瑞士卷) 点云
在实际工程中，我们极少手写 KPCA，而是直接调用 scikit-learn 库。Swiss Roll (瑞士卷) 是三维点云流形学习（Manifold Learning）中最经典的非线性测试用例：数据在 3D 空间中卷曲，但其内在（Intrinsic）结构其实是一张 2D 的平面。
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_swiss_roll
from sklearn.decomposition import PCA, KernelPCA
# 1. 生成 3D 瑞士卷点云
X, color = make_swiss_roll(n_samples=1000, noise=0.05, random_state=42)
# 2. 线性 PCA 降维 (试图降到 2D)
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X)
# 3. Kernel PCA 降维 (使用 RBF 高斯核)
kpca = KernelPCA(n_components=2, kernel="rbf", gamma=0.04, fit_inverse_transform=True)
X_kpca = kpca.fit_transform(X)
# 4. 可视化对比 (省略绘图细节代码)
# 结果会显示：
# - 线性 PCA 只是把 3D 瑞士卷“压扁”在平面上，颜色发生混叠，无法展开。
# - Kernel PCA 能够成功将卷曲的流形“铺平”，不同颜色的点在 2D 平面上呈现清晰的渐变。参数调优提示：在使用
RBF核时，gamma参数（即公式中的 ）至关重要。gamma太大容易过拟合（每个点成为孤岛），太小则退化为线性 PCA。通常需要通过网格搜索（Grid Search）结合交叉验证来寻找最佳gamma。
七、Kernel PCA 在点云处理中的应用
7.1 非线性点云去噪
对于分布在弯曲曲面上的点云，Kernel PCA 可以学习到"沿曲面展开"的低维流形表示。将高维坐标投影到该流形上再重构，可以在保留曲面结构的同时去除噪声。
7.2 非线性特征提取与模式识别
对于具有强烈非线性运动模式或复杂几何拓扑的三维点云数据，Kernel PCA 能比线性 PCA 更紧凑、更本质地表达其内在自由度：
- 三维人体骨架序列分析：人体关节在三维空间中的运动轨迹高度非线性。KPCA 可以将复杂的 3D 骨架动作（如挥手、跑步）降维到 2D/3D 的潜在流形空间，使得动作分类（Action Recognition）变得简单且线性可分。
- 气动/流体点云分析：在分析风洞试验或 CFD（计算流体力学）生成的复杂流场三维点云时，涡流的分布是高度非线性的。KPCA 能够有效提取出主导流场的非线性模态。
7.3 与深度自编码器的比较
| 方法 | 优点 | 缺点 |
|---|---|---|
| Kernel PCA | 有闭式解、数学可证、不需训练 | 核矩阵 内存、σ 需人工选择 |
| 深度自编码器 | 可处理极大样本、自动学习特征 | 需 GPU 训练、超参数多、解释性差 |
八、总结
| 概念 | 标准 PCA | Kernel PCA |
|---|---|---|
| 核心思想 | 寻找方差最大的线性方向 | 在高维特征空间中寻找线性方向 |
| 输入 | 核矩阵 | |
| 计算复杂度 | ||
| 适用数据 | 线性分布 | 非线性流形分布 |
| 核函数 | 不需要 | RBF / 多项式 / 自定义正定核 |
| 新点投影 |
实践建议：对于三维点云处理，如果数据规模较小（）且呈现非线性分布特征（如弯曲管道、螺旋结构），Kernel PCA 是一个数学优美且实用的工具。对于大规模点云，可以考虑使用随机 Fourier 特征近似（Random Kitchen Sinks）来加速核计算。
下一章我们将转向 PCA 在点云处理中的另一个重要应用：基于 PCA 的点云噪声滤波。