STATISTICS GRADUATE LEVEL • ADVANCED TOPICS IN GRADUATE STATISTICS

Principal Component Analysis (PCA)

A foundational dimensionality-reduction technique that reveals the latent structure of high-dimensional data through orthogonal linear transformations.

Historical Context & Motivation

The challenge of analyzing data with many correlated variables has confronted statisticians for well over a century. As experimental sciences matured, researchers found themselves collecting dozens or even hundreds of measurements on each observational unit, and the question naturally arose: is there a principled way to reduce this overwhelming complexity while preserving the essential patterns? Principal Component Analysis (PCA) emerged as the definitive answer to this question, providing a mathematically rigorous framework for identifying the directions of maximum variability in multivariate data. The method transforms a set of possibly correlated variables into a smaller set of uncorrelated variables — the principal components — that capture most of the information contained in the original measurements.

PCA's intellectual roots trace back to the geometry of multidimensional ellipsoids and the algebraic theory of quadratic forms, but its modern formulation crystallized through the work of several key figures across mathematics, psychology, and statistics. Understanding this history illuminates why PCA occupies such a central position in multivariate analysis and why its influence extends far beyond statistics into machine learning, signal processing, genomics, and finance.

1901
Pearson's Line of Closest Fit
Karl Pearson introduced the concept of fitting lines and planes of closest fit to a system of points in multidimensional space. His geometric formulation — minimizing the sum of squared perpendicular distances — laid the conceptual foundation for PCA, framing it as an optimization problem in Euclidean geometry.
1933
Hotelling's Formal Development
Harold Hotelling independently derived PCA from a statistical perspective, defining principal components as linear combinations of variables that successively maximize variance subject to orthogonality constraints. His algebraic treatment connected PCA to eigenvalue decomposition of the covariance matrix, establishing the method's modern mathematical framework.
1936
Eckart–Young Theorem
Carl Eckart and Gale Young proved that the best rank-k approximation of a matrix (in the Frobenius norm) is given by its truncated singular value decomposition. This result provided a deep theoretical justification for PCA as the optimal linear dimensionality reduction method.
1965
Computational Advances
The development of efficient numerical algorithms for eigenvalue decomposition and the singular value decomposition (SVD), combined with increasing computational power, made PCA practically feasible for large datasets. Gene Golub and colleagues refined algorithms that remain foundational to modern linear algebra software.
2000s
Modern Extensions
The explosion of high-dimensional data in genomics, imaging, and text analysis spurred extensions such as sparse PCA, kernel PCA, robust PCA, and randomized algorithms capable of handling matrices with millions of rows and columns. PCA became a standard preprocessing step in modern machine learning pipelines.

The central question PCA addresses can be stated concisely: given a data matrix with p correlated variables observed on n units, can we find a new coordinate system — a rotation of the original axes — in which the data's variance structure is maximally simplified? This question connects classical geometry, linear algebra, and statistical inference in an elegant synthesis that continues to motivate research more than a century after Pearson's original insight.

Core Principles & Definitions

PCA rests on a remarkably elegant set of foundational ideas that connect variance, linear algebra, and geometric intuition. Before diving into the mathematical machinery, it is essential to understand the conceptual pillars that support the entire framework. Each principle below addresses a specific aspect of how PCA transforms data, why that transformation is optimal, and what the resulting components represent.

1

Variance as Information

PCA operates on the premise that directions of greatest variance carry the most information. If a variable barely varies across observations, it contributes little to distinguishing among data points. PCA seeks new axes — linear combinations of the original variables — that sequentially maximize variance.
2

Orthogonality Constraint

Each successive principal component must be orthogonal (uncorrelated) to all preceding components. This constraint ensures that each component captures a unique, non-redundant source of variation. The result is a decorrelated representation of the original data.
3

Eigenvalue Decomposition

The principal components are the eigenvectors of the data's covariance (or correlation) matrix. The associated eigenvalues quantify the variance explained by each component. The spectral theorem for symmetric matrices guarantees a complete orthonormal basis of eigenvectors.
4

Dimensionality Reduction

By retaining only the first k components (where k ≪ p), we approximate the data in a lower-dimensional subspace that captures the majority of total variance. The proportion of variance explained provides a natural criterion for choosing k.
5

Optimal Reconstruction

Among all possible rank-k linear approximations, the PCA solution minimizes the mean squared reconstruction error (equivalently, maximizes preserved variance). This optimality is guaranteed by the Eckart–Young–Mirsky theorem.
KEY TAKEAWAY
Think of PCA as rotating a camera to find the viewing angle from which a cloud of data points looks most spread out. The first principal component is the direction of greatest elongation; the second is the direction of greatest remaining spread perpendicular to the first; and so on. Just as a skilled photographer chooses an angle that reveals the most structure in a landscape, PCA identifies the coordinate system that makes the data's variance structure most transparent — and crucially, it does so in the uniquely optimal way guaranteed by the spectral theorem.

Visual Explanation

The geometric intuition behind PCA is best conveyed through a visualization of a two-dimensional data cloud and its principal component axes. In the diagram below, we see how PCA identifies the directions of maximum variance in a scatter of data points and projects them onto a new coordinate system defined by these directions.

The data cloud (blue points) exhibits a strong positive correlation between X₁ and X₂. The dashed ellipse approximates the concentration of the data. PC1 (cyan arrow) aligns with the major axis of the ellipse — the direction of maximum variance — while PC2 (violet arrow) aligns with the minor axis, capturing the remaining variance orthogonally.

Several critical observations emerge from this diagram. First, the principal component axes are rotated relative to the original X₁ and X₂ axes by an angle that depends entirely on the covariance structure of the data. Second, PC1 captures the majority of the total variance — the data is far more spread out along PC1 than along PC2. If we were to project all points onto the PC1 axis alone, we would lose relatively little information, achieving a reduction from two dimensions to one. Third, the orthogonality of PC1 and PC2 ensures that the scores along these new axes are uncorrelated, simplifying subsequent analyses. This geometric picture generalizes directly to p dimensions: the first k principal components define the k-dimensional hyperplane that best approximates the data cloud in the least-squares sense.

Mathematical Framework

The formal derivation of PCA proceeds from a constrained optimization problem and arrives at an eigenvalue decomposition of the covariance matrix. Let X be an n × p centered data matrix (each column has mean zero), where n is the number of observations and p is the number of variables. The sample covariance matrix is S = (1/(n−1)) Xᵀ X, a p × p symmetric positive semi-definite matrix.

COVARIANCE MATRIX
S = (1 / (n − 1)) Xᵀ X
where X is the n × p centered data matrix, Xᵀ is its transpose, and S is the p × p sample covariance matrix.

The first principal component is the unit vector w₁ ∈ ℝᵖ that maximizes the variance of the projected data. The variance of the projection of X onto a unit vector w is wᵀ S w. We therefore seek to maximize this quadratic form subject to the constraint ‖w‖ = 1.

OPTIMIZATION PROBLEM
w₁ = arg max‖w‖=1 wᵀ S w
The first principal component direction maximizes the projected variance. Using a Lagrange multiplier λ, the stationarity condition yields Sw = λw, i.e., w must be an eigenvector of S.

Introducing a Lagrange multiplier to enforce the unit-norm constraint, we form the Lagrangian L(w, λ) = wᵀSw − λ(wᵀw − 1). Setting the gradient to zero gives Sw = λw, which is precisely the eigenvalue equation for S. The objective value at this solution is wᵀSw = λ, so to maximize variance, we choose the eigenvector corresponding to the largest eigenvalue λ₁. The k-th principal component is the eigenvector associated with the k-th largest eigenvalue, and it is automatically orthogonal to all preceding components because S is symmetric and the spectral theorem guarantees orthogonal eigenvectors.

SPECTRAL DECOMPOSITION
S = W Λ Wᵀ
where W = [w₁ | w₂ | … | wₚ] is the p × p orthogonal matrix of eigenvectors (loadings), and Λ = diag(λ₁, λ₂, …, λₚ) is the diagonal matrix of eigenvalues with λ₁ ≥ λ₂ ≥ … ≥ λₚ ≥ 0.
PROPORTION OF VARIANCE EXPLAINED
PVE(k) = Σᵢ₌₁ᵏ λᵢ / Σᵢ₌₁ᵖ λᵢ
The proportion of total variance captured by the first k principal components. This quantity is central to deciding how many components to retain: a common heuristic is to choose k such that PVE(k) ≥ 0.80 or 0.90, or to use the scree plot elbow criterion.
📐 SVD Connection
PCA can equivalently be computed via the Singular Value Decomposition (SVD) of the centered data matrix X = UDVᵀ. The columns of V are the principal component directions (loadings), D²/(n−1) gives the eigenvalues of S, and UD gives the principal component scores. In practice, the SVD is often preferred because it is numerically more stable than forming Xᵀ X explicitly.

Choosing the Number of Components

One of the most consequential decisions in any PCA application is selecting how many principal components to retain. Keeping too few components discards meaningful signal; keeping too many reintroduces noise and defeats the purpose of dimensionality reduction. Several complementary criteria have been developed, each with its own rationale and limitations. The scree plot is perhaps the most widely used visual tool: it plots the eigenvalues λ₁ ≥ λ₂ ≥ … ≥ λₚ in descending order, and the analyst looks for an "elbow" — a point after which the eigenvalues level off, suggesting that subsequent components capture mostly noise.

A scree plot for six components. The cyan line shows eigenvalues declining sharply through the first two components, then leveling off — the elbow at k = 2 (pink dashed line) suggests retaining two components. The amber dashed line shows cumulative variance explained reaching 73% at two components and 86% at three.

Formal Criteria for Component Selection

Summary of common criteria for selecting the number of principal components
CriterionRuleStrengths / Limitations
Kaiser's RuleRetain components with eigenvalue > 1 (when using the correlation matrix).Simple and widely used, but tends to over-retain in large p and under-retain in small p.
Scree Plot ElbowIdentify the "elbow" where eigenvalue decline levels off; retain components before the elbow.Intuitive but subjective; ambiguous when no clear elbow exists.
Cumulative Variance ThresholdRetain the fewest components capturing ≥ 80–90% of total variance.Practical and interpretable, but the threshold is arbitrary.
Parallel AnalysisCompare observed eigenvalues to those from random data of the same dimensions; retain components whose eigenvalues exceed the random threshold.Statistically principled and performs well empirically; requires simulation.
Cross-ValidationChoose k that minimizes out-of-sample reconstruction error via held-out data.Data-driven and avoids arbitrary thresholds; computationally expensive for large datasets.

In practice, analysts often combine multiple criteria rather than relying on any single rule. A conservative approach might retain components that satisfy both Kaiser's rule and account for at least 80% of variance, while also confirming the choice against parallel analysis. The appropriate strategy also depends on whether PCA is being used for exploratory purposes (where interpretability matters) or as a preprocessing step for a downstream model (where predictive performance on held-out data is the relevant metric).

Worked Example

Consider a dataset with n = 5 observations on p = 3 variables (X₁, X₂, X₃). We will walk through the complete PCA procedure: centering the data, computing the covariance matrix, finding eigenvalues and eigenvectors, and projecting the data onto the principal component space.

Raw data matrix (5 observations × 3 variables)
ObservationX₁X₂X₃
1426
28510
3648
410712
5224
PCA on a 5 × 3 Dataset
1
Step 1 — Center the DataCompute the column means: X̄₁ = (4+8+6+10+2)/5 = 6, X̄₂ = (2+5+4+7+2)/5 = 4, X̄₃ = (6+10+8+12+4)/5 = 8. Subtract the means from each observation to obtain the centered data matrix. For example, observation 1 becomes (4−6, 2−4, 6−8) = (−2, −2, −2).
Centered matrix: [−2,−2,−2 | 2,1,2 | 0,0,0 | 4,3,4 | −4,−2,−4]
2
Step 2 — Compute the Covariance Matrix SS = (1/4) Xᵀ X, where X is the centered matrix and n−1 = 4. Computing each entry: S₁₁ = (4+4+0+16+16)/4 = 10.0, S₂₂ = (4+1+0+9+4)/4 = 4.5, S₃₃ = (4+4+0+16+16)/4 = 10.0, S₁₂ = (4+2+0+12+8)/4 = 6.5, S₁₃ = (4+4+0+16+16)/4 = 10.0, S₂₃ = (4+2+0+12+8)/4 = 6.5.
S = [[10.0, 6.5, 10.0], [6.5, 4.5, 6.5], [10.0, 6.5, 10.0]]
3
Step 3 — Find EigenvaluesSolve det(S − λI) = 0. Notice that columns 1 and 3 of S are identical, which means S is singular (one eigenvalue is zero). The characteristic polynomial simplifies. Computing numerically: λ₁ ≈ 23.70, λ₂ ≈ 0.80, λ₃ = 0.00. The total variance is 10.0 + 4.5 + 10.0 = 24.50, and indeed λ₁ + λ₂ + λ₃ = 23.70 + 0.80 + 0.00 = 24.50.
λ₁ ≈ 23.70 (96.7%), λ₂ ≈ 0.80 (3.3%), λ₃ = 0.00 (0%)
4
Step 4 — Find Eigenvectors (Loadings)For λ₁ ≈ 23.70, solve (S − 23.70I)w = 0 and normalize to unit length. The first eigenvector is approximately w₁ ≈ (0.621, 0.406, 0.621)ᵀ, indicating that PC1 loads roughly equally on X₁ and X₃ with a moderate loading on X₂. This makes sense given the strong correlation structure. For λ₂ ≈ 0.80, w₂ ≈ (−0.345, 0.870, −0.345)ᵀ, which contrasts X₂ against X₁ and X₃.
w₁ ≈ (0.621, 0.406, 0.621)ᵀ, w₂ ≈ (−0.345, 0.870, −0.345)ᵀ
5
Step 5 — Project Data onto PC SpaceThe PC scores are obtained by Z = X_centered × W. For observation 1, the centered values are (−2, −2, −2), so z₁₁ = (−2)(0.621) + (−2)(0.406) + (−2)(0.621) = −3.296. Since PC1 alone explains 96.7% of total variance, a one-dimensional representation using only PC1 scores provides an excellent summary of the three-variable dataset.
PC1 captures 96.7% of total variance — a single component summarizes nearly all information in 3 variables.
🔍 Why λ₃ = 0?
Because X₃ = X₁ in every observation (after centering, their columns are identical), the three variables actually span only a two-dimensional subspace. PCA detects this perfect collinearity automatically: the zero eigenvalue signals that the third component carries no information. In practice, near-zero eigenvalues often indicate approximate collinearities worth investigating.

Strengths, Limitations, and Practical Considerations

PCA is one of the most versatile tools in the statistician's repertoire, but like any method, it has both notable strengths and important limitations. Understanding these is essential for deciding when PCA is appropriate and how to interpret its results. The table below provides a systematic comparison.

Strengths and limitations of PCA
StrengthsLimitations
Optimal linear dimensionality reduction (minimizes reconstruction error among all linear projections of the same rank).Captures only linear relationships; nonlinear structure in the data will be missed entirely.
Unsupervised — requires no class labels; works purely from the data's covariance structure.Components that maximize variance may not maximize class separability (use LDA if class discrimination is the goal).
Computationally efficient via SVD; scales to large datasets with randomized algorithms.Sensitive to the relative scaling of variables; results change dramatically if variables are measured in different units and not standardized.
Produces uncorrelated components, simplifying subsequent analysis (e.g., avoiding multicollinearity in regression).Loadings are dense — every original variable contributes to every component — making components difficult to interpret without further rotation or sparsification.
No tuning parameters beyond the number of components; the mathematics fully determines the solution.Sensitive to outliers: extreme observations can disproportionately influence the covariance matrix and hence the principal directions.

Covariance vs. Correlation Matrix

A perennial practical question is whether to perform PCA on the covariance matrix or the correlation matrix (equivalently, whether to standardize the variables to unit variance before computing PCA). When variables are measured in commensurate units and their natural variances are meaningful — for instance, returns on multiple stocks — the covariance matrix preserves scale information that may be substantively important. However, when variables have heterogeneous scales (e.g., height in centimeters and weight in kilograms), using the covariance matrix allows the highest-variance variable to dominate the first principal component purely because of its measurement unit. In such cases, standardizing to the correlation matrix is almost always preferable.

KEY TAKEAWAY
PCA is the Swiss Army knife of multivariate analysis — remarkably useful for exploration, visualization, noise reduction, and preprocessing. But remember that it is fundamentally a tool for capturing linear variance. If you suspect that the meaningful structure in your data lives along curved manifolds (as in many image and text datasets), consider nonlinear extensions such as kernel PCA or autoencoders. Similarly, if your goal is supervised prediction rather than unsupervised variance decomposition, methods like Partial Least Squares (PLS) or sufficient dimension reduction may be more appropriate.

Connection to Advanced Theory & Extensions

Classical PCA opens the door to a rich family of advanced methods, each addressing a specific limitation of the original technique. Understanding these extensions contextualizes PCA within the broader landscape of modern statistical and machine learning methods. The table below summarizes the key relationships.

Classical PCA and its modern extensions
ExtensionKey InnovationWhen to Use
Kernel PCAApplies PCA in a high-dimensional feature space induced by a kernel function (e.g., RBF), enabling capture of nonlinear structure.Data lies on or near a nonlinear manifold; linear PCA fails to separate meaningful clusters.
Sparse PCAAdds an L₁ (lasso) penalty to the loadings, producing components with many zero weights for interpretability.High-dimensional data (p >> n) where interpretable, feature-selecting components are required.
Robust PCADecomposes the data matrix as a sum of low-rank and sparse components, isolating the effect of outliers or corrupted entries.Data contaminated by gross outliers or missing values; video surveillance background/foreground separation.
Probabilistic PCAEmbeds PCA within a latent variable generative model: x = Wz + μ + ε, where z ~ N(0, I). Enables missing data imputation and Bayesian model selection.Need to handle missing data, perform model comparison, or connect PCA to the EM algorithm and factor analysis.
Functional PCAExtends PCA to infinite-dimensional function spaces, decomposing sample curves into principal modes of variation via eigenanalysis of the covariance operator.Data are curves, spectra, or time series observed over a continuous domain (e.g., growth curves, fMRI time series).

From a theoretical perspective, PCA can be understood as the maximum-likelihood solution under the probabilistic PCA model when the noise variance tends to zero (Tipping & Bishop, 1999). This connection reveals PCA as a special case of factor analysis and opens the door to Bayesian treatments where the number of components can be inferred from the data using automatic relevance determination. For high-dimensional settings where p >> n, the random matrix theory of Marchenko and Pastur provides critical insights into how sample eigenvalues behave under the null hypothesis of no structure, informing principled thresholding rules that outperform classical criteria like Kaiser's rule.

🚀 Looking Ahead
If you are continuing into nonlinear dimensionality reduction, note that t-SNE and UMAP have largely superseded PCA for visualization of high-dimensional data in two or three dimensions. However, PCA remains the standard first step in pipelines that require a linear, invertible, and variance-preserving transformation — for example, as a denoising preprocessor before feeding data into a neural network or Gaussian process.

Practice Problems

PROBLEM 1CONCEPTUAL
Explain why the principal components are necessarily uncorrelated with each other. In your answer, reference the properties of the covariance matrix S and the spectral theorem for symmetric matrices.
PROBLEM 2BASIC CALCULATION
A dataset of 4 variables yields eigenvalues λ₁ = 6.2, λ₂ = 2.8, λ₃ = 0.7, λ₄ = 0.3. Compute the proportion of variance explained by the first two principal components and determine how many components should be retained under Kaiser's rule (assuming the correlation matrix was used).
PROBLEM 3INTERMEDIATE
Suppose the covariance matrix of a bivariate dataset is S = [[4, 3], [3, 9]]. Find the eigenvalues and eigenvectors of S, and express the first principal component as a linear combination of the original variables X₁ and X₂.
PROBLEM 4APPLIED
A genomics researcher measures expression levels for p = 20,000 genes across n = 100 tumor samples. She plans to use PCA for visualization and preprocessing before a classification algorithm. Discuss (a) whether she should use the covariance or correlation matrix, (b) the computational strategy she should employ, and (c) how random matrix theory might inform her choice of k.
PROBLEM 5CRITICAL THINKING
Prove that the rank-k matrix that minimizes the Frobenius-norm distance to an arbitrary matrix A ∈ ℝⁿˣᵖ is given by the truncated SVD Aₖ = UₖDₖVₖᵀ, and explain how this result establishes the optimality of PCA for dimensionality reduction.

Summary

Principal Component Analysis (PCA) is a foundational technique for dimensionality reduction that finds the orthogonal linear transformation maximizing preserved variance. Its solution emerges from the eigenvalue decomposition of the covariance matrix S = WΛWᵀ, where the eigenvectors (loadings) define the new coordinate axes and the eigenvalues quantify the variance along each axis. The Eckart–Young–Mirsky theorem guarantees that PCA yields the optimal rank-k approximation in the Frobenius norm.

The number of retained components is determined by criteria such as the scree plot, cumulative variance threshold, Kaiser's rule, or parallel analysis. While PCA is limited to linear relationships and is sensitive to scaling, modern extensions — including kernel PCA, sparse PCA, and probabilistic PCA — address these limitations and connect classical PCA to the broader landscape of latent variable modeling and manifold learning.

Varsity Tutors • Statistics Graduate Level • Principal Component Analysis (PCA)