Linear Discriminant Analysis

Linear discriminant analysis (LDA) is a classifier that models the images of each class as draws from a Gaussian distribution with its own mean and a covariance matrix shared by all classes, and assigns a new image to the class with the highest posterior probability. Because the covariance is shared, the log-posterior of each class is a linear function of the input, so the decision boundaries are hyperplanes (Murphy, 2022, §9.2.2). The same linear functions arise from a criterion with no probability model at all: Fisher (1936) asked which linear combination of measurements maximizes the squared difference between two class means relative to the variance within the classes, and his answer is the two-class LDA direction. The method also works as supervised dimensionality reduction, projecting the data onto at most $c - 1$ directions for $c$ classes. In computer vision it is the discriminant step of Fisherfaces in Eigenfaces and face recognition, a classifier on texture descriptors, the baseline on HOG features in the GTSRB traffic-sign competition, and, as Gaussian discriminant analysis fitted to network features, the density in DDU’s uncertainty estimate.

The method is also called Fisher’s linear discriminant, or Fisher’s linear discriminant analysis (FLDA) when the projection is meant; Murphy (§9.2) calls the family of Gaussian class-conditional classifiers Gaussian discriminant analysis (GDA). Murphy notes that the name is confusing twice over: LDA is a generative, not a discriminative, classifier, and “LDA” also abbreviates latent Dirichlet allocation, an unrelated model of word counts.

Problem

Given labelled images, each represented by a feature vector $\mathbf{x} \in \mathbb{R}^d$ (pixels of a patch, a descriptor, or a network’s features), predict the class of a new one. LDA solves the problem generatively: it fits a model of how each class produces its images, $p(\mathbf{x} \mid y = k)$, and turns it into posteriors with Bayes’ rule. A second problem it solves is choosing a few linear features that keep the classes apart, where principal component analysis keeps the directions of largest variance whatever their cause.

Inputs and Outputs

• Input: $n$ training vectors $\mathbf{x}_i \in \mathbb{R}^d$ with labels $y_i \in \{1, \dots, c\}$; optionally the class priors.
• Fitted parameters: the class means $\hat{\boldsymbol\mu}_k$, the priors $\hat\pi_k$, and one pooled within-class covariance $\hat\Sigma$, or equivalently one weight vector and one offset per class.
• Output: class posteriors and the most probable class for a new vector, or its projection onto at most $\min(c - 1, d)$ discriminant directions.

Intuition

Mahalanobis distance measures how far a point is from a class mean in units of the spread of the data in each direction (see similarity and distance measures). LDA assigns an image to the class whose mean is nearest in that distance, adjusted by the class priors. The scikit-learn user guide (§1.2.3) puts it another way: LDA is equivalent to first sphering the data, the whitening transform that makes the shared covariance the identity, and then picking the nearest class mean in Euclidean distance. Since the $c$ class means span an affine subspace of at most $c - 1$ dimensions, distances to them can be computed after projecting onto that subspace, which is where the dimensionality reduction comes from.

Fisher’s view gives the same direction for two classes from a different question. Project every image onto a line; the projected classes separate well when their means are far apart and each class is tight along the line. Directions in which the within-class variation is large, such as a texture’s slow brightness fluctuations, are down-weighted, and directions in which the classes differ but the nuisance variation is small are up-weighted. On image patches this is a matched filter that knows the noise: the example below shows that ignoring the covariance and comparing Euclidean distances to the true class means costs more than a third more errors than the Bayes rule on a textured background.

Algorithm

1. Estimate the priors $\hat\pi_k = n_k / n$ (or take given priors).
2. Estimate each class mean $\hat{\boldsymbol\mu}_k$.
3. Estimate the shared covariance from the deviations of each image from its own class mean, pooled over the classes, optionally shrunk toward a simpler matrix when $n$ is not much larger than $d$.
4. For each class compute the weight vector $\mathbf{w}_k = \hat\Sigma^{-1}\hat{\boldsymbol\mu}_k$ and the offset $w_{k0} = -\tfrac{1}{2}\hat{\boldsymbol\mu}_k^\top \hat\Sigma^{-1}\hat{\boldsymbol\mu}_k + \log \hat\pi_k$, by solving linear systems rather than forming the inverse.
5. Classify $\mathbf{x}$ as $\arg\max_k \mathbf{w}_k^\top\mathbf{x} + w_{k0}$; the softmax of these scores gives the posteriors.
6. For dimensionality reduction, solve the generalized eigenproblem $S_B\mathbf{w} = \lambda S_W\mathbf{w}$ for between- and within-class scatter matrices and keep the eigenvectors with the largest eigenvalues, at most $c - 1$ of them.

Mathematical Formulation

Gaussian discriminant analysis, quadratic and linear

With class-conditional densities $p(\mathbf{x} \mid y = k) = \mathcal{N}(\mathbf{x} \mid \boldsymbol\mu_k, \Sigma_k)$ and priors $\pi_k$, the log-posterior is, up to a constant shared by all classes (Murphy, Eq. 9.4),

$$\log p(y = k \mid \mathbf{x}) = \log \pi_k - \tfrac{1}{2}\log\lvert 2\pi\Sigma_k\rvert - \tfrac{1}{2}(\mathbf{x} - \boldsymbol\mu_k)^\top\Sigma_k^{-1}(\mathbf{x} - \boldsymbol\mu_k) + \text{const}.$$

Murphy calls this the discriminant function. Its boundaries are quadratic in $\mathbf{x}$, which gives quadratic discriminant analysis (QDA). If all classes share $\Sigma_k = \Sigma$, the term $\mathbf{x}^\top\Sigma^{-1}\mathbf{x}$ is the same for every class and drops out (Murphy, Eqs. 9.5–9.7):

$$\log p(y = k \mid \mathbf{x}) = \underbrace{\boldsymbol\mu_k^\top\Sigma^{-1}}_{\mathbf{w}_k^\top}\,\mathbf{x} \;\underbrace{-\tfrac{1}{2}\boldsymbol\mu_k^\top\Sigma^{-1}\boldsymbol\mu_k + \log\pi_k}_{w_{k0}} + \text{const},$$

which is linear in $\mathbf{x}$; these $\mathbf{w}_k$ and $w_{k0}$ are what scikit-learn stores as coef_ and intercept_ (user guide §1.2.2.2). The posterior is a softmax of linear scores, the form of multinomial logistic regression, which owns the comparison of the two ways of fitting it. For two classes (Murphy, Eqs. 9.13–9.16),

$$p(y = 1 \mid \mathbf{x}) = \sigma\big(\mathbf{w}^\top(\mathbf{x} - \mathbf{x}_0)\big), \qquad \mathbf{w} = \Sigma^{-1}(\boldsymbol\mu_1 - \boldsymbol\mu_0),$$

and the most probable class is 1 when $\mathbf{w}^\top\mathbf{x} > \tfrac{1}{2}\mathbf{w}^\top(\boldsymbol\mu_0 + \boldsymbol\mu_1) - \log(\pi_1/\pi_0)$. The priors move the threshold along $\mathbf{w}$ but not the direction.

The Bayes error of two Gaussian classes

When the model is true, the two-class LDA rule with the true parameters is the Bayes classifier, the rule with the lowest possible error (see machine learning for the Bayes error). Its error has a closed form, by a short derivation of our own. The projection $z = \mathbf{w}^\top\mathbf{x}$ with $\mathbf{w} = \Sigma^{-1}(\boldsymbol\mu_1 - \boldsymbol\mu_0)$ is Gaussian in each class with variance $\mathbf{w}^\top\Sigma\mathbf{w} = \Delta^2$ and means $\Delta^2$ apart, where

$$\Delta^2 = (\boldsymbol\mu_1 - \boldsymbol\mu_0)^\top\Sigma^{-1}(\boldsymbol\mu_1 - \boldsymbol\mu_0)$$

is the squared Mahalanobis distance between the class means. With equal priors the threshold is halfway, so a point is misclassified when it deviates from its own mean by more than $\Delta/2$ standard deviations in the wrong direction:

$$R^* = \Phi(-\Delta/2).$$

More generally, any linear rule “class 1 if $\mathbf{w}^\top\mathbf{x} > t$” has error $\tfrac{1}{2}\big[\Phi\big((t - \mathbf{w}^\top\boldsymbol\mu_1)/s\big) + \Phi\big((\mathbf{w}^\top\boldsymbol\mu_0 - t)/s\big)\big]$ with $s = \sqrt{\mathbf{w}^\top\Sigma\mathbf{w}}$, which lets the example below score a fitted rule exactly instead of on a test sample. Fisher used the same reasoning in §III of his paper: supposing a plant misclassified if its deviation exceeds half the difference between the species means of his compound, he read the frequency from a table of the normal distribution, about 2.79 per million for his ratio of 5.0018.

Fisher’s criterion

Fisher (1936, §II) took four flower measurements of fifty plants each of Iris setosa and I. versicolor, measured by Dr E. Anderson, and asked which linear function $X = \lambda_1 x_1 + \dots + \lambda_4 x_4$ maximizes the ratio of the difference between the species means to the standard deviation within species. With $d_p$ the differences of the means and $S_{pq}$ the sums of squares and products of deviations from the species means, he maximized $D^2/S$, where $D = \sum_p \lambda_p d_p$ and $S = \sum_p\sum_q \lambda_p\lambda_q S_{pq}$. Setting the derivatives to zero shows that the coefficients are proportional to the solution of $\sum_q S_{pq}\lambda_q = d_p$ (his Eq. 1), that is, $\boldsymbol\lambda \propto S^{-1}\mathbf{d}$. He did not assume normality to obtain the direction; normality enters only in §III, to estimate how often the compound would misclassify a plant.

In matrix form (Murphy, §9.2.6, following Bishop), project with $z = \mathbf{w}^\top\mathbf{x}$ and maximize

$$J(\mathbf{w}) = \frac{\mathbf{w}^\top S_B\mathbf{w}}{\mathbf{w}^\top S_W\mathbf{w}}, \qquad S_B = (\boldsymbol\mu_1 - \boldsymbol\mu_0)(\boldsymbol\mu_1 - \boldsymbol\mu_0)^\top, \qquad S_W = \sum_{k}\sum_{i : y_i = k}(\mathbf{x}_i - \boldsymbol\mu_k)(\mathbf{x}_i - \boldsymbol\mu_k)^\top .$$

The maximizer satisfies the generalized eigenproblem $S_B\mathbf{w} = \lambda S_W\mathbf{w}$, and because $S_B\mathbf{w}$ always points along $\boldsymbol\mu_1 - \boldsymbol\mu_0$, the solution is $\mathbf{w} \propto S_W^{-1}(\boldsymbol\mu_1 - \boldsymbol\mu_0)$ (Murphy, Eqs. 9.36–9.42). $S_W$ is the pooled covariance times a constant, so this is the Gaussian LDA direction. The two derivations differ in what they leave open: Fisher’s criterion fixes only the direction and says nothing about the threshold, while the Gaussian model fixes the threshold from the means and priors. The Gaussian derivation is not Fisher’s; it is the modern reading of his direction (the Cortes–Vapnik paper article discusses an account that credits Fisher with more).

Several classes and the projection

For $c$ classes, Murphy (§9.2.6.2) defines $S_W$ as above and $S_B = \sum_k n_k(\boldsymbol\mu_k - \boldsymbol\mu)(\boldsymbol\mu_k - \boldsymbol\mu)^\top$ around the overall mean, and maximizes $\lvert W^\top S_B W\rvert / \lvert W^\top S_W W\rvert$ over projection matrices $W$. The solution is $W = S_W^{-1/2}U$, where $U$ holds the leading eigenvectors of $S_W^{-1/2}S_B S_W^{-1/2}$; when $S_W$ is singular, Murphy suggests PCA first, which is the Fisherfaces recipe in Eigenfaces. $S_B$ has rank at most $c - 1$, so at most $c - 1$ directions carry any between-class variance; Murphy calls this a rather severe restriction. The directions are $S_W$-orthogonal rather than orthogonal (our reading of the generalized eigenproblem), so, unlike principal components, they are not an orthonormal basis of pixel space.

Fisher’s own three-species analysis (§VI) is not this eigenvector solution. He tested a genetic hypothesis that I. versicolor lies two-thirds of the way from I. setosa to I. virginica: he formed the compound of the species means with coefficients 4, 1 and −5 for I. virginica, I. versicolor and I. setosa, added the three species’ matrices of sums of squares and products with weights 16, 1 and 25 in the same order, and solved for a single discriminating function, whose differences between the species means came out close to the ratio 2 : 1 (§VI, Tables VIII and IX).

Discriminant as regression

Fisher (§IV) also pointed out that the two-class discriminant is a regression. If every plant is given a value of a variate $y$ that depends only on its species, $n_2/(n_1 + n_2)$ for one and $-n_1/(n_1 + n_2)$ for the other, the least-squares regression of $y$ on the measurements has coefficients proportional to the discriminant’s, and the analysis of variance of $y$ reproduces his test of significance. Two-class LDA’s direction can therefore be computed with any least-squares solver.

Which covariance estimate

The scale of the pooled covariance varies between sources and implementations, which changes the printed weights and offsets but not the two-class direction:

• Fisher worked with the sums $S_{pq}$ themselves (98 degrees of freedom for two species of fifty) and, for the variance of his compound, divided by 95, because three adjustable ratios of the coefficients had been fitted (§II–§III, Table V).
• Murphy’s maximum-likelihood estimate divides the pooled sum by $n$ (Eq. 9.21).
• scikit-learn’s lsqr and eigen solvers use $\sum_k \hat\pi_k C_k$, with each $C_k$ the biased covariance of class $k$; with the default priors this is the pooled sum divided by $n$. The svd solver does not use the covariance to fit the rule: it standardizes the within-class deviations and scales them by $1/\sqrt{n - c}$ before the singular value decomposition (discriminant_analysis.py, line 599 in 1.7.2).

With unequal priors the scale matters to the decision, since the log prior ratio is added to scores that scale with the inverse of the covariance estimate; with equal priors, as in the example, it does not.

Parameters

In scikit-learn 1.7.2, LinearDiscriminantAnalysis has these parameters, read from its source and docstring:

• solver ('svd' by default). svd computes the rule from two singular value decompositions without forming the covariance, which the docstring recommends for many features; lsqr solves $\hat\Sigma\mathbf{w}_k = \hat{\boldsymbol\mu}_k$ and does not support transform; eigen solves the generalized eigenproblem of between- and within-class covariance and supports both. Only lsqr and eigen accept shrinkage or a custom covariance estimator; svd with shrinkage raises an error.
• shrinkage (None). 'auto' uses the Ledoit–Wolf formula, and a number in $[0, 1]$ fixes the shrinkage coefficient.
• priors (None, the class proportions). They shift the offsets only.
• n_components (None, giving $\min(c - 1, d)$). It affects only transform.
• tol (1e-4). Singular values below it are discarded in the svd solver, which is how that solver copes with a singular within-class scatter.
• covariance_estimator (None). Any estimator with fit and covariance_, such as sklearn.covariance.OAS.
• store_covariance (False). Makes the svd solver also compute and store the weighted within-class covariance; the other solvers always store it.

QuadraticDiscriminantAnalysis estimates one covariance per class and offers reg_param instead of shrinkage.

Shrinkage

With fewer images than dimensions the pooled covariance is singular, and with not many more it is poorly conditioned, so its inverse amplifies estimation noise in the directions where it is smallest. Murphy (§4.5.2) shows that a MAP estimate under an inverse-Wishart prior is a convex combination $\lambda\Sigma_0 + (1 - \lambda)\hat\Sigma_{\text{mle}}$, and with a prior scatter proportional to the diagonal of the sample covariance this keeps the variances and shrinks the off-diagonal entries toward zero. He calls LDA with this estimate regularized discriminant analysis (Eq. 9.22), and notes that $\lambda$ can be chosen by cross-validation or by a closed-form formula, for which he cites Ledoit and Wolf (and Schaefer and Strimmer) and points to scikit-learn’s LedoitWolf.

scikit-learn’s shrunk estimators use a different target, a multiple of the identity: $(1 - \alpha)S + \alpha\,(\operatorname{tr}S/p)\,I$ (_shrunk_covariance.py). For shrinkage='auto', the 1.7.2 source standardizes each class’s data to unit variance, applies Ledoit–Wolf there, scales the result back by the standard deviations, and weights the classes by their priors (discriminant_analysis.py, lines 67–71, and _class_cov, line 124 on), so the target in the original units is the diagonal matrix of the class’s variances. For a numerical shrinkage, the same source shrinks the unstandardized empirical covariance (line 75, calling shrunk_covariance, lines 151–154 of _shrunk_covariance.py), whose target is the average variance times the identity; on features with unequal variances the two targets differ. The example below checks the 'auto' path; it does not exercise the numerical one, since all its pixels have equal variance. The user guide adds that for Gaussian data the Oracle Approximating Shrinkage estimator (OAS) has a smaller mean squared error than Ledoit and Wolf’s formula and should then classify better.

Complexity

By our count, fitting LDA with an explicit covariance costs $O(nd^2)$ to accumulate the scatter and $O(d^3)$ to factor it, plus $O(cd^2)$ to solve for the $c$ weight vectors; the svd solver’s decomposition of the $n \times d$ matrix of deviations costs $O(nd\min(n, d))$. Prediction is one matrix–vector product per class. QDA stores and factors $c$ covariances, $c$ times the parameters and work. On $8 \times 8$ patches all of this takes milliseconds; on raw images of tens of thousands of pixels the $d \times d$ covariance is the bottleneck, which is one reason Fisherfaces reduce the dimension with PCA first.

Implementation

Error against the exact Bayes error on textured edge patches

The example builds a two-class problem on $8 \times 8$ image patches whose Bayes error is known exactly. Class 0 is a stationary Gaussian texture: a Gaussian-correlated field with unit variance and a correlation length of 1.5 pixels, plus white sensor noise of standard deviation 0.3, so the covariance $\Sigma$ is known in closed form. Class 1 is the same texture plus a faint vertical step edge, a sigmoid of width 0.5 pixels centred between the two middle columns, whose amplitude is set so that $\Delta = 2\Phi^{-1}(0.9)$ and the Bayes error is exactly 0.10. Training sets have $n$ patches, half from each class, for $n$ from 16 to 4,096 around the dimension $d = 64$, with 20 seeds per size.

Each fitted linear rule is scored by its exact error under the true model, the formula above, rather than on a sampled test set. The methods are scikit-learn’s default LDA (svd solver), LDA with Ledoit–Wolf shrinkage (lsqr, shrinkage='auto'), an oracle that knows $\Sigma$ but estimates the class means, and, from 256 patches up, QDA, scored on 20,000 fresh test patches per seed together with the default LDA on the same patches. The code also checks the data generator (C1), the exact-error formula against 100,000 test patches (C2) and at the Bayes rule (C3), and its reading of scikit-learn’s shrinkage path (C4), and it compares Fisher’s direction $S_W^{-1}(\hat{\boldsymbol\mu}_1 - \hat{\boldsymbol\mu}_0)$, computed in NumPy, with scikit-learn’s coef_ and threshold. The hypotheses, fixed before any run:

• H1: every LDA variant approaches the Bayes error, with a mean excess of at most 0.01 at $n = 4096$ and a mean excess that does not increase from $n = 128$ on.
• H2: shrinkage helps most at small $n$: its gain over plain LDA is positive at every $n \le 2d$, largest at some $n \le 2d$, and at $n = 4096$ less than a tenth of its largest value.
• H3: Fisher’s direction and threshold equal scikit-learn’s up to scale, for both the svd and the unshrunk lsqr solver, wherever $S_W$ is invertible ($n \ge 96$).
• H4: QDA, which estimates two covariances where one is true, has a higher error than LDA at every size it runs.

The outputs came from NumPy 2.3.5, SciPy 1.16.3 and scikit-learn 1.7.2, single-threaded; the script takes about 4–5 seconds and prints the same output on repeated runs.

import os
os.environ["OMP_NUM_THREADS"] = os.environ["OPENBLAS_NUM_THREADS"] = os.environ["MKL_NUM_THREADS"] = "1"
import numpy as np
from scipy.stats import norm
from sklearn.discriminant_analysis import (LinearDiscriminantAnalysis as LDA,
                                           QuadraticDiscriminantAnalysis as QDA)

# 8x8 patches of a stationary Gaussian texture (Gaussian correlation, length 1.5 px)
# plus white sensor noise (sd 0.3); class 1 adds a faint vertical step edge.
P, D = 8, 64
yy, xx = np.mgrid[0:P, 0:P]
pos = np.stack([yy.ravel(), xx.ravel()], 1).astype(float)
r2 = ((pos[:, None, :] - pos[None, :, :]) ** 2).sum(-1)
SIGMA = np.exp(-r2 / (2 * 1.5**2)) + 0.3**2 * np.eye(D)
edge = 1 / (1 + np.exp(-(xx.ravel() - 3.5) / 0.5))
A = 2 * norm.ppf(0.9) / np.sqrt(edge @ np.linalg.solve(SIGMA, edge))   # Bayes error 0.10
MU = A * edge                                                          # class means 0 and MU
DELTA = np.sqrt(MU @ np.linalg.solve(SIGMA, MU))
BAYES = norm.cdf(-DELTA / 2)
L = np.linalg.cholesky(SIGMA)

def sample(n, cls, rng):
    return rng.standard_normal((n, D)) @ L.T + (MU if cls else 0.0)

def train_set(n, rng):
    X = np.vstack([sample(n // 2, 0, rng), sample(n // 2, 1, rng)])
    return X, np.r_[np.zeros(n // 2, int), np.ones(n // 2, int)]

def exact_err(w, c):          # rule "class 1 if w.x > c", error under the true model
    s = np.sqrt(w @ SIGMA @ w)
    return 0.5 * (norm.cdf((c - w @ MU) / s) + norm.cdf((0.0 - c) / s))

def rule(m):                  # scikit-learn's binary decision w.x + b > 0
    return m.coef_[0], -m.intercept_[0]

# lw_numpy ports the shrinkage formula of sklearn.covariance.ledoit_wolf_shrinkage
# (scikit-learn 1.7.2, sklearn/covariance/_shrunk_covariance.py), modified (one block,
# no blocking loop). scikit-learn licence:
#
# BSD 3-Clause License
#
# Copyright (c) 2007-2024 The scikit-learn developers.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# * Neither the name of the copyright holder nor the names of its
#   contributors may be used to endorse or promote products derived from
#   this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
def lw_numpy(X):              # Ledoit-Wolf as in sklearn.covariance.ledoit_wolf_shrinkage
    X = X - X.mean(0); n, p = X.shape
    S = X.T @ X / n; mu = np.trace(S) / p; X2 = X**2
    delta_ = ((X.T @ X) ** 2).sum() / n**2
    beta = ((X2.T @ X2).sum() / n - delta_) / (p * n)
    delta = (delta_ - 2 * mu * np.diag(S).sum() + p * mu**2) / p
    sh = 0.0 if beta == 0 else min(beta, delta) / delta
    out = (1 - sh) * S; out[np.diag_indices(p)] += sh * mu
    return out

print(f"edge amplitude a = {A:.6f}, Delta = {DELTA:.6f}, Bayes error = {BAYES:.6f}")
print(f"nearest-mean rule (true means, Sigma ignored): error = "
      f"{norm.cdf(-(MU @ MU) / (2 * np.sqrt(MU @ SIGMA @ MU))):.6f}")

# C1: the generator has the planned covariance and mean
rng = np.random.default_rng(99)
X0, X1 = sample(200_000, 0, rng), sample(200_000, 1, rng)
se = np.sqrt((np.outer(np.diag(SIGMA), np.diag(SIGMA)) + SIGMA**2) / 200_000)
zc = (np.abs(np.cov(X0, rowvar=False, bias=True) - SIGMA) / se).max()
zm = (np.abs(X1.mean(0) - MU) / np.sqrt(np.diag(SIGMA) / 200_000)).max()
print(f"C1 max |z| covariance {zc:.2f}, mean {zm:.2f} (<= 5): {zc <= 5 and zm <= 5}")
# C3: the exact-error formula at the Bayes rule gives Phi(-Delta/2)
wb = np.linalg.solve(SIGMA, MU)
c3 = abs(exact_err(wb, 0.5 * wb @ MU) - BAYES)
print(f"C3 |err(Bayes rule) - Phi(-Delta/2)| = {c3:.1e} (<= 1e-12): {c3 <= 1e-12}")

SIZES = [16, 32, 48, 64, 96, 128, 256, 512, 1024, 4096]
SEEDS = range(20)
err = {k: np.zeros((len(SIZES), 20)) for k in ("plain", "LW", "oracle")}
qda_diff = {}
h3_cos, h3_thr, c2_ok, c4_max = 0.0, 0.0, True, 0.0
print("\nC2 exact vs Monte Carlo (plain LDA, seed 0, 100,000 test patches)")
for i, n in enumerate(SIZES):
    dq = []
    for s in SEEDS:
        rng = np.random.default_rng([s, n])
        X, y = train_set(n, rng)
        plain = LDA().fit(X, y)
        lw = LDA(solver="lsqr", shrinkage="auto").fit(X, y)
        m0, m1 = X[y == 0].mean(0), X[y == 1].mean(0)
        wo = np.linalg.solve(SIGMA, m1 - m0)
        err["plain"][i, s] = exact_err(*rule(plain))
        err["LW"][i, s] = exact_err(*rule(lw))
        err["oracle"][i, s] = exact_err(wo, 0.5 * wo @ (m0 + m1))
        if n >= 96:   # H3: Fisher's direction S_W^-1 (m1 - m0) against scikit-learn
            Sw = sum((X[y == k] - X[y == k].mean(0)).T @ (X[y == k] - X[y == k].mean(0)) for k in (0, 1))
            wf = np.linalg.solve(Sw, m1 - m0)
            cf = 0.5 * wf @ (m0 + m1) / np.linalg.norm(wf)
            for m in (plain, LDA(solver="lsqr").fit(X, y)):
                w, c = rule(m)
                h3_cos = max(h3_cos, 1 - w @ wf / np.linalg.norm(w) / np.linalg.norm(wf))
                h3_thr = max(h3_thr, abs(c / np.linalg.norm(w) - cf))
        if s == 0:    # C2 and C4 on seed 0
            Xt, yt = train_set(100_000, rng)
            w, c = rule(plain)
            mc = np.mean((Xt @ w > c) != yt); ex = err["plain"][i, s]
            tol = 4 * np.sqrt(ex * (1 - ex) / 100_000); c2_ok &= abs(mc - ex) <= tol
            print(f"  n={n:5d}  exact {ex:.5f}  MC {mc:.5f}  |diff| {abs(mc - ex):.5f}  4SE {tol:.5f}")
            ref = sum(0.5 * (np.std(X[y == k], 0)[:, None] * lw_numpy((X[y == k] - X[y == k].mean(0))
                       / np.std(X[y == k], 0)) * np.std(X[y == k], 0)[None, :]) for k in (0, 1))
            c4_max = max(c4_max, np.abs(lw.covariance_ - ref).max())
        if n >= 256:  # QDA: Monte Carlo on 20,000 fresh test patches, paired with plain LDA
            Xt, yt = train_set(20_000, rng)
            q = QDA().fit(X, y)
            w, c = rule(plain)
            dq.append(np.mean(q.predict(Xt) != yt) - np.mean((Xt @ w > c) != yt))
    if dq:
        qda_diff[n] = np.array(dq)
print(f"C2 all within 4 SE: {c2_ok}")
print(f"C4 max |covariance_ - sum_k prior_k * rescaled LW(standardized class k)| = {c4_max:.1e} (<= 1e-10): {c4_max <= 1e-10}")

print(f"\nTest error (exact) minus Bayes error {BAYES:.5f}: mean [min, max] over 20 seeds")
print(f"{'n':>5s}  {'plain LDA (svd)':>27s}  {'Ledoit-Wolf (lsqr)':>27s}  {'oracle Sigma':>27s}")
for i, n in enumerate(SIZES):
    row = "  ".join(f"{(e := err[k][i] - BAYES).mean():.5f} [{e.min():.5f}, {e.max():.5f}]"
                    for k in ("plain", "LW", "oracle"))
    print(f"{n:5d}  {row}")

gap = err["plain"].mean(1) - err["LW"].mean(1)
print("\ngap = mean error plain - mean error LW: " + "  ".join(f"{n}: {g:.5f}" for n, g in zip(SIZES, gap)))
small = [i for i, n in enumerate(SIZES) if n <= 2 * D]
big = [SIZES.index(n) for n in (128, 256, 512, 1024, 4096)]
for k in ("plain", "LW", "oracle"):
    ex = err[k].mean(1) - BAYES
    print(f"H1 {k:6s}: mean excess at 4096 = {ex[-1]:.5f} (<= 0.01): {ex[-1] <= 0.01}; "
          f"non-increasing over 128..4096: {bool(np.all(np.diff(ex[big]) <= 0))}")
print(f"H2a gap > 0 at every n <= 128: {bool(np.all(gap[small] > 0))}")
print(f"H2b largest gap at n = {SIZES[int(np.argmax(gap))]} (<= 128): {SIZES[int(np.argmax(gap))] <= 128}")
print(f"H2c gap(4096) = {gap[-1]:.5f} < 0.1 x max gap = {0.1 * gap.max():.5f}: {gap[-1] < 0.1 * gap.max()}")
print(f"H3 max 1 - cos(Fisher, coef_) = {h3_cos:.1e} (<= 1e-10); max threshold difference = "
      f"{h3_thr:.1e} (<= 1e-8): {h3_cos <= 1e-10 and h3_thr <= 1e-8}")
for n, d in qda_diff.items():
    print(f"H4 n={n:5d}: QDA - LDA test error (MC, paired) mean {d.mean():.5f} [{d.min():.5f}, {d.max():.5f}]: {d.mean() > 0}")

edge amplitude a = 1.260517, Delta = 2.563103, Bayes error = 0.100000
nearest-mean rule (true means, Sigma ignored): error = 0.136562
C1 max |z| covariance 3.88, mean 2.58 (<= 5): True
C3 |err(Bayes rule) - Phi(-Delta/2)| = 0.0e+00 (<= 1e-12): True

C2 exact vs Monte Carlo (plain LDA, seed 0, 100,000 test patches)
  n=   16  exact 0.21224  MC 0.21196  |diff| 0.00028  4SE 0.00517
  n=   32  exact 0.19522  MC 0.19499  |diff| 0.00023  4SE 0.00501
  n=   48  exact 0.23785  MC 0.23632  |diff| 0.00153  4SE 0.00539
  n=   64  exact 0.44037  MC 0.43968  |diff| 0.00069  4SE 0.00628
  n=   96  exact 0.25368  MC 0.25436  |diff| 0.00068  4SE 0.00550
  n=  128  exact 0.19181  MC 0.19033  |diff| 0.00148  4SE 0.00498
  n=  256  exact 0.15290  MC 0.15449  |diff| 0.00159  4SE 0.00455
  n=  512  exact 0.12689  MC 0.12728  |diff| 0.00039  4SE 0.00421
  n= 1024  exact 0.11026  MC 0.10927  |diff| 0.00099  4SE 0.00396
  n= 4096  exact 0.10377  MC 0.10519  |diff| 0.00142  4SE 0.00386
C2 all within 4 SE: True
C4 max |covariance_ - sum_k prior_k * rescaled LW(standardized class k)| = 4.4e-16 (<= 1e-10): True

Test error (exact) minus Bayes error 0.10000: mean [min, max] over 20 seeds
    n              plain LDA (svd)           Ledoit-Wolf (lsqr)                 oracle Sigma
   16  0.12294 [0.06484, 0.20906]  0.10552 [0.04477, 0.22933]  0.15050 [0.06886, 0.22989]
   32  0.13378 [0.08533, 0.24608]  0.07940 [0.04271, 0.12785]  0.10481 [0.06366, 0.18567]
   48  0.16069 [0.10537, 0.29937]  0.05898 [0.03390, 0.13470]  0.07351 [0.05360, 0.11248]
   64  0.30037 [0.27132, 0.37983]  0.04955 [0.02884, 0.09185]  0.05824 [0.03936, 0.10140]
   96  0.18494 [0.13825, 0.23098]  0.04144 [0.02649, 0.06165]  0.04064 [0.02798, 0.05910]
  128  0.11213 [0.08364, 0.14881]  0.03417 [0.02481, 0.05146]  0.03121 [0.02072, 0.04423]
  256  0.04858 [0.03452, 0.06396]  0.02163 [0.01307, 0.03113]  0.01674 [0.01216, 0.02176]
  512  0.02483 [0.01676, 0.03851]  0.01468 [0.00946, 0.02654]  0.00894 [0.00627, 0.01112]
 1024  0.01080 [0.00728, 0.01398]  0.00799 [0.00578, 0.01011]  0.00400 [0.00300, 0.00519]
 4096  0.00298 [0.00217, 0.00379]  0.00270 [0.00198, 0.00342]  0.00108 [0.00080, 0.00146]

gap = mean error plain - mean error LW: 16: 0.01742  32: 0.05438  48: 0.10170  64: 0.25082  96: 0.14350  128: 0.07797  256: 0.02695  512: 0.01015  1024: 0.00281  4096: 0.00028
H1 plain : mean excess at 4096 = 0.00298 (<= 0.01): True; non-increasing over 128..4096: True
H1 LW    : mean excess at 4096 = 0.00270 (<= 0.01): True; non-increasing over 128..4096: True
H1 oracle: mean excess at 4096 = 0.00108 (<= 0.01): True; non-increasing over 128..4096: True
H2a gap > 0 at every n <= 128: True
H2b largest gap at n = 64 (<= 128): True
H2c gap(4096) = 0.00028 < 0.1 x max gap = 0.02508: True
H3 max 1 - cos(Fisher, coef_) = 3.3e-16 (<= 1e-10); max threshold difference = 2.0e-14 (<= 1e-8): True
H4 n=  256: QDA - LDA test error (MC, paired) mean 0.21047 [0.17255, 0.23295]: True
H4 n=  512: QDA - LDA test error (MC, paired) mean 0.15549 [0.13670, 0.18205]: True
H4 n= 1024: QDA - LDA test error (MC, paired) mean 0.09572 [0.07495, 0.11105]: True
H4 n= 4096: QDA - LDA test error (MC, paired) mean 0.03093 [0.02450, 0.03665]: True

All four checks held: the generator has the planned covariance and mean, the exact-error formula reproduces $\Phi(-\Delta/2)$ at the Bayes rule and agrees with 100,000 sampled test patches within four standard errors at every size, and scikit-learn’s 'auto' covariance equals the prior-weighted, rescaled Ledoit–Wolf estimates of the standardized classes to rounding error.

Ignoring the covariance costs over a third more errors. The nearest-mean rule with the true means, which ignores the texture’s covariance, errs on 13.7% of patches against the Bayes rule’s 10%. The edge’s mean pattern and the texture overlap at low spatial frequencies, where the texture has most of its power; $\Sigma^{-1}$ moves the weight to the sharp transition of the edge, where only the sensor noise competes (our reading of the construction).

H1 held for all three LDA variants: the mean excess error at 4,096 patches is 0.0030 for plain LDA, 0.0027 with shrinkage and 0.0011 for the oracle, all below 0.01, and none of the three means rises from 128 patches on. At 4,096 patches the two estimated-covariance variants are still about two to three times as far from the Bayes error as the oracle; by construction the only difference is that they estimate the 2,080 distinct entries of the $64 \times 64$ covariance from the data.

H2 held. Shrinkage lowers the mean error at every size up to 128 patches, its gain is largest at $n = 64$ (0.251), and at 4,096 patches it is 0.0003, less than a tenth of that. No hypothesis was set for the shape of plain LDA’s curve, but it is not monotone: its mean error is worst at $n = d = 64$, about 0.40, close to chance, and lower on both sides.

No hypothesis was set either for the comparison of shrinkage with the oracle, and the output shows a crossing: Ledoit–Wolf LDA has a lower mean error than the rule that knows $\Sigma$ exactly at every size up to 64 patches, nearly the same at 96, and a higher one from 128 on.

H3 held: for every seed and size from 96 patches up, Fisher’s direction and scikit-learn’s agree to $3.3 \times 10^{-16}$ in $1 - \cos$, and their normalized thresholds to $2 \times 10^{-14}$, for both solvers, although the svd solver scales its covariance by $1/(n - c)$ and lsqr by $1/n$; with equal priors the scale cancels.

H4 held: QDA’s error is higher than LDA’s on the same test patches at every size and in every seed, by 0.21 at 256 patches and still 0.031 at 4,096. The two classes share their covariance by construction, so QDA’s second covariance adds 2,080 parameters to estimate and nothing to gain.

The example ran as planned. The plan’s bars were written after planning scripts had computed exact and oracle quantities on other seeds; in particular the 0.01 bar of H1 was fixed after the oracle’s mean excess at 4,096 patches had been printed, so for the oracle H1 is an illustration rather than a test. Plain LDA, shrinkage LDA and QDA were first run in the pilot, whose output the first run reproduced unchanged. The example departs from the batch plan’s description in seven ways, each fixed in the plan before the pilot: the error is computed exactly from the Gaussian model rather than on a test sample, with a Monte Carlo check; the optional QDA row is included, from 256 patches up, with paired Monte Carlo error; the nearest-mean reference value was added; the check of scikit-learn’s shrinkage path (C4) was added; the comparison of Fisher’s direction with scikit-learn was extended to the threshold and to the lsqr solver, and restricted to sizes where $S_W$ is invertible; “Ledoit–Wolf” is scikit-learn’s shrinkage='auto' with the lsqr solver and “oracle” uses the true $\Sigma$ with estimated means; and the texture, noise, edge, Bayes error, sizes and seeds are construction choices made in the plan. After the first run, the comment block carrying scikit-learn’s licence was added above lw_numpy, which follows scikit-learn’s Ledoit–Wolf formula; no executable line changed, and the output is the same.

Properties and Behavior

• Optimal when the model holds. With Gaussian classes of shared covariance and known parameters, LDA is the Bayes classifier and its error is $\Phi(-\Delta/2)$ for two equally likely classes; with estimated parameters it converges toward that error as the training set grows, as H1 shows.
• Invariant to invertible linear transforms. The Mahalanobis distance is unchanged by any invertible linear transform of the data (similarity and distance measures), so, by our reasoning, the decisions of LDA without shrinkage do not depend on the units or the linear colour space of the features, whereas nearest-mean and k-NN rules in Euclidean distance do.
• Closed form. Fitting needs means and one covariance, no iterations, initialization or learning rate; the scikit-learn user guide lists closed-form solutions, inherent multiclass handling and no hyperparameters to tune among the attractions of LDA and QDA.
• Generative and discriminative readings agree on the boundary’s form. The posterior has the same linear-logit form as logistic regression; the logistic regression article sets out when each fit is preferable.
• Direction by least squares. By Fisher’s §IV, the two-class direction is the least-squares regression of coded labels on the features.
• At most $c - 1$ useful projections. The between-class scatter has rank $c - 1$, so transform returns at most $c - 1$ dimensions, and explained_variance_ratio_ measures the share of between-class variance along each, not of total variance (the eigenvalues of $S_W^{-1}S_B$ in eigen, the squared singular values of the scaled class means in svd, in the 1.7.2 source).
• Sample size around the dimension. In the example, plain LDA’s error was worst at $n = d$, and shrinkage LDA’s was far lower there.

Limitations

• Shared covariance. When classes differ in spread or orientation, the linear boundary is a compromise; QDA models the difference at the price of $c$ covariances, which, as H4 shows, costs accuracy when the data cannot support them.
• Gaussian, unimodal classes. A class made of several clusters, such as an object seen from several viewpoints, is summarized by one mean; a mixture per class (Gaussian mixture model) or a nonlinear classifier is needed.
• Sensitivity to outliers. The means and covariance are moment estimates; the logistic regression article reports Hastie and colleagues’ remark that LDA is not robust to gross outliers.
• High-dimensional raw images. With $d$ in the tens of thousands, the covariance cannot be estimated, or even stored comfortably, from realistic numbers of images; Fisherfaces first reduce the dimension with PCA, and even then Eigenfaces shows them magnifying nuisance variation absent from the gallery.
• Linear features only. Whatever separates the classes must be linearly visible in the features; on pixels this rarely holds beyond aligned patches.

Variants

• Quadratic discriminant analysis: one covariance per class, quadratic boundaries (Murphy, §9.2.1).
• Regularized discriminant analysis and shrinkage LDA: a covariance shrunk toward its diagonal (Murphy, Eq. 9.22) or toward a scaled identity with the Ledoit–Wolf or OAS coefficient (scikit-learn).
• Diagonal LDA and naive Bayes: a diagonal shared or per-class covariance, $O(cd)$ parameters instead of $O(cd^2)$ (Murphy, §9.2.4.2).
• Nearest centroid: with uniform priors LDA is the nearest-class-mean rule in Mahalanobis distance; replacing that distance with a learned one gives nearest-class-mean metric learning (Murphy, §9.2.5).
• Fisherfaces: PCA to $N - c$ dimensions, then Fisher’s discriminant to $c - 1$, on face images (Eigenfaces).
• Cross-view quadratic discriminant analysis (XQDA): a discriminant subspace and metric learned from intra- and extra-personal differences for person re-identification.
• Gaussian discriminant analysis on network features: fitted after training to give a feature density for out-of-distribution detection (DDU, in uncertainty estimation).

History

Fisher’s paper appeared in the Annals of Eugenics in 1936, illustrated on measurements of fifty plants each of three species of iris, two of them measured by Dr E. Anderson (§II); we read it in the scan of the University of Adelaide R. A. Fisher Digital Archive, reproduced from his Contributions to Mathematical Statistics (1950), and cite the journal pages. Fisher writes that use had already been made of the idea, at his suggestion, in craniometry, to sex differences in the mandible and to secular trends in cranial measurements (§I). In the same year Mahalanobis published his generalized distance between populations with a common dispersion matrix (similarity and distance measures), the quantity $\Delta$ that sets LDA’s Bayes error. This article takes the probabilistic reading with Gaussian classes, the multiclass eigenvector solution and the shrinkage estimators from Murphy’s textbook and scikit-learn; it does not trace their original papers, and cites Ledoit and Wolf’s only to identify them.