Real data rarely uses its coordinate axes efficiently. Height and weight, say, are two separate numbers for every person, yet they are far from independent: taller people tend to be heavier, so a scatter plot of the two is not a shapeless cloud but a cigar-shaped ellipse tilted across the page. Most of what varies from person to person happens along the long axis of that cigar — a single combined "size" direction — while only a little residual variation is left across its short axis. Principal component analysis (PCA) is the precise recipe for finding that long axis, and every other axis perpendicular to it, directly from the data itself, with no assumptions beyond the covariance structure already present in the numbers.
Maximizing Variance: The First Principal Component
Collect \(n\) observations of \(d\) numeric variables as the rows of a data matrix \(\mathbf{X}\in\mathbb{R}^{n\times d}\), and center each column on its own mean so that every column of \(\mathbf{X}\) sums to zero. A principal component direction is a unit vector \(\mathbf{p}\in\mathbb{R}^d\), \(\|\mathbf{p}\|_2=1\), and projecting every row \(\mathbf{x}_i\) onto it gives one number per observation, \(z_i=\mathbf{x}_i^T\mathbf{p}\), collected as \(\mathbf{z}=\mathbf{Xp}\). The first principal component is the direction \(\mathbf{p}_1\) that makes this one-dimensional summary \(\mathbf{z}\) capture as much of the spread of the original data as possible, i.e. the direction that maximizes the sample variance of \(\mathbf{z}\),
\[\operatorname{Var}(\mathbf{z})=\frac{1}{n-1}\mathbf{z}^T\mathbf{z}=\frac{1}{n-1}\mathbf{p}^T\mathbf{X}^T\mathbf{Xp} =\mathbf{p}^T\mathbf{Cp},\]
where \(\mathbf{C}=\frac{1}{n-1}\mathbf{X}^T\mathbf{X}\) is the \(d\times d\) sample covariance matrix: its diagonal entries are the variances of each individual variable and its off-diagonal entries are the covariances between pairs of variables. Since \(\mathbf{C}\) is symmetric and positive semidefinite, this is exactly the setup for the Rayleigh quotient: maximizing \(\mathbf{p}^T\mathbf{Cp}\) subject to \(\|\mathbf{p}\|_2=1\) is solved by taking \(\mathbf{p}\) to be the eigenvector of \(\mathbf{C}\) belonging to its largest eigenvalue, and the achieved maximum variance is exactly that eigenvalue. The same conclusion also drops out directly from a Lagrange multiplier on the constraint \(\mathbf{p}^T\mathbf{p}=1\): stationary points of \(\mathbf{p}^T\mathbf{Cp}-\lambda(\mathbf{p}^T\mathbf{p}-1)\) satisfy
\[\frac{\partial}{\partial\mathbf{p}}\Big[\mathbf{p}^T\mathbf{Cp}-\lambda(\mathbf{p}^T\mathbf{p}-1)\Big] =2\mathbf{Cp}-2\lambda\mathbf{p}=\mathbf{0}\quad\Longleftrightarrow\quad\mathbf{Cp}=\lambda\mathbf{p},\]
so every candidate direction is forced to be an eigenvector of \(\mathbf{C}\), and among those, \(\mathbf{p}^T\mathbf{Cp}=\lambda\mathbf{p}^T\mathbf{p}=\lambda\) is largest for the largest eigenvalue — the Rayleigh-quotient argument and the Lagrange-multiplier argument are two views of the same fact.
The Remaining Components
Having removed the variance captured by \(\mathbf{p}_1\), the same question is asked again restricted to directions orthogonal to it: which unit vector \(\mathbf{p}_2\perp\mathbf{p}_1\) maximizes \(\mathbf{p}_2^T\mathbf{Cp}_2\)? Because \(\mathbf{C}\) is symmetric, its eigenvectors can always be chosen mutually orthogonal, and restricting the Rayleigh quotient to the subspace orthogonal to \(\mathbf{p}_1\) is again maximized by an eigenvector — specifically the one for the second-largest eigenvalue. Repeating this \(d\) times exhausts the whole space and recovers the full eigendecomposition \(\mathbf{C}=\mathbf{P}\mathbf{\Lambda}\mathbf{P}^T\), where the columns of \(\mathbf{P}\) are the principal components \(\mathbf{p}_1,\ldots,\mathbf{p}_d\) sorted by decreasing eigenvalue \(\lambda_1\ge\lambda_2\ge\cdots\ge\lambda_d\ge0\). The transformed coordinates \(\mathbf{Z}=\mathbf{XP}\), called the scores, are themselves uncorrelated (\(\operatorname{Cov}(\mathbf{Z})=\mathbf{P}^T\mathbf{CP}=\mathbf{\Lambda}\) is diagonal) — PCA is precisely the rotation of the data onto axes along which every cross-correlation vanishes, ordered from most to least variance. The entries of a component \(\mathbf{p}_j\) itself are called its loadings: they are the weights with which each original variable enters that component, so a loading close to zero means the corresponding variable contributes almost nothing to that particular direction, while a large loading (positive or negative) marks a variable that dominates it.
None of this pins down the principal components perfectly uniquely: flipping any \(\mathbf{p}_j\) to \(-\mathbf{p}_j\) leaves \(\mathbf{p}_j^T\mathbf{Cp}_j\) unchanged, and whenever two eigenvalues coincide, any orthonormal basis of their shared eigenspace works equally well — the same sign and rotation ambiguity already discussed for the singular value decomposition's singular vectors. Neither ambiguity affects the variances captured or the quality of any reconstruction below.
An Equivalent View: Minimizing Reconstruction Error
A second, seemingly different way to ask for a good low-dimensional summary is to look for the \(k\)-dimensional subspace, spanned by the orthonormal columns of \(\mathbf{W}\in\mathbb{R}^{d\times k}\) (\(\mathbf{W}^T\mathbf{W}=\mathbf{I}_k\)), onto which projecting and then reconstructing every point loses as little as possible:
\[\min_{\mathbf{W}^T\mathbf{W}=\mathbf{I}_k}\ \sum_{i=1}^n\big\|\mathbf{x}_i-\mathbf{W}\mathbf{W}^T\mathbf{x}_i\big\|_2^2.\]
Here \(\mathbf{W}^T\mathbf{x}_i\) is the projection of \(\mathbf{x}_i\) onto the subspace and \(\mathbf{WW}^T\mathbf{x}_i\) is the reconstruction back in the original \(d\)-dimensional space. For any single point, orthogonality of the projection gives the Pythagorean identity
\[\|\mathbf{x}_i\|_2^2=\big\|\mathbf{W}^T\mathbf{x}_i\big\|_2^2+\big\|\mathbf{x}_i-\mathbf{WW}^T\mathbf{x}_i\big\|_2^2,\]
since the reconstruction \(\mathbf{WW}^T\mathbf{x}_i\) lies in the column space of \(\mathbf{W}\) and the residual \(\mathbf{x}_i-\mathbf{WW}^T\mathbf{x}_i\) is orthogonal to it. Summed over all \(n\) points, the total squared norm \(\sum_i\|\mathbf{x}_i\|_2^2\) is a fixed number that does not depend on \(\mathbf{W}\) at all, so minimizing the reconstruction error on the left is exactly the same as maximizing the captured energy \(\sum_i\|\mathbf{W}^T\mathbf{x}_i\|_2^2=\operatorname{trace}(\mathbf{W}^T\mathbf{X}^T\mathbf{XW})\) on the right — precisely the total variance captured along the \(k\) directions in \(\mathbf{W}\), summed. The variance- maximizing solution from the previous section is therefore also the reconstruction-error-minimizing one: the optimal \(\mathbf{W}\) is the top \(k\) eigenvectors of \(\mathbf{C}\), and the minimum achievable reconstruction error is exactly the sum of the discarded eigenvalues, \(\sum_{j=k+1}^d\lambda_j\).
Choosing the Number of Components: Explained Variance
Since \(\mathbf{C}\) is symmetric, its trace — the sum of the individual variables' own variances — equals the sum of all its eigenvalues, \(\operatorname{trace}(\mathbf{C})=\sum_{j=1}^d\lambda_j\). This total is exactly the quantity that "100% of the variance" refers to, and each eigenvalue's share of it,
\[\text{explained variance of component }j=\frac{\lambda_j}{\sum_{i=1}^d\lambda_i},\]
is the fraction of the data's total variance carried by that one component alone — not an R²-like goodness-of-fit measure and not a probability, just a share of a fixed total that add up to \(100\%\) across all \(d\) components. Keeping the first \(k\) components retains the cumulative fraction \(\sum_{i=1}^k\lambda_i\big/\sum_{i=1}^d\lambda_i\), and a common rule of thumb is to keep the smallest \(k\) for which this cumulative fraction clears a chosen threshold, e.g. 90%. Concretely, for a five-variable dataset with covariance eigenvalues \(\lambda=(11.81,\,5.19,\,1.67,\,0.47,\,0.05)\), the individual shares are \((61.5\%,\,27.1\%,\,8.7\%,\,2.4\%,\,0.3\%)\), so the first two components already explain \(61.5\%+27.1\%=88.6\%\) of the total variance and the first three clear the 90% mark at \(97.3\%\) — the point of diminishing returns is visible directly in how quickly the bars in a plot of these shares taper off:
This variance-based threshold is a useful default, but it is a proxy for how much of the original data is preserved, not a guarantee about how well a downstream task will perform on the reduced representation. When the scores feed into a subsequent classifier, the achieved accuracy typically improves sharply over the first handful of components and then plateaus long before the cumulative explained variance reaches 100%, since the later, low-variance components are increasingly dominated by measurement noise that carries little information relevant to the task — in practice it is worth checking task performance against \(k\) directly rather than assuming a fixed variance threshold such as 90% or 99% is automatically optimal for every use case.
Relation to the Singular Value Decomposition
Directly eigendecomposing \(\mathbf{C}=\frac{1}{n-1}\mathbf{X}^T\mathbf{X}\) is rarely how PCA is computed in practice. As already noted in the singular value decomposition chapter's discussion of this exact application, the eigenvectors of \(\mathbf{X}^T\mathbf{X}\) are precisely the right singular vectors of the centered matrix \(\mathbf{X}\) itself, with eigenvalues \(\sigma_j^2/(n-1)\) — so computing \(\mathbf{X}=\mathbf{U\Sigma V}^T\) directly gives \(\mathbf{P}=\mathbf{V}\) without ever forming \(\mathbf{X}^T\mathbf{X}\) explicitly, avoiding the same squared- condition-number penalty discussed there. Truncating to the top \(k\) singular vectors also connects the reconstruction-error view from above to the Eckart–Young theorem already proved in that chapter: the rank-\(k\) truncated SVD is the best possible rank-\(k\) approximation of \(\mathbf{X}\) in exactly the Frobenius-norm sense used above, so PCA reconstruction and optimal low-rank matrix approximation are the same statement viewed from two different fields.
Computing PCA When There Are Fewer Observations Than Variables
The covariance matrix \(\mathbf{C}\) is \(d\times d\), and eigendecomposing a general \(d\times d\) symmetric matrix costs on the order of \(O(d^3)\) operations — fine while \(d\) is small, but a single \(100\times100\)-pixel image already has \(d=10\,000\) pixel values as its variables, making \(\mathbf{C}\) a \(10\,000\times10\,000\) matrix with \(10^8\) entries. Whenever there are far fewer observations than variables, \(n\ll d\) — exactly the regime a modest collection of high-resolution face images falls into — that cost can be sidestepped entirely, since \(\mathbf{X}\) has at most \(n\) nonzero singular values and \(\mathbf{XX}^T\) is only \(n\times n\). Eigendecomposing this much smaller Gram matrix instead,
\[\mathbf{XX}^T\mathbf{u}_i=\sigma_i^2\mathbf{u}_i,\]
gives the same nonzero singular values \(\sigma_i\) as \(\mathbf{X}\), since \(\mathbf{XX}^T=\mathbf{U\Sigma}^2\mathbf{U}^T\) shares its eigenvalues with \(\mathbf{X}^T\mathbf{X}=\mathbf{C}(n-1)\) by construction of the SVD. The principal component directions themselves are then recovered from the identity \(\mathbf{X}^T\mathbf{u}_i =\sigma_i\mathbf{v}_i\), which rearranges to
\[\boxed{\mathbf{v}_i=\frac{\mathbf{X}^T\mathbf{u}_i}{\sigma_i}}\qquad(\sigma_i>0),\]
at a cost governed by \(n\) rather than \(d\). This is precisely the shortcut that made the Eigenfaces method above computationally practical in the first place: a modest training set of high-resolution face images has \(n\ll d\), so the eigenfaces are computed from the small \(n\times n\) Gram matrix and then lifted back into pixel space with this formula, never touching the enormous \(d\times d\) covariance matrix directly. When neither \(n\) nor \(d\) is small enough for a direct eigendecomposition — as with the activations of a large neural network — principal components can still be recovered without ever forming the covariance matrix by recasting the problem as a multi-player game, the subject of the EigenGame chapter.
Whitening: Making Every Direction Equally Important
PCA by itself only decorrelates the data; the scores along different components can still have very different variances, \(\operatorname{Var}(Z_i)=\lambda_i\). Some downstream methods instead want every direction to carry equal variance, a preprocessing step called whitening (or sphering), since a data cloud with identical variance in every direction looks spherical rather than ellipsoidal. Starting from the scores \(\mathbf{Z}=\mathbf{XP}\), dividing each column \(j\) by \(\sqrt{\lambda_j}\) rescales its variance to exactly one; written as a single linear map applied directly to \(\mathbf{X}\), this is
\[\boxed{\mathbf{X}_{\text{white}}=\mathbf{X}\,\mathbf{P}\mathbf{\Lambda}^{-1/2}},\]
where \(\mathbf{\Lambda}^{-1/2}=\operatorname{diag}(1/\sqrt{\lambda_1},\ldots,1/\sqrt{\lambda_d})\). Its covariance is indeed the identity matrix,
\[\frac{1}{n-1}\mathbf{X}_{\text{white}}^T\mathbf{X}_{\text{white}} =\mathbf{\Lambda}^{-1/2}\mathbf{P}^T\mathbf{C}\mathbf{P}\mathbf{\Lambda}^{-1/2} =\mathbf{\Lambda}^{-1/2}\mathbf{\Lambda}\mathbf{\Lambda}^{-1/2}=\mathbf{I},\]
using \(\mathbf{P}^T\mathbf{CP}=\mathbf{\Lambda}\) from before. This particular choice leaves the whitened data expressed in the rotated principal-component basis; multiplying once more by \(\mathbf{P}^T\) on the right, \(\mathbf{X}\mathbf{P}\mathbf{\Lambda}^{-1/2}\mathbf{P}^T\), rotates it back into the original coordinate axes while keeping the same identity covariance — a symmetric variant of whitening sometimes preferred because it distorts the original axes as little as possible. Either way, whitening only rescales variances and never changes correlations that are already zero, so any two orthogonal directions in the whitened data remain uncorrelated.
Principal Components Are Not a Regression Line
For two variables it is tempting to picture the first principal component as "the regression line", but the two are different objects that only sometimes look similar. Ordinary least-squares linear regression of \(y\) on \(x\) minimizes the sum of squared vertical distances (in \(y\) only) and is not symmetric in \(x\) and \(y\); the first principal component instead minimizes the sum of squared perpendicular distances to the line, treating every coordinate identically, which is exactly the reconstruction error minimized above for \(k=1\). On the correlated example data shown earlier, the regression slope works out to \(0.828\) while the principal-component (orthogonal regression) slope is \(0.839\) — close but not equal, and the gap grows as the noise perpendicular to the main direction grows relative to the spread along it.
Assumptions and Limitations
- Linearity. PCA can only find directions that are straight lines through the data; a dataset whose real structure curves through space (a "manifold") is not summarized well by any small number of linear directions, since PCA is built entirely from second moments and blind to anything else. Nonlinear extensions such as kernel PCA replace the plain dot product used to build \(\mathbf{C}\) with a kernel function, implicitly performing PCA in a higher-dimensional feature space where a curved relationship becomes linear, while dedicated manifold-learning methods aim to preserve local neighborhood distances instead of global variance.
- Sensitivity to scale and outliers. Since \(\mathbf{C}\) is built from squared deviations, a single variable measured on a much larger numeric scale, or a single extreme outlier, can dominate the leading eigenvector regardless of whether it is actually the most meaningful direction — standardizing each column to unit variance first (working with the correlation matrix instead of the covariance matrix) is the usual remedy when variables are not already on comparable units.
- Decorrelation is not independence. The scores \(\mathbf{Z}\) are uncorrelated by construction, but uncorrelated is a strictly weaker property than statistically independent unless the data is jointly Gaussian. When true statistical independence of the extracted components is the actual goal, independent component analysis (ICA) replaces PCA's decorrelation objective with an explicit independence criterion, at the cost of the components no longer being orthogonal or ordered by variance.
- Blind to class labels. PCA is built entirely from the covariance of the data itself, with no notion that some points might belong to one group and others to another; the direction of greatest overall variance and the direction that best separates two known groups can be completely different, or even perpendicular to each other. Linear discriminant analysis replaces variance with the ratio of between-class to within-class scatter whenever class labels are available and the actual goal is classification rather than description.
Applications
- Dimensionality reduction and visualization. Projecting onto the first two or three principal components turns a high-dimensional dataset into a scatter plot that keeps as much of the original spread as any two or three linear directions can, useful both for visual exploration and as a preprocessing step before clustering or classification. A classic illustration is the MNIST dataset of handwritten digit images: each \(28\times28\)-pixel image is a single point in a \(d=784\)-dimensional pixel space, yet projecting onto just its first two principal components already groups images by digit identity into visibly separate regions, since the pixel-level differences between distinct digits typically dominate the variation among different handwriting styles of the same digit.
- Noise reduction and compression. Discarding the low-variance trailing components removes directions that are more likely to be measurement noise than signal, and storing only the leading scores and loadings compresses the data in exactly the same rank-\(k\) sense as SVD-based data compression. In signal and image processing this same construction has an older name of its own: the Hotelling transform (Harold Hotelling, 1933) or, in its continuous formulation, the Karhunen–Loève transform (Kari Karhunen, 1947; Michel Loève, 1955) — both predate the name "principal component analysis" but describe exactly the same eigenvectors-of-the-covariance-matrix construction.
- Eigenfaces. Treating each face image as one long vector of pixel intensities and running PCA over a collection of faces produces a small set of principal directions — the "eigenfaces" — such that every face in the collection is well approximated by a short list of scores along them, an approach introduced by Sirovich and Kirby (1987) and developed into a practical face-recognition system by Turk and Pentland (1991). Recognizing a new face then reduces to comparing its short score vector against the stored score vectors of known faces — nearest neighbor in this low-dimensional space rather than pixel-by-pixel comparison in the original image space — and works well even when only a handful of leading components are kept. The same idea applies beyond raw pixel intensities: running PCA on the \((x,y)\) coordinates of corresponding facial landmarks across many faces instead yields "eigenfaces" of shape rather than texture, an extension used by Lanitis, Taylor, and Cootes (1995) to model how face geometry varies independently of its grey-value appearance.
Schematically, projecting a set of digit images this way and plotting each image at its own \((z_1,z_2)\) coordinate looks roughly like this — images of the same digit fall near each other purely because they share similar pixel patterns, with no digit labels ever entering the computation of the components themselves: