Principal Component Analysis Through the Singular Value Decomposition

Principal component analysis answers a narrow question: which linear subspace of dimension k keeps the most information in a data matrix, where information is measured as squared deviation from the mean. Three different formulations give the same answer, which is why the method shows up under several names. It is the best rank-k approximation to the centred data in the Frobenius norm, it is the projection that maximizes retained variance, and it is the linear transform whose output coordinates are uncorrelated. The link between the three is an orthogonal decomposition of the sample covariance, and the numerically reliable way to obtain it is the singular value decomposition of the data rather than an eigendecomposition of the covariance.

Centring, then one matrix factorization

The input is an n by d matrix X with one row per observation. PCA requires column means to be removed first, because the method maximizes variance and is otherwise steered by the offset of the origin; a constant shift of all data changes the directions completely, which is a useful reminder that the quantity being decomposed is covariance and not raw second moments. After centring, let X = U S V^T be the thin SVD, with S diagonal and non-negative in decreasing order. The columns of V are the principal directions, the columns of U S are the coordinates of the observations in that basis, and the sample covariance matrix is

$$X^T X / (n-1) = V , S^2 V^T / (n - 1)$$

so the sample variance along component i is s_i^2 / (n-1), and the components are mutually uncorrelated in the sample because V is orthogonal. The ratio of s_i^2 to the sum of all squared singular values is the explained variance ratio of that component.

Everything follows from that identity. Retaining the first k columns of V and projecting onto them maximizes the sum of the retained variances over all k-dimensional subspaces, and by the Eckart-Young theorem the same truncation minimizes the sum of squared reconstruction errors among all rank-k matrices, so the discarded singular values squared are exactly the residual error. The ordering is not a convention: it comes from the singular values being sorted.

Why the SVD and not the covariance matrix

The textbook route forms the d by d covariance matrix and diagonalizes it. That step doubles the condition number of the problem, because kappa(X^T X) = kappa(X)^2 for the spectral condition number, and squaring a condition number near 10^8 pushes the computed small eigenvalues into the noise of the arithmetic. Their signs and even their existence can be wrong, and small components are precisely the ones a scree plot is used to interpret. The SVD of the centred data avoids the squaring entirely, and it also gives the scores directly rather than requiring a separate multiplication.

The cost differs by shape. A full SVD of an n by d matrix costs O(n d^2) when n is at least d, with a constant that depends on the implementation; forming the covariance costs O(n d^2) as well and then diagonalizing it costs O(d^3), which dominates for very wide data. When the row count is small relative to the feature count, the useful trick is to decompose the n by n Gram matrix instead, since the non-zero spectra of X X^T and X^T X coincide. The directions are then recovered from the left factors rather than the right, at a cost of O(n^2 d), which is the cheaper side when n is much smaller than d.

When only the leading few components are needed, the randomized route is substantially cheaper. Multiply the data by a random matrix with k + p columns, with p a small oversampling margin, orthonormalize the result with a QR factorization, and compute a small SVD of the data restricted to that basis. The arithmetic is O(n d k) plus the QR work, and the error is bounded probabilistically in terms of the discarded singular values. A few applications of the power iteration (X X^T)^q widen the gap between the retained and discarded spectrum, which matters when the singular values decay slowly and a random sketch would otherwise be poorly aligned with the leading subspace.

What the outputs do and do not mean

A singular vector is defined only up to sign: if v is a principal direction then so is -v, with the scores negated and every other quantity unchanged. The sign is therefore arbitrary in the mathematics, and implementations impose a deterministic convention so that repeated runs on the same data produce identical output; a common rule is to flip a vector so that its entry of largest magnitude is positive. Code that compares components across libraries without accounting for that convention can see apparently different results that are identical up to a sign.

A principal direction is a direction of maximum variance, not a cause, not a cluster, and not a quantity with units. Its entries are weights on the original features, and a large weight means the feature varies a lot along that direction, which may simply reflect the measurement units. Two consequences follow. Analysing the correlation matrix, which is the covariance of standardized features, gives a different set of directions from analysing the covariance matrix, and neither is more correct: the first treats all features as equally important to the extent that their spreads are normalized, the second lets the feature with the largest spread dominate. The same scaling question is the subject of What a Distance Function Decides About Cluster Shape, since a metric and a projection both decide which coordinates matter.

The estimate of the retained variance is upwardly biased when the number of observations is not comfortably larger than the number of features. If n <= d, the centred data has rank at most n - 1, so all components beyond the first n - 1 have sample variance exactly zero, not approximately zero. The leading components absorb a share of what is sampling noise, and a scree plot from such a matrix depicts the rank deficiency as if it were structure. A related sensitivity concerns outliers: because the objective is a sum of squared deviations, a single distant observation can dominate the leading direction, and the ordinary SVD will faithfully report that direction.

Where the projection is useful

Whitening divides each retained component by its standard deviation, s_i / sqrt(n-1), so that the transformed data has identity covariance. That is convenient for subsequent models that assume isotropic noise, and it destroys the variance ordering: after whitening, no direction is more informative than another, and distances in the transformed space are no longer comparable to distances in the original one. A tied full covariance in a Gaussian mixture amounts to the same operation, which is the connection noted in Gaussian Mixtures and the EM Lower Bound.

PCA also supplies the initialization and the sanity check for nonlinear embedding methods. Because it is a linear projection with an exact reconstruction formula, it gives the best linear explanation of the data before any claim is made about curvature, and it is the usual starting point for the optimization described in Why t-SNE Preserves Neighborhoods but Distorts Distances. A nonlinear embedding that shows structure which no pair of linear components shows is evidence of nonlinearity; an embedding that reproduces the leading linear plane is not.

import numpy as np

def pca(X, k):
    Xc = X - X.mean(0, keepdims=True)          # centring is required
    U, s, Vt = np.linalg.svd(Xc, full_matrices=False)
    # deterministic sign convention: largest-magnitude entry of each row positive
    idx = np.argmax(np.abs(Vt), axis=1)
    signs = np.sign(Vt[np.arange(len(Vt)), idx])
    signs[signs == 0] = 1.0
    Vt = Vt * signs[:, None]
    scores = U[:, :k] * s[:k]                  # coordinates in the new basis
    var = (s ** 2) / (len(Xc) - 1)
    return scores, Vt[:k], var[:k] / var.sum()

def reconstruction_error(X, k):
    Xc = X - X.mean(0, keepdims=True)
    U, s, Vt = np.linalg.svd(Xc, full_matrices=False)
    approx = (U[:, :k] * s[:k]) @ Vt[:k]
    return np.linalg.norm(Xc - approx)         # discarded singular values squared

Choosing the number of components

Reconstruction error decreases monotonically as components are added, so it cannot by itself select a truncation; any value of k equal to the full rank reconstructs the centred data exactly, and the question is only how much error is acceptable. The usual diagnostics are geometric or predictive. A cumulative explained-variance curve with a chosen threshold is the crudest and most common. Parallel analysis, in which the observed singular values are compared against those of shuffled or random reference data of the same shape, gives a cutoff with a stated null model. Cross-validated reconstruction holds out rows, fits the projection on the rest, and measures the error of the held-out rows after projecting them onto the retained subspace; that procedure can identify a truncation below which held-out error rises.

The transform is fitted state, so the centring vector and the directions have to be stored and reapplied to any new data. Projecting uncentred new observations onto directions fitted on centred training data shifts every score by a constant and makes the scores incomparable across batches. When a preprocessing chain includes both a projection and a distance-based clustering, the projection changes the geometry before any neighbourhood or centroid decision is made, which is why the choice of retained dimension and the choice of clustering model interact rather than compose independently.