Linear Discriminant Analysis (LDA)

September 2026 · 17 min read

LDA uses class labels to find a direction that separates groups. It favors a large gap between class centers and a small spread of points within each class. This can reduce the number of measurements while keeping information useful for classification.

This builds on the PCA guide. PCA looks for the direction with the most overall spread, even if that direction does not separate the labels. LDA uses the labels to choose its direction. We will compare both methods on eight points, then work out how LDA finds its direction.

The widest direction can hide the classes

There are four points in each class. Class A is (−4, −2), (−4, −1), (4, −2), and (4, −1). Class B is (−4, 1), (−4, 2), (4, 1), and (4, 2). Both classes have the same first-feature values, −4 and 4. So that feature varies a lot, but it does not tell us which class a point belongs to.

PCA picks the horizontal direction because it contains about 86.5% of the total variation. But along that line, A and B overlap completely. The vertical direction separates all eight training points. Use the controls to rotate the line and watch how the class scores change:

This is a small training example, not proof that the result will work on new data. It shows why keeping most of the overall variation does not always keep the information that separates classes.

Measure separation against the spread inside each class

Choose a direction using weights $w$. A point $x_i$ becomes one score, $z_i=w^\top x_i$: multiply each measurement by its weight and add the results. Let $\bar z_A$ and $\bar z_B$ be the average scores for the two classes.

We want those averages far apart, but that gap is only useful if points within each class stay close together. A gap of 3 is not impressive if each class spreads from −100 to 100. Fisher's rule compares the squared gap between class averages with the total squared distance of each point's score from its own class average. If $y_i$ is point $i$'s label, the two-class criterion is:

For the horizontal direction, both class averages are 0, so the ratio is 0. For the vertical direction, their averages are −1.5 and 1.5, giving a squared gap of 9. Each of the eight points is 0.5 from its own class average, so the total within-class spread is $8\times0.5^2=2$. The ratio is $9/2=4.5$.

If we multiply every weight by 10, all scores grow by 10 and both parts of the ratio grow by 100. The ratio stays the same, so stretching the score axis cannot improve it. This number measures separation; it is not a probability or a test of prediction accuracy.

Find the direction without trying every angle

To find the best direction, summarize the spread inside the classes with the within-class scatter matrix, $S_W$. For each point, subtract its own class average. The diagonal entries add the squared differences for each measurement. The off-diagonal entries add products of differences from two measurements, showing whether those differences move together.

Let $\mu_A$ and $\mu_B$ be the class averages, and let $d=\mu_B-\mu_A$ be the gap between them. If $S_W$ can be inverted, the direction that maximizes Fisher's ratio is found by solving $S_Ww=d$. This accounts for how much spread there is inside each class. In our example:

The first equation gives $w_1=0$ and the second gives $w_2=1.5$. Divide by its length to get the easier-to-read direction $w=(0,1)$. Here LDA keeps the second feature because that feature separates the class averages. In other datasets, the off-diagonal entries also matter: they account for measurements that vary together inside the classes.

Fit with labels; transform new points with the saved direction

After fitting the direction, we can score a new point using its measurements and the saved weights. We also use the same score origin as during training. Here, zero is halfway between the two class averages, at (0, 0), so the new point (8, 1.2) gets score 1.2.

import numpy as np
from sklearn.decomposition import PCA
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis

X = np.array([[-4, -2], [-4, -1], [4, -2], [4, -1],
              [-4,  1], [-4,  2], [4,  1], [4,  2]], dtype=float)
y = np.repeat([0, 1], 4)  # A = 0, B = 1.
means = np.array([X[y == c].mean(axis=0) for c in [0, 1]])
residuals = X - means[y]  # Center each row at its own class mean.
SW = residuals.T @ residuals
w = np.linalg.solve(SW, means[1] - means[0])
w /= np.linalg.norm(w)  # Fix a convenient score scale.
origin = (means[0] + means[1]) / 2

def transform(new):
    return (np.asarray(new) - origin) @ w

print(SW)                      # [[128, 0], [0, 2]]
print(w)                       # [0, 1]
print(transform([[8, 1.2]]))    # [1.2]

pca = PCA(n_components=1).fit(X)
print(pca.explained_variance_ratio_)  # About [0.8649].
lda = LinearDiscriminantAnalysis(solver="eigen", n_components=1).fit(X, y)
reference_direction = lda.scalings_[:, 0].copy()
reference_direction /= np.linalg.norm(reference_direction)
# Compare axes, allowing an arbitrary sign.
np.testing.assert_allclose(np.outer(w, w),
                           np.outer(reference_direction, reference_direction), atol=1e-10)
# Coordinate scale and origin can differ across implementations.
np.testing.assert_array_equal(lda.predict(X), y)  # Training check, not a test score.

A projection still needs a classification rule

The projection turns each point into a score, but a score alone does not assign a class. We need a threshold or another rule. In this symmetric example, a threshold of zero separates the observed A and B points.

The probabilistic LDA classifier makes an additional assumption: each class forms a Gaussian, or bell-shaped, cloud. The clouds can have different centers, but share the same covariance. LDA combines these distributions with class priors to make predictions, giving a straight-line decision boundary. Fisher's separation calculation does not require Gaussian data. Either method still needs to be checked on new data.

Why C classes give at most C − 1 directions

With several classes, LDA still looks for directions that separate class averages while accounting for spread inside each class. Two class averages can be separated along one line; three can spread across at most a plane. In general, $C$ class averages give at most $C-1$ independent directions. With $D$ input features, the output can have at most $\min(D,C-1)$ directions, and sometimes fewer.

That is why two classes give at most one LDA direction. Asking for more cannot create new ways to separate their averages. An LDA explained-variance percentage describes its separating directions; it is not the same as the percentage of total data variation reported by PCA.

Where to be careful

• Class averages can hide a pattern. If A occupies two opposite corners of a square and B occupies the other two, their averages are the same. Fisher's mean-gap score is zero even though the classes form different shapes.
• There may not be enough data to estimate spread. With few rows, many measurements, or duplicate columns, the within-class scatter can be singular or unstable. Regularization or a solver that handles low-rank data may help. A tiny denominator does not prove that the classes are reliably separated.
• Keep preprocessing inside training. Fit missing-value handling, feature changes, and regularization using the training data only. In exact arithmetic, unregularized LDA can adjust to a reversible rescaling of features. In practice, poor scaling can cause numerical problems, and regularization can change the result.

Optional implementation: two features, two classes

This optional C++ example fits the Fisher direction and uses it to score new points. It solves a 2 × 2 system, scales the direction to length 1, and puts score zero halfway between the class averages. It stops with an error if the scatter matrix cannot be solved; adding regularization would be a separate change.

#include <array>
#include <cmath>
#include <stdexcept>
#include <vector>
using Point = std::array<double, 2>;
struct FisherProjection { Point origin, direction; };

FisherProjection fit_fisher(const std::vector<Point>& X,
                            const std::vector<int>& y) {
    if (X.size() != y.size() || X.empty())
        throw std::invalid_argument("matching nonempty X and y required");
    std::array<std::size_t, 2> count{0, 0};
    for (std::size_t i = 0; i < X.size(); ++i) {
        if (y[i] != 0 && y[i] != 1)
            throw std::invalid_argument("labels must be 0 or 1");
        for (double v : X[i]) if (!std::isfinite(v))
            throw std::invalid_argument("finite features required");
        ++count[y[i]];
    }
    if (count[0] < 2 || count[1] < 2)
        throw std::invalid_argument("need at least two rows in each class");
    std::array<Point, 2> mean{};
    for (std::size_t i = 0; i < X.size(); ++i)
        for (int j = 0; j < 2; ++j) mean[y[i]][j] += X[i][j] / count[y[i]];
    double a = 0, b = 0, c = 0;
    for (std::size_t i = 0; i < X.size(); ++i) {
        double dx = X[i][0] - mean[y[i]][0], dy = X[i][1] - mean[y[i]][1];
        a += dx*dx; b += dx*dy; c += dy*dy;
    }
    if (!std::isfinite(a) || !std::isfinite(b) || !std::isfinite(c))
        throw std::overflow_error("scatter overflow");
    // Rescale the 2x2 system; a common scalar does not change its direction.
    double scale = std::fmax(a, c);
    if (!(scale > 0)) throw std::invalid_argument("zero within-class scatter");
    a /= scale; b /= scale; c /= scale;
    double det = a*c - b*b;
    if (!(det > 1e-12))
        throw std::invalid_argument("singular or ill-conditioned scatter");
    Point d{mean[1][0] - mean[0][0], mean[1][1] - mean[0][1]};
    Point u{(c*d[0] - b*d[1])/det, (a*d[1] - b*d[0])/det};
    double length = std::hypot(u[0], u[1]);
    if (!std::isfinite(length)) throw std::overflow_error("direction overflow");
    if (!(length > 0)) throw std::invalid_argument("identical class means");
    return {{mean[0][0]/2 + mean[1][0]/2, mean[0][1]/2 + mean[1][1]/2},
            {u[0]/length, u[1]/length}};
}

double transform(const FisherProjection& model, const Point& x) {
    for (double v : x) if (!std::isfinite(v))
        throw std::invalid_argument("finite new point required");
    double z = (x[0] - model.origin[0])*model.direction[0]
             + (x[1] - model.origin[1])*model.direction[1];
    if (!std::isfinite(z)) throw std::overflow_error("score overflow");
    return z;
}

References and further reading

• Fisher (1936): The Use of Multiple Measurements in Taxonomic Problems — the original separation criterion and the linear equations for its maximizing coefficients.
• Scikit-learn: Linear and Quadratic Discriminant Analysis — the probabilistic classifier, dimensionality reduction and covariance regularization.
• Scikit-learn: LinearDiscriminantAnalysis API — transform dimensions, fitted attributes and solver options.