Matrix decompositions and PCA

Eigenvalues, the SVD, low-rank approximation and principal component analysis derived from scratch rather than called as a black box.

Eigenvalues and eigenvectors

A matrix is a linear map. Most vectors are rotated by it, but a few special directions are only stretched. Those are the eigenvectors, and the stretch factors are the eigenvalues: A v = lambda v. In ML the directions with the largest eigenvalues are the directions along which the data varies most.

import numpy as np

A = np.array([[2.0, 1.0],
              [1.0, 2.0]])

values, vectors = np.linalg.eigh(A)   # eigh assumes symmetric: use it for covariances
values                                # [1.0, 3.0]  ascending
vectors                               # columns are the eigenvectors, orthonormal

v = vectors[:, 0]
np.allclose(A @ v, values[0] * v)     # True by definition

• Use eigh for symmetric matrices (covariances, Hessians) — it is faster and returns real, orthonormal eigenvectors. Use eig for the general, possibly complex case.
• Eigenvectors are only defined up to sign: v and -v are both valid. Never compare eigenvector signs across runs without aligning them.
• A positive semi-definite matrix (all eigenvalues >= 0) is exactly what a covariance matrix is, and it is what guarantees a distance is a real distance.
• Eigenvalues of a covariance matrix tell you the variance captured along each principal direction.

The SVD and low-rank approximation

The singular value decomposition writes any matrix as A = U S V^T: U and V are rotations with orthonormal columns, S is a diagonal of non-negative singular values in descending order. Unlike an eigendecomposition it always exists, for any shape.

A = np.random.default_rng(0).normal(size=(200, 50))
U, S, Vt = np.linalg.svd(A, full_matrices=False)

S.shape          # (50,)  singular values, descending
U.shape, Vt.shape  # (200, 50), (50, 50)

# best rank-k approximation, guaranteed by the Eckart-Young theorem
k = 10
A_k = (U[:, :k] * S[:k]) @ Vt[:k, :]
np.linalg.norm(A - A_k)            # error equals sqrt(sum(S[k:]**2))

# energy retained
(S[:k] ** 2).sum() / (S ** 2).sum()

PCA from first principles

PCA is the SVD of the centred data. Centring is not optional: without subtracting the mean, the first component points at the mean vector instead of the direction of greatest variance.

X = np.random.default_rng(1).normal(size=(500, 12)) @ np.diag(np.linspace(3, 1, 12))

Xc = X - X.mean(axis=0)                  # centre, never skip this
U, S, Vt = np.linalg.svd(Xc, full_matrices=False)

explained = S ** 2 / (S ** 2).sum()
np.cumsum(explained)[:5]                 # how many components for 90%?

Z = Xc @ Vt[:3].T                        # project to 3 dimensions
Xr = Z @ Vt[:3] + X.mean(axis=0)         # reconstruct in original space

from sklearn.decomposition import PCA
pca = PCA(n_components=3).fit(X)
np.allclose(np.abs(Vt[:3]), np.abs(pca.components_), atol=1e-8)   # signs may differ

• Scale before PCA when features have different units; otherwise the column with the biggest numbers dominates every component.
• Use SVD on the centred matrix rather than eigh on the covariance: it is numerically more stable and never squares the condition number.
• PCA is linear and unsupervised. It finds variance, not class separation — for supervised projection use LDA, and for nonlinear structure use t-SNE or UMAP.
• PCA is a rotation, not a feature-selection method: every component is a dense combination of all original columns.

FAQ

Related

Vectors and matrices for ML Numerical stability and floating point

Last refreshed 2026-09-18.