2.6. Linear Discriminant Analysis

Linear Discriminant Analysis (LDA) does 2 jobs from a single fit. First, it works as a generative classifier. It models each class as a Gaussian with its own mean and a covariance shared across all classes. It then applies Bayes’ rule to pick the most probable label. Second, it works as a supervised dimensionality reducer. It projects the features onto the directions that separate the classes best. RustyML’s LDA gives both results from one fit call. predict, decision_function, and predict_proba serve the classifier. transform and fit_transform serve the projection. This page covers the API and the in-house linear algebra behind it. It also covers the Gaussian assumptions the model relies on, and what happens when the covariance matrix goes singular.

2.6.1. What the crate implements

RustyML implements LDA under rustyml::machine_learning::{LDA, Shrinkage, DiscriminantSolver}. The prelude also re-exports these types (see the prelude). The model takes an f64 feature matrix and i32 labels. LDA is a classifier, so labels are integer class ids, not floating-point values. The table below lists the API surface:

The crate does not implement QDA (quadratic discriminant analysis). It offers only the shared-covariance, linear-boundary variant. The getters get_classes, get_priors, get_means, get_overall_mean, and get_projection return None before fit runs. After fit runs, each getter returns its value wrapped in Some. Checking any of these getters tells you whether the model is trained.

The solver enum used to carry the plain name Solver. That name collided with the linear models’ own Solver enum. After a prelude glob import, Solver::GradientDescent resolved to the wrong enum and failed to compile. Both enums now carry names that match what they select. This crate uses DiscriminantSolver here and LeastSquaresSolver in the linear models.

2.6.2. A first classifier

Start with 3 well-separated clusters in 2D. LDA::default() resolves the component count automatically. You do not need to set a component count when you only want labels.

use ndarray::array;
use rustyml::machine_learning::LDA;

fn main() {
    // 3 tight, well-separated clusters, 3 samples each.
    let x = array![
        [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],
    ];
    let y = array![0, 0, 0, 1, 1, 1, 2, 2, 2];

    let mut lda = LDA::default();
    lda.fit(&x, &y).unwrap();

    let preds = lda.predict(&x).unwrap();
    println!("predictions: {:?}", preds);

    // Per-class discriminant scores (n_samples, n_classes) and softmax posteriors.
    let scores = lda.decision_function(&x).unwrap();
    let proba = lda.predict_proba(&x).unwrap();
    println!("scores shape:   {:?}", scores.shape());
    println!("proba shape:    {:?}", proba.shape());
    println!("class order:    {:?}", lda.get_classes().unwrap());
}

predict, decision_function, and predict_proba share a single core computation. For a sample x, LDA computes a linear discriminant score for each class. The formula is score_j(x) = x * Sigma^-1 * mu_j - 0.5 * mu_j * Sigma^-1 * mu_j + ln(prior_j). Here Sigma is the shared covariance, mu_j is the class mean, and prior_j is the class frequency in the training set. decision_function returns these raw scores as an (n_samples, n_classes) matrix. predict takes the per-row argmax of the scores and maps it back to the original label. predict_proba applies a numerically stable row-wise softmax to the same scores, so each row sums to 1. The column order of decision_function and predict_proba follows get_classes(). get_classes() sorts the labels in ascending order. Do not assume this order matches the order labels first appeared in y.

The priors are pure class frequencies (n_class / n_samples). For this balanced set, each prior is 1/3. Class imbalance changes this. The ln(prior_j) term shifts the boundary toward the rarer class, the way a Bayes classifier should. To get equal priors, balance the training set. The LDA constructor has no priors= override.

2.6.3. Supervised projection and the n_classes - 1 ceiling

The projection side is what distinguishes LDA from PCA. LDA finds directions that maximize between-class scatter relative to within-class scatter. Only n_classes - 1 such directions carry non-zero separability, because the between-class scatter matrix has rank at most n_classes - 1. The crate enforces this limit. The usable component count is min(n_classes - 1, n_features), resolved at fit time. LDA::default() (or n_components = None) takes that maximum. LDA::new(k) requests exactly k components and fails at fit if k exceeds the cap.

A 2D scatter plot therefore needs at least 3 classes. 2 classes give only 1 discriminant axis. fit_transform trains the model and projects the data in a single call:

use ndarray::array;
use rustyml::machine_learning::{DiscriminantSolver, LDA, Shrinkage};

fn main() {
    let x = array![
        [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],
    ];
    let y = array![0, 0, 0, 1, 1, 1, 2, 2, 2];

    // 3 classes -> at most 2 discriminant axes. Project straight to 2D.
    let mut lda = LDA::new(2)
        .unwrap()
        .with_solver(DiscriminantSolver::Eigen)
        .with_shrinkage(Shrinkage::Auto)
        .unwrap();

    let coords = lda.fit_transform(&x, &y).unwrap(); // (9, 2)
    println!("embedding shape: {:?}", coords.shape());
    for (row, &label) in coords.outer_iter().zip(y.iter()) {
        println!("class {label}: ({:.3}, {:.3})", row[0], row[1]);
    }

    // The projection matrix maps n_features -> n_components.
    let w = lda.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]

The projection matrix W has shape (n_features, n_components). transform computes (x - xbar) * W. Here xbar is the mean of the whole training matrix. This value matches scikit-learn’s xbar_, and get_overall_mean() reads it back. RustyML matches scikit-learn’s (X - xbar_) @ scalings_ formula exactly. That match brings 2 details, and both changed recently. First, the projection is centered by the training mean. A transformed training set sits around the origin, not around the raw feature values. Second, the axes keep the scale the whitening produced. RustyML does not renormalize the axes to unit L2 norm. That scale is not cosmetic. It gives the projected data unit within-class covariance, the property the whole construction targets. A projection read from an older, unit-norm axis now shows a different offset and a different scale.

RustyML orders the columns by class separability, with the largest generalized eigenvalue first. The first component is therefore the single most discriminative axis. This ordering helps if you later want a 1D reduction from the same fit. A discriminant axis and its negation separate the classes equally well. The sign of any column is therefore arbitrary. Treat the orientation of a plotted axis as meaningless. Do not read directionality into it.

transform and fit_transform perform the same reduction step. fit_transform(&x, &y) runs fit(&x, &y) followed by transform(&x) in a single call.

2.6.4. Solvers and the in-house linear algebra

DiscriminantSolver selects how RustyML inverts the shared covariance for the classifier’s scoring coefficients Sigma^-1 * mu_c. This is an important and easy-to-miss fact. The solver choice only affects predict, decision_function, and predict_proba. The transform projection does not depend on the solver. RustyML always derives it from a symmetric eigendecomposition of the within-class covariance. Switching solvers changes the decision boundary’s numerics. It never changes the 2D embedding.

All 3 solvers reach identical labels on well-separated data. They differ only in how they handle an ill-conditioned Sigma. SVD and Eigen both build a pseudo-inverse. This pseudo-inverse discards any direction whose singular value or eigenvalue falls below a relative tolerance. A rank-deficient covariance therefore degrades gracefully instead of failing. LSQR avoids the inverse entirely. It solves each Sigma * coef_c = mu_c system iteratively. This choice saves memory when n_features is large.

The linear algebra is entirely in-house pure Rust. RustyML carries no LAPACK or nalgebra dependency (see prefer in-house numerics). The symmetric eigendecomposition uses Householder tridiagonalization followed by the implicit-shift QL iteration, the classic EISPACK tred2 and tql2 pair. The SVD uses a one-sided Jacobi method that computes small singular values to high relative accuracy. LSQR uses the Golub-Kahan bidiagonalization with running Givens rotations. Each method is deterministic to machine precision on the small-to-mid covariance matrices LDA produces. That determinism is why LDA needs no random seed.

2.6.5. Shrinkage and covariance regularization

With few samples per feature, the raw sample covariance is a noisy estimate. Its inverse amplifies that noise. Shrinkage pulls Sigma toward a scaled identity matrix. This trades a small amount of bias for a large reduction in variance. Shrinkage has 2 forms:

• Shrinkage::Auto: the Ledoit-Wolf closed-form optimal intensity. RustyML computes it from the data, with no tuning knob to set.
• Shrinkage::Manual(alpha): an explicit alpha in [0, 1]. A value of 0 applies no shrinkage. A value of 1 shrinks fully to the identity target. with_shrinkage validates alpha and returns Error::InvalidParameter if it falls outside [0, 1] or is not finite. This validation is why with_shrinkage returns a Result.

2 subtleties matter here. First, RustyML always adds a tiny fixed diagonal ridge to the covariance, about 1e-6 times the average variance. It adds this ridge regardless of the shrinkage choice, only to keep the factorizations stable. Shrinkage::Manual(0.0) and no shrinkage therefore produce identical results. Second, shrinkage here is independent of the solver, and this is a real divergence from scikit-learn. scikit-learn refuses shrinkage with its svd solver. RustyML applies shrinkage to the covariance before any solver runs. DiscriminantSolver::SVD with Shrinkage::Auto is therefore valid, and often the best combination. Choose Shrinkage::Auto as the default whenever you have collinear features or a thin sample-to-feature ratio.

2.6.6. Assumptions, and the symptoms when they break

LDA’s optimality rests on 2 assumptions. Each class follows a Gaussian distribution. All classes share one covariance matrix. The shared covariance is what makes the decision boundary linear, because the quadratic terms cancel. When these assumptions hold, LDA is data-efficient. It estimates a single covariance from all classes pooled together, instead of a separate covariance per class.

When the assumptions break, the symptoms are specific. Say the classes have genuinely different covariances, a condition called heteroscedastic data. One class forms a tight blob, and another forms a broad cloud. The true Bayes boundary is then quadratic. LDA’s forced-linear boundary systematically misclassifies samples in the region where the spreads differ. The wide class tends to dominate territory near the boundary. The crate has no QDA as an alternative. For strongly heteroscedastic data, use a non-linear classifier instead, such as a decision tree, an SVM with an RBF kernel, or a KNN classifier. A class that is multimodal or clearly non-Gaussian can also break the projection, even when the classifier still produces usable labels. Check a transform scatter plot before you trust the reduction.

2.6.7. Singular covariance and high dimensions

This section covers what happens when features are collinear, or when n_features exceeds n_samples, so the covariance becomes singular. RustyML does not raise a dedicated singular-matrix error. The crate places no guard on n_features against n_samples. The only sample-count checks are that n_samples must exceed n_classes, and each class needs at least 2 samples. A rank-deficient covariance gets absorbed in 2 ways. The always-on diagonal ridge nudges the covariance toward invertibility. The SVD and Eigen pseudo-inverse, and the whitening step inside the projection, zero out any direction whose eigenvalue falls below a relative tolerance. This projects onto the covariance’s non-null subspace instead of dividing by a value near 0.

use ndarray::array;
use rustyml::machine_learning::LDA;

fn main() {
    // The third column equals the sum of the first two, so the feature
    // matrix is rank-deficient. The shared covariance is therefore
    // singular.
    let x = array![
        [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],
    ];
    let y = array![0, 0, 0, 1, 1, 1];

    // No singular-matrix error: the ridge and pseudo-inverse absorb the
    // rank deficiency.
    let mut lda = LDA::default();
    lda.fit(&x, &y).unwrap();

    let coords = lda.transform(&x).unwrap(); // (6, 1): 2 classes -> 1 component
    assert!(coords.iter().all(|v| v.is_finite()));
    println!("reduced to {:?}, all finite", coords.shape());
}

The projection has 1 remaining failure mode: Error::Computation, with context set to “Discriminant direction norm too small for stable projection”. RustyML raises this error only when a selected discriminant axis collapses to a near-zero vector, relative to the largest axis. The guard became relative once RustyML stopped renormalizing the axes to unit length. This failure happens when you request more components than the data can support in its degenerate subspace. Treat it as a hard numerical failure, not a validation slip. Reduce n_components, or regularize harder, to fix it. For n_features greater than n_samples, or for heavily collinear data, do not rely on the fixed 1e-6 ridge alone. Turn on Shrinkage::Auto instead. It conditions the covariance far better than the minimal stabilizer, and it makes the difference between a usable projection and one dominated by noise directions.

2.6.8. LDA versus PCA versus logistic regression

3 related methods are worth a direct comparison, because choosing among them is usually the real question.

PCA maximizes variance without ever looking at y. Its top components can therefore align with the direction that separates classes least, if that is where the spread lives. LDA maximizes separation directly, but it is capped at n_classes - 1 axes. Use LDA for a supervised low-dimensional view, or for a fast generative classifier, when the classes are roughly Gaussian. Use PCA for unsupervised compression, or when you need more than n_classes - 1 dimensions.

Against logistic regression, the split is generative versus discriminative. LDA models the class-conditional densities and applies Bayes’ rule. Logistic regression models P(y | x) directly and makes no assumption about the shape of each class. When the Gaussian and equal-covariance assumptions genuinely hold, LDA converges to a good boundary from fewer samples. When they do not hold, logistic regression is the more dependable linear classifier, because it never assumed a class shape. A common approach keeps both models and lets held-out accuracy decide, measured with the classification metrics.

2.6.9. Errors, persistence, and reproducibility

The error surface follows the crate-wide error handling conventions. LDA::new(0) returns Error::InvalidParameter, because components must be positive. An out-of-range Shrinkage::Manual also returns InvalidParameter, caught at build time. Structural problems surface at fit. An empty feature matrix, or a matrix with 0 feature columns, returns Error::EmptyInput. Fewer than 2 classes, n_samples not greater than n_classes, a class with only 1 sample, and an over-large n_components all return Error::InvalidInput. A non-finite entry in x returns Error::NonFinite. An x/y length mismatch returns Error::DimensionMismatch. At predict or transform time, an empty input matrix also returns Error::EmptyInput, a wrong feature count returns Error::DimensionMismatch, and calling an unfitted model returns Error::NotFitted("LDA").

use ndarray::array;
use rustyml::error::Error;
use rustyml::machine_learning::{LDA, Shrinkage};

fn main() {
    // Manual shrinkage must lie in [0, 1]. RustyML rejects 1.5 when you
    // build the model.
    match LDA::new(1).unwrap().with_shrinkage(Shrinkage::Manual(1.5)) {
        Err(Error::InvalidParameter { name, .. }) => println!("rejected parameter: {name}"),
        _ => unreachable!(),
    }

    // n_components is only bounded once the class count is known, so an over-large
    // request survives the constructor and fails at fit.
    let x = array![
        [1.0, 1.0], [1.5, 0.8], [0.8, 1.2],
        [5.0, 5.0], [5.2, 4.8], [4.8, 5.2],
    ];
    let y = array![0, 0, 0, 1, 1, 1]; // 2 classes -> max_components = 1
    let mut lda = LDA::new(2).unwrap();
    match lda.fit(&x, &y) {
        Err(Error::InvalidInput(msg)) => println!("fit rejected: {msg}"),
        _ => unreachable!(),
    }
}

A fitted LDA serializes with save_to_path and load_from_path. Both methods use the postcard binary format. The .bin extension is only a filename convention. The format stays binary regardless of the extension you choose. The round trip restores the classes, the priors, the means, the overall training mean, the projection, and the cached scoring coefficients. A loaded model predicts and transforms identically, with no re-fit needed. The overall mean is a field the format gained when transform started centering the projection. Blobs saved by an older version of RustyML will not load. Re-fit and re-save the model with the current version. See model persistence in depth for the format details.

use ndarray::array;
use rustyml::machine_learning::LDA;
use std::fs;

fn main() {
    let x = array![
        [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],
    ];
    let y = array![0, 0, 0, 1, 1, 1, 2, 2, 2];

    let mut lda = LDA::new(2).unwrap();
    lda.fit(&x, &y).unwrap();
    let before = lda.predict(&x).unwrap();

    let path = "lda_model.bin";
    lda.save_to_path(path).unwrap();
    let loaded = LDA::load_from_path(path).unwrap();
    let after = loaded.predict(&x).unwrap();

    assert_eq!(before, after);
    fs::remove_file(path).unwrap();
    println!("round-trip predictions identical");
}

LDA is fully deterministic. It draws no random numbers, so the same data yields the same projection, the same coefficients, and the same labels on every run. It needs no seed to set, unlike the stochastic models in reproducibility and random seeds. The only run-to-run freedom is the arbitrary sign of each discriminant axis. That sign is a mathematical property of eigenvectors, not a sign of nondeterminism. Fitting parallelizes the per-class statistics and the final label scan, once the work clears the crate’s calibrated threshold. Parallelism never changes the result. It only changes how fast the result arrives. See performance tuning and parallelism for more on this topic.