Book contents
Contents
raw Math
RAW Book Machine Learning Classification

Introduction to Linear Discriminant Analysis

Robert Eisele

Principal component analysis finds the axis along which a cloud of points spreads out the most, using nothing but the coordinates of the points themselves. If those points also carry class labels and the actual goal is to tell the classes apart, chasing raw spread can be actively counterproductive: the direction of greatest overall variance and the direction that best separates two groups are, in general, two completely different lines. Linear discriminant analysis (LDA) replaces "maximize spread" with "maximize the distance between class means relative to how much each class spreads out around its own mean," a criterion first proposed by R. A. Fisher (1936) for two groups and later generalized to any number of classes. The two ideas are easy to confuse because both produce a linear projection of the data, but they answer different questions: PCA asks where the data varies the most, while LDA asks where the classes differ the most.

variance axis separating axis

Projecting every point onto the dashed axis of greatest variance mixes the two classes together completely, since that axis runs along the direction in which both clusters are elongated. Projecting onto the solid axis instead — a direction PCA has no way of preferring, because it carries almost none of the total variance — separates the two classes perfectly. Section "A Worked Numerical Example" below reproduces this exact picture with the underlying numbers.

Between-Class and Within-Class Scatter

Suppose the data falls into \(C\) classes, with class \(c\) containing \(N_c\) points and having mean \(\boldsymbol{\mu}_c\), and let \(\boldsymbol{\mu}\) be the mean of all \(N=\sum_c N_c\) points together. Two scatter matrices summarize, respectively, how spread out each class is around its own mean and how far apart the class means are from the overall mean:

\[\mathbf{S}_W=\sum_{c=1}^{C}\sum_{i\in\text{class }c}(\mathbf{x}_i-\boldsymbol{\mu}_c)(\mathbf{x}_i-\boldsymbol{\mu}_c)^T, \qquad \mathbf{S}_B=\sum_{c=1}^{C}N_c(\boldsymbol{\mu}_c-\boldsymbol{\mu})(\boldsymbol{\mu}_c-\boldsymbol{\mu})^T.\]

\(\mathbf{S}_W\) (within-class scatter) is exactly the sum of each class's own covariance-like scatter matrix; a small \(\mathbf{S}_W\) means the classes are each tightly clustered around their own center. \(\mathbf{S}_B\) (between-class scatter) is built the same way but treats the \(C\) class means themselves as the data, weighted by how many points each class contributes; a large \(\mathbf{S}_B\) means the class centers are far apart. Adding a term \((\mathbf{x}_i-\boldsymbol{\mu}_c)+(\boldsymbol{\mu}_c-\boldsymbol{\mu})=\mathbf{x}_i-\boldsymbol{\mu}\) and expanding the total scatter matrix \(\mathbf{S}_T=\sum_i(\mathbf{x}_i-\boldsymbol{\mu})(\mathbf{x}_i-\boldsymbol{\mu})^T\) shows the cross terms vanish, because \(\sum_{i\in c}(\mathbf{x}_i-\boldsymbol{\mu}_c)=\mathbf{0}\) by definition of \(\boldsymbol{\mu}_c\), leaving the clean decomposition

\[\mathbf{S}_T=\mathbf{S}_W+\mathbf{S}_B,\]

the same total-equals-within-plus-between split used in one-way analysis of variance, just written for vectors instead of scalars.

Fisher's Criterion

A one-dimensional projection is fixed by choosing a direction \(\mathbf{w}\) and mapping each point to the scalar \(y=\mathbf{w}^T\mathbf{x}\). Projected onto \(\mathbf{w}\), the between-class scatter becomes \(\mathbf{w}^T\mathbf{S}_B\mathbf{w}\) and the within-class scatter becomes \(\mathbf{w}^T\mathbf{S}_W\mathbf{w}\), so Fisher's criterion asks for the direction that makes the former as large as possible relative to the latter:

\[J(\mathbf{w})=\frac{\mathbf{w}^T\mathbf{S}_B\mathbf{w}}{\mathbf{w}^T\mathbf{S}_W\mathbf{w}}.\]

This is a generalized Rayleigh quotient: it behaves exactly like the ordinary Rayleigh quotient \(\mathbf{w}^T\mathbf{S}_B\mathbf{w}/\mathbf{w}^T\mathbf{w}\) used to characterize eigenvectors, except the plain inner product \(\mathbf{w}^T\mathbf{w}\) in the denominator is replaced by the quadratic form \(\mathbf{w}^T\mathbf{S}_W\mathbf{w}\). \(J\) is invariant under rescaling \(\mathbf{w}\to\alpha\mathbf{w}\) for any \(\alpha\neq0\), since both numerator and denominator scale by \(\alpha^2\), so only the direction of \(\mathbf{w}\) matters, not its length. Differentiating with the quotient rule and setting the result to zero,

\[\frac{\partial J}{\partial\mathbf{w}}=\frac{2\mathbf{S}_B\mathbf{w}\,(\mathbf{w}^T\mathbf{S}_W\mathbf{w})-2\mathbf{S}_W\mathbf{w}\,(\mathbf{w}^T\mathbf{S}_B\mathbf{w})}{(\mathbf{w}^T\mathbf{S}_W\mathbf{w})^2}=\mathbf{0},\]

and dividing through by \(\mathbf{w}^T\mathbf{S}_W\mathbf{w}\) shows every stationary point satisfies

\[\boxed{\mathbf{S}_B\mathbf{w}=\lambda\,\mathbf{S}_W\mathbf{w}}\qquad\text{with }\lambda=J(\mathbf{w}),\]

a generalized eigenvalue problem. When \(\mathbf{S}_W\) is invertible this is equivalent to the ordinary eigenvalue problem \(\mathbf{S}_W^{-1}\mathbf{S}_B\mathbf{w}=\lambda\mathbf{w}\), and since \(J\) equals \(\lambda\) at any stationary direction, the maximum of Fisher's criterion is exactly the largest eigenvalue of \(\mathbf{S}_W^{-1}\mathbf{S}_B\), attained at its corresponding eigenvector.

The Two-Class Closed Form

With exactly two classes, \(\mathbf{S}_B\) has a special structure that avoids solving an eigenvalue problem at all. Writing \(N_1,N_2\) for the class sizes and \(\boldsymbol{\mu}_1,\boldsymbol{\mu}_2\) for the class means,

\[\mathbf{S}_B=\frac{N_1N_2}{N}(\boldsymbol{\mu}_1-\boldsymbol{\mu}_2)(\boldsymbol{\mu}_1-\boldsymbol{\mu}_2)^T,\]

which has rank one: for any vector \(\mathbf{w}\), the product \(\mathbf{S}_B\mathbf{w}\) is the vector \(\boldsymbol{\mu}_1-\boldsymbol{\mu}_2\) rescaled by the scalar \(\frac{N_1N_2}{N}(\boldsymbol{\mu}_1-\boldsymbol{\mu}_2)^T\mathbf{w}\), so \(\mathbf{S}_B\mathbf{w}\) always points along \(\boldsymbol{\mu}_1-\boldsymbol{\mu}_2\) no matter what \(\mathbf{w}\) is. Substituting this into \(\mathbf{S}_B\mathbf{w}=\lambda\mathbf{S}_W\mathbf{w}\) and absorbing every scalar factor into \(\lambda\) shows that the solution direction must satisfy \(\mathbf{S}_W\mathbf{w}\propto\boldsymbol{\mu}_1-\boldsymbol{\mu}_2\), giving the closed form Fisher originally derived:

\[\boxed{\mathbf{w}\propto\mathbf{S}_W^{-1}(\boldsymbol{\mu}_1-\boldsymbol{\mu}_2)}.\]

No eigenvalue problem needs to be solved in the two-class case — a single linear solve against \(\mathbf{S}_W\) is enough, which is also why this special case was the one Fisher (1936) originally worked out, years before the general multi-class version became standard.

A Worked Numerical Example

The two clusters plotted above consist of seven points each, built from a shared direction \(\mathbf{u}=(1,1)/\sqrt2\) (the direction both clusters are elongated along) and a perpendicular direction \(\mathbf{v}=(1,-1)/\sqrt2\) (the direction the two cluster centers are offset along):

\[ \begin{array}{r|rrrrrrr} \text{class A} & (-5.36,-5.71) & (-4.53,-4.18) & (-3.00,-3.36) & (-2.18,-1.82) & (-0.64,-1.00) & (0.18,0.53) & (1.71,1.36) \\ \text{class B} & (-3.52,-7.55) & (-2.69,-6.02) & (-1.17,-5.20) & (-0.34,-3.66) & (1.20,-2.83) & (2.02,-1.31) & (3.55,-0.48) \end{array} \]

The class means come out to \(\boldsymbol{\mu}_A=(-1.97,-2.03)\) and \(\boldsymbol{\mu}_B=(-0.14,-3.86)\). Computing the scatter matrices from these fourteen points directly from their definitions gives

\[ \mathbf{S}_W=\begin{pmatrix}78.18 & 77.33\\77.33 & 78.18\end{pmatrix},\qquad \mathbf{S}_B=\begin{pmatrix}11.83 & -11.83\\-11.83 & 11.83\end{pmatrix}, \]

and \(\mathbf{S}_T=\mathbf{S}_W+\mathbf{S}_B\) confirms the total-scatter identity exactly. Solving \(\mathbf{w}\propto\mathbf{S}_W^{-1}(\boldsymbol{\mu}_A-\boldsymbol{\mu}_B)\) gives, once normalized to unit length, \(\mathbf{w}_{\text{LDA}}=(-0.707,0.707)\) — exactly the direction \(\mathbf{v}\) the cluster centers were offset along. The largest eigenvalue of \(\mathbf{S}_W^{-1}\mathbf{S}_B\) is \(27.60\), matching \(J(\mathbf{w}_{\text{LDA}})\) exactly, while the second eigenvalue is \(0\), consistent with \(\mathbf{S}_B\) having rank one for two classes. The direction of greatest total variance, by contrast — the leading eigenvector of \(\mathbf{S}_T\), which is what PCA would return — comes out to \((0.707,0.707)\), exactly \(\mathbf{u}\), and evaluating Fisher's criterion there gives \(J(\mathbf{w}_{\text{PCA}})\approx0\): projecting onto it, the two classes overlap almost perfectly, since \(\mathbf{u}\) is precisely the direction the between-class offset has zero component along. Both endpoints of the maximization problem are visible in a single dataset: the axis of maximum variance and the axis of maximum class separation turn out to be exactly perpendicular to each other here.

More Than Two Classes

With \(C>2\) classes, \(\mathbf{S}_B\) is built from \(C\) rank-one terms \(N_c(\boldsymbol{\mu}_c-\boldsymbol{\mu})(\boldsymbol{\mu}_c-\boldsymbol{\mu})^T\), so its rank is at most \(C\) — but one degree of freedom is used up because the \(N_c\)-weighted class-mean deviations always sum to the zero vector, \(\sum_c N_c(\boldsymbol{\mu}_c-\boldsymbol{\mu})=\sum_c N_c\boldsymbol{\mu}_c-N\boldsymbol{\mu}=\mathbf{0}\) by definition of \(\boldsymbol{\mu}\). That leaves at most \(C-1\) linearly independent deviation vectors, so \(\operatorname{rank}(\mathbf{S}_B)\le C-1\): no more than \(C-1\) directions can ever have a nonzero Fisher criterion, however many dimensions the original data lives in. Rather than a single direction, the general problem asks for a whole projection matrix \(\mathbf{W}\) whose columns \(\mathbf{w}_1,\dots,\mathbf{w}_k\) jointly maximize the matrix generalization of Fisher's criterion, \(\operatorname{trace}\!\big((\mathbf{W}^T\mathbf{S}_W\mathbf{W})^{-1}\mathbf{W}^T\mathbf{S}_B\mathbf{W}\big)\), and the solution is exactly the top \(k\le C-1\) eigenvectors of \(\mathbf{S}_W^{-1}\mathbf{S}_B\), ordered by eigenvalue — the same generalized eigenvalue problem as before, just keeping more than one solution.

\(\mathbf{S}_W^{-1}\mathbf{S}_B\) is not symmetric in general, even though both factors are, which is numerically unpleasant: symmetric eigenvalue solvers are faster and more stable than general ones. The Cholesky factorization \(\mathbf{S}_W=\mathbf{L}\mathbf{L}^T\) fixes this. Substituting \(\mathbf{v}=\mathbf{L}^T\mathbf{w}\), so that \(\mathbf{w}=\mathbf{L}^{-T}\mathbf{v}\), turns the generalized problem into an ordinary symmetric one:

\[\mathbf{S}_B\mathbf{L}^{-T}\mathbf{v}=\lambda\mathbf{L}\mathbf{L}^T\mathbf{L}^{-T}\mathbf{v}=\lambda\mathbf{L}\mathbf{v} \quad\Longrightarrow\quad \big(\mathbf{L}^{-1}\mathbf{S}_B\mathbf{L}^{-T}\big)\mathbf{v}=\lambda\mathbf{v},\]

where \(\mathbf{M}=\mathbf{L}^{-1}\mathbf{S}_B\mathbf{L}^{-T}\) is symmetric because \(\mathbf{S}_B\) is. Solving the ordinary symmetric eigenproblem \(\mathbf{M}\mathbf{v}=\lambda\mathbf{v}\) and mapping each eigenvector back with \(\mathbf{w}=\mathbf{L}^{-T}\mathbf{v}\) recovers exactly the same directions and eigenvalues as solving \(\mathbf{S}_W^{-1}\mathbf{S}_B\mathbf{w}=\lambda\mathbf{w}\) directly. Running both routes on the worked example above confirms this: the symmetric route returns \(\mathbf{v}\) mapping back to \(\mathbf{w}=(-0.707,0.707)\) with the same top eigenvalue \(27.60\), matching the direct solution exactly.

Relation to Least-Squares Regression

The two-class closed form has a second derivation that starts from ordinary least-squares regression instead of a scatter-matrix ratio, which is a useful sanity check since the two derivations look nothing alike at first glance. Assign every point in class 1 the numeric target \(t_i=N/N_1\) and every point in class 2 the target \(t_i=-N/N_2\) — targets chosen so that, weighted by class size, they average to zero over the whole dataset. Fitting an ordinary linear regression \(t\approx\mathbf{w}^T\mathbf{x}+b\) by least squares and reading off the coefficient vector \(\mathbf{w}\) (ignoring the intercept \(b\)) reproduces, up to an overall scale factor, exactly \(\mathbf{S}_W^{-1}(\boldsymbol{\mu}_1-\boldsymbol{\mu}_2)\): running this regression on the fourteen points from the worked example above gives a coefficient direction that agrees with \(\mathbf{w}_{\text{LDA}}\) to within floating-point rounding. The reason is structural rather than coincidental: the least-squares normal equations for this design reduce, after expanding the sums, to exactly the same linear system that defines the Fisher direction, with the arbitrary-looking target values chosen precisely so the cross terms between the two classes cancel out. This equivalence only holds for two classes with this specific target encoding; the general multi-class problem does not reduce to plain regression in the same way.

PCA versus LDA: When More Assumptions Beat More Data

Fisher's criterion is provably optimal once \(\mathbf{S}_W\) and \(\mathbf{S}_B\) are known exactly, but in practice both are estimated from a finite sample, and \(\mathbf{S}_W^{-1}\) can amplify estimation noise badly when the sample is small relative to the number of dimensions. Martinez and Kak (2001) studied exactly this failure mode and showed that PCA can outperform LDA at classification precisely in the small-sample regime, despite LDA optimizing a criterion tailored to classification and PCA optimizing one that ignores the labels entirely.

LDA axis PCA axis

The eight points plotted here (solid circles for one class, open circles for the other) are drawn from two populations whose true means are well separated along a diagonal and whose true covariance is strongly elongated along that very same diagonal — the hardest configuration for a small sample, since a handful of points is not enough to reveal the true shape of that elongation. On these eight points alone, computing \(\mathbf{S}_W^{-1}(\boldsymbol{\mu}_1-\boldsymbol{\mu}_2)\) gives an axis that, projected on, separates these particular eight training points quite well by construction — that is, after all, exactly what the LDA criterion was built to optimize on whatever sample it is given. But tested against 20,000 fresh points drawn from the same two true populations, that same axis performs almost no better than chance, since it happens to run nearly perpendicular to the true separating direction: a small, noisy estimate of \(\mathbf{S}_W\) pointed the closed-form solution the wrong way. The direction of greatest total variance, which does not depend on any per-class scatter estimate at all, lands much closer to the true separating direction and keeps working well on the fresh 20,000-point sample. This is not an argument against LDA in general — with enough data per class relative to the number of dimensions, its estimate of \(\mathbf{S}_W\) becomes reliable and it reliably outperforms PCA at the classification task it is built for — it is a warning that the closed-form optimum of a criterion is only as trustworthy as the scatter estimates that go into it.

Applications

Fisher directions remain linear combinations of the original features, so classes that are only separable by a curved boundary defeat LDA the same way a curved manifold defeats PCA. Mika et al. (1999) extended the same scatter-ratio criterion into a reproducing-kernel feature space, mirroring how kernel PCA generalizes ordinary PCA, to recover a nonlinear version of the method for exactly those cases.

References

Related Stuff

raw / Linear Algebra
02 Aug
·
Robert Eisele

Introduction to Principal Component Analysis

This derives principal component analysis (PCA) from first principles, starting with a visual scatter-plot illustration of why correlated data has a natural long axis, then deriving the first principal component as the direction maximizing projected variance via both the Rayleigh quotient and a Lagrange-multiplier argument, showing both lead to the covariance matrix's dominant eigenvector. It extends this to the full eigendecomposition for all remaining orthogonal components, proves the equivalent reconstruction-error-minimization view via a Pythagorean argument and shows it yields the identical optimal subspace, explains precisely what 'percentage of variance explained' means as a share of the trace of the covariance matrix (with a worked five-component example and an explained-variance bar chart), connects PCA to the singular value decomposition and the Eckart-Young theorem for optimal low-rank reconstruction, contrasts the first principal component with an ordinary least-squares regression line via their different (perpendicular vs. vertical) error criteria, illustrated with a diagram overlaying both fitted lines and both residual directions on the same data, discusses assumptions and limitations (linearity, sensitivity to scale and outliers, decorrelation versus true independence) with kernel PCA and independent component analysis mentioned as extensions, and closes with applications in dimensionality reduction and visualization (illustrated with a schematic MNIST digit-clustering diagram), noise reduction and compression (including the older Hotelling and Karhunen-Loeve transform names for the same construction), and Eigenfaces (Sirovich and Kirby 1987; Turk and Pentland 1991; Lanitis, Taylor and Cootes 1995 for shape eigenfaces) with nearest-neighbor recognition in score space. It also notes the sign- and degeneracy-related non-uniqueness of the principal components, mirroring the same ambiguity already discussed for the singular value decomposition, defines loadings alongside scores, derives the dual-PCA Gram-matrix trick for computing components when there are far fewer observations than variables (cross-linked to a dedicated chapter on EigenGame, a game-theoretic reformulation for computing components at even larger scale), derives whitening (sphering) as a rescaling of the scores to unit variance, and cautions that variance-based component selection is only a proxy for downstream task performance.
raw / Linear Algebra
02 Aug
·
Robert Eisele

Introduction to Singular Value Decomposition

This derives the singular value decomposition (SVD) of a real matrix of any shape directly from the eigendecomposition of the symmetric positive semidefinite matrix A^T A, motivated visually by how any matrix maps the unit circle to an ellipse and by the general rotate-stretch-rotate geometric picture this gives for every matrix. It works through a full numerical 2x2 example, shows how the right and left singular vectors are eigenvectors of A^T A and A A^T and how this completes the four fundamental subspaces picture, illustrated with a diagram of the Strang diagram itself, distinguishes the full and compact SVD, answers precisely when the eigendecomposition of a symmetric matrix coincides with its SVD (positive semidefiniteness) and how singular values differ from eigenvalues in general (including for defective and rotation matrices), derives the uniqueness of the SVD up to sign flips and rotations within repeated singular values with a worked cautionary example of what goes wrong when both sides are eigendecomposed independently, explains the numerically robust way to derive one side from the other, proves the Eckart-Young theorem for optimal low-rank approximation, and lists its properties (including the determinant as a product of singular values) and applications in least squares (including the explicit pseudo-inverse sum formula, exact-solvability criterion, and affine solution space), homogeneous least squares, data compression, and principal component analysis.
raw / Linear Algebra
08 May
·
Robert Eisele

Introduction to Eigenvalues and Eigenvectors

This introduces eigenvalues and eigenvectors as invariant directions of linear transformations. It derives the characteristic polynomial (including the alternate det(lambda I - A) sign convention), works through a full numerical example of a symmetric matrix with a repeated eigenvalue to show the characteristic polynomial, row reduction, and eigenspace computation concretely, explains eigenspaces and multiplicity including a defective, non-diagonalisable counter-example and why real symmetric matrices can never be defective, proves the trace/determinant formulas and the linear independence of eigenvectors for distinct eigenvalues, states Cauchy's interlacing theorem for the eigenvalues of a principal submatrix and the Perron-Frobenius theorem for nonnegative matrices, studies orthogonal matrices (proving every eigenvalue has absolute value 1), rotation matrices with complex eigenvalues and the degrees-of-freedom argument behind them, reflection matrices with real +-1 eigenvalues and a diagram of the mirror line and an eigenvector's invariant direction, and three-dimensional rotations via Euler's rotation theorem and the trace-angle relation, presents the eigenvalue decomposition, extends it to arbitrary real matrix powers and the symmetric square root, derives Gram factors, whitening transforms, and their use in generating correlated random vectors, derives the Rayleigh quotient and its variational characterization of the extreme eigenvalues, shows how repeated application of a matrix is governed by its dominant eigenvalue via power iteration, illustrated with an interactive diagram of successive normalized directions spiraling onto the dominant eigenvector and with the Fibonacci sequence and Binet's formula, proves the Cayley-Hamilton theorem (with a common false-proof pitfall explained) and uses it to derive the determinant of a matrix directly from traces of its powers via Newton's identities, and lists typical applications in linear algebra and engineering, including spectral graph theory and PageRank.
All Math