2.11. 核主成分分析KernelPCA在一个隐式特征空间里运行 PCA。它不分解你数据的协方差而是分解一个中心化后的核Gram矩阵。这样就能捕捉到线性投影看不到的非线性结构。核 PCA 处理的是样本两两之间的核函数取值。它的操作对象是一个n x n矩阵而不是d x d矩阵。仅这一点就决定了内存开销、求解器的选择以及最主要的限制核 PCA 没有inverse_transform。本节内容基于src/machine_learning/decomposition/kernel_pca.rs、src/machine_learning/types.rs中共享的核类型以及tests/machine_learning/kernel_pca.rs里的集成测试。核相关的这套机制KernelType枚举和Gamma系数和驱动支持向量机的是同一份代码。核 PCA 原样复用了这份代码没有另外定义自己的核类型。2.11.1. 核技巧以及为什么普通 PCA 看不见圆环PCA 找的是输入空间里方差最大的方向。当结构是线性的时候这个办法很好用。当结构不是线性的时候它就不管用了。最经典的失败案例是 2 个同心圆环。内环半径为0.5。外环半径为3.0。仅凭半径就能把这 2 个类别完全分开。没有一条直线能把它们分开所以任何线性投影也做不到。把这些点交给 PCA主成分只会抓住外环在角度方向上的分布。真正区分类别的半径信息始终不会出现。核技巧把每个点x通过一个非线性特征映射phi(x)映射出去。这会把x映射到一个维度高得多的空间。普通 PCA 随后在那个空间里运行。在这个抬升后的空间里圆环可以变得线性可分。核 PCA 从不直接构造phi(x)。这就是这个技巧的关键。抬升空间里的 PCA 只需要形如phi(xi) * phi(xj)的内积。核函数K(xi, xj)直接算出这个内积不需要构造phi。对 RBF 核来说K(x, y) exp(-gamma * ||x - y||^2)。它的隐式特征空间是无穷维的。每个核函数值仍然只是一个标量你可以直接算出来。RBF 的取值只取决于两点之间的距离。这就编码了普通 PCA 会丢弃的半径结构。这正是 RBF 核能劈开这 2 个圆环的原因。2.11.7 有一个可运行的示例。2.11.2. 双重中心化不那么显然的核心抬升空间里的 PCA 需要中心化的特征phi_centered(xi) phi(xi) - (1/n) * sum_k phi(xk)。你没法直接减掉这个均值因为你手上根本没有phi。核 PCA 真正分解的对象是中心化特征的内积矩阵。这些内积完全可以用原始核矩阵K表达出来。展开phi_centered(xi) * phi_centered(xj)就得到双重中心化恒等式Kc[i, j] K[i, j] - row_means[i] - row_means[j] overall_meanrow_means[i]是K第i行的均值。overall_mean是整个矩阵的均值。写成矩阵形式就是Kc H * K * H。H是中心化矩阵H I - (1/n) * J。I是单位矩阵J是n x n的全 1 矩阵。之所以叫双重中心化是因为H在K的两侧同时作用。它先减掉行均值和列均值再把总体均值加回来这样公式就不会把它减掉两次。自己实现核 PCA 时常常会漏掉那个 overall_mean项。漏掉它会让每一次投影都出现偏差。fit实现的就是这套流程。它先算出训练核矩阵的每行均值和整体均值kernel_means然后就地把每个元素改写为K[i,j] - row_mean[i] - row_mean[j] overall_meancenter_kernel_matrix。H * K * H有一个直接推论Kc的每一行求和都为零于是投影结果的每一列均值也为零。测试test_centering_training_output_has_near_zero_column_means验证了每个投影分量的均值都落在零附近1e-9以内。新样本需要另一套不对称的公式。很多简单粗暴的实现就在这里放弃拒绝支持样本外变换。投影一个新样本时它相对训练集的那一行交叉核需要用训练集的统计量来中心化而不是它自己的统计量。center_cross_kernel_matrix减去训练集的行均值和这条新行自身的均值再加上训练集的整体均值。RustyML 实现了这一整套逻辑所以对未见过的数据调用transform也能得到正确结果。参见 2.11.6。2.11.3. 构造估计器构造函数接收核函数和主成分数量。它会校验两者然后返回一个Resultpub fn new(kernel: KernelType, n_components: usize) - ResultSelf, Error pub fn with_eigen_solver(self, eigen_solver: EigenSolver) - Self参数类型含义kernelKernelType核函数及其参数。RustyML 会提前校验它。参见 2.11.4。n_componentsusize保留多少个主成分。必须 0。在拟合时还必须满足 n_samples。n_components 0会立刻被拒绝返回Error::InvalidParameter并带上字段名。n_components n_samples这个关系没法在构造阶段检查因为那时候还不知道样本数。fit会负责这项检查同样返回Error::InvalidParameter。失败的fit不会改动模型。测试test_fit_n_components_greater_than_n_samples_returns_invalid_parameter确认了拟合失败之后那些读取拟合状态的 getter 仍然是None。你不会得到一个改了一半的估计器。特征值求解器默认是EigenSolver::Dense。用链式构建方法with_eigen_solver来设置它。Default实现给你的是一个gamma 0.1的 RBF 核、n_components 2加上 dense 求解器userustyml::machine_learning::decomposition::kernel_pca::{EigenSolver,KernelPCA};userustyml::machine_learning::{Gamma,KernelType};usendarray::array;fnmain(){// 二维空间里的 6 个点。RBF 核保留 2 个主成分。letxarray![[1.0,0.0],[0.0,1.0],[-1.0,0.0],[0.0,-1.0],[2.0,0.5],[-0.5,2.0],];letmutkpcaKernelPCA::new(KernelType::RBF{gamma:Gamma::Value(0.5)},2).unwrap().with_eigen_solver(EigenSolver::Dense);letprojectedkpca.fit_transform(x).unwrap();assert_eq!(projected.nrows(),6);assert_eq!(projected.ncols(),2);// 拟合后的状态通过 getter 暴露出来。println!(kept {} components,kpca.get_n_components());println!(training samples: {:?},kpca.get_n_samples());// Some(6)leteigenvalueskpca.get_eigenvalues().unwrap();println!(leading eigenvalue: {},eigenvalues[0]);}这些 getter 如实反映内部状态。get_kernel、get_n_components、get_eigen_solver按值返回。get_n_samples和get_n_features返回Optionusize拟合前是None。get_eigenvalues和get_eigenvectors分别返回OptionArray1f64和OptionArray2f64。存下来的特征向量是中心化核矩阵特征分解得到的那些列。它们是每个样本对应的系数习惯上写作alpha形状为n_samples x n_components。它们不是 PCA 会给你的那种输入空间方向。核 PCA 没有什么有意义的载荷向量可看。这是在隐式空间里工作的另一面。2.11.4. 核函数与 gamma 的选择KernelType有 5 个变体和 SVC 共享变体公式参数LinearK(x, y) x*y无Poly { degree, gamma, coef0 }(gamma*x*y coef0)^degreedegree: u32 0、gamma: Gamma、coef0: f64RBF { gamma }exp(-gamma*Sigmoid { gamma, coef0 }tanh(gamma*x*y coef0)gamma: Gamma、coef0: f64Cosine(x*y) / (Linear把核 PCA 退化回普通 PCA差别只在中心化约定上。只把它当基线用。RBF是默认选项。当你怀疑数据里有非线性的、基于距离的结构比如圆环时用RBF。Poly捕捉多项式交互。Cosine把模长归一化掉只保留方向。这对高维稀疏数据很有用。Sigmoid不是真正的Mercer核。它中心化后的 Gram 矩阵可能是不定的。这一点会牵涉到 2.11.6 里描述的特征值处理。构造时的校验很严格而且针对每种核各不相同。Poly要求degree 0、gamma为正的有限值、coef0为有限值。RBF要求gamma为正的有限值。Sigmoid只要求它的参数是有限的。Sigmoid接受gamma 0测试test_new_sigmoid_gamma_zero_accepted因为零系数虽然退化仍然是一个合法的 sigmoid。任何违规都会返回Error::InvalidParameter并点名出错的字段。gamma系数的类型是Gamma。它要么是一个显式的值要么是一条在拟合时才解析的、依赖数据的规则Gamma变体解析为何时使用Gamma::Value(v)v你心里已经有具体的带宽。Gamma::Scale1 / (n_features * Var(X))随特征分散程度自适应的默认值scikit-learn 的scale。Gamma::Auto1 / n_features更简单的1/d规则scikit-learn 的auto。fit只解析一次Scale和Auto用的是训练数据的方差和特征数。它会存下解析出来的值这样训练矩阵和之后每一次transform调用用的都是同一个系数。如果数据方差为零所有特征都是常量Gamma::Scale会失败返回Error::InvalidInput因为公式里要除以它。对 RBF 核来说gamma是带宽平方的倒数gamma 1 / (2 * sigma^2)。gamma大带宽小时核只看得见近邻。Gram 矩阵趋近于单位矩阵每个点看起来都极度独特投影会把噪声也拟合进去。gamma小带宽大时任意两点看起来都很相似。Gram 矩阵趋近于常数矩阵主成分什么也抓不到。真正有用的区间落在近邻点的gamma * ||x - y||^2典型值接近 1 的地方。一个实用的起点是Gamma::Scale。从这里出发如果投影看起来像噪声就调小gamma。如果不同的簇挤成了一团就调大gamma。核 PCA 是无监督的所以没有内建的交叉验证可以用来选这个值。测试套件里那个同心圆环的可分性指标class_separability就是可以拿来调参的那类下游信号。2.11.5. 特征值求解器精确与迭代EigenSolver决定 RustyML 如何从中心化核矩阵里提取前n_components个特征对。3 个求解器都是纯 Rust 的自研实现。这一步不依赖任何第二个线性代数库nalgebra只是一个开发依赖只用于测试时的交叉验证。变体策略适用场景Dense默认完整的对称特征分解先做 Householder 三对角化再做隐式位移 QL 迭代也就是经典的 EISPACK/JAMA 算法组合。完整分解之后取靠前的特征对。中小规模的核矩阵你付得起完整O(n^3)分解的开销。Lanczos带完全再正交化的 Krylov 子空间迭代。把问题化简为一个小的三对角问题再用同一个 dense 求解器精确求解。大核矩阵里靠前的少数几个主成分。PowerIteration幂迭代加 Hotelling 收缩一次求一个主成分。最简单的迭代选项。当 Lanczos 显得多余时用它。求解器的选择只影响速度和数值路径不影响结果。Dense、Lanczos、PowerIteration在靠前的特征值上是一致的。它们产出的投影也是一致的至多每列相差一个符号。测试直接验证了这一点test_eigensolver_dense_vs_lanczos_agree和test_eigensolver_dense_vs_power_iteration_agree分别以1e-5和1e-4的容差比较列范数不看符号和最大特征值。符号的不确定性来自特征向量本身不是 bug。如果你需要一个固定的符号自己在下游处理好。迭代求解器省不了内存。3 个求解器操作的都是同一个n x n中心化核矩阵它必须先完整存在任何分解才能开始。当n_components远小于n_samples时Lanczos和PowerIteration能省掉完整分解的O(n^3)开销。但两者都无法避免O(n^2)的矩阵本身。2.11.8 讨论了这个限制。// 从大核矩阵里取少数几个主成分跳过完整的 O(n^3) 分解。 let kpca KernelPCA::new(KernelType::RBF { gamma: Gamma::Scale }, 3) .unwrap() .with_eigen_solver(EigenSolver::Lanczos);核 PCA 的任何求解器都不带随机性。没有种子要设。同一份数据在同一台机器上跑 2 次输出逐位相同。测试test_determinism_dense_solver断言的是完全相等assert_allclose(..., 0.0)。这一点和 t-SNE 不同它的随机初始化确实需要设种子。2.11.6. 拟合、变换与样本外投影这 3 个入口都是固有方法。KernelPCA也实现了 crate 的Fit、Transform、FitTransformtrait。这些 trait 只是转发到固有方法所以泛型代码可以像对待PCA一样对待KernelPCApub fn fitS(mut self, x: ArrayBaseS, Ix2) - Resultmut Self, Error pub fn transformS(self, x: ArrayBaseS, Ix2) - ResultArray2f64, Error pub fn fit_transformS(mut self, x: ArrayBaseS, Ix2) - ResultArray2f64, Errorfit至少需要 2 个样本。只有 1 行会返回Error::InvalidInput。0 行会返回Error::EmptyInput。非有限的输入会返回Error::NonFinite。fit会解析gamma构建并中心化训练核矩阵提取特征对然后把后续 transform 需要的一切都存起来。这包括完整训练矩阵的一份副本每次transform调用都会复用它。transform能投影任何与训练数据特征数相同的矩阵包括模型从未见过的数据。它会在新样本和存下来的训练样本之间构建交叉核矩阵。它用训练集的统计量把这个矩阵中心化参见 2.11.2再把它投影到存下来的特征向量上。第k个主成分上的投影坐标是(Kc * v_k) / sqrt(lambda_k)。这个1/sqrt(lambda)缩放把原始特征向量变成了归一化恰当的主成分。在fit之前调用transform会返回Error::NotFitted。特征数不匹配会返回Error::DimensionMismatch。userustyml::machine_learning::decomposition::kernel_pca::KernelPCA;userustyml::machine_learning::{Gamma,KernelType};usendarray::array;fnmain(){letx_trainarray![[1.0,0.0],[0.0,1.0],[-1.0,0.0],[0.0,-1.0],[2.0,0.5],[-0.5,2.0],[1.5,-1.5],[-2.0,1.0],];// 默认求解器就是 Dense。不用调用构建方法。letmutkpcaKernelPCA::new(KernelType::RBF{gamma:Gamma::Value(0.5)},2).unwrap();kpca.fit(x_train).unwrap();// 模型从没见过的点特征数相同。letx_newarray![[0.3,0.3],[1.8,-1.2]];letprojected_newkpca.transform(x_new).unwrap();assert_eq!(projected_new.nrows(),2);assert_eq!(projected_new.ncols(),2);println!({projected_new:?});}fit_transform是一个快捷方式。它先拟合模型再对同一份数据做变换。测试test_fit_transform_equals_fit_then_transform确认它和分 2 步调用的结果吻合到1e-10。fit_transform只是写法上的便利不是性能优化。它内部就是先调用fit再对同一个矩阵调用transform。transform无论如何都会从头重建一次核矩阵所以开销和分开调用完全一样。当你需要把另一个矩阵投影过一个已经拟合好的模型时用分 2 步的写法。对于作用在互不相同的点上的合规 Mercer 核中心化后的 Gram 矩阵是半正定的。每个保留下来的特征值都严格为正测试test_eigenvalues_are_positive_after_fit也确认了它们按降序排列。不过核 PCA 不会因为出现非正特征值就直接判定失败。中心化 Gram 矩阵的半正定性只在舍入误差范围内成立。像Sigmoid这样的非 Mercer 核会产生真正为负的末尾特征值。fit只拒绝非有限的特征值NaN 或 Inf 会映射为Error::Computation而不会拒绝整次拟合。任何特征值算不上有意义的正数低于相对阈值1e-12 * lambda_max的主成分都会得到0.0的投影缩放。这会把那一列清零而不是产生Inf或NaN。这样既保住了你要求的n_components维度又悄悄丢掉了退化方向。测试test_fit_indefinite_kernel_negative_eigenvalue_is_tolerated用一个 Sigmoid 核驱动这条路径。它验证了出问题的那一列全为零其余部分保持有限。2.11.7. 实战示例分离同心圆环这正是普通 PCA 搞不定的情形。有 2 个圆环按半径可分却在线性意义上纠缠在一起。跑一遍 RBF 核这 2 个类别就会落进主成分空间里可以区分的区域userustyml::machine_learning::decomposition::kernel_pca::KernelPCA;userustyml::machine_learning::{Gamma,KernelType};usendarray::Array2;usestd::f64::consts::PI;fnmain(){// 内环 r 0.5外环 r 3.0。没有直线能把它们分开。letn12;letmutdata:Vecf64Vec::new();foriin0..n{leta2.0*PI*iasf64/nasf64;data.push(0.5*a.cos());data.push(0.5*a.sin());}foriin0..n{leta2.0*PI*iasf64/nasf64;data.push(3.0*a.cos());data.push(3.0*a.sin());}letxArray2::from_shape_vec((2*n,2),data).unwrap();letmutkpcaKernelPCA::new(KernelType::RBF{gamma:Gamma::Value(0.5)},2).unwrap();letprojkpca.fit_transform(x).unwrap();// RBF 核编码了半径距离所以圆环会沿着某个主成分分开。letinner_mean:f64(0..n).map(|i|proj[[i,0]]).sum::f64()/nasf64;letouter_mean:f64(n..2*n).map(|i|proj[[i,0]]).sum::f64()/nasf64;println!(inner-ring mean of component 0: {inner_mean:.4});println!(outer-ring mean of component 0: {outer_mean:.4});println!(gap between ring means: {:.4},(inner_mean-outer_mean).abs());}把KernelType::RBF { .. }换成KernelType::Linear2 个环的均值之间的间隔就会塌陷。线性投影会被外环的角度变化主导。它从来不会编码半径。测试test_rbf_separates_radial_clusters_better_than_linear用一个 Fisher 风格的可分性分数把这一点量化了。它断言 RBF 投影以一个可观的优势胜过线性投影。2.11.8. 核 PCA 做不到的事核 PCA 没有inverse_transform。普通 PCA 有。你可以把一个低维编码映射回输入空间因为它的投影是一个带有干净转置的线性映射。核 PCA 做不到这一点。这不是疏漏。这是原像问题pre-image problem。投影后的点活在隐式特征空间里。要把它反过来你需要找到一个输入x让它的特征映射phi(x)恰好落在那个位置。对多数核来说尤其是 RBF特征映射是非线性的、无穷维的而且不满射。特征空间里任取一点通常没有精确的原像。它只有近似的原像需要靠另一套非线性优化才能求出来。RustyML 没有附带这套近似算法。核 PCA 严格来说是一个前向的、单向的投影。把它用于可视化、用投影做降噪或者当作喂给下游分类器的非线性特征环节。不要用它来做重建。Gram 矩阵是O(n^2)的这才是真正的天花板。fit会构建一个n x n的f64矩阵。内存按8 * n^2字节增长和特征数无关。这大约是n 10,000时的 800 MB以及n 20,000时的 3.2 GB。时间开销更糟。构建矩阵是O(n^2 * d)通过一次并行 GEMM 完成。dense 特征分解是O(n^3)。当你只需要少数几个主成分时换成Lanczos或PowerIteration能削减分解的开销。但没有任何办法能去掉O(n^2)的矩阵本身。实际用起来核 PCA 在几千个样本这个量级上还比较从容。到了几万这个量级就开始吃力了。超过这个规模之后抽一个有代表性的子集去拟合再对其余数据调用transform每个新 batch 都要为此付出O(m * n * d)。或者换一种从不构造完整核矩阵的方法。每一次transform调用都拖着整个训练集。投影是相对存下来的训练样本定义的所以对m个新点transform会重建一个m x n的交叉核矩阵。PCA 的变换开销和训练集规模无关。核 PCA 的变换开销却会随n永远增长下去。做预算时要把这一点算进去。并行机制会在超过内部尺寸门限时自动启动门限是按核矩阵的元素个数设定的。粗略地说n^2越过几十万个元素之后中心化的扫描会开始并行。越过几百万个元素之后逐元素的中心化也会并行。核 GEMM 有它自己的 FLOPs 门限。你不需要逐次调用去配置这些。可调的门限见 7.3. 性能调优与并行。2.11.9. 持久化KernelPCA派生了Serialize和Deserialize并暴露出标准的那一对方法pub fn save_to_path(self, path: str) - Result(), Error pub fn load_from_path(path: str) - ResultSelf, Error序列化用的是紧凑的 postcard 二进制格式。路径里的.bin、.dat或者任何其他扩展名都只是文件名的一部分。字节内容永远是二进制的。往返一趟之后的模型能原样重现transform的输出。测试test_save_load_round_trip断言相等到1e-12。留意一下什么内容被序列化了。一个拟合好的核 PCA会连同特征向量和中心化统计量一起存下整个训练矩阵因为transform全都要用到它们。保存下来的文件会随着训练集一起变大。这是同一套O(n^2)、存样本设计的又一个后果。在你持久化一个在大语料上拟合的模型之前记住这一点。userustyml::machine_learning::decomposition::kernel_pca::KernelPCA;userustyml::machine_learning::{Gamma,KernelType};usendarray::array;usestd::fs;fnmain(){letxarray![[1.0,0.0],[0.0,1.0],[-1.0,0.0],[0.0,-1.0],[2.0,0.5],[-0.5,2.0],[1.5,-1.5],[-2.0,1.0],];letmutkpcaKernelPCA::new(KernelType::RBF{gamma:Gamma::Value(0.5)},2).unwrap();kpca.fit(x).unwrap();letbeforekpca.transform(x).unwrap();letpathkpca_model.bin;kpca.save_to_path(path).unwrap();letloadedKernelPCA::load_from_path(path).unwrap();letafterloaded.transform(x).unwrap();assert_eq!(before.shape(),after.shape());fs::remove_file(path).unwrap();}文件不存在时会以Error::Io的形式出现测试test_load_from_nonexistent_path_returns_io_error。关于错误的整体分类见 1.6. 错误处理。关于整个 crate 的持久化模式见 7.2. 深入模型持久化。