Principal Component Analysis

Let $X \in \mathbb{R}^{n \times p}$ be a data matrix of $n$ observations in $p$ dimensions, with rows $x_i^\top$. Assume the data is centered, meaning the column means have been subtracted so $\bar{x} = \frac{1}{n}\sum_i x_i = 0$. The empirical covariance matrix is

$S$ is symmetric and positive semi-definite. PCA seeks a unit vector $w_1 \in \mathbb{R}^p$, $\|w_1\|=1$, such that projecting each data point onto $w_1$ — computing the scalar score $z_i = x_i^\top w_1$ — maximizes the variance of the resulting scores $\{z_i\}$. This is the first principal component direction. Subsequent components $w_2, w_3, \ldots$ are chosen the same way, subject to being orthogonal to all previous directions, so that each new component captures the maximal remaining variance. The key terminology: $w_k$ is the $k$-th principal component (also called a loading vector or eigenvector); the projected value $z_i = x_i^\top w_k$ is the score; the variance captured, $\lambda_k$, is the eigenvalue associated with $w_k$; and $\lambda_k / \sum_j \lambda_j$ is the fraction of explained variance. Collecting the top $k$ components as columns of a matrix $W_k \in \mathbb{R}^{p \times k}$, the transformation $Z = X W_k$ maps each observation to a $k$-dimensional code, and $X W_k W_k^\top$ reconstructs an approximation back in the original $p$-dimensional space. Crucially, PCA is a linear, unsupervised, variance-based method — related to but distinct from methods built on Distance Metrics and Similarity Measures like multidimensional scaling, and it makes no use of labels, unlike supervised dimensionality reduction such as linear discriminant analysis. It assumes that variance is a good proxy for "information" and that the interesting structure is linear; when that assumption fails, nonlinear techniques discussed under Manifold Learning and Nonlinear Dimensionality Reduction become necessary.

Step 1 — Center (and optionally standardize) the data. Subtract the mean of each feature so the cloud of points is centered at the origin. If features are on very different scales (e.g., income in dollars vs. age in years), divide each feature by its standard deviation as well; this variant is PCA on the correlation matrix rather than the covariance matrix.

Step 2 — Compute the covariance matrix. $S = \frac{1}{n-1}X^\top X$, a $p \times p$ symmetric matrix whose $(j,k)$ entry is the sample covariance between feature $j$ and feature $k$.

Step 3 — Eigendecompose $S$. Since $S$ is symmetric positive semi-definite, the spectral theorem guarantees a full set of real, non-negative eigenvalues $\lambda_1 \ge \lambda_2 \ge \cdots \ge \lambda_p \ge 0$ with a corresponding orthonormal set of eigenvectors $w_1, \ldots, w_p$. These eigenvectors are exactly the principal component directions, and the eigenvalue $\lambda_k$ equals the variance of the data when projected onto $w_k$.

Step 4 — Truncate. Keep the top $k$ eigenvectors (largest eigenvalues) forming $W_k$, chosen either as a fixed number or by requiring cumulative explained variance $\sum_{i\le k}\lambda_i / \sum_i \lambda_i$ to exceed a threshold (e.g., 95%).

Step 5 — Project / reconstruct. Compute scores $Z = X W_k$ for compression or visualization, and reconstruct $\hat{X} = Z W_k^\top$ as a denoised low-rank approximation. Why does maximizing variance under orthogonality yield eigenvectors of $S$? Consider the first component: maximize $\mathrm{Var}(Xw) = w^\top S w$ subject to $\|w\|=1$. Using a Lagrange multiplier $\lambda$,

and setting the gradient to zero gives

So any stationary point $w$ must be an eigenvector of $S$, and the variance achieved there is $w^\top S w = w^\top(\lambda w) = \lambda$. To maximize variance, pick the eigenvector with the largest eigenvalue. For the second component, maximize $w^\top S w$ subject to $\|w\|=1$ and $w \perp w_1$; the same Lagrangian argument (with an added orthogonality multiplier) shows the optimum is the eigenvector with the next-largest eigenvalue, and inductively for all $k$. This is precisely why the eigenvectors of $S$, sorted by eigenvalue, solve the sequential variance-maximization problem exactly — no iterative search is needed in principle, though in practice numerical eigensolvers are themselves iterative. ``` ALGORITHM PCA(X, k) Input: X, an n × p data matrix (n samples, p features) k, target number of components (k <= p) Output: W_k (p × k matrix of principal directions) Z (n × k matrix of projected scores) explained_variance_ratio (length-k vector) // Step 1: center the data (mean-subtraction is mandatory) mu <- column_mean(X) // length-p vector Xc <- X - broadcast(mu) // n × p, mean-zero columns // (optional) standardize if feature scales differ wildly if standardize: sigma <- column_std(Xc) Xc <- Xc / broadcast(sigma) // Step 2: covariance matrix S <- (1 / (n - 1)) transpose(Xc) @ Xc // p × p, symmetric PSD // Step 3: eigendecomposition S = V diag(lambda) * V^T (lambda, V) <- eigendecompose(S) // lambda: length-p, V: p × p sort (lambda, V) by lambda descending // lambda_1 >= lambda_2 >= ... // Step 4: truncate to top k eigenvectors W_k <- V[:, 1..k] // p × k // Step 5: project Z <- Xc @ W_k // n × k total_var <- sum(lambda) explained_variance_ratio <- lambda[1..k] / total_var return W_k, Z, explained_variance_ratio

ALGORITHM RECONSTRUCT(Z, W_k, mu) // approximate original data from k-dimensional codes Xhat <- Z @ transpose(W_k) + broadcast(mu) return Xhat ```

Time complexity. Forming the covariance matrix $S = X^\top X / (n-1)$ costs $O(np^2)$ (a $p \times n$ times $n \times p$ matrix product). Eigendecomposing the resulting $p \times p$ symmetric matrix costs $O(p^3)$ using standard dense methods (e.g., the symmetric QR algorithm after tridiagonalization). Total time is

When $n \gg p$ (many samples, few features), the $np^2$ term dominates and PCA scales linearly in $n$. When $p$ is very large (e.g., $p \sim 10^4$ features, $n$ moderate), the $p^3$ eigendecomposition becomes the bottleneck, and it is standard to instead perform the Singular Value Decomposition of the centered data matrix directly, $X = U\Sigma V^\top$, which costs $O(\min(np^2, n^2p))$ and avoids explicitly forming $S$ (better numerically conditioned, since squaring $X$ to get $S$ squares its condition number). The right singular vectors $V$ equal the PCA loadings $W$, and singular values relate to eigenvalues via $\lambda_i = \sigma_i^2/(n-1)$. If only the top $k \ll p$ components are needed, randomized or Lanczos-based iterative methods compute them in roughly $O(npk)$ time, avoiding the full $p^3$ or $np^2$ cost entirely.

Space complexity. Storing $X$ costs $O(np)$; storing $S$ costs $O(p^2)$; the output $W_k$ and $Z$ cost $O(pk)$ and $O(nk)$. Total space is $O(np + p^2)$, dominated by the covariance matrix when $p$ is large relative to $n$.

Correctness (why top-$k$ eigenvectors are optimal, not just a heuristic). Two equivalent formal justifications exist: Variance-maximization view. As derived above, the sequential constrained-optimization argument shows each successive principal direction is exactly the eigenvector of the covariance matrix with the next-largest eigenvalue — this is not approximate, it follows directly from setting the Lagrangian gradient to zero and using the spectral theorem's guarantee of an orthonormal eigenbasis for any real symmetric matrix. Best low-rank approximation view (Eckart–Young). PCA is equivalent to finding the rank-$k$ subspace minimizing reconstruction error:

Expanding this objective and using $\|x_i\|^2 = \|WW^\top x_i\|^2 + \|x_i - WW^\top x_i\|^2$ (Pythagoras, since $WW^\top x_i$ is an orthogonal projection), minimizing reconstruction error is equivalent to maximizing $\sum_i \|W^\top x_i\|^2 = \mathrm{trace}(W^\top S W)(n-1)$, the total projected variance. So variance-maximization and squared-error-minimization are the same objective viewed from two sides — this is why PCA is simultaneously "the projection that preserves the most variance" and "the best possible linear compression under squared error," a fact formalized generally by the Eckart–Young theorem for the Singular Value Decomposition.

Honest tradeoffs.

• PCA only captures linear structure; data on a curved manifold (e.g., a Swiss roll) will not be well summarized by a few linear directions, motivating Manifold Learning and Nonlinear Dimensionality Reduction.
• PCA is sensitive to the scale of features — a feature measured in millimeters will dominate variance over one measured in kilometers unless standardized.
• PCA is sensitive to outliers, since variance (a squared quantity) is heavily influenced by extreme points; robust variants (e.g., PCA on ranks, or robust covariance estimators) exist for contaminated data.
• Interpretability of components can be poor: each principal component is a linear combination of all original features, unlike sparse or interpretable feature-selection methods.

Consider $n=5$ two-dimensional points measuring (hours studied, exam score deviation from mean), already roughly linearly related:

Step 1 — center. Column means: $\bar{x}_1 = 5$, $\bar{x}_2 = 5.6$.

Step 2 — covariance matrix. Compute entries with $n-1 = 4$:

So

Step 3 — eigendecomposition. Solve $\det(S - \lambda I) = 0$:

Quadratic formula: $\lambda = \dfrac{18.3 \pm \sqrt{18.3^2 - 8}}{2} = \dfrac{18.3 \pm \sqrt{334.89 - 8}}{2} = \dfrac{18.3 \pm \sqrt{326.89}}{2} = \dfrac{18.3 \pm 18.08}{2}$.

Total variance $= \lambda_1 + \lambda_2 = 18.3$ (matches trace of $S = 10 + 8.3$, confirming the computation). Component 1 explains $18.19/18.3 \approx 99.4\%$ of the variance — expected, since the data lies almost exactly on a line. For $\lambda_1 \approx 18.19$: solve $(S - \lambda_1 I)w = 0$:

Normalizing $(1, 0.910)$ to unit length ($\sqrt{1^2+0.910^2}=1.348$):

The second eigenvector is orthogonal: $w_2 \approx (-0.675,\ 0.742)$.

Step 4/5 — project. Score for point A, centered $(-4,-3.6)$:

A single number per point (the PC1 score) now nearly reproduces the full 2D dataset, since PC1 captures 99.4% of the variance — exactly the compression PCA promises when features are strongly correlated.

PCA underlies preprocessing in nearly every machine learning pipeline that touches high-dimensional numeric data. In computer vision, Eigenfaces and related PCA-based methods compressed face images for recognition before deep learning dominated the field. In genomics, PCA reduces gene-expression or SNP data with tens of thousands of dimensions to a handful of components used for visualization and population-structure detection. In finance, PCA on asset return covariance matrices extracts common risk factors (level, slope, and curvature of yield curves are classic PCA-derived factors). In signal processing and sensor networks, PCA denoises data by discarding low-variance components assumed to be noise. Recommendation systems and NLP pipelines use PCA (or its close relative, Singular Value Decomposition, applied directly as in latent semantic analysis) to compress sparse high-dimensional representations before feeding them to downstream models. Almost every exploratory-data-analysis workflow begins with a PCA scatterplot to visually inspect clusters in two or three dimensions before deploying more expensive supervised or clustering methods.

Ask the tutor: How does PCA change if I compute it via the SVD of the data matrix instead of eigendecomposing the covariance matrix directly, and why is that usually the better choice numerically?