先别急着把半监督判别分析SDA归类成“又一个矩阵运算技巧”。它解决的场景特别具体——标签很贵但没标签的数据不值钱。我最早是在一个分类项目里接触到SDA每个类别只有二十来个标注样本硬盘里却躺着上万条无标签的同类特征。用传统LDA做降维投影方向只要换一次训练子集就飘分类精度忽高忽低。后来改用半监督判别分析SDA把未标注数据的流形结构当作约束项加进目标函数投影方向一下子稳住了MATLAB里实现起来也不过几十行。这篇文章就围绕SDA的目标函数、几何意义和完整实现展开适合已经熟悉LDA、想往半监督方向迈一步的读者。1. LDA在少标注样本下的“失速点”与SDA的解决思路1.1 类内散度不再可靠判别方向开始漂传统线性判别分析LDA最大化类间距离同时最小化类内距离本质是求一组投影方向让投影后的数据“类间尽量分开、类内尽量聚拢”。它的关键输入是类内散度矩阵 S_w 和类间散度矩阵 S_b。问题就出在这类内散度矩阵需要每个类别内部的样本共同估计协方差结构标注样本一旦少到接近特征维度这个估计就非常不可靠。我举一个常见例子。假设数据是64维的图像特征每类只有10个标注样本。此时单个类别的协方差矩阵本来就是病态的几类样本叠在一起估计出来的 S_w 几乎沦为一堆噪声的加权和。LDA再厉害输入是坏的输出自然也是坏的。投影方向会被少数几个偏离中心的标注样本牵着走换一批标注样本反复实验时投影方向的变化幅度肉眼可见。更麻烦的是在很多真实场景里还有类别不平衡。标注样本少的类别在 S_b 中的贡献权重被压低LDA会更多照顾大类别小类别经常被“牺牲”掉。这种问题上单纯增加伪标签或者重采样都容易引入额外偏差治标不治本。1.2 SDA的核心思路把未标注数据当成“结构约束器”SDASemi-supervised Discriminant Analysis最初的动机非常朴素标注数据不够但未标注数据里藏着大量关于数据分布的信息。如果我们相信一个经典假设——相似的数据点在投影之后也应该保持相似那么未标注数据就能充当一把尺子去惩罚那些破坏局部结构的投影方向。这把尺子通常通过图拉普拉斯来实现。先把所有样本标注的加未标注的看作图的节点用k近邻关系或热核距离连接“相似”节点得到邻接矩阵 W。再计算拉普拉斯矩阵 L D - W。SDA在传统LDA的类间散度上加入一项 βXLXᵀ然后求解广义特征值问题。这里 β 控制未标注数据对判别方向的约束强度β越大投影方向就越需要尊重数据原有的局部流形结构。这个思路之所以有效是因为它在判别目标和几何先验之间搭了一个平衡。标注样本负责告诉算法“哪些类需要分开”未标注样本负责告诉算法“哪些方向走起来不自然”。两者合作时投影方向很难再被少量标注样本的噪声带偏。1.3 从LDA到SDA的关系不是替代而是扩展很多人第一次看到SDA公式会问这跟“先用LDA再加正则化”有什么区别从数学上说确实等价于加了一个流形正则化项但实际意义完全不同。LDA回答的是“类间与类内的比值最优”SDA回答的是“判别边界最优且投影后数据结构不发生扭曲”。下面这张表可以帮你快速定位SDA在算法家族里的位置方法是否用标注是否用未标注目标方向PCA否否保留全局方差LPP局部保持投影否是保留局部流形结构LDA是否最大化类间可分性SDA是是类间可分 流形结构保持从表中能看出SDA占据了一个很特殊的位置它同时吸收LDA和LPP的优点。若 β0它退化为LDA若去掉 S_b 项只剩正则项又等价于拉普拉斯特征映射的线性版本。理解这层关系后面调参和扩展就顺手多了。2. 目标函数与拉普拉斯正则项的几何含义2.1 从LDA的目标函数原地扩展LDA的目标函数有很多等价写法最方便扩展到半监督的是这个形式max tr((UᵀS_tU)⁻¹ UᵀS_bU)其中 S_t 是总体散度矩阵S_b 是类间散度矩阵。为什么用它因为总体散度 S_t 可以由全部样本计算哪怕没有标注也能算而类间散度 S_b 只需要每个类的均值向量和类内样本数。这个形式天然适合半监督改造。SDA在分子里加一项 βXLXᵀ变成max tr((UᵀS_tU)⁻¹ Uᵀ(S_b βXLXᵀ)U)这里 XLXᵀ 是 d×d 矩阵d是原始特征维度。整体含义是投影方向不仅要让各类均值充分拉开还要避免在未标注数据的局部近邻结构上“横冲直撞”。对应于广义特征值问题就是(S_b βXLXᵀ) u λS_t u需要注意这里用的 S_t 一般是基于全部样本计算的总体散度这一点与原始LDA的定义略有差异但对实际使用影响不大反而让未标注数据的分布信息也能进入分母。2.2 图拉普拉斯惩罚做了什么拉普拉斯项 XLXᵀ 的几何含义值得多说几句。对于任意投影向量 u样本 x_i 和 x_j 投影后得到标量 uᵀx_i 和 uᵀx_j。可以验证uᵀX L Xᵀu Σᵢⱼ Wᵢⱼ (uᵀxᵢ - uᵀxⱼ)²这个式子非常直观如果两个原始样本在邻接图中被 Wᵢⱼ 连起来了那么投影后它们的距离应该尽量小如果这两个样本根本不相似Wᵢⱼ0惩罚项就直接忽略它们。所以正则项本质是一个“投影光滑性”惩罚——在原始空间相邻的点投影后不许被强行拆散。用生活化的类比LDA像是你在地图上把几个城市中心尽量隔开而SDA的拉普拉斯项像修路规则——你规划的高速路不能把原本紧挨着的村庄一分为二。判别方向照样要拉但必须顺着人群原本的聚集方式走。这个惩罚项特别适合处理“标注样本稀疏但数据本身有清晰流形”的情况。例如手写数字、人脸图像、语音特征往往都分布在低维流形上未标注样本可以帮算法把流形方向大致框出来。2.3 广义特征值问题与求解条件SDA最后归结到广义特征值问题。理论上需要保证 S_t 可逆但实际数据里 S_t 经常奇异尤其是特征维度大于样本数的时候。这种情况下直接用 eig 函数很可能得到 Inf 或 NaN解决方式放在后面实现章节细说。另一个实际问题是特征向量的排序。求解广义特征值问题后我们要取特征值最大的 r 个方向。这里的特征值衡量的是“该方向上的判别性与流形保持性之和”与“总体散布”的比值取前r个就是最大化这一比值的方向组合。对多分类问题传统LDA最多只能取 c-1 个有效方向c是类别数。但SDA因为加入了 XLXᵀ 这项有效方向数有潜力超过 c-1这一点在实践里经常被忽略我后面单独展开。3. MATLAB实现步骤与函数封装3.1 数据准备与预处理先约定输入约定X是 d×n 的矩阵每一列是一个样本y是 n×1 的标签向量标注样本的标签是正整数未标注样本统一设为0。这种约定在实际工程里最顺因为MATLAB的 find、逻辑索引都能直接处理。预处理方面我的建议是先做 zscore 标准化。原因很简单邻接图构建用的是欧氏距离如果某个特征天然数值很大比如像素值0到255而另一个特征天然很小比如归一化后的角度距离会被大数值特征主导k近邻图实际只反映了单个特征的邻居关系。标准化之后每个特征在距离计算中才有“投票权”。X zscore(X, 0, 2); % 按行标准化每列是一个样本调用时注意zscore的维度参数MATLAB里按行运算时第二参是0、第三参是2指沿列方向标准化。3.2 邻接图与拉普拉斯矩阵构建邻接图构建是SDA实现里最容易出问题的地方。常见方案有k近邻和ε-球两种SDA文献里用得最多的是k近邻加热核权重。代码可以这样写function L compute_laplacian(X, k, t) n size(X, 2); D2 pdist2(X, X).^2; [~, idx] sort(D2, 2); mask false(n); for i 1:n mask(i, idx(i, 2:k1)) true; % 排除自身 end W exp(-D2 / (2 * t^2)); W(~mask) 0; W (W W) / 2; % 强制对称允许“互不为近邻”退化 D diag(sum(W, 2)); L D - W; L (L L) / 2; % 数值对称保险 end这里有个细节排序后 idx(i, 2:k1) 从第二列开始取是因为第一列是自身距离一定是0。但少数情况数据里存在完全重复的样本自身距离并不排第一这时取到的邻居会包含自身导致W对角线上出现非零值。稳妥做法是在排序前把对角线距离设成Inf或者在取邻居时显式排除 iD2 pdist2(X, X); D2(eye(n) 1) Inf; [~, idx] sort(D2, 2); mask false(n); for i 1:n mask(i, idx(i, 1:k)) true; % 现在第1列才是真正的最近邻居 end热核宽度 t 的确定我推荐一个省事的经验法取所有样本到其第k近邻距离的中位数。这样自动适配数据尺度不会因为整体数据放缩而需要手调kth_dist zeros(n,1); for i 1:n kth_dist(i) D2(i, idx(i, k)); end t median(kth_dist);3.3 散度矩阵计算与数值稳定处理S_b 和 S_t 的实现比较直接但要注意均值总体的选择。我这里基于全部样本计算总体均值标注样本计算各类均值function [Sb, St] compute_scatter(X, y, mu_all) [d, n] size(X); labels unique(y(y 0)); c length(labels); Sb zeros(d, d); if c 0 for ci 1:c idx_c (y labels(ci)); n_c sum(idx_c); mu_c mean(X(:, idx_c), 2); diff mu_c - mu_all; Sb Sb n_c * (diff * diff); end end Xc X - mu_all; St Xc * Xc; endS_t 奇异的问题必须在这里就处理掉。直接做法是加一个微小平移项也就是在原始理解上给它一个小正则本质是把不可逆问题近似成可逆问题。实际效果很好traceSt trace(St); St_reg St 1e-6 * (traceSt / d) * eye(d);如果你希望更严格可以先用PCA把数据压到不超过 n-1 维再做SDA。这个策略在高维小样本场景下是最稳妥的既降噪又保证非奇异。3.4 主函数与投影应用核心求解部分不长function [Proj, eigVal] sda(X, y, beta, k, t, r) [d, n] size(X); mu_all mean(X, 2); [Sb, St] compute_scatter(X, y, mu_all); St_reg St 1e-6 * (trace(St) / d) * eye(d); L compute_laplacian(X, k, t); M Sb beta * (X * L * X); M (M M) / 2; [V, D] eig(M, St_reg); eigVal real(diag(D)); [~, ord] sort(eigVal, descend); V real(V(:, ord)); Proj V(:, 1:r); for j 1:r v Proj(:, j); nv sqrt(v * St_reg * v); if nv 1e-12 Proj(:, j) v / nv; end end end这里的投影方向我做了 St_reg 正交归一化是为了让不同投影向量的尺度一致后续做KNN分类时避免某些方向天然就很大、某些方向天然就很小。如果你只关心降维而不要求向量本身可解释归一化可以省。投影的应用就是普通矩阵乘法Ztrain Proj * Xtrain; Ztest Proj * Xtest;用KNN分类器即可。自写一个简单的欧氏距离KNN就够function pred knn_classify(Ztrain, ytrain, Ztest, kk) D pdist2(Ztest, Ztrain); [~, idx] sort(D, 2); pred zeros(size(Ztest, 2), 1); for i 1:size(Ztest, 2) nbr_labels ytrain(idx(i, 1:kk)); pred(i) mode(nbr_labels); end end整个主流程从数据到分类结果其实不到60行MATLAB代码。这就是SDA在工程上的友好之处——原理听起来很绕代码却非常紧凑。4. 在合成数据上验证SDA如何救活“标注样本稀缺”的场景4.1 实验设置为了能看到SDA在少标注样本下的优势我建议用三团高斯簇合成数据做实验。每簇200个样本、10维特征簇均值人为拉开协方差矩阵引入一定程度的各向异性使得简单的欧氏距离分类不能很好工作。rng(42); d 10; nPerClass 200; mu1 zeros(d, 1); mu2 3 * ones(d, 1); mu3 [-2; 2; zeros(d-2, 1)]; X [mvnrnd(mu1, eye(d), nPerClass); mvnrnd(mu2, eye(d), nPerClass); mvnrnd(mu3, eye(d), nPerClass)]; y_true [ones(nPerClass, 1); 2 * ones(nPerClass, 1); 3 * ones(nPerClass, 1)];然后让每类只有2个样本带标签其余作为未标注样本参与SDA的邻接图构建。把带标签的这部分单独留出来做LDA基线未标注数据不参与LDA计算。这样能公平对比“只看少量标注”和“少量标注大量未标注”的差距。测试集单独采样与训练/未标注数据完全分开避免信息泄漏。4.2 解读投影与分类精度我跑了多组随机初始化观察到趋势非常一致。只靠每类2个标注样本的LDA投影方向严重偏向两个被标出的样本连线的方向第三类几乎没被照顾到测试集KNN精度在0.6左右徘徊。SDA因为加入了未标注样本构建的拉普拉斯约束投影方向会沿三个簇围成的整体分布走分类精度能稳定提到0.8以上有时接近0.9。典型的一组结果是LDA大约0.62SDA约为0.85。具体数字会因生成的数据不同而波动不用太纠结绝对值重点在于趋势标注样本越少SDA相对LDA的增益越明显。当每类标注样本从2个增加到40个LDA逐渐赶上SDA两者差距缩小到几个百分点。这说明SDA真正发挥作用的前提就是“标注稀缺”。投影可视化也有趣。LDA降维后的二维散点图上三类数据经常挤成两团因为LDA最多只能提供 c-12 个方向且那一点点标注样本限定了方向。SDA的二维投影里第三类往往能露出来。这不是魔法而是未标注数据通过拉普拉斯项提供了标注数据没有的全局结构信息。4.3 关键参数β、k、t的调节经验SDA有三个直接影响结果的参数β、k、t。其中β是注意力最集中的地方。我的建议是网格搜索时在 log 空间取点betas logspace(-2, 2, 5); % 0.01, 0.1, 1, 10, 100原因很简单目标函数是“判别项 β乘正则项”的线性叠加参数每增大10倍对解的影响其实很大线性网格根本没法覆盖合理范围。对合成数据β1附近效果最好对真实数据β0.1到10之间往往能找到甜点。k近邻个数一般取5到10就够了。样本量特别大时可以取 sqrt(n) 附近的整数但别超过50。k太小近邻图容易断裂k太大又把全局结构强行当成局部结构流形假设的优势会被稀释。t 的决定用我前面说的“第k近邻距离中位数”几乎从不出问题。如果邻接矩阵给得太密很多权重都趋近于1可以适当调小t反之如果权重几乎全趋近于0就把t调大。还有一点需要提醒β并不是越大越好。β过大的时候SDA会无限逼近无监督的LPP类别信息被“淹掉”投影方向主要反映数据分布而不是类别可分性。β过小又退化成不稳定LDA。合理的判断标准不是“β取多少精度最高”这一个孤立指标而是看是否稳定——在训练子集变化时投影方向是否仍然稳定。5. 实际使用中的几个坑与扩展思考5.1 矩阵对称性与奇异问题的处理代码层面最容易犯的错是忘记对称化。邻接矩阵W如果不强制对称计算出来的拉普拉斯矩阵可能带有符号特征值eig求解时出现复数导致排序结果完全乱掉。我在封装上面的 compute_laplacian 时故意做了两次对称化一次针对 W一次针对 L充分保证数值上符合理论前提。S_t 奇异是另一个高频坑。特征维度大于样本数时S_t必然奇异直接调用 eig 解广义特征问题通常返回一堆Inf。我建议的处理顺序是先做PCA降维到不超过样本数-1的维度再做SDA。这比单纯加正则项更稳因为PCA能把零空间直接消掉不让无效方向参与后续计算。如果不用PCA那就必须在 S_t 上加上 1e-6 量级的单位阵让矩阵严格可逆。还有一个容易忽略的坑eig(A, B) 对矩阵A、B的对称性非常敏感哪怕其中一个有1e-10的数值不对称特征值就可能出现虚部。求解前用 A(AA)/2 这类操作“洗一次”很有必要。5.2 有效方向数与多分类问题传统LDA能提取的判别方向最多是 c-1 个。SDA在分子上加了 βXLXᵀ因此广义特征问题中的 M 矩阵秩可以显著大于 c-1。换句话说SDA允许你提取超过类别数减一的投影方向这在很多情况下是实用优势。比如一个三分类问题传统LDA只有2个有效方向第三类信息很容易被前两类的均值差异淹没。SDA加上未标注数据的多样化结构后第三个方向也能承载判别信息分类器在更高维度的投影空间里能切出更复杂的边界。这也是为什么我在主函数里保留了参数 r而很少有人讨论SDA的这个特性。5.3 从SDA再往前走一步SDA在2011年前后陆续衍生出一堆变体比如把柔性标签纳入迭代过程的方法用核技巧扩展的非线性版本还有加上稀疏约束的稀疏SDA。如果你在真实项目中体会到“标注太少”的痛苦SDA可以作为第一根拐杖先把方向和流程跑通再考虑要不要上更复杂的深度半监督模型。我个人现在的习惯是凡是遇到标注样本稀缺、未标注数据又明显存在聚类或流形结构的任务先把SDA当作基线降维方法配合KNN跑一遍完整参数区间。这一套下来通常能快速判断“半监督这条路值不值得走”。在MATLAB里从无到有实现全部流程包括邻接图构建、散度矩阵计算、广义特征值求解和KNN评估普通笔记本上处理几千个样本也就是几秒钟的事。这种低成本试错才是SDA这类传统方法在今天依然值得掌握的理由。