Principal component analysis

Problem

A dataset with 300 features is hard to plot, slow to search and full of redundancy: survey questions that measure the same attitude, pixels that move together, sensor channels that track each other. You want fewer numbers per row that keep most of the information, to visualize the data, to speed up a downstream model, to remove noise, or to make clustering and nearest-neighbour search work in fewer dimensions.

"Most of the information" needs a definition. PCA's answer: keep the k-dimensional linear subspace that preserves the most variance, which is the same as the subspace that reconstructs the points with the smallest squared error. That choice is cheap to compute and easy to invert, but it is linear and it is sensitive to the units of each feature, and both limits come up below.

Intuition

Picture a cloud of points shaped like a tilted cigar. Its long axis is the direction where points differ most; knowing only where a point sits along that axis already tells you most of where it is. That axis is the first principal component. The second is the longest axis at right angles to the first, and so on. PCA is a rotation of the coordinate system so that the axes line up with the cloud, sorted from widest to thinnest. Dropping the thin axes is the compression: every point is projected onto the wide ones.

Projecting onto PC1 turns each pair of numbers into one and loses only the short dashed offsets, the part of each point that lies along PC2.

Two things the picture hides. First, PCA only sees straight lines: if the data lies on a curve, no rotation lines an axis up with it, and points far apart along the curve can land on top of each other. Second, "widest" depends on units: measure one feature in grams and another as a fraction, and the grams axis wins whatever the data says. Both are worked through below.

Mechanics

Variance-maximizing directions

Stack n points of d features as rows of X and subtract the column means, giving the centred matrix X_c. The sample covariance is

C = X_cᵀ X_c / (n − 1)

The variance of the data projected onto a unit vector w is wᵀ C w. Maximizing it subject to ‖w‖ = 1 with a Lagrange multiplier λ gives C w = λ w: the best direction is an eigenvector of C, and the variance along it is its eigenvalue, so the first principal component is the eigenvector with the largest eigenvalue. Requiring the next direction to be orthogonal to the first gives the second eigenvector, and so on. Because C is symmetric, its eigenvectors are orthogonal and the eigenvalues are real and non-negative (see the linear algebra toolkit).

The same subspace also minimizes the squared reconstruction error. For any unit w, ‖x‖² = (wᵀx)² + ‖x − w wᵀx‖², so for centred data, maximizing the variance kept is the same as minimizing what is thrown away. That is why "PCA" can be introduced either way.

PCA via the SVD

Forming C squares the condition number, so libraries usually skip it. Take the singular value decomposition X_c = U S Vᵀ. Then

C = V S² Vᵀ / (n − 1)

so the rows of Vᵀ are the principal directions (scikit-learn's components_), the eigenvalues are λ_j = S_j² / (n − 1) (explained_variance_), and the coordinates of the points in the new basis, the scores, are X_c V = U S. Eigenvectors are only defined up to sign, so the same data could return +v or −v. scikit-learn's svd_flip fixes the sign so that the largest-magnitude entry of each row of Vᵀ is positive, which makes results reproducible across runs and LAPACK builds.

scikit-learn's PCA picks among four solvers when svd_solver="auto": "covariance_eigh" (form the d × d covariance and diagonalize it) when there are at most 1,000 features and at least ten times as many rows as features; otherwise "full" (an exact LAPACK SVD) when the larger dimension is at most 500 or n_components="mle"; otherwise "randomized" when 1 ≤ k < 0.8 · min(n, d); and "full" in the remaining cases. The covariance route is fast for tall, thin data and is safe there because the d × d matrix is small; it clips tiny negative eigenvalues caused by rounding to zero.

The randomized solver (Halko, Martinsson and Tropp) multiplies X_c by a random d × (k + p) Gaussian matrix, with oversampling p (n_oversamples=10), sharpens the result with a few power iterations X_cX_cᵀ, orthonormalizes it into a basis Q, and takes the exact SVD of the small matrix Qᵀ X_c. When the spectrum decays, the top k components come out accurate at a fraction of the cost. PyTorch's torch.pca_lowrank does the same: it centres the columns, uses q = min(6, m, n) by default and 2 power iterations, and returns U, S, V, leaving the S² / (m − 1) step to you.

Explained variance and choosing k

The share of the total variance that component j accounts for is

ratio_j = λ_j / Σ_i λ_i

A scree plot shows these ratios in order; the cumulative sum tells you how much a k-component projection keeps. Common ways to pick k:

• A variance threshold. Pass a float such as n_components=0.9, and scikit-learn keeps the smallest k whose cumulative ratio is strictly greater than 0.9 (a searchsorted(..., side="right") + 1 on the cumulative sum).
• The elbow of the scree plot: the point where the ratios flatten into a floor of noise.
• Minka's MLE (n_components="mle"): treat the data as a low-rank signal plus isotropic Gaussian noise and pick the rank with the highest approximate Bayesian evidence. It needs n ≥ d.
• The downstream task. If a classifier on top is what matters, cross-validate k like any other hyperparameter.

After fitting with k components, scikit-learn stores noise_variance_, the mean of the dropped eigenvalues. It is the noise estimate of probabilistic PCA (Tipping and Bishop), which treats each point as W z + μ + ε with z ~ N(0, I) and ε ~ N(0, σ² I); that model is what score and get_covariance use.

Projection and reconstruction

Keeping the first k rows V_kᵀ of Vᵀ, a point x maps to k scores and back:

z = V_kᵀ (x − μ)     x̂ = μ + V_k z

Over the training set, the total squared reconstruction error Σ‖x_i − x̂_i‖² is exactly the sum of the dropped squared singular values, S_k+1² + … + S_r², or (n − 1) times the dropped eigenvalues: Eckart–Young says no other rank-k projection does better. Reconstruction is also how PCA denoises: noise spreads its variance over every direction, while signal concentrates in a few, so projecting onto the top components keeps most of the signal and removes most of the noise. A point whose reconstruction error is far above the training points' is unlike them, a cheap anomaly score.

Whitening

Whitening divides each score by the standard deviation along its component, so that the new features are uncorrelated and all have variance 1:

z_j = v_jᵀ (x − μ) / √λ_j

The covariance of whitened data is the identity. That helps models that assume isotropic, equally scaled inputs, such as RBF-kernel SVMs, k-means or ICA. The cost: low-variance directions, which are often mostly noise, get blown up to the same size as the signal, and inverse_transform has to undo the scaling. In scikit-learn, whiten=True changes transform but not components_; the division uses √λ clipped below at machine epsilon so a zero-variance direction does not divide by zero.

Scaling comes first

Variance has units. A feature in grams with a spread of tens has a variance thousands of times larger than a fraction that varies in the second decimal, and PC1 will simply be "the grams feature". Unless features already share a unit and a meaningful scale (pixel intensities, say), standardize them first with a StandardScaler. PCA on standardized data is PCA of the correlation matrix rather than the covariance matrix.

Nonlinear alternatives: t-SNE and UMAP

When the data lies on a curved surface, any linear projection flattens it by squashing, which can put far-apart points on top of each other. Neighbour-embedding methods drop the goal of a global linear map and instead try to keep each point's nearest neighbours near it in 2 or 3 dimensions.

t-SNE (van der Maaten and Hinton) turns distances in the original space into neighbour probabilities with a Gaussian around each point, whose width is tuned so that every point has the same perplexity (roughly, an effective number of neighbours, 30 by default), and symmetrizes them into p_ij. In the embedding it uses a heavy-tailed Student-t kernel, q_ij ∝ (1 + ‖y_i − y_j‖²)⁻¹ for a 2-D output, and moves the points y by gradient descent to minimize KL(P ‖ Q). The heavy tail is the trick: a 2-D map has too little room for all the moderately distant points of a high-dimensional space (the "crowding problem"), and the t kernel lets them sit much farther apart than a Gaussian would allow. scikit-learn's TSNE runs the first 250 iterations with the p_ij multiplied by an early exaggeration of 12, which pulls clusters together; uses the Barnes-Hut approximation with about 3 × perplexity neighbours per point by default; starts from a PCA projection scaled down so PC1 has standard deviation 10⁻⁴; and sets the learning rate to max(n / 12 / 4, 50).

UMAP (McInnes, Healy and Melville) builds a weighted k-nearest-neighbour graph (n_neighbors, 15 by default in the reference implementation) whose edge weights fall off with distance past each point's nearest neighbour, then lays out a low-dimensional graph to minimize the cross-entropy between the two, with negative sampling. min_dist sets how tightly points may pack. It is usually faster than t-SNE, keeps somewhat more of the global arrangement, and can embed new points with transform; it is a separate package (umap-learn), not part of scikit-learn.

Caveats for both, which matter more than the differences:

• Cluster sizes mean nothing. Each point's neighbourhood is rescaled to the same effective size, so dense and sparse clusters come out about equally wide.
• Distances between clusters mean little. The loss barely penalizes moving far-apart groups, so the gap between two islands is not a measure of how different they are.
• The hyperparameters change the picture. Different perplexities or n_neighbors values give different-looking maps, and random noise at low perplexity can show apparent clusters.
• They are stochastic and, for t-SNE, non-parametric. Different seeds give different maps; scikit-learn's TSNE has no transform for new points.
• They are for looking, not for features. Use the coordinates as model inputs only with care; PCA or a learned embedding is usually the better pipeline step, and running PCA to about 50 dimensions before t-SNE is standard to cut noise and cost.

Worked example

Two big components, then a flat floor of noise: the elbow and the 0.9 threshold agree on k = 2 here. On real data the floor is rarely this flat.

Points 0 and 4 are a quarter of the circle apart along the arc, yet PCA gives them the same coordinate; so do 1 and 3, 9 and 11, and 8 and 12.

Implementation

import numpy as np
from sklearn.decomposition import PCA
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

pca = make_pipeline(StandardScaler(), PCA(n_components=0.9))
Z = pca.fit_transform(X)                       # (n, k), k chosen by the 0.9 threshold
p = pca[-1]
p.components_, p.explained_variance_ratio_     # (k, d) rows of Vᵀ, shares of variance
X_back = pca.inverse_transform(Z)              # reconstruction in the original units

# The same thing from the SVD:
def pca_svd(X, k):
    mu = X.mean(0)
    U, S, Vt = np.linalg.svd(X - mu, full_matrices=False)
    flip = np.sign(Vt[np.arange(len(Vt)), np.abs(Vt).argmax(1)])
    U, Vt = U * flip, Vt * flip[:, None]       # svd_flip: largest entry of each row positive
    var = S**2 / (len(X) - 1)                  # explained_variance_
    return mu, Vt[:k], var[:k], (U * S)[:, :k]  # mean, components, variances, scores

# PyTorch, randomized, on a GPU tensor:
#   U, S, V = torch.pca_lowrank(A, q=k, center=True, niter=2)
#   scores = A_centred @ V[:, :k];  variances = S**2 / (A.shape[0] - 1)

Fit the scaler and PCA on the training split only and reuse them for validation and test data, or information leaks across the split; a pipeline does that automatically. For data that does not fit in memory, IncrementalPCA updates the components batch by batch; for sparse matrices such as TF-IDF, TruncatedSVD skips the centring, which would make the matrix dense. For a 2-D picture of high-dimensional data, reduce to about 50 dimensions with PCA first, then run t-SNE or UMAP on that, and try more than one perplexity or n_neighbors and seed before trusting any structure you see.

Tradeoffs

• PCA is the default because it is deterministic, fast, invertible and has one knob. Its components are hard to name, since each mixes every feature; for interpretable parts, sparse PCA or non-negative matrix factorization trade some variance for sparser loadings.
• Variance is not importance. The directions that separate classes can have small variance; PCA ignores labels, so a supervised method (LDA, or simply the model's own features) can beat it as preprocessing.
• Outliers dominate variance. One extreme point can rotate PC1 toward itself; robust PCA variants or clipping help.
• t-SNE and UMAP make persuasive pictures. Treat clusters in them as hypotheses to check in the original space, not as findings.

Connections

The linear algebra toolkit supplies the eigen-decomposition, the SVD and power iteration this page is built from; feature preprocessing supplies the standardization PCA needs first; and PCA is the usual step before clustering or nearest-neighbour search in many dimensions. An autoencoder with linear layers and squared loss learns the same subspace as PCA, and with nonlinear layers it is a nonlinear generalization.

References

• Pearson. On Lines and Planes of Closest Fit to Systems of Points in Space. Philosophical Magazine 1901. doi:10.1080/14786440109462720
• Hotelling. Analysis of a Complex of Statistical Variables into Principal Components. Journal of Educational Psychology 1933. doi:10.1037/h0071325
• Eckart and Young. The Approximation of One Matrix by Another of Lower Rank. Psychometrika 1936. doi:10.1007/BF02288367
• Tipping and Bishop. Probabilistic Principal Component Analysis. Journal of the Royal Statistical Society B 1999.
• Minka. Automatic Choice of Dimensionality for PCA. NIPS 2000.
• Halko, Martinsson and Tropp. Finding Structure with Randomness: Probabilistic Algorithms for Constructing Approximate Matrix Decompositions. SIAM Review 2011. arXiv:0909.4061
• van der Maaten and Hinton. Visualizing Data using t-SNE. Journal of Machine Learning Research 2008.
• McInnes, Healy and Melville. UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction. 2018. arXiv:1802.03426
• Wattenberg, Viégas and Johnson. How to Use t-SNE Effectively. Distill 2016. doi:10.23915/distill.00002
• scikit-learn 1.9.1: PCA, transform and inverse_transform, TSNE. PyTorch v2.15.0: pca_lowrank.