写KPCA的Matlab实现最让我头疼的不是核函数怎么写而是训练集和测试集之间那套“看着一样、实际完全不同”的预处理逻辑。做故障诊断和人脸识别的人应该都有同感PCA的流程闭着眼都能背换成KPCA之后稍不留神就把测试样本单独做了中心化结果模型输出一眼假。这篇文章我就把整套Train与Test分离的细节掰开讲清楚代码可以直接抄注释里全是实测过的坑。1. KPCA和PCA到底差在哪核函数把问题搬到高维空间1.1 从PCA的“线性假设”说起主成分分析的核心动作是对中心化后的数据矩阵做协方差分解找到方差最大的投影方向。用到的无非是 XX 的特征分解得到线性载荷向量 W然后计算得分 T XW。这整套逻辑非常自洽前提是数据在原始空间里大致服从线性结构。但现实中的数据很少这么听话。工业过程变量之间的耦合、图像像素之间的关联、生物信号里的非线性漂移几乎都是弯弯曲曲的结构。强行套PCA结果就是主成分上保留的“方差”很大可这些方差对应的方向并不能反映真实的流形结构。这时候就需要核主成分分析Kernel PCAKPCA出场了。1.2 核技巧不需要知道φ(x)长什么样KPCA的思路很直接先用一个非线性映射 φ 把原始样本 x 映射到高维特征空间 F在 F 里再做PCA。可是 φ 的显式表达式通常写不出来特征空间动辄几千维直接计算映射后的坐标既不现实也没必要。核技巧在这里给了个近乎作弊的解决方案。因为我们真正需要的不是 φ(x) 本身而是高维空间里的内积 φ(x_i), φ(x_j)。只要这个内积能用某个核函数 k(x_i, x_j) 直接算出来那么PCA过程中的所有计算——协方差分解、投影、重构——都可以只用内积完成。实践里最常用的是高斯核RBF核k(x_i, x_j) exp( -||x_i - x_j||^2 / (2σ^2) )这个核函数隐含了一个无穷维的特征空间但计算代价不过是一个欧氏距离加一次指数运算。我第一次接触的时候觉得这简直像魔术后来才意识到它只是把“显式映射”变成了“内积等价替换”数学上没有任何魔法工程上却大大简化了实现。1.3 中心化其实是在特征空间里做的PCA里先对原始数据减均值KPCA里也一样要中心化区别在于减的是 φ(x) 在特征空间里的均值。因为 φ 的显式形式未知我们不能像PCA那样先算均值再减而是直接对核矩阵做变换。训练集核矩阵 K 的定义是K(i,j) k(x_i, x_j)中心化后的核矩阵记为 Kc计算公式为Kc K - (1/n)·1 1 K - (1/n)·K 1 1 (1/n²)·(1 K 1)·1 1其中 1 是 n 维全1列向量。用Matlab实现时我不建议写成逐个元素的循环直接用矩阵运算效率高得多function Kc center_kernel_matrix(K) n size(K, 1); OneN ones(n, n) / n; Kc K - OneN * K - K * OneN OneN * K * OneN; end这段代码里 OneN 就是 (1/n)·1 1整体对应公式里的 H I - (1/n)·1 1 作用于 K 的两侧。注意一个细节中心化必须作用在核矩阵上而不是先把原始数据标准化再算核矩阵否则特征空间里的均值偏移没有被真正修正后续所有投影都会有系统性误差。还有一个容易忽略的点对称性。理论上 K 和 Kc 都是对称矩阵但浮点运算后可能出现微小不对称。后续做特征分解时建议先 (KcKc)/2 做个对称化省得特征向量出现虚部。2. Train与Test为什么要分开先说数据泄漏的后果2.1 数据泄漏的代价任何降维方法只要涉及“利用数据本身计算统计量”就存在数据泄漏的风险。PCA里泄漏的是均值向量和协方差结构KPCA里泄漏的东西更多训练核矩阵的均值、特征向量、特征值、核参数σ每一个都是在训练集上估计得到的。如果测试阶段把这些统计量重新算一遍会出现什么后果直观地说就是测试样本被映射到一套“属于它自己”的坐标系里而不是训练阶段确定的坐标系。训练时学到的分类器边界、监控阈值全部失效。更隐蔽的问题在于评估指标虚高你用测试集参与了坐标系构建相当于把测试集信息偷渡进了模型CV分数再好看都不可信。我见过一个真实的案例。某人做工业过程故障诊断训练集和测试集来自同一个反应釜的不同批次数据量都不小。他把KPCA的测试过程写成了“独立中心化独立特征分解”结果T²和SPE统计量的正常波动范围比训练阶段小了两个数量级报警阈值怎么调都失灵最后排查下来就是泄漏问题。2.2 训练阶段必须锁定的四项参数为了避免这类事故我的做法是建立一个明确的“模型对象”把训练阶段得到的全部统计量固化进去。以下四项是必须在训练集上计算并保存的项目含义测试阶段如何处理训练核矩阵列均值每个训练样本在特征空间中与其他训练样本内积的平均作为常量不再重新估计训练核矩阵全局均值所有训练样本内积的总体平均作为常量用于测试样本中心化归一化后的特征向量 α高维空间中的主方向系数直接用矩阵乘法投影特征值 λ 与核参数 σ特征空间中的方差信息与核宽度保持不变其中核参数 σ 在测试阶段尤其容易被疏忽。有些人训练时用网格搜索选出了最优σ测试时却又重新算了一遍核矩阵甚至用了不同的σ。这相当于训练和测试使用了两把不同的尺子去量同一个物体后端的分类器必然失效。2.3 一个反面例子的代价下面这段代码是我故意写的错误版本用来演示“测试集独立中心化”的典型错误% 错误写法测试集单独中心化、单独特征分解 K_test kernel_gaussian(X_test, X_test, sigma); % 注意测试与测试 Kc_test center_kernel_matrix(K_test); % 错误self-centering [alpha_test, lambda_test] eig(Kc_test / size(K_test,1));这样做的直接后果有两个。第一测试集的投影方向源于测试集自身的方差结构和训练集的主方向毫无对应关系第二由于两次中心化使用了不同的均值实验里测出来“降维后分类准确率反而更高”时你该警惕——那几乎可以断定是数据泄漏制造的假象。机器学习里有一条朴素原则测试阶段唯一允许做的操作是用训练阶段确定好的变换去“套用”到新样本上。KPCA的测试投影正是这个原则的典型案例。3. 训练阶段落地核矩阵、中心化、特征分解一套带走3.1 高斯核函数的矢量计算先把基础工具写好。计算两个样本集合之间的高斯核矩阵不要用双重循环直接用矩阵运算展开距离公式function K kernel_gaussian(X, Y, sigma) % X: nx-by-m, Y: ny-by-m, 返回 K: nx-by-ny nx size(X, 1); ny size(Y, 1); G X * Y; % 线性内积 X2 sum(X.^2, 2); % nx-by-1 Y2 sum(Y.^2, 2); % ny-by-1 dist2 repmat(X2, 1, ny) repmat(Y2, nx, 1) - 2 * G; dist2(dist2 0) 0; % 防止数值误差导致负距离 K exp(-dist2 / (2 * sigma^2)); end这里有个经常被忽略的数值细节当样本维度较高时X2 Y2 - 2G 在浮点运算下可能出现极小的负值虽然理论上距离平方不可能为负。负值直接喂给 exp 会导致核矩阵出现 NaN。我习惯加一行截断省得排查半天找不到原因。3.2 训练函数的完整实现接下来是核心的训练函数。我会把所有训练阶段需要保存的内容打包进一个结构体 model后续测试投影只需要这一个变量。function model kpca_train(X_train, sigma, n_components) % KPCA训练只使用训练集 % 输入 % X_train: n-by-m 训练样本矩阵 % sigma: 高斯核宽度 % n_components: 保留的主成分个数 % 输出 % model: 结构体包含测试阶段所需的全部统计量 [n, ~] size(X_train); % 1. 计算训练核矩阵并中心化 K kernel_gaussian(X_train, X_train, sigma); Kc center_kernel_matrix(K); Kc (Kc Kc) / 2; % 强制对称 % 2. 特征分解直接分解 Kc/n 或 Kc 均可注意特征值含义的差异 % 这里分解 Kc/n便于与PCA的“方差解释率”对齐 [V, D] eig(Kc / n); lambda diag(D); % 3. 按特征值降序排列 [lambda_sort, idx] sort(lambda, descend); V_sort V(:, idx); % 4. 归一化让特征向量满足 alpha * Kc * alpha 1 % 这样训练得分的方差恰好等于对应特征值语义和PCA一致 alpha zeros(n, n); for i 1:n norm_factor sqrt(V_sort(:, i) * Kc * V_sort(:, i)); if norm_factor eps alpha(:, i) V_sort(:, i) / norm_factor; else alpha(:, i) V_sort(:, i); % 特征值近乎零时跳过硬归一化 end end % 5. 截断到指定主成分数 model.alpha alpha(:, 1:n_components); model.lambda lambda_sort(1:n_components); model.sigma sigma; model.X_train X_train; model.n_components n_components; % 6. 保存训练核矩阵的均值统计量测试阶段中心化必需 model.meanK_all mean(K(:)); model.meanK_col mean(K, 1); % 1-by-n 列均值 end这段代码里有几个设计决策值得解释。特征分解对象我选择了 Kc/n 而不是 Kc这样得到的特征值 λ 直接就是对应主成分解释的方差占比类比量方便后续计算累计贡献率和确定主成分个数。归一化条件用了 alpha * Kc * alpha 1而不是 alpha * alpha 1这背后的原因是经过归一化后训练得分 T_train Kc * alpha 的协方差矩阵正好是对角阵 diag(lambda)每一个主成分的方差可解释性完全对齐PCA。3.3 主成分个数的选择n_components 的选择有硬性上限就是训练样本数 n。因为核矩阵是 n×n 的特征向量最多也只有 n 个。实际使用中我习惯先保留全部特征值再按累计贡献率自动截断explained cumsum(model.lambda) / sum(model.lambda); n_auto find(explained 0.85, 1, first); % 85%累计贡献率对于故障检测场景85%的阈值经常不够因为T²和SPE统计量对少量小特征值也很敏感我会把阈值放到95%甚至更高。分类场景则可以适当放松。这类决策没有绝对标准但至少比拍脑袋定个数要靠谱。4. 测试阶段落地用训练参数投影新样本的完整代码4.1 测试投影的公式拆解测试样本 x_new 在第 k 个主成分上的得分理论公式是t_k(x_new) Σ_i α_i^k · k(x_new, x_i)这里的 α_i^k 是训练阶段得到的归一化特征向量第 i 个分量x_i 是训练样本。所以在Matlab里投影的骨架就是“测试核向量”与“特征向量矩阵”的乘积。难点在于这个“测试核向量”必须先做中心化而中心化必须使用训练集的均值统计量。设训练集核矩阵 K 的列均值为 meanK_col一个 1×n 行向量全局均值为 meanK_all测试样本与所有训练样本的核向量记作 kt中心化后的向量为kt_c(j) kt(j) - meanK_col(j) - mean(kt) meanK_all其中 mean(kt) 是该测试样本核向量的平均值。这个公式是训练集中心化公式在单样本视角下的精确对应减去训练集列均值再减去测试核向量自身的行均值最后加回全局均值。4.2 测试投影函数实现function T_test kpca_project(model, X_test) % KPCA测试投影使用训练阶段保存的全部统计量 % 输入 % model: kpca_train 输出的结构体 % X_test: nt-by-m 测试样本矩阵 % 输出 % T_test: nt-by-p 测试得分矩阵 % 1. 计算测试样本与训练样本的核矩阵 (nt-by-n) Kt kernel_gaussian(X_test, model.X_train, model.sigma); % 2. 用训练集均值统计量做中心化 nt size(X_test, 1); n size(model.X_train, 1); Kt_c Kt - repmat(model.meanK_col, nt, 1) ... - repmat(mean(Kt, 2), 1, n) ... repmat(model.meanK_all, nt, n); % 3. 投影得分 中心化核矩阵 * 归一化特征向量 T_test Kt_c * model.alpha; end注意这里三个 repmat 分别实现了公式里的三项列均值修正、测试样本核向量自身均值修正、全局均值回加。如果哪一项漏了测试得分的数值范围会明显偏离训练得分在T²控制图上表现为所有新样本集体偏移或集体收缩。4.3 一句话理解训练得分和测试得分的关系训练得分在函数内部可以这样算T_train model.Kc * model.alpha; % 虽然model里没存Kc但效果等价如果我想在模型里缓存训练得分建议在kpca_train里额外保存 T_train这样测试阶段可以直接对比训练得分与测试得分的分布是否一致这是一个非常有效的诊断手段。具体来说如果某个测试样本的得分向量与训练得分的总体中心出现系统性偏移往往意味着数据分布发生了漂移或者核参数σ选择不合理。5. 我踩过的坑核宽度、归一化约定与批量处理建议5.1 高斯核宽度σ的选择是个绕不开的坎σ直接决定核矩阵的形状。σ过小核矩阵趋向单位阵每个样本都只和自己“相似”特征空间里样本变成孤岛降维等于没降σ过大核矩阵趋向全1矩阵中心化后数值几乎为零有用信息全部被抹平。我的经验是从“数据两两距离的中位数”入手。计算步骤很简单% 抽样计算两两距离的中位数避免n大时内存爆掉 sample_idx randperm(n, min(n, 500)); Dsq pdist2(X_train(sample_idx, :), X_train(sample_idx, :)).^2; sigma_init sqrt(median(Dsq(:)) / 2);以这个为初始值再用网格搜索或者交叉验证微调。对于分类任务可以配合准确率搜索对于故障检测可以看T²和SPE统计量的漏报率。不要指望一个σ走天下换数据集就得重调。5.2 特征分解时eig与eigs的选择当训练样本数不大比如 n 3000直接用 eig(full(Kc)) 最省事一次性拿到全部特征向量。如果只想要前 p 个主成分可以用 eigs(Kc, p, largestabs)计算量会小很多但我遇到过小特征值接近零时 eigs 收敛不稳定的情况。一个额外提醒如果样本量大到核矩阵存不下n 5000 时 double 矩阵就是 200MB可以考虑Nyström近似或者随机特征近似。这部分属于进阶话题简单的思路是用随机采样出的锚点计算低秩核近似再在近似的特征空间里做投影。工程上够用就行不必强行上大杀器。5.3 归一化约定不同是代码移植时最大的坑不同论文、不同工具箱对特征向量的归一化约定不完全一致有的要求 αα 1有的要求 αKcα 1。我用的版本是后者理由在前面说过能让投影得分方差和特征值对齐。迁移代码时一定要确认对方用的是哪种约定否则训练得分看着正常测出来的T²统计量却可能整体差出一个常数倍。另外有些代码会把 Kc 直接除以 n 再做特征分解有的不除。这不影响特征向量的方向但特征值相差 n 倍。如果你要拿特征值计算累计贡献率必须和核矩阵分解方式保持一致不然解释方差比例会彻底乱套。5.4 批量测试与在线单样本预测的差异在线诊断场景里测试样本往往是一个一个来的而不是一次给一整批。这时候 kpca_project 依然适用因为 Kt 的大小是 1×n完全没有问题。但有一点要注意mean(Kt, 2) 在单样本时就是该样本核向量的平均值批量时是每个样本各自的行均值公式统一不需要特殊分支处理。5.5 如果要做重构的补充说明KPCA 的重构即从低维得分回到原始空间不像PCA那样直接乘载荷矩阵转置因为高斯核隐含的特征空间是无穷维的。常用的做法是求解一个 pre-image 问题用迭代优化找到原始空间里某个样本使它的核映射最接近目标得分。这个过程的细节不少如果只是做特征提取、分类或监控异常投影得分本身已经足够不必额外做重构。在我自己维护的代码里kpca_train 和 kpca_project 这两个函数加起来不过几十行但每次调参和排查问题十有八九都集中在中心化统计量的传递和归一化约定上。把这几个点理清楚之后KPCA在Matlab里落地就是一件非常规整的事。后续如果想扩展到多工况自适应或者增量更新也只需要在这套Train/Test分离的骨架上去改训练阶段的统计量更新策略测试阶段的投影逻辑完全不需要动。
