# 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](../Chapter-01/1.5._The_Prelude_and_Imports.md)). 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:
| Method | Signature (abridged) | Returns |
| --- | --- | --- |
| `LDA::new` | `new(n_components: usize)` | `Result<LDA, Error>` |
| `LDA::default` | `default()` | `LDA` (auto components) |
| `with_solver` | `with_solver(DiscriminantSolver)` | `LDA` |
| `with_shrinkage` | `with_shrinkage(Shrinkage)` | `Result<LDA, Error>` |
| `fit` | `fit(&x, &y)` where `y: &Array1<i32>` | `Result<&mut LDA, Error>` |
| `predict` | `predict(&x)` | `Result<Array1<i32>, Error>` |
| `decision_function` | `decision_function(&x)` | `Result<Array2<f64>, Error>` |
| `predict_proba` | `predict_proba(&x)` | `Result<Array2<f64>, Error>` |
| `transform` | `transform(&x)` | `Result<Array2<f64>, Error>` |
| `fit_transform` | `fit_transform(&x, &y)` | `Result<Array2<f64>, Error>` |
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`](./2.1._Linear_Regression.md) 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.
```rust
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:
```rust
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]
}
```
```text
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.
| `DiscriminantSolver` | How it scores | When to use it |
| --- | --- | --- |
| `SVD` (default) | SVD pseudo-inverse of `Sigma`, small singular values truncated | General default. Handles rank deficiency well |
| `Eigen` | Symmetric eigendecomposition inverse, tiny eigenvalues zeroed | Well-conditioned covariance. Comparable to `SVD` |
| `LSQR` | Iterative Paige-Saunders least squares, forms no explicit inverse | Very high dimension, where an explicit inverse wastes memory |
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](../Chapter-06/6.0._Math_Utilities.md)). 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](./2.4._Decision_Trees.md), an [SVM with an RBF kernel](./2.5._Support_Vector_Machines.md), or a [KNN classifier](./2.3._K_Nearest_Neighbors.md). 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.
```rust
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.
| | Supervised? | Boundary / axes | Assumptions | Component ceiling |
| --- | --- | --- | --- | --- |
| **LDA** | Yes (uses labels) | Directions of max class separation | Gaussian classes, shared covariance | `n_classes - 1` |
| [**PCA**](./2.10._Principal_Component_Analysis.md) | No (ignores labels) | Directions of max total variance | None (just second moments) | `n_features` |
| [**Logistic regression**](./2.2._Logistic_Regression.md) | Yes | Linear decision boundary only | None on class shape | (not a reducer) |
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](../Chapter-05/5.2._Classification_Metrics.md).
## 2.6.9. Errors, persistence, and reproducibility
The error surface follows the crate-wide [error handling](../Chapter-01/1.6._Error_Handling.md) 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")`.
```rust
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](../Chapter-07/7.2._Model_Persistence_in_Depth.md) for the format details.
```rust
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](../Chapter-07/7.1._Reproducibility_and_Random_Seeds.md). 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](../Chapter-07/7.3._Performance_Tuning_and_Parallelism.md) for more on this topic.