Principal Component Analysis

Principal Component Analysis (PCA) is an orthogonal linear transformation and unsupervised dimensionality-reduction technique in computation and algorithms and multivariate statistics that transforms a set of correlated variables into an ordered set of linearly uncorrelated variables called principal components. Geometrically, PCA rotates the coordinate system of a p-dimensional feature space so that the first principal axis aligns with the direction of greatest sample variance, and each subsequent axis maximizes residual variance subject to being orthogonal to all preceding axes. Originating in two dual formulations—Karl Pearson's 1901 derivation of best-fitting orthogonal subspaces minimizing squared reconstruction error, and Harold Hotelling's 1933 algebraic derivation maximizing projected variance—PCA serves as a primary analytical engine across population genetics, single-cell transcriptomics, climate teleconnection analysis, and conformational dynamics.

Data Centering and the Sample Covariance Matrix

Consider an n × p data matrix X representing n observations across p numerical features. Prior to transformation, each column j ∈ {1, ..., p} is mean-centered by subtracting its empirical sample mean μ_j = (1/n) ∑_i=1ⁿ X_ij, yielding the centered matrix X_c whose columns have zero mean. When features are measured in incommensurable physical units (for example, temperature in kelvins alongside pressure in pascals), each column is additionally divided by its sample standard deviation s_j to produce standardized z-scores (correlation-matrix PCA).

The p × p unbiased sample covariance matrix S is given by:

S = [1 / (n − 1)] X_c^T X_c

Because S is real, symmetric (S^T = S), and positive semi-definite (w^TSw ≥ 0 for all w ∈ R^p), the spectral theorem ensures that S admits a complete set of p real, non-negative eigenvalues λ₁ ≥ λ₂ ≥ ... ≥ λ_p ≥ 0 associated with mutually orthonormal eigenvectors v₁, v₂, ..., v_p satisfying v_j^Tv_k = δ_jk.

Hotelling's Variance Maximization vs. Pearson's Reconstruction Error

Hotelling (1933) sought a unit loading vector v₁ ∈ R^p (‖v₁‖₂ = 1) such that the linear projection (score vector) z₁ = X_cv₁ has maximum empirical variance:

Var(z₁) = [1 / (n − 1)] v₁^T X_c^T X_c v₁ = v₁^T S v₁

Maximizing the quadratic form v₁^TS v₁ subject to the unit-norm constraint v₁^Tv₁ = 1 via the Lagrangian L(v₁, λ₁) = v₁^TS v₁ − λ₁(v₁^Tv₁ − 1) and setting ∂L/∂v₁ = 0 yields the eigenvalue equation:

S v₁ = λ₁ v₁

Multiplying on the left by v₁^T shows that Var(z₁) = λ₁. By the Courant–Fischer min-max theorem for the Rayleigh quotient, the maximum variance equals the largest eigenvalue λ₁, achieved when v₁ is the leading eigenvector. Inductively, the k-th principal component loading vector v_k maximizes v_k^TS v_k subject to ‖v_k‖₂ = 1 and orthogonality v_k^Tv_j = 0 for all j < k, yielding the k-th eigenvector with variance λ_k. Furthermore, distinct principal component scores are mutually uncorrelated: Cov(z_j, z_k) = v_j^TS v_k = λ_k v_j^Tv_k = 0 for j ≠ k.

In Pearson's (1901) dual geometric formulation, one seeks a K-dimensional affine subspace (K < p) spanned by orthonormal basis vectors V_K = [v₁, ..., v_K] that minimizes the total squared orthogonal Euclidean distance (Frobenius norm reconstruction error) between the centered data points and their orthogonal projections X̂_K = X_cV_KV_K^T:

‖X_c − X_cV_KV_K^T‖_F² = (n − 1) [∑_j=1^p λ_j − ∑_j=1^K v_j^T S v_j] = (n − 1) ∑_j=K+1^p λ_j

By the Pythagorean theorem in Euclidean space, minimizing the discarded residual variance ∑_j=K+1^p λ_j is algebraically identical to maximizing the retained projected variance ∑_j=1^K λ_j.

Equivalence to SVD and the Eckart–Young–Mirsky Theorem

In numerical linear algebra, explicitly forming the covariance matrix X_c^TX_c squares the condition number κ(X_c)², which degrades floating-point precision for ill-conditioned datasets. Instead, numerical implementations compute the Singular Value Decomposition (SVD) of the centered n × p matrix X_c directly:

X_c = U Σ V^T = ∑_k=1^r σ_k u_k v_k^T

where r = rank(X_c) ≤ min(n, p), U is an n × r matrix of orthonormal left singular vectors (U^TU = I_r), V is a p × r matrix of orthonormal right singular vectors (V^TV = I_r), and Σ = diag(σ₁, σ₂, ..., σ_r) contains the ordered non-negative singular values σ₁ ≥ σ₂ ≥ ... ≥ σ_r > 0. Substituting the SVD into the covariance definition reveals:

S = [1 / (n − 1)] V Σ² V^T

Thus, the right singular vectors v_k are the principal component loadings (eigenvectors of S), the principal component scores are Z = X_cV = UΣ, and the covariance eigenvalues are related to the singular values by:

λ_k = σ_k² / (n − 1)

By the Eckart–Young–Mirsky theorem (1936), truncating the singular value expansion after K < r terms yields X_K = ∑_k=1^K σ_ku_kv_k^T, which is the best rank-K approximation to X_c in both the spectral (L₂ operator) norm and the Frobenius norm. When p ≫ n (such as transcriptomic or genomic matrices with p ≈ 10⁵ features and n ≈ 10³ samples), computing the n × n Gram matrix K = X_cX_c^T or applying Randomized SVD (Halko, Martinsson, & Tropp, 2011) reduces time complexity from O(np²) to O(npK).

Proportion of Variance Explained and Marchenko–Pastur / Tracy–Widom Bounds

Because the trace of a matrix is invariant under orthogonal similarity transformations, total sample variance satisfies ∑_j=1^p Var(X_j) = Tr(S) = ∑_j=1^p λ_j. The Proportion of Variance Explained (PVE) by the k-th component is PVE_k = λ_k / ∑_j=1^p λ_j.

In high-dimensional settings where p and n grow proportionally (p/n → γ > 0), random matrix theory proves that even for pure isotropic Gaussian noise (S_true = I_p), sample eigenvalues spread across the Marchenko–Pastur bulk distribution supported on [(1 − √γ)², (1 + √γ)²], and the largest noise eigenvalue fluctuates according to the Tracy–Widom distribution (TW₁) of order n^−2/3. Only sample eigenvalues exceeding the Baik–Ben Arous–Péché phase-transition threshold λ_signal > √γ emerge from the random noise bulk as statistically significant components.

Population Genetics and Geographic Clines (EIGENSTRAT)

In human and evolutionary genomics, Nick Patterson, Alkes L. Price, and David Reich (2006) formalized the application of PCA to biallelic Single Nucleotide Polymorphism (SNP) genotype matrices (X_ij ∈ {0, 1, 2} copies of the derived allele) along with Tracy–Widom significance testing (EIGENSTRAT). Under spatial isolation-by-distance models (stepping-stone coalescent processes), the genetic covariance matrix acts as an empirical inverse of the spatial graph Laplacian governing gene flow across a geographic domain. Consequently, the top two principal components (PC1 and PC2) of genome-wide SNP variation reconstruct two-dimensional geographic maps of continental and sub-continental ancestry—such as the northwest-to-southeast and north-to-south genetic clines of Europe—and are routinely included as covariates in Genome-Wide Association Studies (GWAS) to prevent false-positive associations caused by population stratification.

Single-Cell Transcriptomics and Biomolecular Essential Dynamics

In single-cell RNA sequencing (scRNA-seq), raw UMI count matrices spanning p ≈ 20,000–30,000 genes across n ≈ 10⁴–10⁶ cells are log-normalized, filtered to the top 2,000–5,000 variable genes, and projected onto their top 30–50 principal components via truncated SVD. This low-rank PCA subspace denoises technical dropout variance and serves as the metric input for k-nearest-neighbor graph construction, Leiden community clustering, and nonlinear manifold visualizations (such as UMAP, t-SNE, and diffusion maps based on the graph Laplacian).

In computational structural biology, applying PCA to the 3N × 3N mass-weighted Cartesian coordinate covariance matrix of α-carbon atoms harvested from molecular dynamics trajectories—after removing overall translational and rotational motion via Kabsch least-squares superposition—is known as Essential Dynamics (Amadei, Linssen, & Berendsen, 1993) or Quasi-Harmonic Analysis. Whereas high-index components with tiny eigenvalues correspond to fast, local harmonic bond vibrations, the top 2 to 5 principal modes capture collective, large-amplitude hinge-bending and domain-closure transitions that govern allosteric enzyme catalysis.

Why is it necessary to standardize features before running PCA?

PCA maximizes variance in the numerical units of the input columns. If one variable is measured in meters (variance ≈ 1) and an identical physical measurement is recorded in millimeters (variance ≈ 10⁶), covariance-matrix PCA will align PC1 almost parallel to the millimeter column simply because of the unit scale factor. Standardizing each feature to unit variance (z-scores, equivalent to eigendecomposition of the Pearson correlation matrix R) gives every variable equal prior weight.

Are principal components guaranteed to be statistically independent?

No. PCA diagonalizes the second-order covariance matrix, which ensures that principal component scores have zero linear correlation (Cov(z_j, z_k) = 0 for j ≠ k). Statistical independence equals zero correlation only when the underlying data follow a multivariate Gaussian distribution. Separating non-Gaussian higher-order dependencies (such as unmixing overlapping audio or electroencephalography signals) requires Independent Component Analysis (ICA), which optimizes fourth-order kurtosis or negentropy.

How does Linear PCA differ from Kernel PCA?

Standard PCA is restricted to orthogonal linear projections in the original p-dimensional Euclidean space and cannot unfold curved or concentric manifolds. Kernel PCA (Schölkopf, Smola, & Müller, 1998) maps data implicitly into a high-dimensional reproducing kernel Hilbert space via a positive-definite kernel function k(x_i, x_j) (such as a Gaussian radial basis function) and diagonalizes the centered n × n kernel Gram matrix, extracting nonlinear principal manifolds.

1. Pearson, K. (1901). On lines and planes of closest fit to systems of points in space. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 2(11), 559–572. https://doi.org/10.1080/14786440109462720
2. Hotelling, H. (1933). Analysis of a complex of statistical variables into principal components. Journal of Educational Psychology, 24(6), 417–441. https://doi.org/10.1037/h0071325
3. Eckart, C., & Young, G. (1936). The approximation of one matrix by another of lower rank. Psychometrika, 1(3), 211–218. https://doi.org/10.1007/BF02288367
4. Patterson, N., Price, A. L., & Reich, D. (2006). Population structure and eigenanalysis. PLoS Genetics, 2(12), e190. https://doi.org/10.1371/journal.pgen.0020190
5. Amadei, A., Linssen, A. B. M., & Berendsen, H. J. C. (1993). Essential dynamics of proteins. Proteins: Structure, Function, and Bioinformatics, 17(4), 412–425. https://doi.org/10.1002/prot.340170408