Linear, Quadratic, and Regularized Discriminant Analysis

The probabilistic interpretation

Each class \(k \in \{1, \ldots, K\}\) is assigned a prior \(\hat{\pi}_k\) such that that \(\sum_{i=1}^k \hat{\pi}_k = 1\). According to Bayes’ rule, the posterior probability is

\[\rm{Pr}(G = k |X = x) = \frac{f_k(x) \pi_k}{\sum_{l=1}^K f_l(x) \pi_l} \]

where \(f_k(x)\) is the density of \(X\) conditioned on \(k\). The maximum-a-posteriori (MAP) estimator simplifies to

\[G(x) = \arg \max_k \rm{Pr}(G = k | X = x) = \arg \max_k f_k(x) \pi_k \]

because the denominator is identical for all classes.

LDA assumes that the density is Gaussian:

\[f_k(x) = {|2 \pi \Sigma_k|}^{-1/2} \exp\left(-\frac{1}{2}(x - \mu_k)^T\Sigma_k^{-1}(x - \mu_k)\right)\]

where \(\Sigma_k\) is the covariance matrix for the samples from class \(k\) and \(|.|\) is the determinant. LDA assumes that all classes have the same covariance matrix, i.e. \(\Sigma_k = \Sigma\,, \forall k\).

Plugging \(f_k\) into the classification function, we arrive at the classification function:

\[G(x) = \arg \max_k \delta_k(x)\]

where

\[\delta_k(x) = x^T \Sigma^{-1} \mu_k -\frac{1}{2} \mu_k^T \Sigma^{-1} \mu_k + \log \pi_k \]

is the discriminant function for class \(k\). So, now that we have a classifier, how we can compute it?

To find the covariance matrix, we simply compute

\[\hat{\Sigma} = \sum_{k=1}^K \frac{1}{N - K} \sum_{g_i = k} (x_i - \hat{\mu}_k) (x_i - \hat{\mu}_k)^T\,.\]

Note that the deviation from the means is divided by \(N-K\), the degrees of freedom, to obtain an unbiased estimator.

The means of the classes, which are also called centroids, are defined by

\[\hat{\mu}_k = \frac{1}{N_k} \sum_{g_i = k} x_i\,.\]

The priors \(\pi_k\) are set to the prevalence ratio of the class-specific observations:

\[\hat{\pi}_k = \frac{N_k}{N}\,.\]

With this, we have defined all parameters required for the classifier.

The dimensionality reduction procedure of LDA involves the within-class variance, \(W = \hat{\Sigma}\), and the between-class variance, \(B\). The between-class variance indicates the deviation of centroids from the overall mean, \(\hat{\mu} = \sum_{k=1}^K \hat{\pi}_k \hat{\mu}_k\), and is defined as:

\[B = \sum_{k=1}^K \hat{\pi}_k(\hat{\mu}_k - \hat{\mu}) (\hat{\mu}_k - \hat{\mu})^T\,.\]

Finding a sequence of optimal substeps involves three steps:

1. Compute the \(K \times p\) matrix \(M\) containing the centroids, \(\mu_k\), and determine the common covariance matrix \(W\).
2. Compute \(M^{\ast} = MW^{-\frac{1}{2}}\) using the eigen-decomposition of W.
3. Compute \(B^{\ast}\) (the between-class covariance) and its eigen-decomposition \(B^{\ast} = V^{\ast} D_B {V^{\ast}}^T\). The columns \(v^{\ast}_l\) of \(V^{\ast}\) define the coordinates of the reduced subspace.

The \(l\)-th discriminant variable (one of the \(K-1\) new dimensions) is determined by \(Z_l = v_l^T X\) with \(v_l = W^{-\frac{1}{2}} v^{\ast}_l\).

The interpretation of Fisher

Fisher’s LDA optimization criterion states that the centroids of the groups should be spread out as far as possible. This amounts to finding a linear combination \(Z = a^T X\) such that \(a^T\) maximizes the between-class variance relative to the within-class variance.

As before, the within-class variance is \(W\) is the pooled covariance matrix, \(\hat{\Sigma}\), which indicates the deviation of all observations from their class centroids. The between-class variance is defined according to the deviation of the centroids from the overall mean, as defined earlier. For \(Z\), the between class variance is \(a^T B a\) and the within-class variance is \(A^T W a\). Thus, LDA can be optimized through the Rayleigh quotient

\[\max_a \frac{a^T B a}{a^T W a}\,, \]

which defines an optimal mapping of \(X\) to the new space \(Z\). Note that \(Z \in \mathbb{R}^{1 \times p}\), that is, the observations are mapped to a single dimension. To obtain additional dimensions, we need to solve the optimization problem for \(a_1, \ldots, a_{K-1}\) where each successive \(a_k\) is constructed to be orthogonal in \(W\) to the previous discriminant coordinates. This leads to the linear transformation \(G = (Z_1^T, Z_2^T, \ldots, Z_{K-1}^T) \in \mathbb{R}^{p \times q}\) with which we can map from \(p\) to \(q\) dimension via \(X G\). Why do we consider \(K-1\) projections? This is because the affine subspace that is spanned by the \(K\) centroids has a rank of at most \(K-1\).

Using Fisher’s formulation of LDA, classification involves two steps:

1. Sphere the data using the common covariance matrix \(\hat{\Sigma} = UDU^T\) (eigendecomposition) such that \(X^{\ast} = D^{-\frac{1}{2}} U^T X\). In this way, the covariance of \(X^{\ast}\) becomes the identity matrix. By eliminating the covariance between the variables in this way, it becomes much easier to separate the classes in the transformed space.
2. Classify observations \(x_i\) to the closest class centroid in the transformed space, taking into account the class priors \(\pi_k\). Here, the intuition is that an observation with equal distance to two centroids would be assigned to the class with the greater prior.

Reduced-rank LDA

LDA performs classification in a reduced subspace. When performing classification, we do not need to use all \(K-1\) dimensions and instead can choose a smaller subspace \(H_l\) with \(l < K-1\). When \(l < K-1\) is used, this is called reduced-rank LDA. The motivation for reduced-rank LDA is that classification basd on a reduced number of discriminant variables can improve performance on the test set when the model is overfitted.

Complexity of the LDA model

The number of effective parameters of LDA can be derived in the following way. There are \(K\) means, \(\hat{\mu}_k\) that are estimated. The covariance matrix does not require additional parameters because it is already defined by the centroids. Since we need to estimate \(K\) discriminant functions (to obtain the decision boundaries), this gives rise to \(K\) calculations involving the \(p\) elements. Additionally, we have \(K-1\) free parameters for the \(K\) priors. Thus, the number of effective LDA parameters is \(Kp + (K-1)\).

Summary of LDA

Here, I summarize the two perspectives on LDA and give a summary of the main properties of the model.

Probabilistic view

LDA uses Bayes’ rule to determine the posterior probability that an observation \(x\) belongs to class \(k\). Due to the normal assumption of LDA, the posterior is defined by a multivariate Gaussian whose covariance matrix is assumed to be identical for all classes. New points are classified by computing the discriminant function \(\delta_k\) (the enumerator of the posterior probability) and returning the class \(k\) with maximal \(\delta_k\). The discriminant variables can be obtained through eigen-decompositions of the within-class and between-class variance.

The phoneme data set

To exemplify linear discriminant analysis, we will use the phoneme speech recognition data set. This data set is useful for showcasing discriminant analysis because it involves five different outcomes.

library(RCurl)
f <- getURL('https://www.datascienceblog.net/data-sets/phoneme.csv')
df <- read.csv(textConnection(f), header=T)
print(dim(df))

## [1] 4509  259

The data set contains samples of digitized speech for five phonemes: aa (as the vowel in dark), ao (as the first vowel in water), dcl (as in dark), iy (as the vowel in she), and sh (as in she). In total, 4509 speech frames of 32 msec were selected. For each speech frame, a log-periodogram of length 256 was computed, on whose basis we want to perform speech recognition. The 256 columns labeled x.1 to x.256 identify the speech features, while the columns g and speaker indicate the phonemes (labels) and speakers, respectively.

For evaluating models later, we will assign each sample either into the training or the test set:

#logical vector: TRUE if entry belongs to train set, FALSE else
train <- grepl("^train", df$speaker)
# remove non-feature columns
to.exclude <- c("row.names", "speaker", "g")
feature.df <- df[, !colnames(df) %in% to.exclude]
test.set <- subset(feature.df, !train)
train.set <- subset(feature.df, train)
train.responses <- subset(df, train)$g
test.responses <- subset(df, !train)$g

Fitting an LDA model in R

We can fit an LDA model in the following way:

library(MASS)
lda.model <- lda(train.set, grouping = train.responses)

Let us take a moment to investigate the relevant components of the model:

• prior contains the group priors \(\pi_k\) and counts the corresponding group counts, \(N_k\).
• means is the centroid matrix \(M \in \mathbb{R}^{K \times p}\) whose components are the mean vectors, \(\mu_k\).
• scaling is the \(N \times (K-1)\) matrix that transforms samples to the space defined by the \(K-1\) discriminant variables.
• svd are the singular values, which indicate the ratio of between- and within-group standard deviations on the linear discriminant variables.

# 1. Manual transformation
# center data around weighted means & transform
means <- colSums(lda.model$prior * lda.model$means)
train.mod <- scale(train.set, center = means, scale = FALSE) %*% 
             lda.model$scaling
# 2. Use the predict function to transform:
lda.prediction.train <- predict(lda.model, train.set)
all.equal(lda.prediction.train$x, train.mod)

## [1] TRUE

# visualize the features in the two LDA dimensions
plot.df <- data.frame(train.mod, "Outcome" = train.responses)
library(ggplot2)
ggplot(plot.df, aes(x = LD1, y = LD2, color = Outcome)) + geom_point()

library(RColorBrewer)
colors <- brewer.pal(8, "Accent")
my.cols <- colors[match(lda.prediction.train$class, levels(df$g))]
plot(lda.model, dimen = 4, col = my.cols)

plot_lda_centroids <- function(lda.model, train.set, response) {
    centroids <- predict(lda.model, lda.model$means)$x
    library(RColorBrewer)
    colors <- brewer.pal(8, "Accent")
    my.cols <- colors[match(lda.prediction.train$class, levels(df$g))]
    my.points <- predict(lda.model, train.set)$x
    no.classes <- length(unique(response))
    par(mfrow=c(no.classes -1, no.classes -1), mar=c(1,1,1,1), oma=c(1,1,1,10))
    for (i in 1:(no.classes - 1)) {
        for (j in 1:(no.classes - 1)) {
            y <- my.points[, i]
            x <- my.points[, j]
            cen <- cbind(centroids[, j], centroids[, i])
            if (i == j) {
                plot(x, y, type="n") 
                max.y <- max(my.points[, i])
                max.x <- max(my.points[, j])
                min.y <- min(my.points[, i])
                min.x <- min(my.points[, j])
                max.both <- max(c(max.x, max.y))
                min.both <- max(c(min.x, min.y))
                center <- min.both + ((max.both - min.both) / 2)
                text(center, center, colnames(my.points)[i], cex = 3)}
            else {
                plot(x, y, col = my.cols, pch = as.character(response), xlab ="", ylab="")
                points(cen[,1], cen[,2], pch = 21, col = "black", bg = colors, cex = 3)
            }
        }
    }
    par(xpd = NA)
    legend(x=par("usr")[2] + 1, y = mean(par("usr")[3:4]) + 20, 
            legend = rownames(centroids), col = colors, pch = rep(20, length(colors)), cex = 3)
}
plot_lda_centroids(lda.model, train.set, train.responses)

posteriors <- lda.prediction.train$posterior # N x K matrix
# MAP classification for sample 1:
pred.class <- names(which.max(posteriors[1,])) # <=> lda.prediction.train$class[1]
print(paste0("Posterior of predicted class '", pred.class, 
    "' is: ", round(posteriors[1,pred.class], 2)))

## [1] "Posterior of predicted class 'sh' is: 1"

# what are the mean posteriors for individual groups?
classes <- colnames(posteriors)
res <- do.call(rbind, (lapply(classes, function(x) apply(posteriors[train.responses == x, ], 2, mean))))
rownames(res) <- classes
print(round(res, 3)) 

##        aa    ao   dcl    iy    sh
## aa  0.797 0.203 0.000 0.000 0.000
## ao  0.123 0.877 0.000 0.000 0.000
## dcl 0.000 0.000 0.985 0.014 0.002
## iy  0.000 0.000 0.001 0.999 0.000
## sh  0.000 0.000 0.000 0.000 1.000

dims <- 1:4 # number of canonical variables to use
accuracies <- rep(NA, length(dims))
for (i in seq_along(dims)) {
    lda.pred <- predict(lda.model, test.set, dim = dims[i])
    acc <- length(which(lda.pred$class == test.responses))/length(test.responses)
    accuracies[i] <- acc
}
reduced.df <- data.frame("Rank" = dims, "Accuracy" = round(accuracies, 2))
print(reduced.df)

##   Rank Accuracy
## 1    1     0.51
## 2    2     0.71
## 3    3     0.86
## 4    4     0.92

lda.pred <- predict(lda.model, test.set)
plot.df <- data.frame(lda.pred$x[,1:2], "Outcome" = test.responses, 
                    "Prediction" = lda.pred$class)
ggplot(plot.df, aes(x = LD1, y = LD2, color = Outcome, shape = Prediction)) +
        geom_point()

QDA in R

We can perform QDA in the following way:

qda.model <- qda(train.set, grouping = train.responses)

The main difference between the QDA and the LDA object is that the QDA has a \(p \times p\) transformation matrix for every class \(k \in \{1, \ldots, K\}\). These matrices ensure that the within-group covariance matrix is spherical but do not induce a reduced subspace. Thus, QDA cannot be used as a visualization technique.

Let us determine whether QDA outperforms LDA on the phoneme data set:

qda.preds <- predict(qda.model, test.set)
acc <- length(which(qda.preds$class == test.responses))/length(test.responses)
print(paste0("Accuracy of QDA is: " , round(acc, 2)))

## [1] "Accuracy of QDA is: 0.84"

The accuracy of QDA is slightly below that of full-rank LDA. This could indicate that the assumption of common covariance is suitable for this data set.

RDA in R

To perform RDA, we need to load the rda package. Note that rda was removed from CRAN, so you need to install it manually now (2020-07). A regularized discriminant analysis model can be fit using the rda function, which has two main parameters: \(\alpha\) as introduced before and \(\delta\), which defines the threshold for values.

library(rda)
# note: optimization is time-intensive, especially if many alpha/gammas are used
alphas <- c(0.1, 0.25, 0.5, 0.75, 0.9)
rda.model <- rda(t(as.matrix(train.set)), y = train.responses,
                alpha = alphas)

Having obtained the fitted model, let us determine the performance of RDA on the test set. The following code is a bit more complicated than the code for LDA/QDA because we need to iterate over the results for both \(\alpha\) and \(\delta\):

rda.preds <- predict(rda.model, t(train.set), train.responses, t(test.set))
# determine performance for each alpha:
rda.perf <- vector("list", dim(rda.preds)[1])
for (i in seq(dim(rda.preds)[1])) {
    # performance for each gamma
    res <- apply(rda.preds[i,,], 1, function(x) length(which(x == as.numeric(test.responses))) / length(test.responses))
    rda.perf[[i]] <- res
}
rda.perf <- do.call(rbind, rda.perf)
rownames(rda.perf) <- alphas

The calculated accuracies suggest that the RDA classifier with \(\delta = 0\) and \(\alpha = 0.25\) (92.2% accuracy) slightly outperforms full-rank LDA (91.9% accuracy).