Linear Discriminant Analysis: Complete Guide to Linear Discriminant Analysis

Summary

Linear Discriminant Analysis (LDA) is a fundamental supervised classification and dimensionality reduction algorithm. Unlike PCA which simply looks for the directions of greatest variance, LDA actively seeks to maximize separation between classes while minimizing the dispersion within each class. Developed by Ronald Fisher in 1936 under the name “discriminant analysis,” this mathematically elegant method remains one of the most widely used classifiers in practice, particularly in the fields of face recognition, bioinformatics, and signal processing.

In this guide, we will explore in depth the mathematical principle behind LDA, its geometric intuition, a complete Python implementation from scratch, as well as its advantages, limitations, and concrete use cases.

Mathematical Principle

LDA is based on a central idea: project the data into a lower-dimensional space where the classes are as separated as possible. To formalize this intuition, Fisher introduced two fundamental scatter matrices.

Inter-class scatter matrix (S_B)

The inter-class scatter matrix measures the separation between the centers of each class and the overall center of the data. It is computed as follows:

S_B = Σ_k n_k · (μ_k − μ)(μ_k − μ)^T

where:
– n_k is the number of samples in class C_k
– μ_k is the mean vector of class C_k
– μ is the global mean vector of all the data
– T denotes the transpose

The larger the eigenvalues of this matrix, the farther apart the class centers are from each other.

Intra-class scatter matrix (S_W)

The intra-class scatter matrix measures the spread of points within each class. It is computed as follows:

S_W = Σk Σ{i ∈ C_k} (x_i − μ_k)(x_i − μ_k)^T

where:
– x_i is the i-th sample belonging to class C_k
– μ_k is the mean of class C_k

A small S_W matrix indicates that the points of each class are grouped compactly around their center.

Fisher’s Criterion

The goal of LDA is to find a projection vector w that maximizes the Fisher ratio:

J(w) = (w^T · S_B · w) / (w^T · S_W · w)

This ratio precisely measures the desired trade-off: a large numerator means the classes are well separated after projection, while a small denominator means the points of each class remain grouped.

Maximizing the criterion: to find the optimal w, we must solve the constrained optimization problem. Using the Lagrange multiplier, we derive that the optimal vector satisfies the generalized eigenvalue equation:

S_B · w = λ · S_W · w

which is equivalent to:

S_W^{-1} · S_B · w = λ · w

Solution for Two Classes

In the special case of two classes (binary classification), the solution has a remarkably simple analytical form:

w = S_W^{-1} · (μ_1 − μ_2)

This vector directly gives the optimal projection direction. The data is then projected onto this axis, and a threshold (usually the midpoint between the projected means) is used to classify new points.

Multi-class Case

For K classes, LDA can extract at most K − 1 discriminant components. We diagonalize the matrix S_W^{-1} · S_B and retain the K − 1 eigenvectors associated with the largest eigenvalues. These vectors form the basis of the optimal discriminant space.

Geometric Intuition

Imagine you have two crowds of people in a room — say, fans of two different football teams. You are standing on a balcony looking at the scene from above (2D projection). The two groups are a bit mixed; it’s hard to distinguish them clearly.

Now imagine you could change your viewing angle. If you look from a particular angle, the two groups appear clearly separated: Team A on the left, Team B on the right, with an empty space in between. This is exactly what LDA does — it finds the optimal projection axis that pushes class centers as far apart as possible while tightening the points within each class.

Why is this approach so effective? Because it leverages class label information, unlike PCA which completely ignores this information. PCA simply looks for the direction where the data varies most, but that direction is not necessarily the one that best separates the classes. LDA, on the other hand, is guided by supervision: it knows which data belongs to which class and takes advantage of that.

Important assumption: LDA assumes that each class follows a multivariate Gaussian (normal) distribution and that all classes share the same covariance matrix. This is a strong assumption that works well in practice when classes have similar shapes and sizes in feature space.

Python Implementation

1. LDA From Scratch with S_B and S_W computation

import numpy as np
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt


class LDAClassifier:
    """Linear Discriminant Analysis implemented from scratch."""

    def __init__(self):
        self.w = None          # Projection vector(s)
        self.means = None      # Means per class
        self.classes = None    # Unique labels
        self.threshold = None  # Decision threshold (binary case)

    def _compute_S_W(self, X, y):
        """Compute the within-class scatter matrix."""
        n_features = X.shape[1]
        S_W = np.zeros((n_features, n_features))
        for c in self.classes:
            X_c = X[y == c]
            mean_c = np.mean(X_c, axis=0)
            diff = X_c - mean_c
            S_W += diff.T @ diff
        return S_W

    def _compute_S_B(self, X, y):
        """Compute the between-class scatter matrix."""
        n_features = X.shape[1]
        overall_mean = np.mean(X, axis=0)
        S_B = np.zeros((n_features, n_features))
        for c in self.classes:
            X_c = X[y == c]
            n_c = X_c.shape[0]
            mean_c = np.mean(X_c, axis=0)
            diff = (mean_c - overall_mean).reshape(-1, 1)
            S_B += n_c * (diff @ diff.T)
        return S_B

    def fit(self, X, y):
        """LDA training."""
        self.classes = np.unique(y)
        self.means = {}
        for c in self.classes:
            self.means = np.mean(X[y == c], axis=0)

        S_W = self._compute_S_W(X, y)
        S_B = self._compute_S_B(X, y)

        # Regularizing S_W to ensure invertibility
        S_W += np.eye(S_W.shape[0]) * 1e-6

        # Solving the generalized eigenvalue problem
        S_W_inv = np.linalg.inv(S_W)
        A = S_W_inv @ S_B
        eigenvalues, eigenvectors = np.linalg.eigh(A)

        # Sorting eigenvalues in descending order
        idx = np.argsort(eigenvalues)[::-1]
        eigenvalues = eigenvalues[idx]
        eigenvectors = eigenvectors[:, idx]

        # For two classes: take the first eigenvector
        n_components = min(len(self.classes) - 1, X.shape[1])
        self.w = eigenvectors[:, :n_components]

        # Computing the threshold (binary case)
        if len(self.classes) == 2:
            proj0 = X[y == self.classes[0]] @ self.w
            proj1 = X[y == self.classes[1]] @ self.w
            self.threshold = (np.mean(proj0) + np.mean(proj1)) / 2

        return self

    def transform(self, X):
        """Project data into discriminant space."""
        return X @ self.w

    def predict(self, X):
        """Predict classes."""
        projected = self.transform(X)

        if len(self.classes) == 2:
            predictions = np.where(projected >= self.threshold,
                                    self.classes[1], self.classes[0])
        else:
            # Multi-class case: assign to nearest center
            predictions = np.zeros(len(X), dtype=int)
            for i in range(len(X)):
                dists = [np.linalg.norm(projected[i] - self.means @ self.w)
                         for c in self.classes]
                predictions[i] = self.classes[np.argmin(dists)]
        return predictions

    def score(self, X, y):
        """Classification accuracy."""
        return np.mean(self.predict(X) == y)


# --- Test on synthetic data ---
X, y = make_classification(n_samples=600, n_features=2, n_informative=2,
                           n_redundant=0, n_clusters_per_class=1,
                           class_sep=1.5, random_state=42)

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.25, random_state=42
)

lda = LDAClassifier()
lda.fit(X_train, y_train)
accuracy = lda.score(X_test, y_test)
print(f"LDA From Scratch — Accuracy: {accuracy:.4f}")

2. LDA with scikit-learn

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.metrics import accuracy_score, classification_report

# Training with sklearn
lda_sklearn = LinearDiscriminantAnalysis()
lda_sklearn.fit(X_train, y_train)

predictions = lda_sklearn.predict(X_test)
precision = accuracy_score(y_test, predictions)

print(f"LDA sklearn — Accuracy: {precision:.4f}")
print("\nClassification report:")
print(classification_report(y_test, predictions))

# Coefficients and threshold
print(f"\nProjection vectors: {lda_sklearn.coef_}")
print(f"Class means: {lda_sklearn.means_}")

3. 2D → 1D Visualization

def visualize_lda_projection(X, y, lda_model, title="LDA Projection"):
    """Visualizes the LDA projection from 2D to 1D."""
    X_proj = lda_model.transform(X)

    fig, axes = plt.subplots(1, 2, figsize=(14, 5))

    # Original 2D space
    ax1 = axes[0]
    for c in np.unique(y):
        mask = y == c
        ax1.scatter(X[mask, 0], X[mask, 1],
                    label=f"Class {c}", alpha=0.7, s=40)

    # LDA direction vector
    w = lda_model.w.ravel()
    w_norm = w / np.linalg.norm(w)
    ax1.quiver(np.mean(X[:, 0]), np.mean(X[:, 1]),
               w_norm[0] * 3, w_norm[1] * 3,
               color='red', scale=5, width=0.01, label='LDA direction')
    ax1.set_title("Original space (2D)")
    ax1.set_xlabel("Feature 1")
    ax1.set_ylabel("Feature 2")
    ax1.legend()
    ax1.grid(True, alpha=0.3)

    # 1D projection
    ax2 = axes[1]
    for c in np.unique(y):
        mask = y == c
        ax2.scatter(X_proj[mask], np.zeros_like(X_proj[mask]),
                    label=f"Class {c}", alpha=0.5, s=40)
    # Adding slight vertical noise for readability
    for c in np.unique(y):
        mask = y == c
        jitter = np.random.normal(0, 0.1, size=X_proj[mask].shape)
        ax2.scatter(X_proj[mask], jitter,
                    label=f"Class {c}", alpha=0.3, s=20)

    ax2.set_title("LDA Projection (1D)")
    ax2.set_xlabel("Discriminant axis")
    ax2.set_yticks([])
    ax2.legend()
    ax2.grid(True, alpha=0.3)

    plt.suptitle(title, fontsize=14, fontweight='bold')
    plt.tight_layout()
    plt.show()


visualize_lda_projection(X_test, y_test, lda, "LDA: From 2D to 1D")

4. LDA vs PCA Comparison

from sklearn.decomposition import PCA

def compare_lda_pca(X, y):
    """Visually compares LDA and PCA projections."""
    lda = LDAClassifier()
    lda.fit(X, y)
    X_lda = lda.transform(X)

    pca = PCA(n_components=1)
    X_pca = pca.fit_transform(X)

    fig, axes = plt.subplots(1, 2, figsize=(14, 5))

    ax1 = axes[0]
    for c in np.unique(y):
        mask = y == c
        jitter = np.random.normal(0, 0.1, size=X_lda[mask].shape)
        ax1.scatter(X_lda[mask] + np.random.normal(0, 0.05, size=jitter.shape),
                    jitter, label=f"Class {c}", alpha=0.4, s=15)
    ax1.set_title("LDA Projection (supervised)")
    ax1.set_xlabel("Discriminant axis")
    ax1.set_yticks([])
    ax1.legend()
    ax1.grid(True, alpha=0.3)

    ax2 = axes[1]
    for c in np.unique(y):
        mask = y == c
        jitter = np.random.normal(0, 0.1, size=X_pca[mask].shape)
        ax2.scatter(X_pca[mask], jitter,
                    label=f"Class {c}", alpha=0.4, s=15)
    ax2.set_title("PCA Projection (unsupervised)")
    ax2.set_xlabel("Principal component")
    ax2.set_yticks([])
    ax2.legend()
    ax2.grid(True, alpha=0.3)

    plt.suptitle("LDA vs PCA — Class separation quality",
                 fontsize=14, fontweight='bold')
    plt.tight_layout()
    plt.show()


compare_lda_pca(X, y)

Key observation: in the LDA vs PCA visualization, you typically observe that LDA separates the classes much better than PCA, because it uses label information when searching for the projection axis. PCA, by contrast, only maximizes total variance, which is not synonymous with good class separation.

Hyperparameters

The scikit-learn implementation exposes several important hyperparameters that are essential to understand for optimizing performance.

n_components

Number of discriminant components to extract. The maximum value is K − 1 where K is the number of classes. For a binary problem, only one component can be extracted. For a 4-class problem, up to 3 can be extracted. Reducing this parameter allows dimensionality reduction before applying another classifier.

lda = LinearDiscriminantAnalysis(n_components=2)  # For 3+ classes

solver

Solver algorithm used to compute the discriminant vectors:

• “svd” (default): uses singular value decomposition. Fast, does not compute the covariance matrix. Incompatible with shrinkage.
• “lsqr”: uses the least squares method. More efficient on large datasets. Supports shrinkage.
• “eigen”: directly solves the generalized eigenvalue problem. Supports both shrinkage and optimization of class separation.

lda = LinearDiscriminantAnalysis(solver='lsqr', shrinkage='auto')

shrinkage

Regularization technique that blends the empirical covariance matrix with an identity matrix. Particularly useful when the number of features is close to the number of samples (small sample size problem).

• None: no shrinkage (default)
• “auto”: automatic shrinkage using the Ledoit-Wolf lemma
• Float between 0 and 1: manual shrinkage coefficient (0 = none, 1 = pure identity)

lda = LinearDiscriminantAnalysis(solver='lsqr', shrinkage=0.3)

priors

Prior probabilities of the classes. By default, estimated from frequencies in the training data. They can be manually specified to correct for class imbalance or to incorporate domain knowledge.

lda = LinearDiscriminantAnalysis(priors=[0.3, 0.7])  # Class 0: 30%, Class 1: 70%

tol

Convergence tolerance for iterative solvers (lsqr, eigen). A stricter (smaller) value guarantees higher precision but increases computation time.

lda = LinearDiscriminantAnalysis(solver='eigen', tol=1e-6)

Advantages and Limitations

Advantages

1. Computational efficiency: very fast at training and inference, since the solution is analytical (no iterative optimization like neural networks).
2. No sensitive hyperparameters: unlike random forests or SVMs, LDA works well with its default parameters.
3. Built-in dimensionality reduction: LDA naturally projects data into a K−1 dimensional space, facilitating visualization and speeding up pipelines.
4. Interpretability: the coefficients w are directly interpretable — they indicate the importance of each feature for discriminating between classes.
5. Well-calibrated probabilities: LDA provides well-calibrated posterior probabilities (thanks to the Gaussian assumption), unlike other classifiers that require recalibration (Platt scaling, isotonic regression).

Limitations

1. Equal covariance assumption: LDA assumes that all classes share the same covariance matrix. If this assumption is strongly violated, performance drops. In that case, QDA (Quadratic Discriminant Analysis) is more appropriate.
2. Gaussian assumption: each class is assumed to follow a multivariate normal distribution. For highly skewed or multimodal data, performance can be poor.
3. Linearity: decision boundaries are linear. If the classes are not linearly separable (e.g., XOR), LDA cannot distinguish them.
4. Sensitivity to outliers: the S_B and S_W matrices use means, which are highly sensitive to extreme values. Robust preprocessing (robust scaling, trimming) is often necessary.
5. Limited number of components: at most K − 1 components can be extracted. For a binary problem, only one discriminant dimension is available.

4 Concrete Use Cases

1. Face Recognition (Eigenfaces & Fisherfaces)

LDA is used in the famous Fisherfaces algorithm for face recognition. After a preliminary PCA step to reduce dimensionality, LDA is applied to maximize the separation between the faces of different individuals. Each person forms a class, and the discriminant projection produces highly discriminative face signatures.

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.decomposition import PCA

# Typical Fisherfaces pipeline
pipeline = make_pipeline(PCA(n_components=50), LinearDiscriminantAnalysis(n_components=9))
pipeline.fit(training_faces, identities)

2. Medical Diagnosis

In bioinformatics, LDA is commonly used for tumor classification from gene expression data. The thousands of genes (features) are reduced to a few discriminant axes that effectively separate healthy tissue from cancerous tissue. Its interpretability is a major asset: the w coefficients identify the most discriminating genes, offering concrete biological leads.

3. Anti-Spam Filtering

An LDA classifier can be used as a basic anti-spam filter. Each email is represented by a feature vector (word frequency, presence of links, length, etc.) and LDA learns to separate spam from legitimate emails. Although less powerful than modern methods, it remains fast, interpretable, and effective on well-structured datasets.

4. Sentiment Analysis

LDA serves as a baseline classifier for sentiment analysis (positive/negative) on vectorized texts (TF-IDF, embeddings). Its training speed makes it an excellent comparison point: if a complex method (BERT, neural networks) doesn’t significantly outperform LDA, the problem may already be sufficiently solved by a simple approach.

See Also

• Résoudre le Problème d’Entretien Maximum Subarray avec Python
• Implémenter l’Algorithme de Dekker en Python : Synchronisation des Threads Efficace