Covariance Matrices & Eigen Decomposition — Covariance matrices, correlation, and eigen decomposition - Statistics Graduate Level

SECTION 1

Historical Context & Motivation

The study of how multiple random variables move together is one of the oldest and most consequential problems in statistics. Long before the age of machine learning or big data, scientists needed to understand the joint behavior of measurements—crop yields and rainfall, astronomical observations across multiple instruments, or anthropometric measurements of human bodies. The covariance matrix emerged as the central algebraic object encoding these pairwise relationships, while eigen decomposition provided the mathematical key to unlocking its geometric meaning. Together, these tools form the backbone of multivariate statistics, principal component analysis, factor analysis, discriminant analysis, and a host of modern dimensionality-reduction techniques.

1846 • Bravais and Bivariate Correlation

Auguste Bravais formalized the mathematical treatment of the bivariate normal distribution, establishing early foundations for measuring linear association between two variables.

1896 • Pearson's Correlation Coefficient

Karl Pearson published the product-moment correlation coefficient, providing a standardized measure of linear dependence that would become the building block of the correlation matrix.

1901 • Pearson's Line of Best Fit

Pearson introduced the concept of fitting lines and planes of closest fit to systems of points in space, foreshadowing principal component analysis and the geometric interpretation of eigen decomposition.

1933 • Hotelling and PCA

Harold Hotelling formalized principal component analysis by showing that the eigenvectors of the covariance matrix define orthogonal axes of maximal variance, completing the bridge between linear algebra and multivariate statistics.

1960s–present • Computational Era

Efficient numerical algorithms such as QR iteration and the singular value decomposition made large-scale eigen decomposition practical, enabling modern applications in genomics, finance, signal processing, and deep learning.

The fundamental question that drove all of this development remains strikingly relevant today: given a cloud of data points in high-dimensional space, what are the natural axes along which the data varies most, and how can we separate meaningful structure from noise? Answering this question requires us to understand the covariance matrix as a geometric object and to deploy eigen decomposition as the tool that reveals its principal directions.

SECTION 2

Core Principles & Definitions

Before diving into computation, it is essential to establish the foundational objects and their interpretations. We work with a random vector X = (X₁, X₂, …, Xₚ)ᵀ taking values in ℝᵖ. The covariance matrix captures all second-order information about the joint distribution of these p variables, and eigen decomposition reveals the intrinsic coordinate system aligned with the data's variability.

Covariance Matrix (Σ)

The p × p symmetric positive semi-definite matrix whose (i, j)-entry is Cov(Xᵢ, Xⱼ) = E[(Xᵢ − μᵢ)(Xⱼ − μⱼ)]. Diagonal entries are variances; off-diagonal entries are covariances.

Correlation Matrix (R)

The standardized version of Σ obtained by dividing each entry by the product of the corresponding standard deviations: Rᵢⱼ = Σᵢⱼ / (σᵢ σⱼ). All diagonal entries equal 1, and off-diagonal entries lie in [−1, 1].

Eigen Decomposition

The factorization Σ = QΛQᵀ where Q is an orthogonal matrix of eigenvectors and Λ = diag(λ₁, …, λₚ) contains the eigenvalues. Because Σ is real symmetric, this decomposition always exists and is real-valued.

Positive Semi-Definiteness

For any non-zero vector a ∈ ℝᵖ, aᵀΣa ≥ 0. This guarantees all eigenvalues are non-negative and ensures the covariance matrix defines a valid ellipsoidal geometry in data space.

Spectral Theorem

Every real symmetric matrix admits an orthonormal basis of eigenvectors with real eigenvalues. This theorem guarantees that the eigen decomposition of any covariance matrix yields geometrically interpretable principal axes.

✦

KEY TAKEAWAY

Think of a covariance matrix as an instruction manual for drawing an ellipse (or ellipsoid in higher dimensions) around your data. The eigenvectors tell you which direction each axis of the ellipse points, and the eigenvalues tell you how long each axis is. If you were stretching a rubber sheet to cover your data cloud, the eigenvectors are the directions you stretch, and the square roots of the eigenvalues are the stretch factors.

SECTION 3

Visual Explanation — The Covariance Ellipse

The most illuminating way to understand what a covariance matrix encodes is to visualize it as a covariance ellipse in two dimensions. A bivariate normal distribution with mean μ and covariance matrix Σ has contours of constant probability density that form ellipses. The eigen decomposition of Σ directly determines the orientation and shape of these ellipses. The following diagram illustrates a bivariate data cloud, its covariance ellipse, and the eigenvector directions.

The pink arrow (e₁) points along the direction of maximum variance (first eigenvector), while the cyan arrow (e₂) is orthogonal and corresponds to the direction of minimum variance (second eigenvector). The dashed ellipse is the 1-standard-deviation contour of the bivariate normal, tilted at the angle determined by the eigenvectors. The amber dot marks the sample mean.

Notice that the data cloud is not aligned with the original X₁ and X₂ axes but rather with the eigenvector directions. This is precisely the insight behind principal component analysis: the eigenvectors of the covariance matrix define a rotated coordinate system in which the variables are uncorrelated. In this new coordinate system, the variance along each axis is given by the corresponding eigenvalue, and the off-diagonal covariances vanish. The length of each semi-axis of the ellipse is √λ, making large eigenvalues correspond to directions of high spread and small eigenvalues to directions where the data is tightly concentrated.

SECTION 4

Mathematical Framework

We now formalize the three central objects—covariance matrix, correlation matrix, and eigen decomposition—and establish the key relationships among them. Let X = (X₁, …, Xₚ)ᵀ be a random vector with mean vector μ = E[X] and suppose that Var(Xᵢ) < ∞ for all i.

COVARIANCE MATRIX

Σ = E[(X − μ)(X − μ)ᵀ] ⟹ Σᵢⱼ = Cov(Xᵢ, Xⱼ) = E[(Xᵢ − μᵢ)(Xⱼ − μⱼ)]

Σ is a p × p real symmetric matrix. Diagonal entries Σᵢᵢ = Var(Xᵢ) are the marginal variances. The matrix is positive semi-definite: for all a ∈ ℝᵖ, aᵀΣa = Var(aᵀX) ≥ 0.

SAMPLE COVARIANCE MATRIX

S = (1/(n − 1)) Σᵢ₌₁ⁿ (xᵢ − x̄)(xᵢ − x̄)ᵀ = (1/(n − 1)) X̃ᵀX̃

Here xᵢ ∈ ℝᵖ are observed data vectors, x̄ is the sample mean, and X̃ is the (n × p) mean-centered data matrix. S is the unbiased estimator of Σ. The factor 1/(n − 1) applies Bessel's correction.

CORRELATION MATRIX

R = D⁻¹ΣD⁻¹ where D = diag(σ₁, σ₂, …, σₚ)

D is the diagonal matrix of standard deviations. The (i, j)-entry of R is ρᵢⱼ = Σᵢⱼ/(σᵢσⱼ) ∈ [−1, 1]. Note that R is also positive semi-definite and symmetric, so it too admits an eigen decomposition.

SPECTRAL (EIGEN) DECOMPOSITION

Σ = QΛQᵀ = Σₖ₌₁ᵖ λₖ qₖ qₖᵀ

Q = [q₁ | q₂ | … | qₚ] is orthogonal (QᵀQ = Iₚ), with columns qₖ being the eigenvectors. Λ = diag(λ₁, …, λₚ) with λ₁ ≥ λ₂ ≥ … ≥ λₚ ≥ 0. The outer-product form Σₖ λₖqₖqₖᵀ shows that Σ is a weighted sum of rank-1 projection matrices, each weighted by the corresponding eigenvalue.

ℹ

Key Relationship

The total variance of the system is trace(Σ) = Σᵢ Var(Xᵢ) = Σₖ λₖ. Therefore the fraction of total variance captured by the k-th principal direction is λₖ / trace(Σ). This identity is the foundation for choosing how many principal components to retain in dimensionality reduction.

SECTION 5

Detailed Breakdown — From Covariance to Principal Axes

The passage from a raw covariance matrix to its geometric interpretation involves several conceptual steps. First, we recognize that Σ encodes an inner product on ℝᵖ via the Mahalanobis metric: for any two vectors u and v, the form uᵀΣ⁻¹v measures their similarity after accounting for the covariance structure. Second, the eigen decomposition diagonalizes this metric, yielding an orthonormal basis in which the metric becomes simply Euclidean but with scaled axes. Third, the connection between covariance and correlation matrices clarifies when to use each: eigen decomposition of Σ is appropriate when the original units are meaningful and comparable; eigen decomposition of R is appropriate when variables are on different scales and we want to analyze standardized structure.

The spectral decomposition pipeline: raw data is mean-centered, then the covariance matrix is computed, and finally eigen decomposition extracts the principal directions (eigenvectors) and variance magnitudes (eigenvalues). The total variance is preserved through the decomposition.

An important subtlety arises when variables are measured on vastly different scales—for instance, height in centimeters and weight in kilograms. In such cases, the covariance matrix will be dominated by the variable with the largest variance (height, due to its larger numerical range), and the first eigenvector will essentially point in the direction of that variable. Standardizing to the correlation matrix places all variables on an equal footing. However, if the original units are meaningful and comparable (e.g., returns on different stocks measured in the same currency), working with the covariance matrix preserves the information about relative magnitudes.

Feature	Covariance Matrix (Σ)	Correlation Matrix (R)
Diagonal entries	Variances (σᵢ²)	All equal to 1
Off-diagonal range	(−∞, +∞)	[−1, +1]
Scale dependence	Yes — changes under rescaling	No — invariant to rescaling
When to use for PCA	Variables in same/comparable units	Variables in different units or vastly different scales
Relationship	R = D⁻¹ΣD⁻¹	Σ = DRD

Comparison of covariance and correlation matrices

SECTION 6

Worked Example — Eigen Decomposition of a 2 × 2 Covariance Matrix

Consider a dataset of n = 100 observations on two variables. After centering, the sample covariance matrix is:

GIVEN COVARIANCE MATRIX

S = [ 5 3 ; 3 2 ]

Var(X₁) = 5, Var(X₂) = 2, Cov(X₁, X₂) = 3. Total variance = trace(S) = 7.

Finding Eigenvalues and Eigenvectors of S

1
Step 1 — Set Up the Characteristic Equation

We solve det(S − λI) = 0. The characteristic polynomial is:(5 − λ)(2 − λ) − 9 = 0. Expanding: λ² − 7λ + 10 − 9 = λ² − 7λ + 1 = 0.

Characteristic equation: λ² − 7λ + 1 = 0
2
Step 2 — Solve for Eigenvalues

Applying the quadratic formula: λ = (7 ± √(49 − 4)) / 2 = (7 ± √45) / 2 = (7 ± 3√5) / 2. Numerically, √5 ≈ 2.2361.

λ₁ ≈ 6.854, λ₂ ≈ 0.146
3
Step 3 — Verify Total Variance

λ₁ + λ₂ = 6.854 + 0.146 = 7.000 = trace(S) ✓. The first principal component captures λ₁/7 ≈ 97.9% of the total variance, confirming that the data is essentially one-dimensional.

97.9% of variance in first component
4
Step 4 — Find the First Eigenvector

For λ₁ ≈ 6.854, solve (S − λ₁I)q₁ = 0. The system becomes (5 − 6.854)q₁₁ + 3q₁₂ = 0, so −1.854q₁₁ + 3q₁₂ = 0, giving q₁₂/q₁₁ = 1.854/3 = 0.618. Normalizing: q₁ = (0.851, 0.526)ᵀ (approximately).

q₁ ≈ (0.851, 0.526)ᵀ
5
Step 5 — Find the Second Eigenvector

By orthogonality, q₂ must satisfy q₁ᵀq₂ = 0. Since q₁ ≈ (0.851, 0.526)ᵀ, we have q₂ ≈ (−0.526, 0.851)ᵀ. We can verify: (S − λ₂I)q₂ = 0 holds with λ₂ ≈ 0.146.

q₂ ≈ (−0.526, 0.851)ᵀ
6
Step 6 — Reconstruct and Interpret

The decomposition is S = QΛQᵀ where Q ≈ [0.851, −0.526; 0.526, 0.851] and Λ = diag(6.854, 0.146). The correlation between X₁ and X₂ is r = 3/√(5 × 2) = 3/√10 ≈ 0.949, indicating a very strong positive linear relationship. The first eigenvector points roughly 31.7° from the X₁-axis (arctan(0.526/0.851)).

S = QΛQᵀ verified; correlation r ≈ 0.949

SECTION 7

Strengths, Limitations & Practical Considerations

Eigen decomposition of covariance matrices is a workhorse of multivariate statistics, but it carries assumptions and limitations that practitioners must respect. Understanding the boundaries of the method is as important as understanding the method itself.

Strengths	Limitations
Captures all second-order (linear) dependencies in one compact matrix; computationally efficient (O(p³) for dense matrices)	Only captures linear associations; completely blind to nonlinear relationships (e.g., quadratic, circular dependencies)
Spectral theorem guarantees real eigenvalues and orthogonal eigenvectors; decomposition always exists for real symmetric matrices	Highly sensitive to outliers; a single extreme observation can dramatically distort the covariance matrix and its eigenstructure
Provides an optimal low-rank approximation (Eckart–Young theorem); basis of PCA, the most widely used dimensionality reduction method	When p > n (high-dimensional setting), the sample covariance matrix is singular with at most n − 1 nonzero eigenvalues; regularization is needed
Eigenvalues directly quantify explained variance; provides a principled criterion for dimensionality reduction (scree plot, cumulative variance threshold)	Eigenvalues of the sample covariance are biased estimators of population eigenvalues; largest eigenvalues are overestimated and smallest are underestimated (Marchenko–Pastur law)
Framework extends naturally to kernel PCA, sparse PCA, and robust PCA for more complex settings	Results depend on whether Σ or R is decomposed; the choice is often ambiguous and can lead to different conclusions

Strengths and limitations of covariance-based eigen decomposition

⚙️

PRACTICAL WISDOM

Think of eigen decomposition of the covariance matrix like using a GPS to find the principal highways through a city. It does an excellent job revealing the major thoroughfares (linear trends), but it completely misses winding side streets (nonlinear structure) and can be thrown off by construction detours (outliers). In high-dimensional genomics or finance, where p ≫ n, the sample covariance is like a blurry map—you need regularization techniques (shrinkage estimators, random matrix corrections) to sharpen the picture before extracting meaningful eigenvectors.

SECTION 8

Connection to Advanced Theory

The eigen decomposition of covariance matrices sits at a crossroads connecting numerous areas of advanced statistics and machine learning. Understanding these connections transforms the topic from a self-contained technique into a gateway to much of modern data science.

Core Concept	Advanced Extension	Key Relationship
PCA via eigen decomposition of Σ	Singular Value Decomposition (SVD)	If X̃ = UDVᵀ is the SVD of the centered data matrix, then S = VD²Vᵀ/(n−1), so the right singular vectors V are the eigenvectors of S.
Linear covariance structure	Kernel PCA	Replacing the linear inner product with a kernel function κ(xᵢ, xⱼ) allows eigen decomposition in an implicit infinite-dimensional feature space, capturing nonlinear structure.
Sample covariance S	Shrinkage Estimation	Ledoit–Wolf shrinkage: Σ̂ = αI + (1−α)S regularizes eigenvalues toward a common value, reducing estimation error when p is large relative to n.
Eigen decomposition of Σ	Random Matrix Theory	The Marchenko–Pastur law describes the limiting distribution of eigenvalues of sample covariance matrices when p/n → γ ∈ (0, ∞), providing a null model for distinguishing signal from noise eigenvalues.
Full-rank eigen decomposition	Sparse PCA	Imposing ℓ₁ penalties on the eigenvectors produces sparse loadings for interpretability, at the cost of losing orthogonality and requiring iterative optimization (e.g., LASSO-type algorithms).

Connections from covariance eigen decomposition to advanced topics

Perhaps the most profound theoretical connection is between eigen decomposition and the Gaussian likelihood. For a multivariate normal distribution N(μ, Σ), the log-likelihood involves log det(Σ) = Σₖ log λₖ and the quadratic form (x − μ)ᵀΣ⁻¹(x − μ) = Σₖ (qₖᵀ(x − μ))²/λₖ. Both expressions decompose naturally in the eigenbasis, making maximum likelihood estimation, hypothesis testing (e.g., likelihood ratio tests for covariance structure), and Bayesian inference on covariance matrices intimately connected to the spectral properties of Σ. The Wishart distribution, which governs the sampling distribution of S under normality, has its own eigenvalue distribution (the joint distribution of ordered eigenvalues), providing the theoretical underpinning for inference about population eigenvalues from sample eigenvalues.

SECTION 9

Practice Problems

Question 1: CONCEPTUAL

Explain why the covariance matrix Σ is always positive semi-definite. What does this property imply about its eigenvalues, and what would it mean geometrically if one of the eigenvalues were exactly zero?

Show answer

For any vector a ∈ ℝᵖ, aᵀΣa = Var(aᵀX) ≥ 0 because variance is non-negative. Since all eigenvalues λₖ satisfy λₖ = qₖᵀΣqₖ for the corresponding eigenvector qₖ, every eigenvalue must be non-negative. If an eigenvalue is exactly zero, the data has zero variance in the corresponding eigenvector direction, meaning the data lies in a lower-dimensional affine subspace of ℝᵖ. Geometrically, the covariance ellipsoid degenerates into a lower-dimensional ellipsoid (e.g., an ellipse instead of an ellipsoid in 3D). This occurs, for instance, when one variable is a perfect linear combination of the others, creating a rank-deficient covariance matrix.

Question 2: BASIC CALCULATION

Given the covariance matrix S = [4, −2; −2, 1], find the eigenvalues and eigenvectors. Compute the correlation matrix R and verify that the correlation coefficient lies in [−1, 1].

Show answer

The characteristic equation is (4 − λ)(1 − λ) − 4 = 0, giving λ² − 5λ = 0, so λ₁ = 5 and λ₂ = 0. Since λ₂ = 0, the matrix is singular and the data is perfectly collinear. For λ₁ = 5: (4 − 5)q₁ − 2q₂ = 0 gives q₁ = −2q₂, yielding q₁ = (−2/√5, 1/√5)ᵀ ≈ (−0.894, 0.447)ᵀ. For λ₂ = 0: 4q₁ − 2q₂ = 0 gives q₂ = 2q₁, yielding q₂ = (1/√5, 2/√5)ᵀ ≈ (0.447, 0.894)ᵀ. The correlation is r = −2/√(4 × 1) = −1, confirming perfect negative linear dependence, consistent with the zero eigenvalue.

Question 3: INTERMEDIATE

A 3 × 3 covariance matrix has eigenvalues λ₁ = 12, λ₂ = 3, λ₃ = 0.5. (a) What fraction of total variance is explained by the first two principal components? (b) If we project onto the first two eigenvectors, what is the covariance matrix of the projected data? (c) Write the rank-2 approximation of Σ.

Show answer

(a) Total variance = 12 + 3 + 0.5 = 15.5. The first two components explain (12 + 3)/15.5 = 15/15.5 ≈ 96.8%. (b) When we project X onto the subspace spanned by q₁ and q₂, the projected scores Z = [q₁, q₂]ᵀX have covariance matrix diag(12, 3), a 2 × 2 diagonal matrix, because the eigenvectors decorrelate the data. (c) The rank-2 approximation is Σ₂ = λ₁q₁q₁ᵀ + λ₂q₂q₂ᵀ = 12q₁q₁ᵀ + 3q₂q₂ᵀ. By the Eckart–Young theorem, this is the closest rank-2 positive semi-definite matrix to Σ in the Frobenius norm, with approximation error ‖Σ − Σ₂‖_F = λ₃ = 0.5.

Question 4: APPLIED

A financial analyst computes the daily return covariance matrix for three assets and obtains S = [0.04, 0.02, 0.01; 0.02, 0.09, 0.03; 0.01, 0.03, 0.01]. She wants to find the minimum-variance portfolio weights w (where 1ᵀw = 1). Express the solution in terms of S⁻¹ and explain how the eigenvalues of S influence portfolio risk.

Show answer

The minimum-variance portfolio has weights w* = S⁻¹1 / (1ᵀS⁻¹1), where 1 is the vector of ones. The portfolio variance is σ²ₚ = 1/(1ᵀS⁻¹1). In the eigenbasis, S⁻¹ = Σₖ (1/λₖ)qₖqₖᵀ, so the portfolio variance is dominated by the smallest eigenvalues of S. If the smallest eigenvalue is very small (near-collinear returns), S⁻¹ amplifies noise in that direction, leading to extreme and unstable portfolio weights. This is why practitioners use shrinkage estimators (e.g., Ledoit–Wolf) or regularized covariance matrices in portfolio optimization—they inflate the small eigenvalues, preventing S⁻¹ from blowing up and producing more robust weight estimates.

Question 5: CRITICAL THINKING

Consider a setting where p = 500 variables and n = 100 observations. (a) What is the rank of the sample covariance matrix S, and how many eigenvalues are exactly zero? (b) Why does this create problems for techniques that require S⁻¹? (c) Describe two distinct approaches to address this issue, and explain how each modifies the eigenvalue spectrum of the estimated covariance matrix.

Show answer

(a) The rank of S is at most min(n − 1, p) = 99, so at least 401 eigenvalues are exactly zero. (b) S is singular and cannot be inverted; this makes Mahalanobis distance, discriminant analysis, and the Wishart likelihood undefined. Even the nonzero eigenvalues are poorly estimated: the Marchenko–Pastur law shows they spread out more than the population eigenvalues, with the largest overestimated and the smallest underestimated. (c) Approach 1: Ledoit–Wolf shrinkage estimates Σ̂ = αI + (1−α)S, which adds α to every eigenvalue, making all 500 strictly positive and shrinking the spectrum toward a common value. Approach 2: Truncated PCA retains only the top k eigenvectors (say k = 20), effectively setting the remaining 480 eigenvalues to zero and working in the reduced k-dimensional subspace. This avoids inverting small eigenvalues but discards potential signal in components k+1 through 99. A hybrid approach uses random matrix theory to identify the bulk noise eigenvalues (predicted by Marchenko–Pastur) and only retains eigenvalues that exceed the theoretical upper edge of the bulk distribution.

SUMMARY

Summary

The covariance matrix Σ encodes all pairwise linear dependencies among p random variables, with diagonal entries representing variances and off-diagonal entries representing covariances. The correlation matrix R = D⁻¹ΣD⁻¹ standardizes these associations to the [−1, 1] scale, making it the appropriate choice when variables have incomparable units. The key algebraic property—positive semi-definiteness—guarantees non-negative eigenvalues and a well-defined ellipsoidal geometry.

The eigen decomposition Σ = QΛQᵀ decomposes the covariance structure into eigenvectors (principal directions of variation) and eigenvalues (variance magnitudes along those directions). The total variance is preserved as trace(Σ) = Σλₖ, enabling principled dimensionality reduction by retaining only the top eigenvectors. This framework underpins PCA, connects to the SVD of the data matrix, and extends through shrinkage estimation and random matrix theory to high-dimensional settings where p ≫ n.