Linear and Quadratic Discriminant Analysis

In the figures below, modeling the same observations with a shared covariance gives a straight separator, while allowing separate covariances gives a curved one. Suppose we need to classify continuous observations into two or more groups, and we have a modest number of labeled examples. K-nearest neighbors classified observations by local proximity; naive Bayes modeled class-conditional feature distributions under conditional independence; logistic regression gave us a probabilistic model of $p(y \mid x)$ via a link function. Discriminant analysis takes a third route. Instead of modeling the label given the features, we model the features given the label: each class is assumed to be a multivariate Gaussian cloud, and we classify by asking which cloud most plausibly produced the observation.

When the class distributions are Gaussian and their means, covariances, and priors are known, this approach gives the Bayes-optimal rule under zero-one misclassification loss. In practice, LDA and QDA estimate those parameters, so a fitted classifier does not inherit an exact finite-sample optimality guarantee. Linear discriminant analysis (LDA) and quadratic discriminant analysis (QDA) differ in one covariance assumption: LDA uses a shared covariance and has linear score differences, while QDA allows class-specific covariances and quadratic score differences. Those quadratic differences can produce curved boundaries, including degenerate cases discussed below.

A name disclaimer is worth making early. In modern machine learning, "LDA" is also used as the acronym for latent Dirichlet allocation, a text-topic model in unsupervised learning. That model is unrelated to the discriminant analysis we develop here. In this chapter, LDA always means linear discriminant analysis. Fisher (1936) introduced a supervised projection criterion. For two Gaussian classes with distinct means and a shared positive-definite covariance, a maximizing Fisher direction is normal to the Bayes decision hyperplane.

Naive Bayes models class-conditional distributions under conditional independence. Here we allow features to be correlated within each class through a full covariance matrix. Whether those matrices are shared or class specific determines whether the score differences are linear or quadratic.

We will work throughout in the two-class setting for the core derivations and extend to $K$ classes when it makes the mathematics clearer. Notation is standard: $\mathbf{x} \in \mathbb{R}^d$ is the feature vector, $y \in \{1, 2, \ldots, K\}$ the class label, $\pi_k$ the prior probability of class $k$, $\boldsymbol{\mu}_k$ the class-conditional mean, and $\boldsymbol{\Sigma}_k$ the class-conditional covariance matrix. We assume $\pi_k > 0$ and symmetric positive-definite covariance matrices throughout the density-based derivations, so the logarithms, determinants, and inverses are defined. In the two-class derivations and examples, we label the classes 0 and 1.

With the true class means, priors, and shared positive-definite covariance, LDA is the Bayes-optimal classifier for Gaussian classes under zero-one loss. Here, "optimal" means minimizing the expected probability of misclassification under that generative model: at each feature vector, predicting the class with the largest true posterior minimizes conditional error. The Gaussian discriminant derivation explains this population rule. Estimating its parameters introduces uncertainty, so fitted LDA need not attain the Bayes error in a finite sample.

LDA shares covariance information across classes. It estimates one $d \times d$ covariance matrix, $K$ mean vectors, and $K - 1$ priors. For $d = 100$ features and two classes ($K = 2$), the covariance contributes $\frac{100(101)}{2} = 5{,}050$ distinct entries, the means contribute $2 \cdot 100 = 200$, and the priors contribute $K - 1 = 1$, for exactly $5{,}251$ free generative parameters. Ordinary binary logistic regression has $d + 1 = 101$ coefficients including its intercept; fitting two separate one-vs-rest models would store $2(d + 1) = 202$ coefficients, with redundant representations for two complementary classes. The larger generative parameter count does not by itself determine classification efficiency. Efron (1975) studied the asymptotic relative classification efficiency of logistic regression and normal discriminant analysis for two multivariate normal classes with different means and a common covariance. That comparison concerns asymptotic classification performance, rather than a guarantee of more accurate finite-sample coefficient estimates.

QDA gives the corresponding population Bayes rule for Gaussian classes with their own positive-definite covariances, provided the true parameters and priors are used. Fitting it requires $K$ full covariance estimates, each with $\frac{d(d+1)}{2}$ free parameters, instead of one. When class covariances differ, its quadratic score differences can represent curved decision boundaries that a linear classifier cannot express. For independent Gaussian samples with a truly common covariance, the pooled covariance has lower entrywise sampling variance than an unbiased class estimate when that class has at least two observations and the other classes add residual degrees of freedom.

For dense full-covariance matrix kernels, computing within-class scatter from $N$ observations once class subsets are available costs $O(Nd^2)$. QDA then factors $K$ covariance matrices at a cost of $O(Kd^3)$ and scores one new observation in $O(Kd^2)$. LDA factors a single shared covariance and can precompute its class weights: the additional work is $O(d^3 + Kd^2)$, after which a prediction costs $O(Kd)$. These are matrix-kernel bounds; solver choices can change them. The primitive fits below additionally scan all labels once per class, adding O(NK) class-mask comparisons, plus separate input validation and np.unique work. Large feature counts also impose quadratic covariance storage, so an adequate sample size alone does not make either dense method practical.

Skew, heavy tails, or multiple modes within a class can make a Gaussian density a poor approximation and can harm probability calibration or classification. LDA and QDA do not adapt their density family to represent arbitrary class shapes. A misspecified density can still give the correct hard-label boundary, so a failure of Gaussianity is a reason to examine held-out performance rather than proof that every prediction is wrong. These methods remain useful when their assumptions are plausible or as baselines for comparison.

Covariance estimation is a major difficulty when the feature dimension is large relative to the available observations per class. In high-dimensional settings, $d$ can rival the number of observations per class. A centered class sample has rank at most $N_k-1$. For $N_k>1$, its empirical covariance is singular when $N_k\le d$; with only one observation, the unbiased covariance is undefined. The ordinary unregularized QDA density formula is then undefined because the covariance inverse does not exist. Even when $N_k > d$, dependent features or estimation noise can leave the covariance singular or poorly conditioned. LDA pools within-class scatter across classes, but its pooled covariance can have rank no greater than $N - K$ and is subject to the same rank and conditioning checks.

For two classes, a nondegenerate LDA boundary is one hyperplane regardless of sample size. It cannot represent a bounded winning region surrounded by the other class or a winning region split across two separated tails. In the multiclass case, each pairwise tie equation is linear and each class winning region is an intersection of half-spaces. QDA can represent additional geometries, but it estimates $K\,\frac{d(d+1)}{2}$ covariance parameters rather than LDA's $\frac{d(d+1)}{2}$. Estimating these separate matrices can add uncertainty in small samples; both methods remain sensitive to outliers, and there is no universal ordering of their fitted errors.

Unregularized LDA and QDA are equivariant under invertible affine transformations when their means and covariances are transformed consistently and calculations are exact. This algebraic property does not require the observed data to be exactly Gaussian. Feature standardization can improve numerical conditioning when units differ substantially, but it does not remove collinearity or guarantee an accurate covariance estimate. Regularization toward a fixed identity target can also change the rule under rescaling, so preprocessing and shrinkage must be considered together.

We first combine Gaussian class densities with Bayes' rule, then examine which score terms cancel under a shared covariance. This separates the population decision rule from the covariance estimates and numerical calculations used to fit it.

A discriminant classifier begins with a generative thought experiment. Imagine that nature first chooses a class according to its prior probability, and then draws a feature vector from a probability density that belongs to that class. To classify a new observation, we reverse the order of that story: instead of asking "what does class $k$ tend to look like?", we ask "how plausible is this particular $\mathbf{x}$ under each class, after accounting for how common each class is?" That reversal is Bayes' rule, and the entire geometry of LDA and QDA follows from the density we choose.

The first modeling decision is therefore to describe each class by a full density rather than by a single prototype. An ordinary Euclidean nearest-centroid rule scores a point by its Euclidean distance from a class mean, so equal-length displacements have equal distance scores regardless of direction. A class density can instead describe one direction as tightly concentrated and another as broad or tilted. The multivariate Gaussian is the canonical choice for this job. It is specified by a center and a covariance, its level sets are ellipsoids, and its log-density is quadratic, which is exactly what makes the later decision rules tractable. We assume

$$p(\mathbf{x} \mid y = k) = \mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k),$$

where $\mathcal{N}$ denotes the multivariate normal density

$$\mathcal{N}(\mathbf{x} \mid \boldsymbol{\mu}, \boldsymbol{\Sigma}) = (2\pi)^{-d/2} |\boldsymbol{\Sigma}|^{-1/2} \exp\!\left(-\tfrac{1}{2} (\mathbf{x} - \boldsymbol{\mu})^\top \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu})\right).$$

where:

• $\boldsymbol{\mu}$: the mean vector of the Gaussian, which locates the center of the ellipsoidal cloud
• $\boldsymbol{\Sigma}$: the covariance matrix, which sets the shape and orientation of the cloud
• $|\boldsymbol{\Sigma}|$: the covariance determinant; its square root is proportional to the volume of an ellipsoid with fixed Mahalanobis radius
• $(\mathbf{x} - \boldsymbol{\mu})^\top \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu})$: the squared Mahalanobis distance from $\mathbf{x}$ to $\boldsymbol{\mu}$
• $d$: the number of feature dimensions

The mean $\boldsymbol{\mu}_k$ locates the center of a class cloud, and the covariance $\boldsymbol{\Sigma}_k$ determines its shape and orientation. For ellipsoids with the same Mahalanobis radius, volume is proportional to $\sqrt{|\boldsymbol{\Sigma}|}$. A larger determinant therefore means a larger such ellipsoid and a lower Gaussian peak density. The quadratic form $(\mathbf{x} - \boldsymbol{\mu})^\top \boldsymbol{\Sigma}^{-1}(\mathbf{x} - \boldsymbol{\mu})$ is the squared Mahalanobis distance. It measures distance relative to the class covariance, so a large Euclidean displacement along a broad direction can still have a small standardized distance. The density decreases continuously with this distance. There is no single hard neighborhood boundary built into the Gaussian model.

A density alone is not yet a classifier, because classes may occur with different frequencies. Bayes' rule combines the class-conditional density with the prior:

$$p(y = k \mid \mathbf{x}) = \frac{p(\mathbf{x} \mid y = k)\, \pi_k}{p(\mathbf{x})},$$

where:

• $p(y = k \mid \mathbf{x})$: the posterior probability that the observation belongs to class $k$
• $p(\mathbf{x} \mid y = k)$: the class-conditional density of $\mathbf{x}$ given class $k$
• $\pi_k$: the prior probability of class $k$
• $p(\mathbf{x}) = \sum_{j=1}^K p(\mathbf{x} \mid y = j)\, \pi_j$: the marginal density of $\mathbf{x}$ under the generative model, shared by all classes and normalizing the posteriors to sum to one

The denominator is the same for every class, so it cannot change which class has the largest posterior probability. The rule "predict the class with the largest posterior probability" is therefore equivalent to "predict the class with the largest $p(\mathbf{x} \mid y = k)\, \pi_k$." This is a small algebraic observation with a large practical consequence: we never need to compute the marginal density in order to classify.

Comparing products of densities is awkward, and in high dimensions the products can underflow. Taking logarithms preserves the argmax because the logarithm is monotonic. So

$$\log\!\left[p(\mathbf{x} \mid y = k)\, \pi_k\right] = -\tfrac{d}{2} \log(2\pi) - \tfrac{1}{2} \log |\boldsymbol{\Sigma}_k| - \tfrac{1}{2} (\mathbf{x} - \boldsymbol{\mu}_k)^\top \boldsymbol{\Sigma}_k^{-1} (\mathbf{x} - \boldsymbol{\mu}_k) + \log \pi_k.$$

where:

• $\mathbf{x}$: the feature vector being classified
• $\boldsymbol{\mu}_k$: the mean vector of class $k$
• $\boldsymbol{\Sigma}_k$: the covariance matrix of class $k$
• $\pi_k$: the prior probability of class $k$
• $d$: the number of features
• $|\boldsymbol{\Sigma}_k|$: the determinant of the class covariance

The first term $-\tfrac{d}{2} \log(2\pi)$ is identical for every $k$ and drops out of any comparison. Dropping it, we obtain the discriminant function

$$\delta_k(\mathbf{x}) = -\tfrac{1}{2} (\mathbf{x} - \boldsymbol{\mu}_k)^\top \boldsymbol{\Sigma}_k^{-1} (\mathbf{x} - \boldsymbol{\mu}_k) - \tfrac{1}{2} \log |\boldsymbol{\Sigma}_k| + \log \pi_k,$$

where:

• $\delta_k(\mathbf{x})$: the discriminant score for class $k$ at point $\mathbf{x}$
• $\mathbf{x}$: the feature vector
• $\boldsymbol{\mu}_k$: the mean vector of class $k$
• $\boldsymbol{\Sigma}_k$: the covariance matrix of class $k$
• $\pi_k$: the prior probability of class $k$
• $(\mathbf{x} - \boldsymbol{\mu}_k)^\top \boldsymbol{\Sigma}_k^{-1} (\mathbf{x} - \boldsymbol{\mu}_k)$: the squared Mahalanobis distance from $\mathbf{x}$ to class $k$'s mean, in units of class $k$'s own covariance
• $|\boldsymbol{\Sigma}_k|$: the determinant of the class covariance, penalizing diffuse class clouds

The Bayes decision rule is then $\hat{y} = \arg\max_k \delta_k(\mathbf{x})$.

Each term has a clear interpretation, and it is worth reading them as a sequence of questions. The first term rewards closeness to the class mean, measured in the class's own Mahalanobis metric: does $\mathbf{x}$ fall inside the class's ellipsoidal neighborhood? The second term is a normalization penalty: at a fixed feature dimension, a larger covariance determinant lowers the Gaussian peak. The density at a particular point also depends on its Mahalanobis distance from the class mean. The third term rewards frequently observed classes, pulling the decision toward a prior that may be strong. Classification is thus a balance between geometric fit, concentration, and prior belief.

Two facts about this derivation are useful. First, adding the same quantity to every class score, even a quantity depending on $\mathbf{x}$, does not change the predicted class. Second, the difference between two class scores is exactly the corresponding log-posterior ratio, because all common terms cancel:

$$\log \frac{p(y = k \mid \mathbf{x})}{p(y = k' \mid \mathbf{x})} = \delta_k(\mathbf{x}) - \delta_{k'}(\mathbf{x}),$$

where:

• $k, k'$: the two classes being compared
• $\delta_k(\mathbf{x}), \delta_{k'}(\mathbf{x})$: the discriminant scores for classes $k$ and $k'$ at point $\mathbf{x}$, each defined up to the same additive constant
• $p(y = k \mid \mathbf{x}), p(y = k' \mid \mathbf{x})$: the posterior probabilities of classes $k$ and $k'$ given $\mathbf{x}$

The right-hand side is a score difference, and its zero set is the pairwise tie locus. With two classes, a change in the winning label can occur only at a tie; a degenerate tie need not separate two open decision regions. With more than two classes, a tie between $k$ and $k'$ contributes to their shared decision boundary only where their common score is at least as large as every other class score. We now examine how the covariance assumptions determine these score differences.

The LDA assumption is that all classes share the same covariance matrix: $\boldsymbol{\Sigma}_k = \boldsymbol{\Sigma}$ for all $k$. The shared-covariance condition is a model assumption about the data-generating process. It says that all class clouds have the same shape, orientation, and volume, and differ only in where they are centered and how common they are. Under that assumption, the Gaussian ellipsoids are translated copies of one another. The discriminant simplifies because the $-\tfrac{1}{2} \log |\boldsymbol{\Sigma}|$ term is the same for all classes and drops out:

$$\delta_k(\mathbf{x}) = -\tfrac{1}{2} (\mathbf{x} - \boldsymbol{\mu}_k)^\top \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu}_k) + \log \pi_k.$$

where:

• $\delta_k(\mathbf{x})$: the discriminant score for class $k$ at point $\mathbf{x}$
• $\mathbf{x}$: the feature vector being classified
• $\boldsymbol{\mu}_k$: the mean vector of class $k$
• $\boldsymbol{\Sigma}$: the shared covariance matrix common to all classes
• $\pi_k$: the prior probability of class $k$
• $(\mathbf{x} - \boldsymbol{\mu}_k)^\top \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu}_k)$: the squared Mahalanobis distance from $\mathbf{x}$ to class $k$'s mean in the shared metric

We now expand the quadratic form. In matrix notation,

$$(\mathbf{x} - \boldsymbol{\mu}_k)^\top \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu}_k) = \mathbf{x}^\top \boldsymbol{\Sigma}^{-1} \mathbf{x} - 2\, \boldsymbol{\mu}_k^\top \boldsymbol{\Sigma}^{-1} \mathbf{x} + \boldsymbol{\mu}_k^\top \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}_k.$$

After multiplying the expansion by $-\tfrac{1}{2}$, the term $-\tfrac{1}{2}\mathbf{x}^\top\boldsymbol{\Sigma}^{-1}\mathbf{x}$ is identical for every class and can be dropped. The other terms retain their signs and coefficients; multiplying an entire score by a negative number would reverse its ranking. The shared quadratic term is the cancellation that removes curvature from the score differences. What remains is

$$\delta_k(\mathbf{x}) = \boldsymbol{\mu}_k^\top \boldsymbol{\Sigma}^{-1} \mathbf{x} - \tfrac{1}{2}\, \boldsymbol{\mu}_k^\top \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}_k + \log \pi_k.$$

This is a linear function of $\mathbf{x}$. Define

$$\mathbf{w}_k = \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}_k, \qquad b_k = -\tfrac{1}{2}\, \boldsymbol{\mu}_k^\top \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}_k + \log \pi_k,$$

where:

• $\mathbf{w}_k$: the covariance-adjusted weight vector for class $k$
• $b_k$: the class-$k$ bias (intercept), combining the mean-quadratic term and the log prior
• $\boldsymbol{\Sigma}, \boldsymbol{\mu}_k, \pi_k$: the shared covariance, class-$k$ mean, and class-$k$ prior, as defined earlier

so that $\delta_k(\mathbf{x}) = \mathbf{w}_k^\top \mathbf{x} + b_k$ is a linear score. For two classes, the decision boundary between them is $\delta_0(\mathbf{x}) = \delta_1(\mathbf{x})$, or

$$\mathbf{w}^\top \mathbf{x} + b = 0,$$

with

$$\mathbf{w} = \boldsymbol{\Sigma}^{-1}(\boldsymbol{\mu}_0 - \boldsymbol{\mu}_1), \qquad b = -\tfrac{1}{2}(\boldsymbol{\mu}_0 + \boldsymbol{\mu}_1)^\top \boldsymbol{\Sigma}^{-1} (\boldsymbol{\mu}_0 - \boldsymbol{\mu}_1) + \log \frac{\pi_0}{\pi_1}.$$

where:

• $\mathbf{w}$: the covariance-adjusted normal to the boundary, pointing toward increasing class-0-versus-class-1 score
• $b$: the bias term that shifts the boundary based on the class means and priors
• $\boldsymbol{\Sigma}, \boldsymbol{\mu}_0, \boldsymbol{\mu}_1, \pi_0, \pi_1$: the shared covariance, class means, and class priors, as defined earlier

The sign convention matters: $\mathbf{w}^\top\mathbf{x} + b$ equals $\delta_0(\mathbf{x})-\delta_1(\mathbf{x})$, so a positive value favors class 0 after accounting for the priors. In general, $\mathbf{w}$ is not parallel to $\boldsymbol{\mu}_0-\boldsymbol{\mu}_1$. If the means differ and the priors are equal, the boundary is the hyperplane through their midpoint on which their Mahalanobis distances are equal. Increasing the prior of one class expands its winning region and moves this hyperplane toward the other class. If the means are identical, the normal is zero and the priors alone decide, or every point is tied when the priors are also equal.

Fisher's original formulation sought a linear combination $\mathbf{w}^\top\mathbf{x}$ maximizing the ratio of between-class to within-class variance. In the two-class case with nonsingular within-class covariance, the between-class scatter of the projected values is $(\mathbf{w}^\top(\boldsymbol{\mu}_0-\boldsymbol{\mu}_1))^2$ and the within-class scatter is $\mathbf{w}^\top\boldsymbol{\Sigma}\mathbf{w}$, so the ratio is

$$J(\mathbf{w}) = \frac{\bigl(\mathbf{w}^\top(\boldsymbol{\mu}_0-\boldsymbol{\mu}_1)\bigr)^2}{\mathbf{w}^\top\boldsymbol{\Sigma}\mathbf{w}},$$

where:

• $\mathbf{w}$: the direction onto which the data are projected
• $\boldsymbol{\mu}_0-\boldsymbol{\mu}_1$: the difference between the class means, setting the between-class separation
• $\boldsymbol{\Sigma}$: the shared (within-class) covariance, setting the within-class spread

When the class means differ, a maximizing nonzero direction is proportional to $\boldsymbol{\Sigma}^{-1}(\boldsymbol{\mu}_0-\boldsymbol{\mu}_1)$. Its sign and nonzero scale do not change the variance ratio. If the means coincide, $J(\mathbf w)=0$ for every nonzero direction, so this criterion does not select an informative direction. Under the shared-covariance Gaussian model with distinct means, the maximizing direction gives the generative classifier when paired with its prior-dependent threshold. The direction alone does not specify that threshold.

LDA has $Kd$ mean parameters and $\frac{d(d+1)}{2}$ shared symmetric-covariance parameters. Estimating the class priors adds $K-1$ free parameters under $\sum_k \pi_k=1$, giving $Kd+\frac{d(d+1)}{2}+(K-1)$ in total. For fixed $K$, the covariance term grows quadratically in $d$, while the mean term grows linearly. In the multiclass setting, LDA needs at most $\min(d,K-1)$ discriminant directions. The centered class means have at most $K-1$ independent directions, and applying an invertible covariance adjustment does not change this rank. The posterior-relevant directions in the original feature coordinates are the span of $\boldsymbol{\Sigma}^{-1}(\boldsymbol{\mu}_k-\boldsymbol{\mu}_{\mathrm{ref}})$, where $\boldsymbol{\mu}_{\mathrm{ref}}$ is one chosen reference-class mean. Projecting onto these directions retains the Gaussian LDA score differences; projecting onto the raw means need not.

QDA relaxes the shared-covariance assumption: each class has its own $\boldsymbol{\Sigma}_k$. This can describe one class that is narrow and another that is broad or differently oriented. The shared-covariance cancellation is no longer available in general. Some terms can still coincide for particular parameter choices, but the general score must retain the full form

$$\delta_k(\mathbf{x}) = -\tfrac{1}{2} (\mathbf{x} - \boldsymbol{\mu}_k)^\top \boldsymbol{\Sigma}_k^{-1} (\mathbf{x} - \boldsymbol{\mu}_k) - \tfrac{1}{2} \log |\boldsymbol{\Sigma}_k| + \log \pi_k.$$

where:

• $\delta_k(\mathbf{x})$: the discriminant score for class $k$ at point $\mathbf{x}$
• $\mathbf{x}$: the feature vector being classified
• $\boldsymbol{\mu}_k$: the mean vector of class $k$
• $\boldsymbol{\Sigma}_k$: the class-specific covariance matrix for class $k$
• $\pi_k$: the prior probability of class $k$
• $(\mathbf{x} - \boldsymbol{\mu}_k)^\top \boldsymbol{\Sigma}_k^{-1} (\mathbf{x} - \boldsymbol{\mu}_k)$: the squared Mahalanobis distance from $\mathbf{x}$ to class $k$'s mean, in class $k$'s own metric
• $|\boldsymbol{\Sigma}_k|$: the determinant of the class $k$ covariance, penalizing diffuse class clouds

When $\boldsymbol{\Sigma}_k^{-1}$ varies with $k$, the $\mathbf{x}^\top\boldsymbol{\Sigma}_k^{-1}\mathbf{x}$ contribution is class dependent and must be retained. The log-determinant penalty must likewise be retained whenever the covariance determinants differ; different covariance matrices can still have the same determinant. The general discriminant is quadratic in $\mathbf{x}$, with linear and second-order terms in its coordinates. Class-specific covariance therefore permits boundaries beyond a single flat cut, while particular parameter choices can still produce simpler or degenerate boundaries.

For two classes ($k = 0$ and $k = 1$), the boundary $\delta_0(\mathbf{x}) = \delta_1(\mathbf{x})$ is

$$-\tfrac{1}{2} \mathbf{x}^\top (\boldsymbol{\Sigma}_0^{-1} - \boldsymbol{\Sigma}_1^{-1}) \mathbf{x} + (\boldsymbol{\mu}_0^\top \boldsymbol{\Sigma}_0^{-1} - \boldsymbol{\mu}_1^\top \boldsymbol{\Sigma}_1^{-1}) \mathbf{x} + c_0 = 0,$$

which is the equation $\delta_0(\mathbf{x}) - \delta_1(\mathbf{x}) = 0$ written with the class-0 terms leading.

where:

• $\mathbf{x}$: the feature vector at which the boundary is evaluated
• $\boldsymbol{\Sigma}_0, \boldsymbol{\Sigma}_1$: the covariance matrices of classes 0 and 1 (both assumed invertible)
• $\boldsymbol{\mu}_0, \boldsymbol{\mu}_1$: the mean vectors of classes 0 and 1
• $c_0$: the constant term, defined explicitly below

The constant $c_0$ collects the terms that do not depend on $\mathbf{x}$. Writing out $\delta_0(\mathbf{x}) - \delta_1(\mathbf{x}) = 0$ in full (the form matching the displayed boundary equation) and keeping only the terms without $\mathbf{x}$ gives

$$c_0 = -\tfrac{1}{2} \boldsymbol{\mu}_0^\top \boldsymbol{\Sigma}_0^{-1} \boldsymbol{\mu}_0 + \tfrac{1}{2} \boldsymbol{\mu}_1^\top \boldsymbol{\Sigma}_1^{-1} \boldsymbol{\mu}_1 - \tfrac{1}{2} \log |\boldsymbol{\Sigma}_0| + \tfrac{1}{2} \log |\boldsymbol{\Sigma}_1| + \log \frac{\pi_0}{\pi_1}.$$

where:

• $\mathbf{x}$: the feature vector
• $\boldsymbol{\Sigma}_0, \boldsymbol{\Sigma}_1$: the covariance matrices of classes 0 and 1
• $\boldsymbol{\mu}_0, \boldsymbol{\mu}_1$: the class mean vectors
• $\pi_0, \pi_1$: the class prior probabilities
• $c_0$: the constant term collecting all contributions that do not depend on $\mathbf{x}$ (the mean-quadratic and log-determinant terms plus the log prior ratio)

The boundary equation has the form of a general quadratic in $\mathbf{x}$: the symmetric matrix coefficient of the second-order term is $-\tfrac{1}{2}(\boldsymbol{\Sigma}_0^{-1} - \boldsymbol{\Sigma}_1^{-1})$, the linear coefficient is $\boldsymbol{\mu}_0^\top \boldsymbol{\Sigma}_0^{-1} - \boldsymbol{\mu}_1^\top \boldsymbol{\Sigma}_1^{-1}$, and $c_0$ is the scalar constant.

In two dimensions, a quadratic tie equation can describe conic sections such as ellipses, hyperbolas, and parabolas; in higher dimensions, the corresponding objects are quadrics. The eigenstructure of $\boldsymbol{\Sigma}_0^{-1}-\boldsymbol{\Sigma}_1^{-1}$ controls the quadratic curvature, but the linear and constant terms also matter. Together they determine the boundary's location, whether it is degenerate, and whether it exists at all. Some QDA models produce bounded class regions, with a narrow class winning near its center and a broad class winning in the tails. QDA can also reduce to a linear rule when the quadratic difference vanishes.

To see where the constants come from, write $\delta_0(\mathbf{x}) - \delta_1(\mathbf{x})$ term by term, consistent with the $\delta_0 - \delta_1 = 0$ sign convention used above. The $\mathbf{x}$-independent contribution is

$$-\tfrac{1}{2} \boldsymbol{\mu}_0^\top \boldsymbol{\Sigma}_0^{-1} \boldsymbol{\mu}_0 + \tfrac{1}{2} \boldsymbol{\mu}_1^\top \boldsymbol{\Sigma}_1^{-1} \boldsymbol{\mu}_1 - \tfrac{1}{2} \log |\boldsymbol{\Sigma}_0| + \tfrac{1}{2} \log |\boldsymbol{\Sigma}_1| + \log \frac{\pi_0}{\pi_1},$$

matching $c_0$ above, while the $\mathbf{x}$-dependent part yields the quadratic and linear coefficients already stated.

The number of free parameters in QDA is $K \cdot d$ (means) plus $K - 1$ (priors) plus $K \cdot \frac{d(d+1)}{2}$ (covariances), so its covariance term grows $K$ times larger than LDA's. For independent Gaussian samples with a truly common covariance, the pooled covariance has lower entrywise sampling variance than an unbiased class estimate when that class has at least two observations and the other classes add residual degrees of freedom. Whether their flexibility improves held-out prediction depends on the population covariances and the available samples; parameter count alone does not determine the prediction error.

The covariance parameter count grows quadratically in the feature dimension and is multiplied by the number of classes for QDA. The next figure isolates that multiplier for a three-class problem, holding means and priors aside so the covariance cost is the only changing quantity.

The formulas above describe the population rule when the means, covariances, and priors are known. In practice, these parameters must be estimated. The maximum-likelihood class priors are the empirical proportions, $\hat{\pi}_k=N_k/N$, where $N_k$ counts observations in class $k$ and $N=\sum_j N_j$. The maximum-likelihood class means are $\hat{\boldsymbol{\mu}}_k=\frac{1}{N_k}\sum_{i:y_i=k}\mathbf{x}_i$. For covariance estimates, this chapter uses the unbiased sample convention rather than the Gaussian maximum-likelihood convention:

$$\hat{\boldsymbol{\Sigma}}_k = \frac{1}{N_k - 1} \sum_{i: y_i = k} (\mathbf{x}_i - \hat{\boldsymbol{\mu}}_k)(\mathbf{x}_i - \hat{\boldsymbol{\mu}}_k)^\top$$

where:

• $\hat{\boldsymbol{\Sigma}}_k$: the estimated covariance matrix of class $k$
• $N_k$: the number of observations in class $k$
• $\hat{\boldsymbol{\mu}}_k$: the empirical mean of class $k$
• $\mathbf{x}_i$: the $i$-th feature vector, associated with label $y_i$

The $N_k-1$ denominator gives the unbiased class covariance estimator for $N_k>1$. When the centered class scatter is positive definite, the Gaussian covariance maximum-likelihood estimator divides that same scatter by $N_k$. If the scatter is singular, no finite positive-definite covariance maximizes this unconstrained Gaussian density likelihood. The relative difference between the two normalization factors vanishes as $N_k$ grows. For LDA, the unbiased pooled covariance estimator is

$$\hat{\boldsymbol{\Sigma}} = \frac{1}{N - K} \sum_{k=1}^K \sum_{i: y_i = k} (\mathbf{x}_i - \hat{\boldsymbol{\mu}}_k)(\mathbf{x}_i - \hat{\boldsymbol{\mu}}_k)^\top,$$

where:

• $\hat{\boldsymbol{\Sigma}}$: the pooled estimate of the shared covariance matrix used by LDA
• $N = \sum_{j=1}^K N_j$: the total number of observations across all classes
• $K$: the number of classes
• $N_k$: the number of observations in class $k$
• $\hat{\boldsymbol{\mu}}_k$: the empirical mean of class $k$
• $\mathbf{x}_i$: the $i$-th feature vector, associated with label $y_i$

This pools the within-class scatter and divides by the residual degrees of freedom $N-K$. When the pooled scatter is positive definite, the shared-covariance Gaussian maximum-likelihood estimator instead divides that same scatter by $N$. If the pooled scatter is singular, the unconstrained positive-definite Gaussian covariance likelihood has no finite maximum. These are different finite-sample normalization conventions; the worked example and primitive implementations below consistently use the unbiased denominators.

Pooling is the estimation analogue of the shared-covariance assumption. Under independent Gaussian samples with a truly common covariance, the pooled estimate has $N-K$ residual degrees of freedom rather than one class's $N_k-1$. For a class with $N_k>1$, pooling has no larger covariance-entry sampling variance under that model; the reduction is strict when the other classes add residual degrees of freedom, so $N-K>N_k-1$. If class covariances differ, pooling instead estimates a compromise matrix.

For a centered class sample with $N_k>1$, $\operatorname{rank}(\hat{\boldsymbol{\Sigma}}_k)\le\min(d,N_k-1)$. Thus $1<N_k\le d$ forces singularity: for $d=100$ and $N_k=50$, there are at most 49 nonzero eigenvalues and the inverse does not exist. Having $N_k>d$ is necessary, but not sufficient, for nonsingularity. LDA's pooled scatter similarly has rank at most $\min(d,N-K)$. Even when these matrices are full rank, small samples and strong feature dependencies can make covariance estimation unreliable. Rank, conditioning, and estimation uncertainty are separate checks.

Shrinkage estimators regularize the covariance by blending the empirical estimate with a low-dimensional target:

$$\hat{\boldsymbol{\Sigma}}_k(\alpha) = (1 - \alpha) \hat{\boldsymbol{\Sigma}}_k + \alpha \mathbf{T},$$

where:

• $\hat{\boldsymbol{\Sigma}}_k(\alpha)$: the shrunk (regularized) covariance estimate for class $k$ at mixing level $\alpha$
• $\alpha$: the mixing (shrinkage) parameter in $[0, 1]$ controlling the bias-variance tradeoff; $\alpha = 0$ recovers the empirical covariance and $\alpha = 1$ yields the pure target
• $\hat{\boldsymbol{\Sigma}}_k$: the empirical covariance of class $k$
• $\mathbf{T}$: the shrinkage target, for example the identity matrix $\mathbf{I}$ or a scaled version $\sigma^2 \mathbf{I}$

At $\alpha=0$, the estimate is the empirical covariance; at $\alpha=1$, it is the target $\mathbf{T}$. If the empirical covariance is positive semidefinite and the target is positive definite, any $\alpha>0$ makes their convex blend positive definite. A shared identity or positive scaled-identity target is spherical. The Ledoit-Wolf (2004) estimator provides a data-driven shrinkage approach. It should not be conflated with manually choosing a regularization strength. In scikit-learn 1.8, QDA's default 'svd' solver uses reg_param; its 'eigen' solver instead supports shrinkage, including 'auto'. LDA supports shrinkage='auto' with 'lsqr' or 'eigen'. The solver and rank limitations are specified below.

Two further limitations matter. A sufficiently extreme observation can dominate an empirical covariance and move the fitted boundary. The empirical mean and covariance have unbounded influence even when the reference distribution is Gaussian. Moving a contaminating point farther away makes its influence grow without bound. Gaussianity therefore does not make these estimators robust to contamination. Robust covariance estimators address a different estimation problem and are outside this chapter's scope.

Discriminant analysis inherits several structural properties from its Gaussian foundation.

• Invertible affine transformations preserve unregularized LDA and QDA predictions in exact arithmetic when their estimates are transformed consistently. If $\mathbf{z}=A\mathbf{x}+\mathbf{b}$ with invertible $A$, then $\boldsymbol{\mu}'_k=A\boldsymbol{\mu}_k+\mathbf{b}$ and $\boldsymbol{\Sigma}'_k=A\boldsymbol{\Sigma}_kA^\top$. Mahalanobis distances are unchanged, while $|\boldsymbol{\Sigma}'_k|=|\det A|^2|\boldsymbol{\Sigma}_k|$. The full Gaussian discriminants therefore change by $-\log|\det A|$, common to every class. Reduced LDA scores, after dropping their shared quadratic term, can acquire a different common term depending on $\mathbf{x}$, but all score differences are still unchanged. This algebraic equivariance does not require Gaussian data; finite precision and a fixed regularization target can limit the invariance in practice.
• Posterior probabilities are given by the softmax of the discriminant scores (the logits):

$$p(y = k \mid \mathbf{x}) = \frac{\exp(\delta_k(\mathbf{x}))}{\sum_{j=1}^K \exp(\delta_j(\mathbf{x}))}.$$

This means LDA and QDA are not just classifiers but probabilistic models: the same machinery that predicts labels can also predict class probabilities. In practice, these probabilities can be miscalibrated if the Gaussian assumption fails, which motivates the calibration methods discussed in later parts of the book. In the implementation below, we compute the class posteriors directly from the discriminant scores via this softmax.

• Under the shared-covariance Gaussian model, the binary posterior is logistic in $\mathbf{x}$. In the earlier sign convention, $\delta_0-\delta_1=\mathbf{w}^\top\mathbf{x}+b$, so $p(y=1\mid\mathbf{x})=\sigma(-\mathbf{w}^\top\mathbf{x}-b)$, where $\sigma(t)=1/(1+\exp(-t))$. Equivalently, its class-1 logit has coefficient $\boldsymbol{\Sigma}^{-1}(\boldsymbol{\mu}_1-\boldsymbol{\mu}_0)$. Logistic regression fits this conditional linear-logit family directly. LDA obtains coefficients by fitting a joint Gaussian model; this imposes assumptions about feature distributions, rather than a smaller family of linear posterior logits.
• QDA does not generally satisfy a linear-logistic relationship. Its log-odds is quadratic in $\mathbf{x}$, so $p(y = 1 \mid \mathbf{x}) = \sigma(f(\mathbf{x}))$ holds with $f(\mathbf{x}) = \log(\pi_1/\pi_0) - \frac{1}{2}\log(|\boldsymbol{\Sigma}_1|/|\boldsymbol{\Sigma}_0|) - \frac{1}{2}(\mathbf{x}-\boldsymbol{\mu}_1)^\top\boldsymbol{\Sigma}_1^{-1}(\mathbf{x}-\boldsymbol{\mu}_1) + \frac{1}{2}(\mathbf{x}-\boldsymbol{\mu}_0)^\top\boldsymbol{\Sigma}_0^{-1}(\mathbf{x}-\boldsymbol{\mu}_0)$, a quadratic function of $\mathbf{x}$. The quadratic term enters before the logistic transformation, so the gradient of the log-odds generally varies affinely across feature space.

We now illustrate the geometric differences on a synthetic two-class problem with deliberately different covariances per class. Both classes are drawn from two-dimensional Gaussians. Class 1 has larger marginal variances and a stronger positive feature correlation than class 0. This violates LDA's shared-covariance assumption. QDA can represent the population Gaussian Bayes boundary here, although fitting its extra covariance parameters does not guarantee higher training or test accuracy.

The following cells prepare the data, fit both reference classifiers, and compute the decision grid and class probabilities before plotting. For binary classification, decision_function returns the one-dimensional log-posterior ratio $\log p(y=1\mid\mathbf{x})-\log p(y=0\mid\mathbf{x})$. The common normalization cancels exactly in this difference. We place a zero column beside that ratio to make two class scores with the same difference, then apply softmax across the two columns to recover the posterior probabilities.

The first figure shows LDA's straight boundary, whose normal uses the pooled estimate $\hat{\boldsymbol{\Sigma}}$ and the class means, $\hat{\boldsymbol{\Sigma}}^{-1}(\hat{\boldsymbol{\mu}}_0-\hat{\boldsymbol{\mu}}_1)$. Its shared covariance is estimated from the within-class scatter of both classes. It is a compromise, rather than class 0's covariance imposed on class 1. LDA can represent correlated and tilted Gaussian clouds, but it forces both clouds to have the same estimated shape and orientation. The plot compares this shared-shape restriction with QDA's class-specific fitted shapes.

The second figure shows QDA. The boundary is now curved. QDA allows class 1 to "reach out" along its long axis, so points that fall inside the elongated direction of the class-1 cloud are more likely to be attributed to class 1. QDA's distance penalty is smaller along a class's broad covariance directions than along its narrow directions. The log-determinant penalty depends on the full covariance volume and stays constant as the observation moves. The price of this flexibility is that QDA estimates $K = 2$ full $2 \times 2$ covariances. In higher dimensions, the per-class covariance estimates are more sensitive to sample size. This tradeoff is the subject of the next section.

We work through a small one-dimensional example with concrete numbers so that every step is visible. The point of the example is not to show a realistic data analysis, but to make the discriminant scores themselves tangible: we can compute them by hand, compare the LDA and QDA decisions, and see exactly where the quadratic boundary comes from. Consider a two-class problem with the following samples:

$$\text{Class 0}:\ \{0, 1, 2, 3\}, \qquad \text{Class 1}:\ \{4, 6, 8, 10\}.$$

Both classes have four observations, so $\hat{\pi}_0 = \hat{\pi}_1 = 0.5$. Equal priors keep the comparison focused on the densities: any difference in the decision must come from the class means and variances, not from one class being more common than the other.

Step 1: compute means and variances.

For class 0, the mean is $\hat{\mu}_0 = (0 + 1 + 2 + 3)/4 = 1.5$, and the unbiased variance is

$$\hat{\sigma}_0^2 = \frac{1}{3}\sum_{i}\bigl(x_i - 1.5\bigr)^2 = \frac{(-1.5)^2 + (-0.5)^2 + (0.5)^2 + (1.5)^2}{3} = \frac{5}{3} \approx 1.667.$$

For class 1, the mean is $\hat{\mu}_1 = (4 + 6 + 8 + 10)/4 = 7$, and the unbiased variance is

$$\hat{\sigma}_1^2 = \frac{(-3)^2 + (-1)^2 + (1)^2 + (3)^2}{3} = \frac{20}{3} \approx 6.667.$$

The empirical variance of class 1 is four times that of class 0, by construction of this small example. These sample spreads make the two fitted covariance conventions differ visibly. Four observations per class do not establish the covariances of an underlying population, so we treat this as a worked calculation rather than evidence that a real population violates LDA's assumption.

Step 2: pooled variance for LDA.

The pooled variance divides the total within-class squared deviations by the degrees of freedom $N - K = 8 - 2 = 6$. The numerator is the sum of the within-class squared deviations: class 0 contributes $5$ and class 1 contributes $20$, for a total of $5 + 20 = 25$. Thus

$$\hat{\sigma}_{\text{pooled}}^2 = \frac{5 + 20}{6} = \frac{25}{6} \approx 4.167.$$

This pooled value is a compromise between the tight class 0 variance and the broad class 1 variance. It is the single shape that LDA will use for both classes.

Step 3: LDA discriminant scores.

Using the exact pooled estimate $\hat{\sigma}^2=25/6\approx4.167$, the one-dimensional linear discriminant is

$$\delta_k(x) = \frac{\hat{\mu}_k}{\hat{\sigma}^2}\, x - \frac{\hat{\mu}_k^2}{2 \hat{\sigma}^2} + \log \hat{\pi}_k$$

where:

• $\hat{\mu}_k$: the estimated mean of class $k$
• $\hat{\sigma}^2 = 25/6 \approx 4.167$: the pooled variance shared by both classes
• $\hat{\pi}_k = 0.5$: the estimated prior probability of class $k$
• $x$: the (scalar) feature value being classified

Working through the two classes:

$$\delta_0(x) = \frac{1.5}{25/6}x - \frac{2.25}{2 \cdot (25/6)} + \log(0.5) \approx 0.36 x - 0.27 - 0.693 = 0.36 x - 0.963,$$ $$\delta_1(x) = \frac{7}{25/6}x - \frac{49}{2 \cdot (25/6)} + \log(0.5) \approx 1.68 x - 5.88 - 0.693 = 1.68 x - 6.573.$$

At a test point $x = 3$:

$$\delta_0(3) \approx 1.08 - 0.963 = 0.117, \qquad \delta_1(3) \approx 5.04 - 6.573 = -1.533.$$

The larger score is $\delta_0$, so LDA predicts class 0. At a test point $x = 5$:

$$\delta_0(5) \approx 1.80 - 0.963 = 0.837, \qquad \delta_1(5) \approx 8.40 - 6.573 = 1.827.$$

Here LDA predicts class 1. The boundary between the two classes is where $\delta_0(x) = \delta_1(x)$, i.e., $0.36x - 0.963 = 1.68x - 6.573$. Solving: subtract $0.36x$ from both sides to get $-0.963 = 1.32x - 6.573$; add $6.573$ to both sides to get $5.61 = 1.32x$; divide by $1.32$ to get $x \approx 4.25$. As expected, this is the midpoint of the two class means, $\tfrac{1}{2}(1.5 + 7) = 4.25$. The equal-prior assumption removes any prior-driven shift, and the shared variance makes the boundary exactly halfway between the centers.

Step 4: QDA discriminant scores.

QDA keeps the class-specific variances. The discriminant is

$$\delta_k^{\text{QDA}}(x) = -\frac{(x - \hat{\mu}_k)^2}{2 \hat{\sigma}_k^2} - \tfrac{1}{2} \log \hat{\sigma}_k^2 + \log \hat{\pi}_k.$$

For class 0:

$$\delta_0^{\text{QDA}}(x) = -\frac{(x-1.5)^2}{2(5/3)}-\tfrac{1}{2}\log(5/3)+\log(0.5) \approx -0.3(x-1.5)^2-0.949.$$

where:

• $\hat{\mu}_0 = 1.5$: the estimated mean of class 0
• $\hat{\sigma}_0^2 = 5/3 \approx 1.667$: the unbiased estimated variance of class 0
• $\hat{\pi}_0 = 0.5$: the estimated prior probability of class 0
• $x$: the test point at which the class-0 QDA score is evaluated

For class 1:

$$\delta_1^{\text{QDA}}(x) = -\frac{(x-7)^2}{2(20/3)}-\tfrac{1}{2}\log(20/3)+\log(0.5) \approx -0.075(x-7)^2-1.642.$$

At $x = 3$:

$$\begin{aligned} \delta_0^{\text{QDA}}(3) &= -0.675-\tfrac{1}{2}\log(5/3)+\log(0.5) \approx -1.624,\\ \delta_1^{\text{QDA}}(3) &= -1.200-\tfrac{1}{2}\log(20/3)+\log(0.5) \approx -2.842. \end{aligned}$$

The quadratic contributions are exactly $-(3-1.5)^2/[2(5/3)]=-0.675$ for class 0 and $-(3-7)^2/[2(20/3)]=-1.200$ for class 1. The combined log-variance and log-prior constants are approximately $-0.948560$ and $-1.641707$, respectively; retaining them until the final rounding gives the scores above.

In each score, the quadratic contribution is minus one-half of the squared standardized distance. The remaining terms are the log-variance penalty and log prior. Both classes are evaluated at the same test point $x=3$.

QDA also predicts class 0 at $x=3$ because $-1.624>-2.842$. At $x=4.25$, where the exact LDA scores tie:

$$\delta_0^{\text{QDA}}(4.25) \approx -3.217, \qquad \delta_1^{\text{QDA}}(4.25) \approx -2.209.$$

QDA favors class 1 at $x=4.25$, while LDA is indifferent there; a claim about different hard labels at that exact point would require an LDA tie-breaking convention. This point is 2.75 raw units from both means. Relative to the class standard deviations, that is about 2.13 standard deviations for class 0 and 1.07 for class 1. The quadratic contributions are $-2.26875$ and $-0.5671875$. Class 1 pays the larger log-variance penalty, approximately $-0.949$ rather than $-0.255$, but its less negative quadratic contribution more than offsets that penalty. Thus the broader class has the higher QDA score at the midpoint.

Step 5: QDA decision boundary.

Solving $\delta_0^{\text{QDA}}(x)=\delta_1^{\text{QDA}}(x)$ gives a quadratic equation. We keep the exact variances and logarithms until the numerical evaluation:

$$-\frac{(x-1.5)^2}{10/3}-\frac{1}{2}\log(5/3)+\log(0.5) = -\frac{(x-7)^2}{40/3}-\frac{1}{2}\log(20/3)+\log(0.5).$$

The equal log-prior terms cancel. Multiplying by $40/3$ and collecting the log-variance difference gives

$$\begin{aligned} 4(x-1.5)^2-(x-7)^2 &=\frac{20}{3}\log\!\left(\frac{20/3}{5/3}\right)\\ &=\frac{20}{3}\log 4\\ &=\frac{40}{3}\log 2. \end{aligned}$$

The left-hand side expands to $4(x^2-3x+2.25)-(x^2-14x+49)=4x^2-12x+9-x^2+14x-49=3x^2+2x-40$. Define the boundary polynomial

$$P(x)=3x^2+2x-40-\frac{40}{3}\log 2.$$

Here $P(x)$ is the exact quadratic whose zeros are the QDA ties, and $\log$ is the natural logarithm. Since $\frac{40}{3}\log 2\approx9.241962$, the constant term is approximately $-49.241962$. The quadratic formula gives

$$\begin{aligned} x &=\frac{-2\pm\sqrt{\,4+12\left(40+\frac{40}{3}\log 2\right)\,}}{6}\\ &\approx -4.398441 \quad\text{or}\quad 3.731774. \end{aligned}$$

To determine the winning regions, subtract the exact class-0 score from the class-1 score:

$$\frac{40}{3}\left(\delta_1^{\text{QDA}}(x)-\delta_0^{\text{QDA}}(x)\right)=P(x).$$

The multiplier is positive, so the score difference has the same sign as $P(x)$. Its positive leading coefficient makes $P(x)$ negative between the roots and positive outside them. Class 0 therefore wins between the two roots, class 1 wins outside, and the endpoints are ties. The wider class wins in both far tails because its quadratic Mahalanobis contribution grows more slowly, while the class-0 penalty along its own metric dominates far from both centers. This two-tail allocation cannot be expressed by a single nondegenerate LDA threshold.

The worked example produces more than isolated labels: it produces two score curves whose crossings define the decision regions. The figure below evaluates the fitted LDA and QDA scores on a grid whose interval spans both QDA roots and the LDA midpoint. Compare the two classes within one method: each method may discard a different term common to its classes, so score heights across methods are not on a shared probability scale.

We now implement both classifiers from first principles using NumPy, and then compare the results with scikit-learn's reference implementations. The fitting and scoring functions are defined first, then exercised on synthetic data. The from-scratch implementation is not intended to replace the library version; it is intended to make the mechanism auditable. In particular, it shows where the covariance inverse, the log-determinant, and the prior enter the score, and why numerical implementations avoid forming explicit matrix inverses whenever possible.

We compare two noisy two-dimensional Gaussian populations: the original unequal-covariance case and a shared-covariance case. Each has 80 training observations and 400 held-out observations per class. Separate random generators preserve the original unequal-covariance sample when we add the shared case.

The fitting routines accept a finite two-dimensional feature array X of shape $(N,d)$ and a one-dimensional label array y of length $N$. Labels may be noncontiguous numbers or strings; the sorted classes array maps score columns back to those labels. Both routines require at least two classes and finite, positive-definite estimated covariances. Finite feature inputs can still overflow during fitting; both routines reject nonfinite fitted state with FloatingPointError. The sample-count conditions below are necessary, not sufficient, for positive definiteness.

LDA divides pooled within-class scatter by $N-K$. It solves $\hat{\boldsymbol{\Sigma}}W^\top=M^\top$ once, where the rows of $M$ are the $K$ class means. The rows of $W$ are the weights $\mathbf{w}_k^\top$ from the derivation, and the length-$K$ intercept vector contains $-\frac12\hat{\boldsymbol{\mu}}_k^\top\mathbf{w}_k+\log\hat{\pi}_k$. Prediction evaluates X @ coef.T + intercept, an $(n,K)$ score matrix for $n$ new observations, with no covariance solve during prediction. After a finiteness check, a Cholesky step checks numerical positive definiteness before that solve.

QDA estimates one unbiased covariance per class and stores its Cholesky factor $L_k$, satisfying $\hat{\boldsymbol{\Sigma}}_k=L_kL_k^\top$, plus its log-determinant. To evaluate the Mahalanobis distance, the scoring routine solves $L_k Z=(X-\hat{\boldsymbol{\mu}}_k)^\top$ by forward substitution. Here $Z$ has shape $(d,n)$; summing its squared entries down each column yields the $n$ squared distances. This triangular solve avoids an explicit inverse and reuses the factors fitted once per class.