Add LDA for supervised dimensionality reduction (#136)

[ChrisJr404]

Fixes #136

Checklist

• My branch is up-to-date with main branch.
• Everything works and tested on latest stable Rust.
• Coverage and Linting have been applied

Current behaviour

decomposition only offers unsupervised reduction (PCA, SVD). There is no way to reduce the features using the class labels, so the directions that best separate the classes are not available.

New expected behaviour

Adds LDA, linear discriminant analysis for supervised dimensionality reduction, next to PCA in the same module. LDA::fit(x, y, params) learns the discriminant directions and transform projects onto them, so it plugs into the existing Transformer interface.

Implementation notes:

• It solves the generalized eigenproblem Sb w = lambda Sw w by whitening with the within-class scatter matrix. Sw is factored once with a symmetric eigen decomposition to build Sw^(-1/2), then Sw^(-1/2) Sb Sw^(-1/2) is symmetric and a second symmetric eigen decomposition gives the directions. That keeps everything on the existing symmetric evd path, no external numeric crates and no unsafe.
• By default it keeps min(n_classes - 1, n_features) components, matching the number of useful directions. with_n_components lets you ask for fewer.
• Singular within-class scatter (for instance more features than samples per class) returns Failed instead of producing garbage.
• The directions and the projection match scikit-learn 1.9.0 LinearDiscriminantAnalysis(solver="eigen") up to the usual per-axis sign, and the known-answer test pins those reference values (compared on magnitudes, the way the PCA tests do).

I scoped this to the dimensionality-reduction side to keep it close to PCA. Happy to follow up with an LDA classifier (predict) in a separate PR if that is wanted.

Change logs

Added

• decomposition::lda::LDA with LDAParameters / LDASearchParameters, a Transformer implementation, a scalings getter, doctest and unit tests.

[codecov]

Codecov Report

❌ Patch coverage is 72.50000% with 33 lines in your changes missing coverage. Please review.
✅ Project coverage is 63.86%. Comparing base (9eaae9e) to head (c037d14).
⚠️ Report is 184 commits behind head on main.

Files with missing lines	Patch %	Lines
src/decomposition/lda.rs	72.50%	33 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##             main     #456       +/-   ##
===========================================
+ Coverage   43.97%   63.86%   +19.88%     
===========================================
  Files          85       96       +11     
  Lines        7281     8341     +1060     
===========================================
+ Hits         3202     5327     +2125     
+ Misses       4079     3014     -1065     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[Mec-iS]

thanks. please consider adding tests as suggested by the coverage bot

[Mec-iS]

also: v0.6.13 has already been published.
New one is v0.6.14

[Mec-iS]

Thanks for this contribution — the implementation is clean and well-structured. Here are my notes from a full diff review.

---

✅ What works well

• The whitening approach (Sw^(-1/2) Sb Sw^(-1/2)) is mathematically correct and a good fit for the existing symmetric evd path. No new dependencies needed.
• min(n_classes - 1, n_features) default for n_components is correct and matches sklearn.
• The singularity guard (li <= T::epsilon()) is sensible, though see note below.
• The known-answer test against sklearn 1.9.0 with sign-agnostic magnitude comparison is the right approach, consistent with the PCA tests.
• serde round-trip test is appreciated.
• Doc-comment math and references are thorough.

---

🐛 Bugs / Correctness

1. Singularity threshold is too tight

if li <= T::epsilon() {

T::epsilon() is the machine epsilon (~2.2e-16 for f64), which is too small a floor — numerically near-singular matrices will pass this check and produce garbage whitening. A practical threshold like li < T::from(1e-10).unwrap() (or relative to the largest eigenvalue) would be safer. sklearn uses 1e-4 relative to the largest singular value.

2. CHANGELOG version mismatch

The CHANGELOG entry is filed under [0.6.13], but the maintainer noted that v0.6.13 has already been published. This entry should be bumped to [0.6.14].

---

⚠️ Missing coverage (33 uncovered lines per codecov)

The bot flagged 72.5% patch coverage. The paths most likely uncovered are:

• n_components = Some(c) validation branch (the c < 1 || c > max_components error arm)
• The n != y.shape() dimension mismatch error path
• The transform shape-mismatch error path
• The whitening / intermediate evd failure path

Suggested additions:

#[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); // wrong n_features
    assert!(lda.transform(&bad).is_err());
}

---

🔧 Minor / Style

3. PartialEq uses T::epsilon() for floating-point equality

Structural equality comparisons across eigenvectors are sensitive to floating-point order. Consider documenting this limitation or using a relative tolerance, consistent with how approx::relative_eq! is used in tests.

4. scalings field name shadows the getter

The private field scalings and the public method scalings() have the same identifier. This compiles fine in Rust (method vs field), but it may confuse readers. Consider renaming the field to scalings_ or projection_matrix to make the distinction explicit — this is also the convention used in the PCA struct.

5. pub mod lda; ordering in mod.rs

The new pub mod lda; line is placed before the pca entry with no blank line. The comment above it describes LDA but the /// PCA is a popular approach... doc comment for the pca module is now visually detached. A blank line between the two module declarations would keep the module-level docs tidy:

/// LDA is a supervised approach ...
pub mod lda;

/// PCA is a popular approach ...
pub mod pca;
pub mod svd;

6. polyfill.io script tag in doc comment

//! <script src="https://polyfill.io/v3/polyfill.min.js?features=es6"></script>

polyfill.io was acquired and its CDN has had documented supply-chain concerns. The PCA module uses MathJax directly without polyfill.io — recommend removing this line for consistency and safety.

---

💡 Follow-up suggestions (non-blocking)

• A predict method (nearest-centroid classification on the projected space) would complete the LDA story and aligns with the issue's original request. Great to scope that as a separate PR as proposed.
• Consider exposing explained_variance_ratio_ (each eigenvalue divided by the sum) analogous to PCA, useful for scree plots.
• The LDASearchParameters struct and iterator are complete but the default only iterates [None]. A short note in the docs clarifying how to construct a non-trivial grid would help users.

---

Overall this is a solid first pass. The core math is right, the API is idiomatic, and the tests pin reference values well. Addressing the CHANGELOG version, the singularity threshold, and the missing test coverage are the main blockers before merge.

[Mec-iS]

Hi @ChrisJr404 — thanks for this contribution, the core math and API shape are solid. I reviewed the full diff and have the following notes (ordered roughly by severity).

---

🐛 Bugs / Correctness

1. Singularity threshold is too tight

// current
if li <= T::epsilon() {

T::epsilon() is the machine epsilon (~2.2e-16 for f64), which is far too tight. Numerically near-singular Sw matrices will sail right past this guard and produce garbage whitening vectors. A safer floor is a small absolute threshold or a value relative to the largest eigenvalue:

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-6).unwrap() * l_max);
for &li in &l {
    if li <= tol {
        return Err(Failed::fit(
            "Within-class scatter matrix is singular …",
        ));
    }
}

scikit-learn uses 1e-4 relative to the largest singular value for the same purpose.

---

2. CHANGELOG entry is under the wrong version

v0.6.13 has already been published (as noted by the maintainer). This entry should be moved to ## [0.6.14].

---

⚠️ Missing test coverage (33 lines, ~27.5 % of the patch)

The uncovered paths are almost certainly:

Missing branch	Suggested test
n_components = Some(c) validation (c < 1)	zero_n_components_is_rejected
n != y.shape() length mismatch	mismatched_x_y_is_rejected
transform wrong-feature-count path	transform_wrong_features_is_rejected
evd failure / intermediate whitening failure	hard to trigger directly; partial coverage OK

Concrete additions you could drop into the existing tests module:

#[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); // wrong n_features
    assert!(lda.transform(&bad).is_err());
}

---

🔧 Minor / Style

3. PartialEq uses T::epsilon() for floating-point equality

Same issue as point 1: comparing eigenvectors with T::epsilon() will be brittle in practice. Consider using a small fixed tolerance (e.g. T::from(1e-10).unwrap()) or documenting the limitation explicitly, consistent with how approx::relative_eq! is used in the tests.

4. Private field scalings shadows the public getter scalings()

Rust resolves this without ambiguity, but the dual use of the identifier confuses readers. The PCA struct uses a distinct field name (components) from its public accessors. Consider renaming the field to projection_matrix (or scalings_ at minimum):

pub struct LDA<T, X> {
    projection_matrix: X,   // was `scalings`
    eigenvalues: Vec<T>,
    n_features: usize,
}

5. Missing blank line between pub mod lda and pub mod pca in mod.rs

The doc comment for the pca module is currently visually attached to the lda declaration. A blank line restores the original formatting:

/// LDA is a supervised approach …
pub mod lda;

/// PCA is a popular approach …
pub mod pca;
pub mod svd;

6. polyfill.io script tag in doc comment

//! <script src="https://polyfill.io/v3/polyfill.min.js?features=es6"></script>

polyfill.io was acquired in 2024 and its CDN has had documented supply-chain concerns. The PCA module already uses MathJax directly without polyfill.io. Recommend removing this line entirely for consistency and safety.

---

💡 Non-blocking suggestions

• explained_variance_ratio_ — exposing each eigenvalue divided by the total (analogous to PCA) would make scree plots trivial for users.
• LDASearchParameters docs — the default grid only iterates [None]. A short doc note explaining how to construct a non-trivial hyperparameter grid would help users adopting this with the existing grid-search infra.
• predict (nearest-centroid) — great call scoping that to a follow-up PR; it would complete the supervised story opened in Add Linear Discriminant Analysis (LDA) support #136.

---

Overall this is a clean first pass. The whitening algebra is correct, the sklearn reference test is the right approach, and the API fits idiomatically into the existing decomposition module. Addressing the CHANGELOG version, singularity threshold, and missing test coverage are the main items before merge.