Linear and Quadratic Discriminant Analysis › Implementation in NumPy and scikit-learn

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.

Both fits check the sample-count prerequisites before estimating covariances. These checks raise ValueError. Computed means and covariances must be finite before factorization, and factors, log-determinants, and LDA coefficients and intercepts must be finite before a model is returned; otherwise the routines raise FloatingPointError. A Cholesky step raises LinAlgError when it detects a failure of numerical positive definiteness. If factorization fails, check numerical range, feature scales, rank, and conditioning before choosing a remedy. Shrinkage toward a positive-definite target, as discussed earlier, can make a singular empirical estimate invertible. The 80 observations per class and two features in this example avoid the sample-count obstruction, but those counts alone do not guarantee positive definiteness. The finiteness and factorization checks test these numerical conditions for the actual estimates. Its quadratic form comes from a triangular-system solve, and its log-determinant comes from the factor's diagonal.

We now generate training and test data, fit both from-scratch classifiers and the scikit-learn equivalents, and report test accuracy and prediction agreement. The comparison has two purposes: it checks that our implementation recovers the same decision rule, and it separates the statistical question (which classifier fits this data better?) from the numerical question (do the implementations agree?).

Each row compares accuracy on 800 held-out observations and label agreement between the primitive and reference implementation. The default reference LDA uses its SVD solver; our LDA precomputes linear coefficients with a solve, while our QDA stores Cholesky factors and uses forward substitution. For these full-rank samples in scikit-learn 1.8.0, the reference LDA SVD path uses the same $N-K$ within-class normalization as our LDA. The default unregularized QDA SVD path uses the same $N_k-1$ class-covariance normalization as our QDA. Other LDA solvers can use a different covariance normalization, so changing solvers requires checking that convention again.

Agreement and accuracy answer different questions. Agreement measures whether the implementations return the same labels on these observations; it does not prove identical behavior for every input. Accuracy compares those labels with the generated test classes. The unequal-covariance population admits a Gaussian QDA Bayes rule, whereas LDA imposes a shared covariance. The shared-covariance population satisfies LDA's covariance assumption, while QDA estimates two covariances without a more appropriate population model. The displayed accuracies describe these two finite samples; neither covariance condition guarantees which fitted classifier will win on another held-out sample.

The main parameters that affect how the fitted LDA and QDA models behave are listed below. These settings belong to the scikit-learn estimators; the primitive fits implement only the unregularized covariance estimates and an optional priors array in sorted classes order. Both primitive fits store class means and covariance arrays. LDA also stores precomputed weights (coef) and intercepts; QDA stores Cholesky factors and log-determinants. These are different return dictionaries.

• priors: Class prior probabilities (default: None, which uses the observed class frequencies). Set them explicitly when the deployment population differs from the training sample, because the prior enters each class score. Changing prior ratios changes score differences and can move, create, or remove decision boundaries; identical class densities can still give a constant winning class.
• shrinkage (LDA, scikit-learn): Covariance regularization (default: None). A value in $[0,1]$ blends the pooled covariance toward a scaled identity target, while 'auto' applies the Ledoit-Wolf estimate. It is available only with the 'lsqr' or 'eigen' solver.
• solver (LDA, scikit-learn): Factorization method (default: 'svd'). With the default store_covariance=False, the 'svd' path avoids forming the covariance explicitly; it cannot be combined with shrinkage. Use 'lsqr' or 'eigen' when shrinkage is needed.
• reg_param (QDA, scikit-learn's default 'svd' solver): regularization strength in $[0,1]$ (default: 0.0). It blends the retained empirical covariance eigenvalues, the squared singular values divided by $N_k-1$, toward one. When the SVD contains a full feature-space basis, this corresponds to blending the class covariance toward identity. With $N_k<d$, the compact SVD lacks feature directions, so this parameter cannot restore full rank; the installed 1.8 implementation raises an error. In version 1.8, solver='eigen' with a suitable shrinkage value supplies a full-covariance alternative.
• store_covariance (LDA, scikit-learn): For the 'svd' solver, whether to compute and retain the weighted within-class covariance (default: False). The 'lsqr' and 'eigen' solvers always compute and store it. Enable it with 'svd' when covariance diagnostics such as eigenvalues or condition numbers are part of the analysis.