Gaussian Discriminant Analysis, machine learning lesson

A machine reports two continuous sensor readings, and an engineer wants the probability that it is running in operating mode A rather than mode B. Binary word features suited the naive Bayes example of the previous lesson, but continuous readings need a density for each mode, and the two readings may vary together. Gaussian discriminant analysis gives each class a multivariate Gaussian density and applies Bayes’ rule. This lesson shows why one shared covariance matrix leads to a straight boundary, why separate covariance matrices allow a curved one, and what the extra flexibility costs in data.

Key terms

Class-conditional density

Definition For continuous features, the class-conditional density $f_c$ of class $c$ is the probability density of the feature vector $\mathbf{X}$ given $Y = c$ , so that $P(\mathbf{X} \in R \mid Y = c) = \int_R f_c(\mathbf{x})\, d\mathbf{x}$ for every region $R$ of feature space. With class priors $\pi_c > 0$ for $c = 1, \dots, K$ , the posterior probability of class $c$ at $\mathbf{x}$ is

\begin{aligned} & P(Y = c \mid \mathbf{X} = \mathbf{x}) \\[0.4em] &= \frac{\pi_c f_c(\mathbf{x})}{\sum_{k=1}^{K} \pi_k f_k(\mathbf{x})} \end{aligned}

wherever the denominator is positive. A density value is not a probability. It can be larger than $1$ , and the probability of any single exact vector is $0$ .

In plain words Each class gets its own description of where its feature vectors tend to fall. Bayes’ rule compares how densely each class populates the neighborhood of $\mathbf{x}$ , weighted by how common the class is. Comparing densities works because the small neighborhood around $\mathbf{x}$ has the same size for every class, so its size cancels from the ratio.

Example Two dimensionless, preprocessed sensor readings describe each machine. Four constructed training machines in mode A read $(-1, 0)$ , $(1, 0)$ , $(0, -1)$ and $(0, 1)$ , and four in mode B read $(0, 0)$ , $(4, 0)$ , $(2, -1)$ and $(2, 1)$ . The Gaussian densities fitted in the cards below give the reading $(0, 0)$ the values $f_A = 1/\pi \approx 0.3183$ and $f_B = e^{-1}/(2\pi) \approx 0.0585$ . With equal priors,

\begin{aligned} P(B \mid \mathbf{x}) &= \frac{0.0585}{0.3183 + 0.0585} \\[0.4em] &\approx 0.1554 . \end{aligned}

Multivariate Gaussian distribution

Definition A random vector $\mathbf{X}$ in $\mathbb{R}^p$ has the multivariate Gaussian, or normal, distribution $N(\boldsymbol{\mu}, \Sigma)$ with mean vector $\boldsymbol{\mu}$ in $\mathbb{R}^p$ and a symmetric positive-definite $p \times p$ covariance matrix $\Sigma$ if it has the density

\begin{aligned} f(\mathbf{x}) &= \frac{\exp\bigl(-\tfrac{1}{2}\Delta^2(\mathbf{x})\bigr)}{(2\pi)^{p/2} (\det \Sigma)^{1/2}} , \\[0.4em] \Delta^2(\mathbf{x}) &= \mathbf{z}^{\top} \Sigma^{-1} \mathbf{z} , \end{aligned}

where $\mathbf{z} = \mathbf{x} - \boldsymbol{\mu}$ , for every $\mathbf{x}$ in $\mathbb{R}^p$ . Then $E[\mathbf{X}] = \boldsymbol{\mu}$ and the covariance matrix of $\mathbf{X}$ is $\Sigma$ . Its diagonal entries are the variances of the individual features, and its off-diagonal entries are the covariances, which say how deviations of two features move together. The number $\Delta^2(\mathbf{x})$ is the squared Mahalanobis distance from $\mathbf{x}$ to $\boldsymbol{\mu}$ . The density is constant on each set where $\Delta^2$ is constant, an ellipse for $p = 2$ and an ellipsoid in general, with axes along the eigenvectors of $\Sigma$ .

In plain words The mean is the center of the cloud of points, and the covariance matrix sets its spread and its tilt. Mahalanobis distance measures how far a point lies from the center in units of the spread in that direction, so a step along a direction of high variance counts for less than the same step along a direction of low variance. Probabilities come from integrating the density over a region. The density itself can exceed $1$ , and for the two-dimensional Gaussian with mean $\mathbf{0}$ and covariance $0.01\,I$ it is about $15.9$ at the center.

Example The Gaussian fitted to mode B has

\begin{gathered} \hat{\boldsymbol{\mu}}_B = \begin{bmatrix} 2 \\ 0 \end{bmatrix} , \\[0.4em] \hat{\Sigma}_B = \begin{bmatrix} 2 & 0 \\ 0 & 0.5 \end{bmatrix} . \end{gathered}

The points $(4, 0)$ and $(2, 2)$ both lie at Euclidean distance $2$ from the mean. Their squared Mahalanobis distances are $2^2/2 = 2$ and $2^2/0.5 = 8$ , so their densities are $e^{-1}/(2\pi) \approx 0.0585$ and $e^{-4}/(2\pi) \approx 0.0029$ , compared with $1/(2\pi) \approx 0.1592$ at the mean. The second point is much less typical of mode B, because the second sensor varies much less.

Positive-definite covariance matrix

Definition A real symmetric $p \times p$ matrix $\Sigma$ is positive definite if $\mathbf{v}^{\top} \Sigma \mathbf{v} > 0$ for every nonzero vector $\mathbf{v}$ in $\mathbb{R}^p$ , which holds exactly when all eigenvalues of $\Sigma$ are positive. Every covariance matrix is symmetric and positive semidefinite, because $\mathbf{v}^{\top} \Sigma \mathbf{v} = \operatorname{Var}(\mathbf{v}^{\top}\mathbf{X}) \ge 0$ . It is positive definite exactly when no nonzero linear combination $\mathbf{v}^{\top}\mathbf{X}$ of the features has variance $0$ . A positive-definite matrix is invertible with $\det \Sigma > 0$ , which the Gaussian density formula needs. A symmetric $2 \times 2$ matrix with entries $s_{11}$ , $s_{12}$ , $s_{12}$ , $s_{22}$ is positive definite exactly when $s_{11} > 0$ and $s_{11}s_{22} - s_{12}^2 > 0$ . The maximum likelihood covariance estimate from $n_c$ points in $\mathbb{R}^p$ has rank at most $n_c - 1$ , so it is singular whenever $n_c \le p$ .

In plain words Positive definite means the data spread out in every direction. If some direction has no spread at all, the cloud is flat, the density formula would divide by zero, and Mahalanobis distance in that direction has no meaning. Software that inverts such a matrix either fails or returns numbers dominated by rounding error.

Example The diagonal matrix with entries $1.25$ and $0.5$ has those two numbers as eigenvalues, so it is positive definite. If the two sensors were exact copies of each other, each with variance $1$ , their covariance matrix would be

\begin{bmatrix} 1 & 1 \\ 1 & 1 \end{bmatrix} ,

with eigenvalues $2$ and $0$ and determinant $0$ . The direction $\mathbf{v} = (1, -1)$ gives $\mathbf{v}^{\top}\Sigma\mathbf{v} = 0$ , since the difference of the two readings is always $0$ . Adding $0.1\,I$ gives eigenvalues $2.1$ and $0.1$ and a positive-definite matrix, but that regularization changes the model, and its size should be chosen on validation data. The symmetric matrix with rows $(1, 2)$ and $(2, 1)$ has eigenvalues $3$ and $-1$ and cannot be a covariance matrix at all.

Gaussian discriminant analysis

Definition Gaussian discriminant analysis is the generative classifier whose class-conditional densities are multivariate Gaussian, $f_c = N(\boldsymbol{\mu}_c, \Sigma_c)$ with each $\Sigma_c$ positive definite, combined with class priors $\pi_c > 0$ . Taking the logarithm of $\pi_c f_c(\mathbf{x})$ and dropping the term $-\tfrac{p}{2}\ln(2\pi)$ , which is the same for every class, gives the discriminant score

\begin{aligned} \delta_c(\mathbf{x}) &= \ln \pi_c \\[0.4em] &\quad - \tfrac{1}{2} \ln \det \Sigma_c \\[0.4em] &\quad - \tfrac{1}{2} \Delta_c^2(\mathbf{x}) , \end{aligned}

where $\Delta_c^2(\mathbf{x})$ is the squared Mahalanobis distance from $\mathbf{x}$ to $\boldsymbol{\mu}_c$ under $\Sigma_c$ . The classifier predicts the class with the largest $\delta_c(\mathbf{x})$ , the posterior probabilities are $e^{\delta_c(\mathbf{x})} \big/ \sum_k e^{\delta_k(\mathbf{x})}$ , and with two classes $P(B \mid \mathbf{x}) = \sigma\bigl(\delta_B(\mathbf{x}) - \delta_A(\mathbf{x})\bigr)$ with the sigmoid $\sigma$ . The maximum likelihood estimates are $\hat{\pi}_c = n_c/n$ , the class mean $\hat{\boldsymbol{\mu}}_c$ and

\hat{\Sigma}_c = \frac{1}{n_c} \sum_{i:\, y_i = c} \mathbf{d}_i \mathbf{d}_i^{\top} ,

where $\mathbf{d}_i = \mathbf{x}_i - \hat{\boldsymbol{\mu}}_c$ . Many libraries divide by $n_c - 1$ instead, which changes the numbers slightly.

In plain words The method fits an elliptical cloud to each class and asks which cloud makes a new point most plausible, after weighting by how common each class is. The score has three parts, the prior, a penalty for a widely spread class, and the Mahalanobis distance to the class center. A widely spread class has a lower density near its center but loses density more slowly as points move away.

Example Both modes have four training machines, so $\hat{\pi}_A = \hat{\pi}_B = 1/2$ . The means are $\hat{\boldsymbol{\mu}}_A = (0, 0)$ and $\hat{\boldsymbol{\mu}}_B = (2, 0)$ . The deviations from the mean in mode B are $(-2, 0)$ , $(2, 0)$ , $(0, -1)$ and $(0, 1)$ , whose outer products add up to the scatter matrix with diagonal entries $8$ and $2$ . The cross-products cancel in both modes, so

\begin{gathered} \hat{\Sigma}_A = \begin{bmatrix} 0.5 & 0 \\ 0 & 0.5 \end{bmatrix} , \\[0.4em] \hat{\Sigma}_B = \begin{bmatrix} 2 & 0 \\ 0 & 0.5 \end{bmatrix} . \end{gathered}

With the divisor $3$ instead of $4$ , the diagonals would be $2/3, 2/3$ and $8/3, 2/3$ .

Linear discriminant analysis

Definition Linear discriminant analysis, LDA in this lesson and not latent Dirichlet allocation, is Gaussian discriminant analysis with a shared covariance matrix. Every class has $f_c = N(\boldsymbol{\mu}_c, \Sigma)$ with the same positive-definite $\Sigma$ . The terms $-\tfrac{1}{2}\ln \det \Sigma$ and $-\tfrac{1}{2}\mathbf{x}^{\top}\Sigma^{-1}\mathbf{x}$ are then common to all classes and drop out, which leaves the linear discriminant

\begin{aligned} \delta_c(\mathbf{x}) &= \mathbf{x}^{\top} \Sigma^{-1} \boldsymbol{\mu}_c \\[0.4em] &\quad - \tfrac{1}{2} \boldsymbol{\mu}_c^{\top} \Sigma^{-1} \boldsymbol{\mu}_c \\[0.4em] &\quad + \ln \pi_c . \end{aligned}

For two classes $P(B \mid \mathbf{x}) = \sigma(\mathbf{w}^{\top}\mathbf{x} + w_0)$ with

\begin{aligned} \mathbf{w} &= \Sigma^{-1}(\boldsymbol{\mu}_B - \boldsymbol{\mu}_A) , \\[0.4em] w_0 &= \ln \tfrac{\pi_B}{\pi_A} - \tfrac{1}{2} \boldsymbol{\mu}_B^{\top} \Sigma^{-1} \boldsymbol{\mu}_B \\[0.4em] &\quad + \tfrac{1}{2} \boldsymbol{\mu}_A^{\top} \Sigma^{-1} \boldsymbol{\mu}_A , \end{aligned}

so when $\boldsymbol{\mu}_A \ne \boldsymbol{\mu}_B$ the boundary $\mathbf{w}^{\top}\mathbf{x} + w_0 = 0$ is a hyperplane, a line when $p = 2$ . The shared matrix is estimated by the pooled within-class covariance $\hat{\Sigma} = \sum_c (n_c/n)\, \hat{\Sigma}_c$ , the sum of all within-class scatter divided by $n$ , which is its maximum likelihood estimate. Dividing the scatter by $n - K$ instead gives the usual unbiased version. However many classes there are, the shared matrix has $p(p + 1)/2$ free entries.

In plain words LDA assumes every class has a cloud of the same shape and orientation, only centered in different places. The quadratic parts of the scores are then identical and cancel, so two classes are compared through a single linear combination of the features. Pooling all classes to estimate one matrix keeps the estimate stable with modest data, but it misses real differences in spread. The posterior has the same logistic form as logistic regression, yet its coefficients come from class means and a pooled covariance instead of from maximizing the conditional likelihood, so the two methods usually give different fits to the same data.

Example Pooling the two fitted matrices gives

\begin{aligned} \hat{\Sigma} &= \tfrac{4}{8}\hat{\Sigma}_A + \tfrac{4}{8}\hat{\Sigma}_B \\[0.4em] &= \begin{bmatrix} 1.25 & 0 \\ 0 & 0.5 \end{bmatrix} . \end{aligned}

Then $\mathbf{w} = (2/1.25, 0) = (1.6, 0)$ and $w_0 = -\tfrac{1}{2} \cdot 2^2/1.25 = -1.6$ , so $\delta_B - \delta_A = 1.6\,x_1 - 1.6$ . The boundary is the vertical line $x_1 = 1$ , halfway between the means because the priors are equal. The reading $(0, 0)$ gets $P(B \mid \mathbf{x}) = \sigma(-1.6) \approx 0.1680$ , and the reading $(3, 0)$ gets $\sigma(3.2) \approx 0.9608$ .

Quadratic discriminant analysis

Definition Quadratic discriminant analysis, QDA, is Gaussian discriminant analysis with a separate positive-definite covariance matrix $\Sigma_c$ for each class. Each discriminant keeps the term $-\tfrac{1}{2}\mathbf{x}^{\top}\Sigma_c^{-1}\mathbf{x}$ , which differs between classes whenever their covariance matrices differ. The difference $\delta_B(\mathbf{x}) - \delta_A(\mathbf{x})$ is then a quadratic function of $\mathbf{x}$ , and the boundary where it equals $0$ is a quadratic surface, in two dimensions a conic section, possibly degenerate, such as an ellipse, a parabola, a hyperbola or a pair of lines. Each $\hat{\Sigma}_c$ has $p(p + 1)/2$ distinct entries estimated from the $n_c$ points of class $c$ alone, and its maximum likelihood estimate is singular when $n_c \le p$ . Common remedies are fewer features and shrinkage, for example $(1 - \alpha)\hat{\Sigma}_c + \alpha \hat{\Sigma}$ toward the pooled matrix with $0 \le \alpha \le 1$ , or adding a positive multiple of the identity, with the amount chosen on validation data.

In plain words QDA lets each class keep its own shape, so a widely spread class can claim points far out on both sides of a tight class, and the boundary can bend around the tight class. The flexibility costs data, because every class must supply enough points to estimate its own covariance, and with few points the estimates are noisy and the boundary unstable. Whether QDA predicts better than LDA is an empirical question, settled by comparing both on data not used for fitting.

Example With $\hat{\Sigma}_A$ and $\hat{\Sigma}_B$ from above and equal priors,

\begin{aligned} \delta_B - \delta_A &= x_1^2 - \tfrac{(x_1 - 2)^2}{4} \\[0.4em] &\quad - \ln 2 . \end{aligned}

The second sensor drops out because both modes have variance $0.5$ there, and $-\ln 2$ comes from the determinants $1$ and $0.25$ . At $(0, 0)$ the difference is $-1 - \ln 2 \approx -1.6931$ , so $P(B \mid \mathbf{x}) \approx 0.1554$ , as in the first card. At $(3, 0)$ it is $9 - 0.25 - \ln 2 \approx 8.0569$ , giving $P(B \mid \mathbf{x}) \approx 0.9997$ against $0.9608$ under LDA, because a large first reading is far less surprising for the widely spread mode B. Setting the difference to $0$ gives $3x_1^2 + 4x_1 - 4 - 4\ln 2 = 0$ , with roots $x_1 \approx 0.9771$ and $x_1 \approx -2.3104$ . QDA therefore assigns mode B to readings with $x_1 > 0.9771$ and also to readings with $x_1 < -2.3104$ , a boundary of two vertical lines that a single LDA line cannot express. This shows what the covariance assumption does. Eight constructed points cannot show which model predicts better.

Check yourself

Keep the fitted LDA model, with $\delta_B - \delta_A = 1.6\,x_1 - 1.6$ under equal priors, but suppose the machines in service run in mode B with probability $0.8$ and in mode A with probability $0.2$ , with the class-conditional densities unchanged. Find the new boundary and say which way it moved.

Show the answer

The prior odds of B are $0.8/0.2 = 4$ , so the score difference gains $\ln 4$ and becomes $1.6\,x_1 - 1.6 + \ln 4$ . Setting it to $0$ gives

x_1 = 1 - \frac{\ln 4}{1.6} \approx 0.1336 .

The boundary moved from $x_1 = 1$ toward the mean of mode A, so the region assigned to the more common mode B grew. The answer assumes that the new prior is right and that the sensor distributions within each mode did not change. A new prior does not correct for drift in the sensors themselves.

A machine reads $(1, 0)$ . Find the probability of mode B under the fitted LDA model and under the fitted QDA model, both with equal priors.

Show the answer

Under LDA the score difference is $1.6 \cdot 1 - 1.6 = 0$ , so the reading lies on the boundary and $P(B \mid \mathbf{x}) = \sigma(0) = 0.5$ . A rule for ties is needed to give it a label. Under QDA the difference is

\begin{aligned} & 1 - \tfrac{(1 - 2)^2}{4} - \ln 2 \\[0.4em] &= 0.75 - \ln 2 \approx 0.0569 , \end{aligned}

so $P(B \mid \mathbf{x}) = \sigma(0.0569) \approx 0.5142$ and the equal-cost decision is mode B. The QDA boundary at $x_1 \approx 0.9771$ lies just to the left of this reading.

A data set has $p = 30$ sensors, with $20$ training machines in mode A and $100$ in mode B. Can QDA be fitted as defined here? Can LDA? What can be done?

Show the answer

QDA needs a positive-definite estimate of each class covariance. The estimate for mode A comes from $20$ points in $\mathbb{R}^{30}$ , so its rank is at most $19$ , and it is singular. Its inverse and its log determinant do not exist, so QDA as defined cannot be fitted. The pooled estimate for LDA uses the deviations of all $120$ machines from their class means and has rank at most $120 - 2 = 118$ , so it can be positive definite, and with sensors that are not exact linear combinations of each other it usually is. LDA estimates $465$ covariance entries, while QDA needs $2 \times 465 = 930$ , of which the $465$ for mode A would come from only $20$ points. The options are to use LDA, to shrink the mode A covariance toward the pooled matrix or toward a positive multiple of the identity, or to reduce the number of sensors, with the choice and the amount of shrinkage made on validation data. If two sensors were exact copies, even the pooled matrix would be singular, and one of them would have to be removed or the matrix regularized.

Gaussian Discriminant Analysis

Key terms

Class-conditional density

Multivariate Gaussian distribution

Positive-definite covariance matrix

Gaussian discriminant analysis

Linear discriminant analysis

Quadratic discriminant analysis

Check yourself

How this lesson connects

Builds on

Used later in