Generalized Eigenvalue Problems
What LDA and CCA actually solve — Av = λBv, not just Av = λv.
On this page
Beginner: section 1.7's eigenvalue problem, Av = λv, asks "which directions does A only stretch, never rotate?" The generalized eigenvalue problem, Av = λBv, asks a related but different question: "which directions does A stretch by exactly λ times as much as B stretches them?" — comparing two matrices' behavior against each other, instead of describing just one matrix alone.
Intermediate: this shows up constantly whenever a method wants to maximize one notion of "spread" relative to another. Linear Discriminant Analysis (LDA) wants to find the projection direction that maximizes between-class separation relative to within-class spread — exactly a generalized eigenvalue problem, with A the between-class scatter matrix and B the within-class scatter matrix.
Advanced: Canonical Correlation Analysis (CCA), which finds the most correlated linear combinations of two different sets of variables (used to relate two different views or modalities of the same data), is also a generalized eigenvalue problem underneath. Both LDA and CCA can technically be reduced to an ordinary eigenvalue problem by multiplying through by B⁻¹, but doing that explicitly is numerically worse than using a solver built specifically for the generalized case.
When B is the identity matrix, this is exactly the ordinary eigenvalue problem from section 1.7 — the generalized version is a strict superset, not a different topic.
When B is symmetric positive definite (the usual ML case), Cholesky-factor it (section 1.15) as B = LLᵀ. Substitute into the generalized problem and insert L⁻ᵀLᵀ = I in a useful spot:
Now define the change of variables y = Lᵀv, so v = L⁻ᵀy, and substitute on the left:
This is now an ordinary eigenvalue problem (section 1.7) for the matrix C = L⁻¹AL⁻ᵀ, with the exact same eigenvalues λ as the original generalized problem — and if A is symmetric, C is symmetric too, so a standard symmetric eigensolver applies directly. The original eigenvectors are recovered via v = L⁻ᵀy. This is precisely the reduction production solvers use internally, which is why calling a dedicated generalized eigensolver is both correct and no more expensive than doing this transformation by hand.
Where this is used: this exact Cholesky-based reduction is what LDA and CCA solvers do internally, and it's also the standard way vibration analysis software solves Kv = λMv for a structure's natural frequencies and mode shapes.
Rotate the projection line — the Fisher score (between-class spread over within-class spread) is exactly what the generalized eigenvalue problem solves for directly.
scipy.linalg.eigh(A, B) solves the generalized problem directly — this is the real, standard way LDA is implemented, not by explicitly forming B⁻¹A. The from-scratch tabs implement exactly the Cholesky-based reduction derived above, by hand.
#include <cmath>
#include <cstdio>
#include <vector>
#include <random>
#include <array>
using Mat2 = std::array<std::array<double, 2>, 2>;
using Vec2 = std::array<double, 2>;
Mat2 matmul2(const Mat2& X, const Mat2& Y) {
Mat2 R{};
for (int i = 0; i < 2; ++i)
for (int j = 0; j < 2; ++j)
for (int k = 0; k < 2; ++k) R[i][j] += X[i][k] * Y[k][j];
return R;
}
Mat2 transpose2(const Mat2& X) {
return {{{X[0][0], X[1][0]}, {X[0][1], X[1][1]}}};
}
Mat2 cholesky2(const Mat2& M) {
double l00 = std::sqrt(M[0][0]);
double l10 = M[1][0] / l00;
double l11 = std::sqrt(M[1][1] - l10 * l10);
return {{{l00, 0.0}, {l10, l11}}};
}
Mat2 invertLowerTriangular(const Mat2& L) {
double l00 = L[0][0], l10 = L[1][0], l11 = L[1][1];
double i00 = 1.0 / l00, i11 = 1.0 / l11;
double i10 = -l10 * i00 * i11;
return {{{i00, 0.0}, {i10, i11}}};
}
// Eigenvalues/eigenvectors of a symmetric 2x2 matrix, closed form.
void eigSymmetric2x2(const Mat2& M, double lambda[2], Vec2 vec[2]) {
double a = M[0][0], b = M[0][1], d = M[1][1];
double tr = a + d, det = a * d - b * b;
double disc = std::sqrt(std::max(tr * tr - 4 * det, 0.0));
lambda[0] = (tr + disc) / 2;
lambda[1] = (tr - disc) / 2;
for (int k = 0; k < 2; ++k) {
Vec2 v = std::fabs(b) > 1e-12 ? Vec2{b, lambda[k] - a} : Vec2{1.0, 0.0};
double norm = std::sqrt(v[0] * v[0] + v[1] * v[1]);
vec[k] = {v[0] / norm, v[1] / norm};
}
}
int main() {
std::mt19937 rng(0);
std::normal_distribution<double> gauss(0.0, 1.0);
auto makeClass = [&](int n, double mx, double my, double corr) {
std::vector<Vec2> pts;
for (int i = 0; i < n; ++i) {
double x = gauss(rng) * 3.0;
double y = corr * x + gauss(rng);
pts.push_back({x + mx, y + my});
}
return pts;
};
auto classA = makeClass(50, -3.0, 0.0, 1.0 / 3.0);
auto classB = makeClass(50, 3.0, 0.0, 1.0 / 3.0);
auto meanVec = [](const std::vector<Vec2>& pts) {
Vec2 m{0.0, 0.0};
for (auto& p : pts) { m[0] += p[0]; m[1] += p[1]; }
m[0] /= pts.size(); m[1] /= pts.size();
return m;
};
auto covariance = [](const std::vector<Vec2>& pts, const Vec2& mean) {
Mat2 cov{};
for (auto& p : pts) {
double d0 = p[0] - mean[0], d1 = p[1] - mean[1];
cov[0][0] += d0 * d0; cov[0][1] += d0 * d1;
cov[1][0] += d1 * d0; cov[1][1] += d1 * d1;
}
double n = static_cast<double>(pts.size()) - 1.0;
for (auto& row : cov) for (auto& v : row) v /= n;
return cov;
};
Vec2 meanA = meanVec(classA), meanB = meanVec(classB);
Vec2 meanDiff = {meanA[0] - meanB[0], meanA[1] - meanB[1]};
Mat2 A = {{{meanDiff[0] * meanDiff[0], meanDiff[0] * meanDiff[1]},
{meanDiff[1] * meanDiff[0], meanDiff[1] * meanDiff[1]}}};
Mat2 covA = covariance(classA, meanA), covB = covariance(classB, meanB);
Mat2 B{};
for (int i = 0; i < 2; ++i)
for (int j = 0; j < 2; ++j) B[i][j] = covA[i][j] + covB[i][j];
Mat2 L = cholesky2(B);
Mat2 Linv = invertLowerTriangular(L);
Mat2 LinvT = transpose2(Linv);
Mat2 C = matmul2(matmul2(Linv, A), LinvT);
double lambda[2];
Vec2 y[2];
eigSymmetric2x2(C, lambda, y);
int best = lambda[0] > lambda[1] ? 0 : 1;
Vec2 v = {LinvT[0][0] * y[best][0] + LinvT[0][1] * y[best][1],
LinvT[1][0] * y[best][0] + LinvT[1][1] * y[best][1]};
double norm = std::sqrt(v[0] * v[0] + v[1] * v[1]);
v[0] /= norm; v[1] /= norm;
std::printf("largest eigenvalue: %.6f\n", lambda[best]);
std::printf("best direction: [%.6f, %.6f]\n", v[0], v[1]);
return 0;
}- LDA is used both as a classifier and as a supervised dimensionality reduction technique — unlike PCA (section 1.7/1.10), it uses class labels to choose directions that separate categories, not just directions of maximum variance.
- CCA underlies multi-view learning — relating text and image embeddings of the same concept, or brain-imaging signals to stimulus features, by finding maximally correlated projections of each.
- Vibration/structural analysis (the mechanical engineering example from section 1.7) is actually a generalized eigenvalue problem in its full form,
Kv = λMv, relating stiffness (K) and mass (M) matrices.
- Explicitly computing
B⁻¹Aand then solving an ordinary eigenvalue problem — this works in theory but is markedly less numerically stable than a dedicated generalized eigensolver, especially when B is close to singular. - Forgetting that a generalized eigenvalue problem needs B to be invertible (or at least positive definite, section 1.15, for the well-behaved real-eigenvalue case) — an ill-conditioned within-class scatter matrix is a common practical failure mode of LDA on small or collinear datasets.
Going deeper
When B is symmetric positive definite (the common case in ML — scatter and covariance matrices are always PSD, section 1.15), the generalized eigenvalue problem has an elegant reduction: Cholesky-factor B = LLᵀ, then solve the ordinary eigenvalue problem for L⁻¹A(L⁻¹)ᵀ — which is exactly what production-grade generalized eigensolvers do internally, tying this lesson directly back to section 1.15.
At the master level: LDA assumes each class's within-class scatter is well-estimated, which fails badly in high dimensions with few samples per class — regularized/shrinkage LDA adds a small multiple of the identity to the within-class scatter matrix before solving, the exact same "ridge" trick from section 1.9 applied here to keep B safely invertible.