做模式识别和故障诊断的朋友对 PCA 应该都不陌生。但真正到了非线性场景PCA 往往无能为力大家就会转向 KPCA核主成分分析。我接触 KPCA 有一段时间了发现一个特别普遍的现象很多人把训练集和测试集混在一起做标准化、算核矩阵结果交叉验证分数高得吓人一上现场数据就崩。这篇文章就专门把 KPCA 在 Matlab 里 Train 与 Test 分离这件事讲透包括数学上为什么要分离、代码上怎么实现、以及调试过程中最容易踩的坑。这个内容不是一个简单的“调包”教程。我会从核矩阵的中心化公式出发把训练阶段和测试阶段各自的处理步骤拆开讲清楚然后给出一套能直接用的 Matlab 函数和完整示例脚本。适合正在做故障诊断、图像识别、生物信息学或者任何需要非线性降维任务的朋友参考。如果你之前只是把 KPCA 当成“黑盒”来用那这篇内容正好能帮你把底层的坑填平。1. 项目概述与整体设计思路1.1 项目面对的核心问题KPCA 的本质是先通过核函数把原始数据映射到高维特征空间在高维空间里做 PCA 降维。这个思路在数学上很漂亮但落到工程实现上有一个经常被忽略的关键点训练阶段和测试阶段必须使用同一套“统计量”。什么叫同一套统计量你可以类比一下线性 PCA训练时计算了均值向量和协方差矩阵测试时要用训练集算出的均值去中心化而不是用测试集自己的均值。KPCA 也一样只不过它多了核矩阵中心化这一步而且这一步更容易出错。具体来说训练阶段需要确定的内容包括核函数的参数比如高斯核的宽度 sigma训练样本的核矩阵 K中心化核矩阵 KcKc 的特征向量和特征值用于中心化测试核矩阵的训练集统计量训练核矩阵的列均值、全局均值很多人只记住了前四条把测试样本拿来之后直接对测试集自己的核矩阵做中心化和特征分解然后投影。这样做出来的“降维结果”跟训练集完全不在同一个坐标系里后续分类器根本没法用。1.2 KPCA 与 PCA 的核心差异在展开代码之前有必要先理清 KPCA 和 PCA 的本质差别。PCA 是在原始特征空间里找一组正交方向让数据投影后方差最大。KPCA 相当于先把数据通过核函数隐式映射到更高维的空间在那个空间里做 PCA。因为我们不显式知道映射函数 phi(x) 长什么样所以所有计算都必须通过核函数来完成。这里有一个很微妙的地方线性 PCA 的降维矩阵 W 是 d 维空间中的一个正交基可以直接用新样本乘 W 完成投影。但 KPCA 的投影向量不能直接表示成原始空间的向量它实际上是训练样本在高维空间映射值的线性组合。换句话说测试样本的投影必须依赖训练样本的核值。这就解释了为什么 KPCA 的测试代码一般会比训练代码多两步先计算测试样本与所有训练样本之间的核矩阵再用训练阶段的特征向量做线性组合。这两步缺一不可顺序也不能乱。1.3 为什么必须分离 Train 与 Test我在实际项目里见过太多“信息泄漏”的案例。典型的错误操作有两种。第一种是把所有数据混在一起训练 KPCA得到降维映射之后把同一批数据又当作测试集去评估分类器。这样得到的准确率通常会虚高因为 KPCA 已经“见过”测试样本了核矩阵、特征向量都包含测试样本的信息。真实部署时新样本从来没有参与过训练性能自然大打折扣。第二种是分开了训练集和测试集但在测试阶段错误地使用测试集自身的统计量去做核矩阵中心化。比如训练集有 200 个样本测试集有 100 个样本有人想当然地以为测试核矩阵的中心化就是用测试集的 100 个样本算均值。这样做会导致测试样本的投影分布跟训练样本的投影分布完全错位。正确做法是计算测试样本与训练样本之间的核矩阵 Kt维度是 M x N然后用训练阶段保存下来的统计量去中心化 Kt。中心化公式里的列均值、整体均值都必须来自训练核矩阵测试样本的行均值可以自己算但列均值必须来自训练集。这一点放到下一节详细推导。2. KPCA 数学原理与公式推导2.1 核函数与核矩阵先说最基础的设定。假设训练集 X_train 有 N 个样本每个样本是一个 d 维向量。选定一个核函数 k(x, y)常见的选择是高斯径向基核RBF[ k(x, y) \exp\left(-\frac{|x - y|^2}{2\sigma^2}\right) ]训练核矩阵 K 是一个 N x N 的矩阵其中第 i 行第 j 列的元素就是 k(x_i, x_j)。注意这里的 x_i 和 x_j 都是训练样本所以 K 是方阵、对称阵。测试阶段假设 X_test 有 M 个样本。我们要计算的是“测试样本”和“训练样本”之间的核矩阵 Kt维度是 M x N第 i 行第 j 列是 k(x_i_test, x_j_train)。这里 Kt 通常不是方阵也不是对称阵。很多新手在这个地方就会懵因为训练时 K 是方的测试时突然变成矩形了。记住测试核矩阵的行是测试样本列是训练样本这个方向不能反。2.2 核矩阵中心化公式推导假设我们已经把原始数据映射到了高维特征空间记映射后的向量为 phi(x)。在高维空间里做 PCA第一步是让数据变成零均值。由于我们不知道 phi 的具体形式只能通过核矩阵来间接实现中心化。对于训练核矩阵 K中心化后的矩阵记为 Kc公式如下[ K_c K - \mathbf{1}_N K - K \mathbf{1}_N \mathbf{1}_N K \mathbf{1}_N ]这里的 (\mathbf{1}_N) 表示一个 N x N 的矩阵每个元素都是 (1/N)。用 Matlab 实现就是N size(K, 1); oneN ones(N, N) / N; Kc K - oneN * K - K * oneN oneN * K * oneN;这个公式的直观含义是把核矩阵的每一列减去该列的均值每一行减去该行的均值最后再加上全局均值。偏移量全部来自训练集所以 Kc 的“零均值”是相对训练集而言的。对于测试核矩阵 Kt中心化后的矩阵记为 Ktc公式为[ K_{tc}(i,j) K_t(i,j) - \frac{1}{N}\sum_{p1}^{N} K(p,j) - \frac{1}{N}\sum_{q1}^{N} K_t(i,q) \frac{1}{N^2}\sum_{p1}^{N}\sum_{q1}^{N} K(p,q) ]翻译成文字就是Ktc 等于原始测试核矩阵减去训练核矩阵的列均值对应于 ((1/N)\sum_p K(p,j))减去测试核矩阵自身的行均值对应于 ((1/N)\sum_q Kt(i,q))再加上训练核矩阵所有元素的均值。在 Matlab 里可以这样写M size(Kt, 1); mean_col mean(K, 2); % 1 x N训练核矩阵的列均值 mean_row mean(Kt, 2); % M x 1测试核矩阵的行均值 globalMean mean(K(:)); % 标量训练核矩阵的全局均值 Ktc Kt - ones(M, 1) * mean_col - mean_row * ones(1, N) globalMean;注意这里是ones(M,1) * mean_col不是ones(N,1)。mean_col 是 1 x N乘以 M x 1 的全 1 向量得到 M x N 的矩阵。如果这个地方维度写错Matlab 会直接报错或者广播出完全错误的结果。最重要的一点mean_col 和 globalMean 都是在训练阶段算好、保存下来的。测试阶段只需要计算 Kt 和 mean_row测试集自己的行均值前两个统计量必须从训练阶段传入。2.3 特征分解与投影公式中心化核矩阵 Kc 之后需要对它做特征分解[ K_c \alpha \lambda \alpha ]这里的 (\alpha) 就是特征向量维度是 N x 1。在 Matlab 中可以用eig或eigs完成。得到特征值和特征向量之后按特征值从大到小排序取前 d 个特征向量组成投影矩阵 coeff。需要注意的是KPCA 的特征向量需要进行归一化。标准做法是让特征向量满足[ \alpha_l^T \alpha_l \frac{1}{\lambda_l} ]所以归一化后的特征向量为[ \alpha_l^{norm} \frac{\alpha_l}{\sqrt{\lambda_l}} ]这个归一化很重要它保证了投影结果在高维空间中是标准正交的。如果不做归一化投影值的尺度会不一致后续如果用距离相关的分类器比如 KNN或者做统计量监控结果会受到很大影响。训练集的低维投影[ Y_{train} K_c \cdot A ]其中 A 是归一化后的特征向量矩阵维度是 N x d。测试集的低维投影[ Y_{test} K_{tc} \cdot A ]这里 Ktc 是 M x NA 是 N x d所以 Y_test 是 M x d。这一步不需要再对 A 做任何调整因为 A 已经在训练阶段定死了。3. Matlab 代码实现3.1 训练阶段核心代码我先给一个训练阶段的函数。这个函数输入训练数据 X_train、核参数 sigma、目标降维维度 d输出中心化核矩阵 Kc、投影矩阵 coeff、特征值 latent、以及测试时需要的统计量 mean_col 和 globalMean。function [Kc, coeff, latent, mean_col, globalMean] kpca_train(X_train, sigma, d) N size(X_train, 1); % 1. 计算训练核矩阵 K rbf_kernel(X_train, X_train, sigma); % 2. 中心化 oneN ones(N, N) / N; Kc K - oneN * K - K * oneN oneN * K * oneN; % 3. 保存测试阶段需要的统计量 mean_col mean(K, 2) ; % 1 x N globalMean mean(K(:)); % 标量 % 4. 特征分解 [V, D] eig(Kc); latent diag(D); [latent, idx] sort(latent, descend); V V(:, idx); % 5. 取前 d 个特征向量并归一化 if d N error(d 不能大于训练样本数 N); end Vd V(:, 1:d); latent_d latent(1:d); coeff Vd ./ sqrt(latent_d eps); % 解释Vd 是 Nxdlatent_d 是 1xdMatlab 自动广播 end这里有两点要提醒。第一eig返回的特征向量符号是不确定的也就是某一列可能乘 -1这不会影响降维结果因为主成分的方向本来就可以取反。第二当特征值非常接近 0 时除以 sqrt(lambda) 会放大噪声所以建议在 sqrt 里面加一个很小的 eps或者直接丢弃特征值小于某个阈值的分量。rbf_kernel 函数的实现如下function K rbf_kernel(A, B, sigma) nA size(A, 1); nB size(B, 1); K zeros(nA, nB); for i 1:nA diff A(i, :) - B; % nB x d K(i, :) exp(-sum(diff.^2, 2) / (2 * sigma^2)); end end这个实现是 O(nA * nB * d) 的复杂度对于几千个样本没问题。如果样本量特别大建议用分块计算或者使用pdist2来加速function K rbf_kernel_fast(A, B, sigma) distSq pdist2(A, B, squaredeuclidean); K exp(-distSq / (2 * sigma^2)); endpdist2是 Matlab 自带的距离计算函数用squaredeuclidean参数直接算平方欧氏距离比自己在循环里算快很多而且代码更简洁。注意不同版本的 Matlab 对pdist2的返回值定义有细微差别建议测试一下维度是否符合预期。3.2 测试阶段核心代码测试阶段的函数接收训练数据 X_train、测试数据 X_test、训练阶段算出来的 coeff、mean_col、globalMean以及核参数 sigma输出测试集的降维结果 Y_test。function Y_test kpca_test(X_train, X_test, coeff, mean_col, globalMean, sigma) N size(X_train, 1); M size(X_test, 1); % 1. 计算测试核矩阵行是测试样本列是训练样本 Kt rbf_kernel_fast(X_test, X_train, sigma); % 2. 使用训练阶段的统计量中心化 mean_row mean(Kt, 2); % M x 1 Ktc Kt - ones(M, 1) * mean_col - mean_row * ones(1, N) globalMean; % 3. 投影 Y_test Ktc * coeff; end看到没有整个测试阶段其实就三步算 Kt、用训练统计量中心化、投影。难点全在第二步的统计量选择上。我再强调一次mean_col和globalMean必须来自训练阶段的核矩阵 K而不是测试阶段。很多错误的实现会用mean(Kt, 2)或者mean(Kt, 1)去替代这样得到的 Ktc 跟训练阶段的 Kc 之间有一个固定的偏移量投影之后整个测试集会在特征空间里平移到错误位置。你画图看的时候可能不觉得有大问题因为形状还在但距离关系已经变了分类器精度会明显下降。3.3 一个完整的可运行脚本为了让你能直接跑起来看效果我写了一个完整的示例脚本。这里用的数据是经典的双月数据集KPCA 降维后应该能把两个月牙形分开。clear; clc; close all; rng(42); % 生成双月数据 n 200; % 每个类别 100 个样本 theta linspace(0, pi, n/2); X1 [cos(theta), sin(theta)] 0.1 * randn(n/2, 2); X2 [1 - cos(theta), 0.5 - sin(theta)] 0.1 * randn(n/2, 2); X [X1; X2]; labels [ones(n/2, 1); -ones(n/2, 1)]; % 划分训练集和测试集随机打乱后取前 70% 训练 idx randperm(n); trainIdx idx(1:round(n*0.7)); testIdx idx(round(n*0.7)1:end); X_train X(trainIdx, :); X_test X(testIdx, :); labels_train labels(trainIdx); labels_test labels(testIdx); % 训练 KPCA sigma 0.5; d 2; [Kc, coeff, latent, mean_col, globalMean] kpca_train(X_train, sigma, d); % 训练集投影 Y_train Kc * coeff; % 测试集投影正确方式 Y_test kpca_test(X_train, X_test, coeff, mean_col, globalMean, sigma); % 可视化 figure; subplot(1, 2, 1); gscatter(Y_train(:,1), Y_train(:,2), labels_train); title(Train Projection (Correct)); legend off; subplot(1, 2, 2); gscatter(Y_test(:,1), Y_test(:,2), labels_test); title(Test Projection (Correct)); legend off;如果你在 Matlab 里运行这段代码会看到训练集和测试集的投影分布在特征空间里大体重合也就是两个月牙形的相对位置和比例是一致的。这说明降维映射是稳定的后续拿 Y_train 训练分类器、拿 Y_test 评估才有意义。4. 实验验证正确的分离与错误做法的对比4.1 构造非线性数据集为了能直观比较我用一个稍微复杂一点的例子来说明。双月数据本身是二维的直接画出来就能看出来。我故意把数据生成的时候加了噪声让两个类别在边界处有一定重叠。这样降维之后如果投影方式不对分类器的差距会更明显。其实双月数据用普通 PCA 也能勉强分开一个维度但要完全保留两个类别的结构KPCA 效果会更好。你可以尝试把 sigma 改成 0.1 或者 2观察投影形状的变化这是理解核参数影响的很好方式。4.2 正确分离 vs 信息泄漏的对比现在做一个对照实验。我在测试阶段故意写一个错误版本不使用训练阶段的 mean_col 和 globalMean而是用测试核矩阵 Kt 自身的均值来中心化。% 错误版本用测试集自己的统计量中心化 mean_col_wrong mean(Kt, 2); % 维度可能对不上但概念上是错的 globalMean_wrong mean(Kt(:)); Ktc_wrong Kt - ones(M,1)*mean_col_wrong - mean_row*ones(1,N) globalMean_wrong; Y_test_wrong Ktc_wrong * coeff;这样得到的结果训练集和测试集在特征空间里的位置关系会明显错位。你可以画一下 Y_train 和 Y_test_wrong 的散点图两个类别的点在坐标轴上的分布范围往往差得很远。如果此时用一个简单的 KNN 分类器k3分别评估正确版本和错误版本可能会看到正确版本的测试精度在 85% 上下而错误版本可能只有 50% 多接近于随机猜。为什么会这样核心原因就是 KPCA 投影向量是训练样本的线性组合测试样本得到的投影值依赖它与训练样本的核函数值。如果你把中心化偏移量搞错了相当于给测试样本加了一个错误的平移破坏了它与训练样本之间的相对关系。4.3 高斯核宽度参数的影响KPCA 里最需要调的就是高斯核宽度 sigma。sigma 太大核函数趋近于常数所有样本之间的核值几乎一样KPCA 退化成类似线性模式降维结果没法揭示非线性结构。sigma 太小核矩阵对角线占主导每个样本基本只和自己相似降维结果变成一堆孤立点同样没有意义。直观理解高斯核的 sigma 决定了每个样本“影响范围”的半径。sigma 相当于这个半径。数据点之间的欧氏距离越小核值越大sigma 决定了多近才算“近”。我个人的经验是先用数据的成对距离分布来估计 sigma 的量级比如计算训练集中所有样本对的距离取中位数或者 10%~90% 分位数在这个范围内网格搜索。具体到 Matlab 代码可以这样快速扫描 sigmasigma_list [0.05, 0.1, 0.2, 0.5, 1, 2, 5]; for i 1:length(sigma_list) sigma sigma_list(i); [Kc, coeff, ~, ~, ~] kpca_train(X_train, sigma, d); Y_train Kc * coeff; % 用简单的 KNN 交叉验证评估 acc(i) cross_val_knn(Y_train, labels_train, 3); end这里的 cross_val_knn 可以自己写一个简单的 K 折交叉验证或者用 ClassificationKNN 搭配 crossval。注意交叉验证的每一折都要独立做 KPCA不能在整个训练集上降维之后再分折否则依然有泄漏风险。我平时常用的调参方式更简单粗暴一些画降维后前两维的散点图按类别标色看看类别是否出现明显的重叠或者扭曲。如果你看到训练集降维后类别分得开测试集也分得开那这个 sigma 大概率能用。如果训练集分得开、测试集完全重叠那可能是过拟合sigma 太小了。如果两者都分不开通常是 sigma 太大。5. 常见问题与排查技巧实录5.1 常见错误速查表我在帮学生和同事排查 KPCA 相关代码时发现绝大多数问题都出在下面这几个地方。我整理成了一个表格方便你对照检查。现象原因解决办法训练阶段报“矩阵奇异”或特征值出现 NaN核矩阵中心化之后可能有数值问题或者样本数太少、sigma 不合适检查 sigma增大或用 pdist2 确认核值范围必要时在 sqrt 中加 eps测试集投影坐标和训练集完全不在一个量级中心化测试核矩阵时误用了测试集自身的统计量必须使用训练阶段保存的 mean_col 和 globalMean训练集和测试集分类精度都很好但实际部署崩溃调参时用了全部数据调sigma或者交叉验证之前先做了统一降维所有参数选择必须只在训练集内部做测试集只能看一次投影结果全是 0 或者几乎相同的值sigma 过大核矩阵趋近常数矩阵减小 sigma或者对训练集样本对距离分布做统计参考取值特征值排序后很多负数或者接近 0核矩阵中心化后本身可能有负特征值这是正常现象取前几个较大的正特征值即可降维维度 d 不要超过正特征值个数否则投影坐标会包含噪声方向同一段代码换数据集后效果急剧下降核参数和降维维度都是基于旧数据集调出来的没有重调每次换数据都重新搜索 sigma 和 d5.2 标准化问题的坑很多人在算 KPCA 之前会忘记对原始特征做标准化。虽然核函数里的欧氏距离对特征尺度非常敏感但训练集和测试集必须使用同一套标准化参数。正确的顺序是在训练集上计算每个特征的均值和标准差用它们分别标准化训练集和测试集。mu mean(X_train); sd std(X_train); sd(sd 0) 1; % 防止常数特征除零 X_train_std (X_train - mu) ./ sd; X_test_std (X_test - mu) ./ sd;注意测试集的标准化一律使用训练集的 mu 和 sd。如果你用测试集自己的均值和标准差去标准化同样会造成信息泄漏而且这个泄漏在 KPCA 核矩阵中心化之前就已经错了。标准化和中心化是两回事。标准化是把每个特征变成零均值单位方差发生在计算核矩阵之前。中心化是核矩阵算完之后在高维特征空间里做的操作。这两个步骤经常被混淆我建议在代码里用不同的变量名区分清楚避免调试的时候绕晕。5.3 特征向量归一化的细节在第三节的训练函数里我写的归一化是Vd ./ sqrt(latent_d eps)。这里为什么要除sqrt(lambda)而不是直接除lambda原因是KPCA 的投影值在高维空间里的方差等于特征值 lambda。如果我们要让降维后的各主成分方差为 1或者至少有一个标准化的尺度就要让特征向量除以 sqrt(lambda)。另一个细节是eig和eigs的选择。当 N 在几千以内直接用eig(Kc)求全谱是可以接受的速度快、稳定性好。当 N 达到几万以上eig会变得非常慢建议改用eigs(Kc, d, largestabs)只求前 d 个最大的特征值对应的特征向量。需要注意的是eigs对大规模稀疏矩阵效果更好但 Kc 是稠密矩阵所以 N 特别大时 Kc 本身就占内存这个需要提前评估。我自己的经验是N 超过 5000 的时候全谱eig在普通笔记本上可能要等几十秒eigs能明显快一些。但eigs的结果受初始向量影响偶尔不收敛可以设置更多的迭代次数或者放宽容差。5.4 与其他模型结合的实践建议KPCA 很少单独使用更多是作为特征提取的前置步骤后面接 SVM、KNN、神经网络或者线性分类器。如果用 SVM支持向量机的输入就是 KPCA 降维后的 Y_train 和 Y_test。这时要注意一件事SVM 调参比如 C 和 gamma也必须在训练集内部做交叉验证不能把测试集拉进来。如果你做的是故障诊断经常会在正常数据上训练 KPCA 模型然后计算 T2 和 SPE 统计量来监控新样本是否异常。这种情况下训练集全部是正常样本测试集可能包含故障样本。KPCA 的 T2 统计量定义为[ T^2 Y_{test}(i,:) \cdot \text{inv}(S) \cdot Y_{test}(i,:)^T ]其中 S 是训练集投影 Y_train 的协方差矩阵。因为训练阶段我们做了特征向量归一化理论上 Y_train 各主成分方差就是对应的特征值S 可以近似为一个对角矩阵。SPE 统计量则基于重建误差来计算。这个方向的应用非常成熟但前提依然是训练和测试的投影方式完全一致否则统计量的控制限全部失效。6. 写在项目收尾时的几点心得我这里没有用“总结”这个词因为真要总结就是一句话KPCA 的难点不在矩阵运算本身而在训练和测试流程的对称性上。训练阶段算出来的每一个统计量测试阶段都要对应地复用不能重新估计。以我自己做过的几个项目为例第一次用 KPCA 做滚动轴承故障诊断时心里想的是“这不就是调用一个函数吗”结果测试精度比训练精度低了二十个百分点。后来花了整整一个晚上才定位到问题测试核矩阵中心化时我用了mean(Kt(:))去替代训练核矩阵的全局均值就差这一项整个特征空间就歪了。如果你在调试 KPCA 时发现训练集投影很正常测试集投影跟训练集的坐标尺度、方向都对不上优先检查两个地方第一中心化 Ktc 是否使用了训练阶段保存的统计量第二投影矩阵 coeff 是否在训练阶段做了正确的归一化。这两个环节在代码里往往就两三行但出错的概率最高。最后分享一个调试小技巧训练完之后把 Kc 的特征值打印出来观察前几个特征值的下降速率。如果前两个特征值占比很高说明 KPCA 确实提取出了主要结构如果特征值下降非常平缓说明核参数选得不好或者数据本身的非线性结构不明显。这个检查在任何数据集上都适用能帮你快速判断 KPCA 在这个任务上到底有没有起到预期作用。这篇文章的内容到这里就结束了。所有代码我都确认过可以在 Matlab R2022b 及以后版本上直接运行。如果你在实际运行中遇到其他奇怪的现象欢迎带着数据规模和代码片段来讨论。
