Linear Discriminant Analysis (LDA)

The Story That Explains LDA

Linear Discriminant Analysis (LDA) is a supervised dimensionality-reduction and classification algorithm invented by Ronald Fisher in 1936. Unlike PCA which maximises total variance, LDA maximises the ratio of between-class scatter to within-class scatter — deliberately finding the projection that makes classes most separable.

The Core Intuition — Two Objectives at Once

3D Interactive Visualisation — LDA in Action

The diagram below shows three classes in 3D space (left) and the same data projected onto the two LDA discriminant axes (right). Drag to rotate the 3D view. Notice how the LDA projection pulls the three clusters apart far more cleanly than random projections would.

Step-by-Step Maths — A Worked Example

Numerical Mini-Example (2 Features, 2 Classes, 4 Points)

LDA Assumptions — The Fine Print

Python Implementation — Wine Dataset

LDA as a Dimensionality Reducer + Visualiser

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.datasets import load_wine
from sklearn.preprocessing import StandardScaler
import numpy as np
import matplotlib.pyplot as plt

# Load the UCI Wine dataset — 178 samples, 13 chemical features, 3 classes
wine = load_wine()
X, y = wine.data, wine.target
labels = wine.target_names  # ['class_0', 'class_1', 'class_2']

# Standardise features (recommended for LDA)
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# Fit LDA — project down to 2 components (C-1 = 3-1 = 2 max)
lda = LinearDiscriminantAnalysis(n_components=2)
X_lda = lda.fit_transform(X_scaled, y)

# Explained variance ratio of each discriminant
print("Explained variance ratio:", lda.explained_variance_ratio_)
print(f"LD1 explains: {lda.explained_variance_ratio_[0]*100:.1f}%")
print(f"LD2 explains: {lda.explained_variance_ratio_[1]*100:.1f}%")

# Plot the 2D projection
colors = ['#f87171', '#60a5fa', '#34d399']
plt.figure(figsize=(8,6))
for i, label in enumerate(labels):
    mask = y == i
    plt.scatter(X_lda[mask, 0], X_lda[mask, 1],
                label=label, color=colors[i], alpha=0.8, edgecolors='white', s=60)
plt.xlabel('LD1'); plt.ylabel('LD2')
plt.title('Wine Dataset — LDA 2D Projection')
plt.legend(); plt.tight_layout(); plt.show()

LDA as a Classifier

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.model_selection import cross_val_score, StratifiedKFold
from sklearn.metrics import classification_report, confusion_matrix
from sklearn.datasets import load_wine
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline

# Build a proper pipeline: scale → LDA classifier
pipe = Pipeline([
    ('scaler', StandardScaler()),
    ('lda',    LinearDiscriminantAnalysis())   # n_components=None → uses all C-1
])

wine = load_wine()
X, y = wine.data, wine.target

# 5-fold stratified cross-validation
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
scores = cross_val_score(pipe, X, y, cv=cv, scoring='accuracy')

print(f"CV Accuracy: {scores.mean():.4f} +/- {scores.std():.4f}")

# Full fit + detailed report
from sklearn.model_selection import train_test_split
X_tr, X_te, y_tr, y_te = train_test_split(
    X, y, test_size=0.25, stratify=y, random_state=42
)
pipe.fit(X_tr, y_tr)
y_pred = pipe.predict(X_te)

print(classification_report(y_te, y_pred, target_names=wine.target_names))
print("Confusion Matrix:")
print(confusion_matrix(y_te, y_pred))

Regularised LDA — Handling Small Samples & Many Features

Standard LDA breaks when features outnumber samples (singular S_W) or when the training set is small. Regularised LDA (RLDA) adds a shrinkage penalty that pulls the within-class covariance estimate toward the identity matrix:

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.datasets import make_classification
from sklearn.model_selection import cross_val_score
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline

# Simulate high-dimensional, small-sample scenario
# 200 samples, 150 features — classic small-n-large-p regime
X, y = make_classification(
    n_samples=200, n_features=150, n_informative=30,
    n_classes=3, n_clusters_per_class=1, random_state=42
)

# Standard LDA — likely to fail or be unstable
std_pipe = Pipeline([
    ('sc', StandardScaler()),
    ('lda', LinearDiscriminantAnalysis(solver='svd'))    # SVD-based, no shrinkage
])

# Regularised LDA — Ledoit-Wolf automatic shrinkage
reg_pipe = Pipeline([
    ('sc', StandardScaler()),
    ('lda', LinearDiscriminantAnalysis(
        solver='lsqr',
        shrinkage='auto'   # Ledoit-Wolf optimal shrinkage
    ))
])

std_scores = cross_val_score(std_pipe, X, y, cv=5, scoring='accuracy')
reg_scores = cross_val_score(reg_pipe, X, y, cv=5, scoring='accuracy')

print(f"Standard LDA CV:    {std_scores.mean():.4f} +/- {std_scores.std():.4f}")
print(f"Regularised LDA CV: {reg_scores.mean():.4f} +/- {reg_scores.std():.4f}")

# Manual shrinkage tuning (grid search over λ)
from sklearn.model_selection import GridSearchCV
param_grid = {'lda__shrinkage': [0.0, 0.1, 0.3, 0.5, 0.7, 0.9, 1.0, 'auto']}
gs = GridSearchCV(reg_pipe, param_grid, cv=5, scoring='accuracy')
gs.fit(X, y)
print(f"Best shrinkage: {gs.best_params_}")
print(f"Best CV accuracy: {gs.best_score_:.4f}")

Real-World Example — Face Recognition (AT&T Dataset)

from sklearn.datasets import fetch_olivetti_faces
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.decomposition import PCA
from sklearn.pipeline import Pipeline
from sklearn.model_selection import cross_val_score, StratifiedKFold
from sklearn.preprocessing import StandardScaler

# Olivetti Faces: 400 images, 4096 features, 40 classes
data = fetch_olivetti_faces(shuffle=True, random_state=42)
X, y = data.data, data.target  # X shape: (400, 4096)

# Fisherfaces pipeline: PCA (Eigenfaces) → LDA (Fisherfaces)
fisherfaces = Pipeline([
    ('scale', StandardScaler()),
    ('pca',   PCA(n_components=150, whiten=True, random_state=42)),
    ('lda',   LinearDiscriminantAnalysis())   # max 39 components (40 classes - 1)
])

cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
scores = cross_val_score(fisherfaces, X, y, cv=cv, scoring='accuracy')

print("=== Fisherfaces (PCA → LDA) ===")
print(f"CV Accuracy: {scores.mean()*100:.1f}% +/- {scores.std()*100:.1f}%")
print(f"Input dims:  4096 → PCA: 150 → LDA: 39")

# Compare: PCA-only nearest centroid
from sklearn.neighbors import NearestCentroid
pca_only = Pipeline([
    ('scale', StandardScaler()),
    ('pca',   PCA(n_components=150, whiten=True, random_state=42)),
    ('clf',   NearestCentroid())
])
pca_scores = cross_val_score(pca_only, X, y, cv=cv, scoring='accuracy')
print("\n=== PCA-only nearest centroid ===")
print(f"CV Accuracy: {pca_scores.mean()*100:.1f}% +/- {pca_scores.std()*100:.1f}%")

LDA vs PCA — When to Use Each

LDA vs QDA — Covariance Assumptions

When LDA's equal-covariance assumption fails, Quadratic Discriminant Analysis (QDA) steps in. QDA estimates a separate covariance matrix Σ_c for each class, resulting in quadratic (curved) decision boundaries rather than linear ones.

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis, QuadraticDiscriminantAnalysis
from sklearn.datasets import make_classification
from sklearn.model_selection import cross_val_score
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline

# Generate data with UNEQUAL class covariances — LDA assumption violated
import numpy as np
np.random.seed(42)

X_c1 = np.random.multivariate_normal([0,0], [[1,0],[0,1]], 200)     # circular
X_c2 = np.random.multivariate_normal([3,3], [[4,2],[2,4]], 200)     # elongated ellipse
X_c3 = np.random.multivariate_normal([6,0], [[0.5,0],[0,3]],200)    # narrow tall

X = np.vstack([X_c1, X_c2, X_c3])
y = np.array([0]*200 + [1]*200 + [2]*200)

lda_pipe = Pipeline([('sc', StandardScaler()), ('clf', LinearDiscriminantAnalysis())])
qda_pipe = Pipeline([('sc', StandardScaler()), ('clf', QuadraticDiscriminantAnalysis())])

lda_cv = cross_val_score(lda_pipe, X, y, cv=5).mean()
qda_cv = cross_val_score(qda_pipe, X, y, cv=5).mean()

print(f"LDA CV Accuracy: {lda_cv:.4f}  (assumption violated)")
print(f"QDA CV Accuracy: {qda_cv:.4f}  (no equal-covariance assumption)")

Medical Diagnosis Example — Breast Cancer LDA

from sklearn.datasets import load_breast_cancer
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.model_selection import cross_val_score, StratifiedKFold
from sklearn.metrics import roc_auc_score, classification_report
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
import numpy as np

bc = load_breast_cancer()
X, y = bc.data, bc.target   # 569 samples, 30 features, 2 classes

pipe = Pipeline([
    ('sc',  StandardScaler()),
    ('lda', LinearDiscriminantAnalysis())
])

# 10-fold stratified CV — important for medical data
cv = StratifiedKFold(n_splits=10, shuffle=True, random_state=42)
acc_scores = cross_val_score(pipe, X, y, cv=cv, scoring='accuracy')
auc_scores = cross_val_score(pipe, X, y, cv=cv, scoring='roc_auc')

print(f"CV Accuracy: {acc_scores.mean():.4f} +/- {acc_scores.std():.4f}")
print(f"CV AUC-ROC:  {auc_scores.mean():.4f} +/- {auc_scores.std():.4f}")

# Full fit for feature weights (discriminant loadings)
pipe.fit(X, y)
lda_model = pipe.named_steps['lda']

# Top contributing features to LD1
coeff = np.abs(lda_model.coef_[0])
top_idx = np.argsort(coeff)[::-1][:5]
print("\nTop 5 most discriminating features:")
for i in top_idx:
    print(f"  {bc.feature_names[i]:35s}: {coeff[i]:.4f}")

LDA Solvers in Scikit-Learn

LDA for Multi-Class — Iris 3D Projection

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.datasets import load_iris
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import cross_val_score
from sklearn.pipeline import Pipeline
import numpy as np

iris = load_iris()
X, y = iris.data, iris.target   # 150 samples, 4 features, 3 classes

pipe = Pipeline([
    ('sc',  StandardScaler()),
    ('lda', LinearDiscriminantAnalysis(n_components=2))
])

# Cross-validated accuracy
scores = cross_val_score(pipe, X, y, cv=10)
print(f"CV Accuracy: {scores.mean():.4f} +/- {scores.std():.4f}")

# Fit and show explained variance per discriminant
pipe.fit(X, y)
lda = pipe.named_steps['lda']
print(f"\nLD1 explains {lda.explained_variance_ratio_[0]*100:.1f}% of between-class variance")
print(f"LD2 explains {lda.explained_variance_ratio_[1]*100:.1f}% of between-class variance")

# Class separation quality: Mahalanobis-like distances between centroids in LDA space
X_lda = lda.transform(pipe.named_steps['sc'].transform(X))
centroids = [X_lda[y==i].mean(axis=0) for i in range(3)]
for i in range(3):
  for j in range(i+1, 3):
    d = np.linalg.norm(np.array(centroids[i])-np.array(centroids[j]))
    print(f"Distance {iris.target_names[i]} ↔ {iris.target_names[j]}: {d:.3f}")

When to Use (and Avoid) LDA

Golden Rules