2.6. 线性判别分析线性判别分析LDA一次fit做两件事。第一它是一个生成式分类器把每个类别建模成一个高斯分布各类均值不同、协方差矩阵在所有类别间共享再用贝叶斯规则挑出最可能的标签。第二它是一个有监督的降维器把特征投影到最能分开各类别的方向上。RustyML 的LDA一次fit调用就同时给出这两种结果predict、decision_function、predict_proba服务于分类器一侧transform和fit_transform服务于投影一侧。本页覆盖这套 API、它背后的自研线性代数、这个模型依赖的高斯假设以及协方差矩阵退化为奇异时会发生什么。2.6.1. crate 实现了什么RustyML 把 LDA 实现在rustyml::machine_learning::{LDA, Shrinkage, DiscriminantSolver}下prelude 也重新导出了这些类型参见 prelude。模型吃一个f64特征矩阵和i32标签。LDA 是分类器所以标签是整数类别 id不是浮点数。下表列出了完整的 API 接口方法签名简化返回值LDA::newnew(n_components: usize)ResultLDA, ErrorLDA::defaultdefault()LDA成分数自动确定with_solverwith_solver(DiscriminantSolver)LDAwith_shrinkagewith_shrinkage(Shrinkage)ResultLDA, Errorfitfit(x, y)其中y: Array1i32Resultmut LDA, Errorpredictpredict(x)ResultArray1i32, Errordecision_functiondecision_function(x)ResultArray2f64, Errorpredict_probapredict_proba(x)ResultArray2f64, Errortransformtransform(x)ResultArray2f64, Errorfit_transformfit_transform(x, y)ResultArray2f64, Errorcrate 没有实现 QDA二次判别分析只提供共享协方差、线性边界这一种变体。getterget_classes、get_priors、get_means、get_overall_mean、get_projection在fit运行前返回None。fit运行后每个 getter 都会把值包在Some里返回。查看任意一个 getter就能知道模型是否已经训练过。求解器枚举过去用的是不加修饰的名字Solver这个名字和线性模型自己的Solver枚举撞了。经过 prelude 的 glob 导入之后Solver::GradientDescent会解析到错误的那个枚举上导致编译失败。现在两个枚举都按各自选择的对象命名这里是DiscriminantSolver线性模型那边是LeastSquaresSolver。2.6.2. 第一个分类器先从 2D 里 3 个彼此分得很开的簇入手。LDA::default()会自动定下成分数只想要标签的时候不需要自己设置成分数。usendarray::array;userustyml::machine_learning::LDA;fnmain(){// 3 个紧凑、彼此分得很开的簇每簇 3 个样本letxarray![[1.0,1.0],[1.5,0.8],[0.8,1.2],[5.0,5.0],[5.2,4.8],[4.8,5.2],[9.0,1.0],[9.2,0.8],[8.8,1.2],];letyarray![0,0,0,1,1,1,2,2,2];letmutldaLDA::default();lda.fit(x,y).unwrap();letpredslda.predict(x).unwrap();println!(predictions: {:?},preds);// 每个类别的判别得分n_samples, n_classes与 softmax 后验概率letscoreslda.decision_function(x).unwrap();letprobalda.predict_proba(x).unwrap();println!(scores shape: {:?},scores.shape());println!(proba shape: {:?},proba.shape());println!(class order: {:?},lda.get_classes().unwrap());}predict、decision_function、predict_proba共用同一个核心计算。对一个样本xLDA 为每个类别算一个线性判别得分公式是score_j(x) x * Sigma^-1 * mu_j - 0.5 * mu_j * Sigma^-1 * mu_j ln(prior_j)。这里Sigma是共享协方差mu_j是类别均值prior_j是该类别在训练集里的频率。decision_function把这些原始得分以(n_samples, n_classes)矩阵的形式返回。predict取每一行得分的 argmax再映射回原始标签。predict_proba对同一批得分做数值稳定的逐行 softmax因此每一行的和为 1。decision_function和predict_proba的列顺序跟着get_classes()走而get_classes()按升序排列标签。不要假设这个顺序跟标签在y里首次出现的顺序一致。这里的先验就是纯粹的类别频率n_class / n_samples。对这个类别均衡的数据集来说每个先验都是1/3。类别不均衡会改变这一点ln(prior_j)这一项会把边界推向更稀有的类别一个贝叶斯分类器本该如此。想要先验相等就把训练集配平——LDA的构造器没有priors这样的覆盖项。2.6.3. 有监督投影与 n_classes - 1 上限投影这一侧才是 LDA 区别于 PCA 的地方。LDA 找的是让类间散度相对类内散度最大化的方向。只有n_classes - 1个方向具备非零区分度因为类间散度矩阵的秩至多为n_classes - 1。crate 强制执行这条上限可用成分数是min(n_classes - 1, n_features)在 fit 时定下来。LDA::default()或n_components None取这个上限LDA::new(k)精确请求k个成分若k超过上限则在fit时失败。所以画一张 2D 散点图至少需要 3 个类别。两个类别只给得出一根判别轴。fit_transform一次调用就完成训练加投影usendarray::array;userustyml::machine_learning::{DiscriminantSolver,LDA,Shrinkage};fnmain(){letxarray![[1.0,1.0],[1.5,0.8],[0.8,1.2],[5.0,5.0],[5.2,4.8],[4.8,5.2],[9.0,1.0],[9.2,0.8],[8.8,1.2],];letyarray![0,0,0,1,1,1,2,2,2];// 3 个类别 - 至多 2 根判别轴直接投到 2DletmutldaLDA::new(2).unwrap().with_solver(DiscriminantSolver::Eigen).with_shrinkage(Shrinkage::Auto).unwrap();letcoordslda.fit_transform(x,y).unwrap();// (9, 2)println!(embedding shape: {:?},coords.shape());for(row,label)incoords.outer_iter().zip(y.iter()){println!(class {label}: ({:.3}, {:.3}),row[0],row[1]);}// 投影矩阵把 n_features - n_componentsletwlda.get_projection().unwrap();println!(projection shape: {:?},w.shape());// [2, 2]}embedding shape: [9, 2] ... one line per sample: class label: (axis 1, axis 2) ... projection shape: [2, 2]投影矩阵W的形状是(n_features, n_components)。transform算的是(x - xbar) * W。这里xbar是整个训练矩阵的均值对应 scikit-learn 的xbar_可以通过get_overall_mean()读出来。RustyML 严格对齐 scikit-learn 的(X - xbar_) scalings_公式。这次对齐带来两处细节而且都是最近才改的。第一投影会按训练均值做中心化因此变换后的训练集落在原点附近而不是落在原始特征值本身的位置。第二各条轴保留白化产生的尺度RustyML 不会把它们重新归一化到单位 L2 范数。这个尺度不是装饰性的它让投影后的数据具有单位类内协方差这正是整套构造要达成的目标。如果你以前是从一个单位范数的旧轴上读坐标的现在会看到不同的偏移和不同的尺度。RustyML 按类别区分度给各列排序广义特征值最大的排在最前所以第一个成分就是那根最具区分力的轴。如果你之后想从同一次拟合里做 1D 降维这个排序会有用。一根判别轴和它的反向对分类的效果完全一样所以任何一列的符号都是任意的。把画出来的坐标轴方向当成没有意义不要从中读出方向性。transform和fit_transform执行的是同一个降维步骤。fit_transform(x, y)用一次调用依次跑fit(x, y)和transform(x)。2.6.4. 求解器与自研线性代数DiscriminantSolver决定 RustyML 怎么为分类器的评分系数Sigma^-1 * mu_c求共享协方差的逆。这是一个重要、也极易被忽略的事实求解器的选择只影响predict、decision_function、predict_proba。transform投影不依赖求解器RustyML 始终从类内协方差的对称特征分解里得到它。换求解器会改变决策边界的数值但永远不会改变 2D 嵌入。DiscriminantSolver怎么评分什么时候用它SVD默认对Sigma做 SVD 伪逆截断小奇异值通用默认能很好地处理秩亏Eigen对称特征分解求逆极小特征值归零协方差良态时使用效果和SVD相当LSQR迭代式 Paige-Saunders 最小二乘不形成显式逆维度极高显式求逆会浪费内存3 种求解器在分得开的数据上得到的标签完全一致差别只在于它们怎么应对病态的Sigma。SVD和Eigen都会构造一个伪逆丢掉那些奇异值或特征值低于相对容差的方向因此秩亏的协方差会优雅地退化而不是直接失败。LSQR完全绕开求逆直接迭代求解每个Sigma * coef_c mu_c方程组这在n_features很大时能省内存。线性代数完全是自研的纯 Rust 实现RustyML 不带 LAPACK 或nalgebra依赖见优先自研数值实现。对称特征分解用的是 Householder 三对角化接上隐式位移 QL 迭代也就是经典的 EISPACKtred2和tql2这对组合。SVD 用的是单边 Jacobi 方法能把小奇异值算到很高的相对精度。LSQR用的是带滚动 Givens 旋转的 Golub-Kahan 双对角化。在 LDA 产生的中小规模协方差矩阵上这些方法都能确定到机器精度这也是 LDA 不需要随机种子的原因。2.6.5. 收缩与协方差正则化每个特征只有很少几个样本时原始样本协方差是个带噪的估计它的逆会把噪声放大。收缩shrinkage把Sigma往一个缩放后的单位矩阵拉用一点点偏差换来方差的大幅下降。Shrinkage有两种形式Shrinkage::AutoLedoit-Wolf 闭式最优强度RustyML 直接从数据算出来没有调参旋钮。Shrinkage::Manual(alpha)显式给一个落在[0, 1]里的alpha。取0就是不收缩取1就完全收缩到单位矩阵目标。with_shrinkage会校验alpha如果它落在[0, 1]之外或者不是有限值就返回Error::InvalidParameter——这正是with_shrinkage返回Result的原因。这里有两个细节值得注意。第一不管你怎么选收缩RustyML 总会给协方差加一条极小的固定对角脊大约是平均方差的1e-6倍纯粹是为了让分解保持稳定。所以Shrinkage::Manual(0.0)和完全不收缩会得到一模一样的结果。第二这里的收缩和求解器是相互独立的这是和 scikit-learn 一处实打实的分歧scikit-learn 在它的svd求解器下拒绝收缩RustyML 则在任何求解器跑之前就把收缩施加到协方差上所以DiscriminantSolver::SVD配Shrinkage::Auto完全合法而且往往是最好的组合。只要遇到共线特征或者样本-特征比很吃紧的情况就把Shrinkage::Auto当默认选择。2.6.6. 假设以及假设失效时的症状LDA 的最优性建立在两条假设上每个类别都服从高斯分布且所有类别共享同一个协方差矩阵。共享协方差正是让决策边界成为线性的原因因为二次项相互抵消了。这些假设成立时LDA 很省数据它是把所有类别汇到一起估一个协方差而不是每个类别各估一个。假设失效时症状很具体。假设各类别的协方差确实不同这种情况叫异方差一个类别是紧致的一团另一个是散开的一片云那么真正的贝叶斯边界就是二次的。LDA 被逼出来的线性边界会在离散程度不同的区域里系统性地误分类那个更宽的类别往往会在边界附近占据更多地盘。crate 里没有 QDA 可以替代所以对强异方差数据改用非线性分类器比如决策树、带 RBF 核的 SVM或者 KNN 分类器。一个类别如果是多峰的或明显非高斯的也可能破坏投影即便分类器还能给出可用的标签。在信任这个降维结果之前先看一眼transform的散点图。2.6.7. 奇异协方差与高维这一节讲的是特征共线或者n_features超过n_samples、协方差因此变得奇异时会发生什么。RustyML不会抛一个专门的奇异矩阵错误。crate 对n_features相对n_samples的关系没有设防。唯一的样本数检查是n_samples必须大于n_classes且每个类别至少要有 2 个样本。秩亏的协方差通过两条途径被吸收常开的对角脊把协方差往可逆方向推SVD和Eigen伪逆以及投影内部的白化步骤会把特征值低于相对容差的方向直接归零转而投影到协方差的非零子空间上而不是去除以一个接近 0 的数。usendarray::array;userustyml::machine_learning::LDA;fnmain(){// 第三列等于前两列之和所以特征矩阵秩亏// 共享协方差因此奇异letxarray![[1.0,2.0,3.0],[1.5,2.5,4.0],[2.0,3.0,5.0],[6.0,5.0,11.0],[6.5,4.5,11.0],[7.0,5.5,12.5],];letyarray![0,0,0,1,1,1];// 没有奇异矩阵错误对角脊和伪逆吸收了秩亏letmutldaLDA::default();lda.fit(x,y).unwrap();letcoordslda.transform(x).unwrap();// (6, 1)2 个类别 - 1 个成分assert!(coords.iter().all(|v|v.is_finite()));println!(reduced to {:?}, all finite,coords.shape());}投影还剩一种失败模式Error::Computation其context字段是 “Discriminant direction norm too small for stable projection”。只有当某根被选中的判别轴相对最大那根轴塌缩成一个近乎零的向量时RustyML 才会抛出这个错误。这道判据变成相对的是因为 RustyML 不再把各条轴归一化到单位长度。这种失败发生在你要求的成分数超过了数据在其退化子空间里能支撑的数量时。把它当成一次硬性的数值失败而不是校验层面的疏漏。想解决它就减少n_components或者加大正则化力度。当n_features大于n_samples或者数据严重共线时不要只依赖那条固定的1e-6对角脊改用Shrinkage::Auto。它对协方差的调理比那个最小稳定器好得多也是能用的投影和被噪声方向主导的投影之间的分水岭。2.6.8. LDA、PCA 与逻辑回归的对比有 3 个相关方法值得直接放在一起比较因为在它们之间做选择往往才是真正的问题。有监督边界 / 轴假设成分上限LDA是用标签类别分离度最大的方向高斯类别共享协方差n_classes - 1PCA否忽略标签总方差最大的方向无只用二阶矩n_features逻辑回归是只有线性决策边界对类别形状无假设不是降维器PCA 最大化方差从头到尾都不看y。所以如果数据的离散主要落在那个最不能分开类别的方向上它的头部成分也会毫不犹豫地和这个方向对齐。LDA 直接最大化分离度但被卡在n_classes - 1根轴上。当你想要一个有监督的低维视图或者一个快速的生成式分类器、并且类别大致高斯时用 LDA。当你想要无监督压缩或者需要多于n_classes - 1个维度时用 PCA。和逻辑回归的分野在于生成式对判别式。LDA 建模类别条件密度再套用贝叶斯规则。逻辑回归直接建模P(y | x)对每个类别的形状不作任何假设。当高斯和等协方差假设确实成立时LDA 能用更少的样本收敛到一个好边界。假设不成立时逻辑回归是更靠得住的线性分类器因为它从未假设过类别的形状。一种常见的做法是两个模型都保留让留出集上的准确率——用分类指标量出来——来决定用哪个。2.6.9. 错误、持久化与可复现性错误面遵循 crate 全局的错误处理约定。LDA::new(0)返回Error::InvalidParameter因为成分数必须为正。越界的Shrinkage::Manual同样返回InvalidParameter在构建时就被抓住。结构性问题在fit时浮现空的特征矩阵或者特征列数为 0 的矩阵返回Error::EmptyInput类别少于 2 个、n_samples不大于n_classes、某个类别只有 1 个样本、n_components过大这些情况全都返回Error::InvalidInputx里出现非有限值返回Error::NonFinitex/y长度不一致返回Error::DimensionMismatch。到 predict 或 transform 时空的输入矩阵同样返回Error::EmptyInput特征数不对返回Error::DimensionMismatch调用未拟合的模型返回Error::NotFitted(LDA)。usendarray::array;userustyml::error::Error;userustyml::machine_learning::{LDA,Shrinkage};fnmain(){// 手动收缩必须落在 [0, 1] 里RustyML 会在你构建模型时拒绝 1.5matchLDA::new(1).unwrap().with_shrinkage(Shrinkage::Manual(1.5)){Err(Error::InvalidParameter{name,..})println!(rejected parameter: {name}),_unreachable!(),}// n_components 只有在类别数已知后才有界所以一个过大的请求// 能通过构造器、却在 fit 时失败letxarray![[1.0,1.0],[1.5,0.8],[0.8,1.2],[5.0,5.0],[5.2,4.8],[4.8,5.2],];letyarray![0,0,0,1,1,1];// 2 个类别 - max_components 1letmutldaLDA::new(2).unwrap();matchlda.fit(x,y){Err(Error::InvalidInput(msg))println!(fit rejected: {msg}),_unreachable!(),}}一个已拟合的LDA用save_to_path和load_from_path序列化两者都走 postcard 二进制格式。.bin扩展名只是个文件名习惯不管你选什么扩展名格式始终是二进制的。这趟往返会还原类别、先验、均值、训练集总体均值、投影以及缓存的评分系数。加载回来的模型不用重新拟合就能给出一模一样的预测和变换结果。那个总体均值是transform开始做中心化时格式新增的字段。旧版本 RustyML 保存的二进制块加载不了需要用当前版本重新拟合、重新保存。格式细节见深入模型持久化。usendarray::array;userustyml::machine_learning::LDA;usestd::fs;fnmain(){letxarray![[1.0,1.0],[1.5,0.8],[0.8,1.2],[5.0,5.0],[5.2,4.8],[4.8,5.2],[9.0,1.0],[9.2,0.8],[8.8,1.2],];letyarray![0,0,0,1,1,1,2,2,2];letmutldaLDA::new(2).unwrap();lda.fit(x,y).unwrap();letbeforelda.predict(x).unwrap();letpathlda_model.bin;lda.save_to_path(path).unwrap();letloadedLDA::load_from_path(path).unwrap();letafterloaded.predict(x).unwrap();assert_eq!(before,after);fs::remove_file(path).unwrap();println!(round-trip predictions identical);}LDA 是完全确定的它不抽任何随机数所以同样的数据每次都产出同样的投影、同样的系数、同样的标签不需要设种子这一点和可复现性与随机种子里那些随机模型不同。唯一的运行间自由度是每根判别轴那个任意的符号这是特征向量的数学性质不是不确定性的表现。拟合会在工作量越过 crate 标定的阈值后把逐类别统计和最后的标签扫描并行化。并行永远不会改变结果它只改变结果送达的速度。更多内容见性能调优与并行。