fix(lda): address review follow-ups from #456 (#457)

* fix(lda): address review follow-ups from #456

- Singularity threshold: T::epsilon() → relative tolerance (1e-4 * l_max), matching sklearn
- PartialEq: T::epsilon() → fixed 1e-10 tolerance for float comparison
- Rename field scalings → projection_matrix to avoid shadowing the public getter
- Add missing test coverage: X/y mismatch, zero n_components, wrong transform features
- Remove polyfill.io script tag (supply-chain concerns)
- Add blank line between lda and pca module declarations in mod.rs
- Bump version to 0.6.14 and move CHANGELOG entry accordingly

* fix(lda): raise PartialEq tolerance from 1e-10 to 1e-6 for f32 safety

1e-10 rounds to 0.0 in f32 (min positive normal ~1.2e-7), making
(a - b).abs() > 0.0 almost always true for distinct values. The 1e-6
floor works correctly for both f32 and f64.

Addresses review feedback on #457.

Author: Mec-iS

CHANGELOG.md

@@ -4,11 +4,11 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-## [0.6.13]
+## [0.6.14]
 ### Added
 - `decomposition/lda.rs`: `LDA`, linear discriminant analysis for supervised dimensionality reduction (#136). It projects the data onto the directions that best separate the classes, keeping `min(n_classes - 1, n_features)` components by default, and implements the `Transformer` interface next to `PCA`. Directions match scikit-learn's `LinearDiscriminantAnalysis(solver="eigen")` up to sign.
 
-## [0.6.12]
+## [0.6.13]
 ### Fixed
 - `xgboost/xgb_regressor.rs`: `XGRegressor::fit` no longer panics when `subsample` is less than 1.0 on a small dataset (#444). The sample for each tree now keeps a minimum of one row, as scikit-learn does for its own `subsample` parameter. Sample sizes of one row or more are unchanged.
 

Cargo.toml

@@ -2,7 +2,7 @@
 name = "smartcore"
 description = "Machine Learning in Rust."
 homepage = "https://smartcorelib.github.io/"
-version = "0.6.13"
+version = "0.6.14"
 authors = ["smartcore Developers"]
 edition = "2024"
 rust-version = "1.85"

src/decomposition/lda.rs

@@ -43,7 +43,6 @@
 //! * ["An Introduction to Statistical Learning", James G., Witten D., Hastie T., Tibshirani R., 4.4 Linear Discriminant Analysis](http://faculty.marshall.usc.edu/gareth-james/ISL/)
 //! * ["Pattern Classification", Duda R.O., Hart P.E., Stork D.G., 2nd ed., 3.8.3 Multiple Discriminant Analysis](https://www.wiley.com/en-us/Pattern+Classification%2C+2nd+Edition-p-9780471056690)
 //!
-//! <script src="https://polyfill.io/v3/polyfill.min.js?features=es6"></script>
 //! <script id="MathJax-script" async src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js"></script>
 use std::cmp::Ordering;
 use std::fmt::Debug;
@@ -63,7 +62,7 @@ use crate::numbers::realnum::RealNumber;
 #[derive(Debug)]
 pub struct LDA<T: Number + RealNumber, X: Array2<T> + EVDDecomposable<T>> {
     // Projection matrix, one column per kept discriminant direction (n_features x n_components).
-    scalings: X,
+    projection_matrix: X,
     // Generalized eigenvalue of every kept direction, largest first.
     eigenvalues: Vec<T>,
     // Number of features expected by `transform`.
@@ -72,18 +71,19 @@ pub struct LDA<T: Number + RealNumber, X: Array2<T> + EVDDecomposable<T>> {
 
 impl<T: Number + RealNumber, X: Array2<T> + EVDDecomposable<T>> PartialEq for LDA<T, X> {
     fn eq(&self, other: &Self) -> bool {
+        let tol = T::from(1e-6).unwrap();
         if self.n_features != other.n_features
             || self.eigenvalues.len() != other.eigenvalues.len()
             || self
-                .scalings
+                .projection_matrix
                 .iterator(0)
-                .zip(other.scalings.iterator(0))
-                .any(|(&a, &b)| (a - b).abs() > T::epsilon())
+                .zip(other.projection_matrix.iterator(0))
+                .any(|(&a, &b)| (a - b).abs() > tol)
         {
             return false;
         }
         for i in 0..self.eigenvalues.len() {
-            if (self.eigenvalues[i] - other.eigenvalues[i]).abs() > T::epsilon() {
+            if (self.eigenvalues[i] - other.eigenvalues[i]).abs() > tol {
                 return false;
             }
         }
@@ -276,8 +276,13 @@ impl<T: Number + RealNumber, X: Array2<T> + EVDDecomposable<T>> LDA<T, X> {
         let sw_evd = sw.evd(true)?;
         let u = sw_evd.V;
         let l = sw_evd.d;
+        let l_max = l
+            .iter()
+            .cloned()
+            .fold(T::zero(), |a, b| if b > a { b } else { a });
+        let tol = T::from(1e-10).unwrap().max(T::from(1e-4).unwrap() * l_max);
         for &li in &l {
-            if li <= T::epsilon() {
+            if li <= tol {
                 return Err(Failed::fit(
                     "Within-class scatter matrix is singular, provide more samples per class or fewer features",
                 ));
@@ -310,17 +315,17 @@ impl<T: Number + RealNumber, X: Array2<T> + EVDDecomposable<T>> LDA<T, X> {
             }
         });
 
-        let mut scalings = X::zeros(m, n_components);
+        let mut projection_matrix = X::zeros(m, n_components);
         let mut eigenvalues = vec![T::zero(); n_components];
         for (col, &src) in order.iter().take(n_components).enumerate() {
             eigenvalues[col] = g[src];
             for row in 0..m {
-                scalings.set((row, col), *directions.get((row, src)));
+                projection_matrix.set((row, col), *directions.get((row, src)));
             }
         }
 
         Ok(LDA {
-            scalings,
+            projection_matrix,
             eigenvalues,
             n_features: m,
         })
@@ -336,12 +341,12 @@ impl<T: Number + RealNumber, X: Array2<T> + EVDDecomposable<T>> LDA<T, X> {
                 ncols, self.n_features
             )));
         }
-        Ok(x.matmul(&self.scalings))
+        Ok(x.matmul(&self.projection_matrix))
     }
 
     /// Get the projection matrix, one column per discriminant direction.
     pub fn scalings(&self) -> &X {
-        &self.scalings
+        &self.projection_matrix
     }
 }
 
@@ -499,6 +504,27 @@ mod tests {
         assert!(result.is_err());
     }
 
+    #[test]
+    fn mismatched_x_y_is_rejected() {
+        let (x, _) = three_class_data();
+        let y_short = vec![0i32; x.shape().0 - 1];
+        assert!(LDA::fit(&x, &y_short, LDAParameters::default()).is_err());
+    }
+
+    #[test]
+    fn zero_n_components_is_rejected() {
+        let (x, y) = three_class_data();
+        assert!(LDA::fit(&x, &y, LDAParameters::default().with_n_components(0)).is_err());
+    }
+
+    #[test]
+    fn transform_wrong_features_is_rejected() {
+        let (x, y) = three_class_data();
+        let lda = LDA::fit(&x, &y, LDAParameters::default()).unwrap();
+        let bad = DenseMatrix::<f64>::zeros(3, 2);
+        assert!(lda.transform(&bad).is_err());
+    }
+
     #[cfg_attr(
         all(target_arch = "wasm32", not(target_os = "wasi")),
         wasm_bindgen_test::wasm_bindgen_test

src/decomposition/mod.rs

@@ -13,6 +13,7 @@
 
 /// LDA is a supervised approach that projects the data onto the directions that best separate the classes.
 pub mod lda;
+
 /// PCA is a popular approach for deriving a low-dimensional set of features from a large set of variables.
 pub mod pca;
 pub mod svd;