Principal Component Analysis

Principal component analysis (PCA) finds the orthogonal directions along which a set of vectors varies most, and re-expresses each vector by its coordinates along the few most important of them. Applied to images, the vectors are whole images, patches or descriptors, and the directions are themselves images: patterns such as an illumination gradient or the average difference between two faces, ranked by how much of the variation in the collection they account for. Keeping the first few gives a compressed, decorrelated description of each image; dividing each coordinate by its spread whitens the data. The method is also called the Karhunen–Loève transform (Bishop, 2006, §12.1).

Vision has used it to represent faces and keypoint descriptors, to shorten and whiten descriptors for retrieval, to whiten image patches before feature learning, to estimate surface normals and to augment colour (see Uses in Computer Vision). Computationally it is one singular value decomposition of a centred data matrix, and statistically it is a statement about the data’s covariance matrix, which Probability and Statistics defines.

Problem

Given $n$ images of $d$ pixels each, stored one per row in $X \in \mathbb{R}^{n \times d}$ as scikit-learn stores them, find a $k$-dimensional subspace, $k \ll d$, that represents the images as well as possible. There are two natural meanings of “as well as possible”, and PCA satisfies both at once:

• Maximum variance. Project the images onto the subspace and keep as much of their variance as possible, so that the projections remain as spread out, and therefore as distinguishable, as the subspace allows (Hotelling, 1933).
• Minimum reconstruction error. Approximate each image by a point in the subspace and make the mean squared distance between images and their approximations as small as possible (Pearson, 1901).

The two criteria were introduced separately, by Pearson for fitting lines and planes and by Hotelling for analysing correlated test scores (see History), and they have the same solution.

Inputs and Outputs

Inputs: a data matrix $X$ with one sample per row (an image flattened to $d$ pixel values, a patch, a descriptor), and the number of components $k$ or a rule for choosing it.

Outputs:

• the mean $\boldsymbol{\mu} \in \mathbb{R}^d$, the average image;
• the principal axes (components, loadings) $\mathbf{v}_1, \dots, \mathbf{v}_k$, orthonormal vectors in $\mathbb{R}^d$ that can be displayed as images, stacked as the columns of $V_k \in \mathbb{R}^{d \times k}$;
• the eigenvalues (explained variances) $\lambda_1 \ge \dots \ge \lambda_k$, the variance of the data along each axis;
• the scores $Z = (X - \mathbf{1}\boldsymbol{\mu}^\top) V_k \in \mathbb{R}^{n \times k}$, the coordinates of each sample in the new basis, which can be computed for new samples with the same formula.

Intuition

Picture 500 photographs of the same object taken under changing light. Every image has hundreds of pixels, but the images differ from each other in only a few ways: the scene gets brighter or darker as a whole, the light comes more from the left or from the right, the object stands out more or less from its background. Each of these changes moves all the pixels together in a fixed pattern. PCA discovers those patterns without being told about them: the first principal axis is the single pixel pattern that accounts for the most variation between the images, the second accounts for the most of what is left, and so on. When the variation really comes from three such patterns plus noise, the variances along the axes drop sharply after the third, and three numbers per image describe it almost completely.

The same picture explains the method’s main limit. PCA looks for patterns that are added to the image in varying amounts; a change that moves the object by a few pixels is not of that kind, and it is spread over many components (the first example below shows both cases).

Algorithm

1. Centre. Compute the mean image $\boldsymbol{\mu}$ and subtract it from every row: $X_c = X - \mathbf{1}\boldsymbol{\mu}^\top$. Without centring, the first axis tends to point towards the mean rather than along the variation.
2. Optionally rescale. If the variables are in different units, divide each column by its standard deviation (see the section on standardization); pixel intensities in the same units are normally left as they are.
3. Decompose. Compute the thin SVD $X_c = U \Sigma V^\top$. The principal axes are the right singular vectors (the rows of NumPy’s Vt), and the eigenvalues are $\lambda_j = \sigma_j^2 / (n - 1)$. Equivalently, compute the eigendecomposition of the covariance matrix $S = X_c^\top X_c / (n - 1)$ and sort the eigenvalues in decreasing order.
4. Choose $k$ from the eigenvalues (see Parameters).
5. Project: $Z = X_c V_k$, which also equals $U_k \Sigma_k$.
6. Reconstruct, if needed: $\hat{X} = Z V_k^\top + \mathbf{1}\boldsymbol{\mu}^\top$.
7. Whiten, if needed: divide each column of $Z$ by $\sqrt{\lambda_j}$, and for ZCA whitening rotate back with $V_k^\top$.

New images are centred with the training mean and projected with the same $V_k$.

Mathematical Formulation

Maximum variance

The covariance matrix of the data, with the denominator scikit-learn uses, is

$$S = \frac{1}{n-1} X_c^\top X_c .$$

Bishop (2006, §12.1.1) and Tipping and Bishop (1999) divide by $n$ instead; the axes are the same and the eigenvalues differ by the factor $n/(n-1)$, as Probability and Statistics explains. The variance of the projections onto a unit vector $\mathbf{v}$ is $\mathbf{v}^\top S \mathbf{v}$. Maximizing it subject to $\mathbf{v}^\top\mathbf{v} = 1$ with a Lagrange multiplier $\lambda$ gives

$$S\mathbf{v} = \lambda \mathbf{v}, \qquad \mathbf{v}^\top S \mathbf{v} = \lambda ,$$

so $\mathbf{v}$ must be an eigenvector of $S$, and the variance is largest for the eigenvector with the largest eigenvalue (Bishop, 2006, §12.1.1). Requiring each further direction to be orthogonal to the previous ones and repeating the argument gives the eigenvectors in order of decreasing eigenvalue: the $k$-dimensional subspace of maximum projected variance is spanned by the $k$ leading eigenvectors of $S$, and the variance it keeps is $\lambda_1 + \dots + \lambda_k$. Since $S$ is symmetric and positive semidefinite, these eigenvalues are real and non-negative and the eigenvectors can be chosen orthonormal (linear algebra). In the eigenbasis the scores are uncorrelated: their covariance is $\operatorname{diag}(\lambda_1, \dots, \lambda_k)$ (Tipping and Bishop, 1999, §1).

Minimum reconstruction error, and why the two views agree

Let $V_k$ have orthonormal columns. The best approximation of a centred sample $\mathbf{x}$ within their span is its orthogonal projection $V_k V_k^\top \mathbf{x}$, and Pythagoras splits its squared norm into a kept and a lost part:

$$\lVert \mathbf{x} \rVert^2 = \lVert V_k^\top \mathbf{x} \rVert^2 + \lVert \mathbf{x} - V_k V_k^\top \mathbf{x} \rVert^2 .$$

Summing over the centred samples and dividing by $n - 1$, the left side is the total variance $\operatorname{tr}(S) = \sum_{j=1}^{d} \lambda_j$, which does not depend on $V_k$; the first term on the right is the variance kept by the projection; the second is the mean squared reconstruction error. Maximizing the first is therefore the same as minimizing the second, which is why Hotelling’s and Pearson’s criteria have the same solution. At the optimum the error is the sum of the discarded eigenvalues,

$$\frac{1}{n-1} \sum_{i=1}^{n} \lVert \mathbf{x}_i - \hat{\mathbf{x}}_i \rVert^2 = \sum_{j=k+1}^{d} \lambda_j ,$$

which Bishop derives directly from the minimum-error formulation, also showing that the optimal subspace passes through the mean (Bishop, 2006, §12.1.2, eq. 12.18).

In matrix form, the reconstruction $\hat{X}_c = X_c V_k V_k^\top = U_k \Sigma_k V_k^\top$ is the truncated SVD of $X_c$, so the minimum-error property is the Eckart–Young theorem applied to the centred data matrix: no rank-$k$ matrix is closer to $X_c$ in the Frobenius norm (singular value decomposition).

The eigendecomposition route and the SVD route

Substituting the thin SVD $X_c = U\Sigma V^\top$ into the covariance gives

$$S = \frac{1}{n-1} V \Sigma^2 V^\top ,$$

an eigendecomposition: the principal axes are the right singular vectors of the centred data matrix, the eigenvalues are $\lambda_j = \sigma_j^2/(n-1)$, and the scores are $Z = X_c V = U\Sigma$ (Murphy, 2022, §20.1.3.3). The two routes agree in exact arithmetic, but forming $X_c^\top X_c$ squares the singular values and loses the smallest of them to rounding, as the singular value decomposition article explains. Murphy (2022, §20.1.3.3) calls the SVD route often preferable for computational reasons, and scikit-learn’s documentation warns that its covariance-based solver is less numerically stable on data with a large range of singular values.

The covariance route has its own advantage. When there are many more samples than dimensions, $S$ is a small $d \times d$ matrix that can be accumulated in one pass over the data, and its eigendecomposition costs $O(d^3)$ regardless of $n$. When there are fewer samples than dimensions, the usual case for whole images, the $n \times n$ matrix $X_c X_c^\top$ has the same non-zero eigenvalues, and its eigenvectors $\mathbf{u}_j$ give the principal axes as $\mathbf{v}_j = X_c^\top \mathbf{u}_j / \sqrt{(n-1)\lambda_j}$, at a cost of $O(n^3)$ instead of $O(d^3)$ (Bishop, 2006, §12.1.4). The thin SVD contains both tricks: it works with whichever dimension is smaller.

With $n < d$, at most $n - 1$ eigenvalues are non-zero, because $n$ centred points span at most $n - 1$ dimensions (Bishop, 2006, §12.1.4); the remaining directions carry no variance, so there is little point in asking for more than $n - 1$ components.

Explained variance

The explained variance ratio of component $j$ is $\lambda_j / \sum_{i} \lambda_i$, the fraction of the total variance it accounts for, and its cumulative sum $F_k = \sum_{j \le k} \lambda_j / \sum_i \lambda_i$ is the fraction kept by $k$ components (Murphy, 2022, eq. 20.27). By the identity above, $1 - F_k$ is the fraction of the total squared deviation from the mean image that the reconstruction misses. A plot of $\lambda_j$ against $j$ is a scree plot, and the cumulative fraction carries the same information in increasing form (Murphy, 2022, §20.1.4.2).

Whitening: PCA and ZCA

A whitening transform maps centred data to data whose covariance is the identity: zero mean, unit variance in every direction, no correlations. Writing $\Lambda = \operatorname{diag}(\lambda_1, \dots, \lambda_d)$,

$$\mathbf{y}_{\text{PCA}} = \Lambda^{-1/2} V^\top (\mathbf{x} - \boldsymbol{\mu})$$

rotates onto the principal axes and divides each coordinate by its standard deviation; its covariance is $\Lambda^{-1/2} V^\top S V \Lambda^{-1/2} = I$ (Bishop, 2006, §12.1.3, eqs. 12.24–12.25, which calls this whitening or sphering).

Whitening matrices are not unique. A matrix $W$ whitens exactly when $W^\top W = S^{-1}$, so any rotation of a whitened vector is still white; Bell and Sejnowski (1997, eq. 5) described the special solutions of this family, two of which matter here. The PCA solution above, $W_P = \Lambda^{-1/2} V^\top$, has orthogonal rows; when image statistics are stationary, they noted, its filters are global Fourier filters. The symmetric solution,

$$W_Z = S^{-1/2} = V \Lambda^{-1/2} V^\top ,$$

rotates back into pixel space after rescaling. They called it ZCA because its filters are zero-phase (symmetric), and described them as local, centre-surround filters, each whitening one pixel while preserving the spatial arrangement of the image and flattening its amplitude spectrum. Krizhevsky (2009, Appendix A) derived the same matrix from the conditions that $W$ be symmetric and whiten, and called it the “Zero Components Analysis” solution; the two sources expand the acronym differently. ZCA output still looks like an image: Krizhevsky (2009, §1.3.3) found that it preserves edge information and sets regions of relatively uniform colour to zero. Both used it as preprocessing for learning image features: Bell and Sejnowski whitened $12 \times 12$ natural-image patches this way before training independent component analysis, and Krizhevsky whitened tiny images before training restricted Boltzmann machines on them.

Two practical points follow from the formula (our reasoning). Whitening divides by $\sqrt{\lambda_j}$, so directions with tiny eigenvalues, which in images are mostly noise, are amplified to the same variance as the strongest structure; a small constant is therefore added to the eigenvalues, or only the leading $k$ components are kept before rescaling. Jégou and Chum observed this harm in retrieval, where whitening too many dimensions hurt (Embeddings reports their analysis). And whitening removes the information about which directions had large variance, so Euclidean distances between whitened vectors weigh all retained directions equally, the Mahalanobis distance of Similarity and Distance Measures restricted to the principal subspace. In scikit-learn, PCA(whiten=True) applies the PCA form: the output has uncorrelated components with unit variance, and the documentation notes that this removes the relative variance scales of the components.

Standardization: covariance or correlation

PCA is not invariant to the units of the variables. Multiplying one variable by a constant stretches the data along that axis and pulls the principal axes towards it. Hotelling saw this in 1933: since changes of units can squeeze and stretch the ellipsoids in any way, the method applies only when each variable has a unit of unique importance, and its results depend on the metric chosen. He worked with standardized scores, each variable divided by its standard deviation, which amounts to the eigendecomposition of the correlation matrix instead of the covariance matrix. Murphy (2022, §20.1.3.1) recommends the correlation matrix, because PCA can otherwise be misled by directions whose variance is high merely because of the measurement scale; Bishop (2006, §12.1.3) standardizes when the variables are measured in different units or differ greatly in variability.

For image data the choice is less simple, because pixels already share a unit. Standardizing every pixel gives a nearly constant background pixel, whose small variance is all noise and illumination, the same weight as a pixel on the object. The third example below shows that this can hand the first component to a minor factor. Standardization repairs variables with incomparable units, such as a feature vector that mixes colour statistics with a shape measure, or one channel accidentally stored on a different scale; for raw pixels in common units, the covariance matrix is usually the right object (our reasoning from the example).

Probabilistic PCA

Standard PCA is not a probability model: it assigns no likelihood to an image, and an image that lies close to the principal subspace but far from all the training images still has a small reconstruction error (Bishop, 2006, §12.2). Tipping and Bishop (1999) obtained PCA as the maximum-likelihood solution of a latent-variable model. A $q$-dimensional latent vector $\mathbf{z} \sim \mathcal{N}(\mathbf{0}, I)$ generates the observation

$$\mathbf{x} = W\mathbf{z} + \boldsymbol{\mu} + \boldsymbol{\epsilon}, \qquad \boldsymbol{\epsilon} \sim \mathcal{N}(\mathbf{0}, \sigma^2 I),$$

so that $\mathbf{x} \sim \mathcal{N}(\boldsymbol{\mu}, WW^\top + \sigma^2 I)$: a Gaussian whose covariance is a rank-$q$ part plus isotropic noise. They showed that the likelihood is maximized by

$$W_{\text{ML}} = U_q (\Lambda_q - \sigma^2 I)^{1/2} R, \qquad \sigma^2_{\text{ML}} = \frac{1}{d - q} \sum_{j = q+1}^{d} \lambda_j ,$$

where $U_q$ holds the $q$ principal eigenvectors of the sample covariance, $\Lambda_q$ their eigenvalues and $R$ is an arbitrary rotation; other choices of eigenvectors are saddle points (their eqs. 7–8, §3.2). The noise variance is the variance lost in the projection, averaged over the discarded dimensions, and the columns of $W$ span the principal subspace of ordinary PCA. Bishop (2006, §12.2.1) attributes a similar result, conjectured without proof, to Roweis (1998). Tipping and Bishop list what the model adds: a likelihood that can be compared with other density models and used for classification or novelty detection, an EM algorithm for the principal subspace, PCA with missing data, and mixtures of PCA models. A sample can be summarized by the posterior mean of $\mathbf{z}$, which is pulled towards the origin relative to the orthogonal projection when $\sigma^2 > 0$ (Bishop, 2006, §12.2.1).

scikit-learn’s PCA stores this noise estimate as noise_variance_ and uses the model in its score and score_samples methods, citing Tipping and Bishop (scikit-learn 1.7.2 documentation). Its documented estimate averages the smallest $\min(n, d) - k$ eigenvalues rather than $d - q$ of them, which is the same when $n \ge d$ but not when there are fewer images than pixels, because most of the zero eigenvalues of the sample covariance are then left out of the average; the first example shows the difference.

Parameters

Number of components $k$. There is no single rule; the sourced ones are:

• A variance threshold. Keep the smallest $k$ whose cumulative explained variance exceeds a fraction such as 0.95. scikit-learn does this when n_components is a float between 0 and 1 (documented for the 'full' solver).
• The scree plot. Look for the “knee” or “elbow” where the eigenvalues stop falling quickly (Murphy, 2022, §20.1.4.2–§20.1.4.3). Murphy describes a profile-likelihood method that locates the elbow automatically by fitting a change-point model to the eigenvalues.
• Held-out reconstruction error does not work as a criterion: unlike a supervised model, PCA’s error on held-out data keeps decreasing as $k$ grows, because more components always approximate new data better (Murphy, 2022, §20.1.4.1).
• A likelihood. Probabilistic PCA turns the choice into model selection. scikit-learn offers n_components='mle', Minka’s automatic choice of dimensionality, and its user guide uses the probabilistic score for cross-validation.
• The downstream task. Ke and Sukthankar (2004) chose the dimension of PCA-SIFT by matching performance, and found that more components eventually hurt (see Uses in Computer Vision).

A threshold on explained variance counts noise as variance to be explained. With strong signal and weak noise, a 95% threshold may stop at the signal’s dimension while a 99% threshold reaches deep into the noise, as the first example shows.

Centring. Always on in PCA. Omitting it gives the uncentred SVD, which scikit-learn offers separately as TruncatedSVD, an alternative for sparse data that does not centre them.

Scaling. Covariance or correlation, as discussed above.

Whitening and its regularizer. Whether to divide by $\sqrt{\lambda_j}$, and how large a constant to add to the eigenvalues first.

Solver. In scikit-learn 1.7.2, svd_solver is one of:

• 'full': the exact thin SVD of the centred data through LAPACK (scipy.linalg.svd), truncated afterwards;
• 'covariance_eigh': the eigendecomposition of the covariance matrix, efficient when $n \gg d$ and $d$ is small, impractical for large $d$, and less stable numerically (added in version 1.5);
• 'arpack': a truncated SVD through ARPACK (scipy.sparse.linalg.svds), which needs $0 < k < \min(n, d)$ and supports sparse input;
• 'randomized': the randomized SVD of Halko, Martinsson and Tropp;
• 'auto' (the default): 'covariance_eigh' when the data have fewer than 1,000 features and more than 10 times as many samples as features; otherwise 'randomized' when the data are larger than $500 \times 500$ and $k$ is below 80% of the smaller dimension; otherwise 'full'.

That is the documented policy. The source code of version 1.7.2 applies the first test with non-strict inequalities (at most 1,000 features and at least 10 times as many samples), reads “larger than $500 \times 500$” as the larger dimension exceeding 500, switches sparse input to 'arpack', and falls back to 'full' when $k$ is given as a fraction or as 'mle'.

Complexity

For a dense $n \times d$ matrix, the exact thin SVD costs on the order of $nd \min(n, d)$ operations. scikit-learn’s user guide states the exact method as $O(n_{\max}^2 n_{\min})$ in time and $n_{\max} n_{\min}$ in memory, and the randomized method as $O(n_{\max}^2 k)$ in time and about $2 n_{\max} k$ in memory, with $n_{\max}$ and $n_{\min}$ the larger and smaller of $n$ and $d$. Forming the covariance costs $O(nd^2)$, and computing only its leading $k$ eigenvectors by the power method costs $O(kd^2)$ instead of $O(d^3)$ (Bishop, 2006, §12.1.1).

Randomized SVD. Halko, Martinsson and Tropp (2011) find the leading singular subspace by multiplying the matrix by a small random Gaussian matrix $\Omega$ with $k$ plus a few columns, orthonormalizing the result $Y = A\Omega$ to get a basis $Q$ that captures most of the matrix’s range, and computing the SVD of the small matrix $Q^\top A$ exactly. When the singular values decay slowly, they apply the matrix several times, $Y = (AA^\top)^q A\Omega$, which raises the singular values to the power $2q + 1$ and makes the dominant directions stand out; in floating point the powers must be interleaved with orthonormalizations (their Algorithms 4.1, 4.3 and 4.4 in the arXiv version). Their abstract gives $O(mn \log k)$ operations for a dense $m \times n$ input, against $O(mnk)$ for classical methods (the faster count uses a structured random matrix, their §4.6), and a constant number of passes over data too large for fast memory. scikit-learn’s randomized solver follows this paper, with a number of power iterations set by iterated_power ('auto' by default) and 10 extra random vectors (n_oversamples=10) by default. It is approximate: in the first example it recovers the three signal components exactly but returns different vectors for the components that only span noise, and the user guide warns that its inverse_transform is not the exact inverse of transform.

Incremental PCA. When the data do not fit in memory, scikit-learn’s IncrementalPCA processes mini-batches and updates the components, with memory proportional to the batch size times $d$; its results almost exactly match PCA, according to the user guide. The documentation cites Ross, Lim, Lin and Yang’s incremental learning for robust visual tracking (International Journal of Computer Vision, 2008).

For $64 \times 64$ face images, $d = 4{,}096$; scikit-learn’s user guide reports computing the first 16 components of the 400 Olivetti face images with the randomized solver in under a second.

Implementation

The examples use NumPy 2.3.5 and SciPy 1.16.3 for the mechanics and scikit-learn 1.7.2 as a cross-check, on synthetic images with exactly known structure. Each took under two seconds with single-threaded linear algebra when checked.

Images that vary along three known factors

The script renders 500 grayscale images of a disk on a background. Three factors change linearly from image to image: overall brightness, a horizontal illumination gradient and the disk’s contrast; independent noise with standard deviation 0.02 is added to every pixel. It computes PCA both ways, checks that the reconstruction error equals the sum of the discarded eigenvalues, compares with scikit-learn’s solvers, and then repeats the analysis after letting the disk also move by up to three pixels.

import numpy as np
from sklearn.decomposition import PCA
np.set_printoptions(suppress=True)

# 500 synthetic 24 x 24 grayscale images (intensities near [0, 1]) of a disk on a
# background. Three known factors change linearly from image to image: global
# brightness, a left-to-right illumination gradient, and the disk's contrast.
# Every pixel also gets independent Gaussian noise with standard deviation 0.02.
rng = np.random.default_rng(0)
n, H, W = 500, 24, 24
yy, xx = np.mgrid[0:H, 0:W]
ramp = (xx - (W - 1) / 2) / (W - 1)                         # -0.5 ... 0.5
disk = ((yy - 12) ** 2 + (xx - 12) ** 2 < 6 ** 2).astype(float)

def render(n, shift=0):
    b = rng.uniform(0.2, 0.6, n)                            # global brightness
    g = rng.uniform(-0.3, 0.3, n)                           # illumination gradient
    c = rng.uniform(0.1, 0.4, n)                            # disk contrast
    dx = rng.integers(-shift, shift + 1, n)                 # optional horizontal shift
    imgs = [b[i] + g[i] * ramp + c[i] * np.roll(disk, dx[i], axis=1) for i in range(n)]
    return np.array(imgs).reshape(n, -1) + rng.normal(0, 0.02, (n, H * W))

X = render(n)                                               # n x d, one image per row
mu = X.mean(axis=0)
Xc = X - mu                                                 # centre every pixel

# Route 1: SVD of the centred data matrix (Vt holds the principal axes as rows).
U, s, Vt = np.linalg.svd(Xc, full_matrices=False)
lam = s ** 2 / (n - 1)                                      # eigenvalues, denominator n - 1
# Route 2: eigendecomposition of the covariance matrix (eigh sorts ascending).
C = Xc.T @ Xc / (n - 1)
lam_eig, E = np.linalg.eigh(C)
lam_eig, E = lam_eig[::-1], E[:, ::-1]

print("top 6 eigenvalues (SVD):", np.round(lam[:6], 5))
print("top 6 eigenvalues (eigh):", np.round(lam_eig[:6], 5))
print("noise variance per pixel:", 0.02 ** 2)
print("|cos| between SVD and eigh axes 1-3:", np.round(np.abs(np.sum(Vt[:3] * E[:, :3].T, axis=1)), 10))

# Reconstruction with k components: mean squared error equals the discarded eigenvalues.
for k in [1, 2, 3, 4, 10]:
    Z = Xc @ Vt[:k].T                                       # scores, n x k
    Xhat = Z @ Vt[:k] + mu                                  # back to pixel space
    err = np.sum((X - Xhat) ** 2) / (n - 1)
    print(f"k={k:2d}  sum ||x - xhat||^2 / (n-1) = {err:.5f}   sum of discarded eigenvalues = {lam[k:].sum():.5f}"
          f"   explained = {lam[:k].sum() / lam.sum():.4f}")

# Cross-check with scikit-learn 1.7.2: eigenvalues, axes (up to sign) and signs.
for solver in ["full", "covariance_eigh", "randomized"]:
    p = PCA(n_components=5, svd_solver=solver, random_state=0).fit(X)
    cos = np.sum(p.components_ * Vt[:5], axis=1)            # +1 or -1 if the axes agree
    print(f"sklearn {solver:16s} top-3 eigenvalues match: {np.allclose(p.explained_variance_[:3], lam[:3])}"
          f"  cos with NumPy axes 1-5: {np.round(cos, 3)}")
big = lambda V: np.sign(V[np.arange(len(V)), np.abs(V).argmax(axis=1)]).astype(int)
print("sign of each axis's largest entry: NumPy", big(Vt[:5]), " sklearn", big(p.components_))
d = H * W
print("noise variance for k=3: Tipping-Bishop eq. (8), mean of d-3 eigenvalues = %.6f;"
      " sklearn noise_variance_ = %.6f" % (lam[3:].sum() / (d - 3), PCA(3, svd_solver="full").fit(X).noise_variance_))
print("PCA(n_components=0.95) keeps", PCA(n_components=0.95, svd_solver="full").fit(X).n_components_, "components")

# The same three factors, but the disk also moves by up to 3 pixels left or right.
Xs = render(n, shift=3)
ls = np.linalg.svd(Xs - Xs.mean(0), compute_uv=False) ** 2 / (n - 1)
print("with shifts, top 8 eigenvalues:", np.round(ls[:8], 5))
cum = np.cumsum(ls) / ls.sum()
print("components for 95% / 99% of the variance: fixed disk",
      np.searchsorted(np.cumsum(lam) / lam.sum(), [0.95, 0.99]) + 1, " shifted disk",
      np.searchsorted(cum, [0.95, 0.99]) + 1)

Output:

top 6 eigenvalues (SVD): [7.85934 1.47899 0.68904 0.00166 0.00165 0.00163]
top 6 eigenvalues (eigh): [7.85934 1.47899 0.68904 0.00166 0.00165 0.00163]
noise variance per pixel: 0.0004
|cos| between SVD and eigh axes 1-3: [1. 1. 1.]
k= 1  sum ||x - xhat||^2 / (n-1) = 2.39643   sum of discarded eigenvalues = 2.39643   explained = 0.7663
k= 2  sum ||x - xhat||^2 / (n-1) = 0.91744   sum of discarded eigenvalues = 0.91744   explained = 0.9105
k= 3  sum ||x - xhat||^2 / (n-1) = 0.22840   sum of discarded eigenvalues = 0.22840   explained = 0.9777
k= 4  sum ||x - xhat||^2 / (n-1) = 0.22674   sum of discarded eigenvalues = 0.22674   explained = 0.9779
k=10  sum ||x - xhat||^2 / (n-1) = 0.21711   sum of discarded eigenvalues = 0.21711   explained = 0.9788
sklearn full             top-3 eigenvalues match: True  cos with NumPy axes 1-5: [ 1. -1. -1.  1.  1.]
sklearn covariance_eigh  top-3 eigenvalues match: True  cos with NumPy axes 1-5: [ 1. -1. -1.  1.  1.]
sklearn randomized       top-3 eigenvalues match: True  cos with NumPy axes 1-5: [ 1.    -1.    -1.     0.594 -0.771]
sign of each axis's largest entry: NumPy [ 1 -1 -1  1  1]  sklearn [1 1 1 1 1]
noise variance for k=3: Tipping-Bishop eq. (8), mean of d-3 eigenvalues = 0.000399; sklearn noise_variance_ = 0.000460
PCA(n_components=0.95) keeps 3 components
with shifts, top 8 eigenvalues: [7.74475 1.83017 0.85062 0.53672 0.21711 0.13903 0.08857 0.06693]
components for 95% / 99% of the variance: fixed disk [  3 114]  shifted disk [ 5 98]

The spectrum has three large eigenvalues and then drops by a factor of more than 400 to a floor, so it reveals the three generating factors without being told about them. The floor sits near $0.0016$ rather than at the noise variance of $0.0004$: with only 500 images of 576 pixels, the sample eigenvalues of pure noise spread out around the true value, the largest of them several times above it (our reading of the numbers; at most 499 eigenvalues are non-zero here). Both routes give the same eigenvalues and axes on this well-conditioned problem. The reconstruction error equals the sum of the discarded eigenvalues for every $k$, the identity derived above; beyond $k = 3$ it barely falls, because what remains is noise, $0.0004$ per pixel in 573 discarded dimensions, about 0.229.

scikit-learn’s exact solvers agree with NumPy up to sign. Two of the first three axes come out flipped, because scikit-learn 1.7.2 makes the largest entry of every component positive, while LAPACK’s signs are arbitrary. The randomized solver finds the three signal axes exactly but returns other vectors for components 4 and 5, which lie in the noise floor where the eigenvalues are nearly equal and the axes are not well determined. The noise estimates show the difference noted under probabilistic PCA: Tipping and Bishop’s formula, averaged over all 573 discarded dimensions, recovers the true $0.0004$, while scikit-learn’s, averaged over the smallest 497 of the 500 eigenvalues it computes, gives $0.00046$. The variance thresholds behave as described under Parameters: 95% picks exactly the three factors, 99% reaches deep into the noise.

When the disk also moves by up to three pixels, the fourth eigenvalue rises from the noise floor to 0.54 and several more follow it, the spreading of a translation over many components described under Intuition.

Whitening image patches

The script cuts 20,000 random $8 \times 8$ patches from a synthetic texture (blurred noise plus bright rectangles and sensor noise of standard deviation 0.01), whitens them with the PCA and ZCA transforms, with and without a regularizing constant, and prints the ZCA filter of one pixel and two principal axes.

import numpy as np
from scipy.ndimage import gaussian_filter
from sklearn.decomposition import PCA
np.set_printoptions(precision=2, suppress=True)

# A 256 x 256 synthetic texture: white noise blurred with a Gaussian (sigma = 2 px),
# so neighbouring pixels are strongly correlated, plus bright rectangles and pixel noise.
rng = np.random.default_rng(1)
img = gaussian_filter(rng.normal(size=(256, 256)), 2.0)
img = (img - img.min()) / (img.max() - img.min())
for _ in range(15):
    r, c = rng.integers(0, 230, 2)
    img[r:r + rng.integers(8, 26), c:c + rng.integers(8, 26)] += 0.5
img += rng.normal(0, 0.01, img.shape)                       # sensor noise, sigma = 0.01

# 20,000 random 8 x 8 patches, one per row (d = 64).
P = 8
rows, cols = rng.integers(0, 256 - P, (2, 20000))
X = np.array([img[r:r + P, c:c + P].ravel() for r, c in zip(rows, cols)])
mu = X.mean(axis=0)
Xc = X - mu
n = len(X)
C = Xc.T @ Xc / (n - 1)
lam, E = np.linalg.eigh(C)
lam, E = lam[::-1], E[:, ::-1]
print("eigenvalues 1, 2, 10, 30, 64:", ", ".join(f"{v:.2g}" for v in lam[[0, 1, 9, 29, 63]]),
      "  ratio largest/smallest: %.0f" % (lam[0] / lam[-1]))
print("correlation of horizontally adjacent pixels: %.3f" % np.corrcoef(Xc[:, 27], Xc[:, 28])[0, 1])

def whiteners(eps):
    W_pca = np.diag(1 / np.sqrt(lam + eps)) @ E.T           # rotate to the principal axes, rescale
    W_zca = E @ W_pca                                       # ... and rotate back: (C + eps I)^(-1/2)
    return W_pca, W_zca

Xs = Xc / Xc.std(axis=0)                                    # standardized patches, for comparison
cosine = lambda A, B: np.mean(np.sum(A * B, 1) / (np.linalg.norm(A, axis=1) * np.linalg.norm(B, axis=1)))
for eps in [0.0, 1e-3]:
    for name, Wm in zip(["PCA", "ZCA"], whiteners(eps)):
        Y = Xc @ Wm.T
        dev = np.abs(np.cov(Y, rowvar=False) - np.eye(64)).max()
        print(f"eps={eps:g} {name} whitening: max |cov - I| = {dev:.3f}   W symmetric: {str(np.allclose(Wm, Wm.T)):5s}"
              f"   mean cosine to the standardized patch: {cosine(Y, Xs):.3f}")

W_pca, W_zca = whiteners(0.0)
k = W_zca[3 * P + 3].reshape(P, P)                          # ZCA filter that outputs pixel (3, 3)
print("ZCA filter for pixel (3, 3), rows/cols 1-5, divided by its centre value:\n", k[1:6, 1:6] / k[3, 3])
print("principal axes 1 and 4 (columns of E), first two rows of the 8 x 8 patch:")
print(E[:, 0].reshape(P, P)[:2], "\n", E[:, 3].reshape(P, P)[:2])

# scikit-learn's whitening keeps the PCA rotation (no rotation back).
Ys = PCA(whiten=True, svd_solver="full").fit_transform(X)
print("sklearn PCA(whiten=True): max |cov - I| = %.4f" % np.abs(np.cov(Ys, rowvar=False) - np.eye(64)).max())

Output:

eigenvalues 1, 2, 10, 30, 64: 1.2, 0.23, 0.016, 0.00051, 0.00011   ratio largest/smallest: 10771
correlation of horizontally adjacent pixels: 0.937
eps=0 PCA whitening: max |cov - I| = 0.000   W symmetric: False   mean cosine to the standardized patch: 0.003
eps=0 ZCA whitening: max |cov - I| = 0.000   W symmetric: True    mean cosine to the standardized patch: 0.470
eps=0.001 PCA whitening: max |cov - I| = 0.903   W symmetric: False   mean cosine to the standardized patch: 0.004
eps=0.001 ZCA whitening: max |cov - I| = 0.633   W symmetric: True    mean cosine to the standardized patch: 0.670
ZCA filter for pixel (3, 3), rows/cols 1-5, divided by its centre value:
 [[ 0.01  0.03 -0.11  0.03  0.01]
 [ 0.02  0.03 -0.24  0.03  0.02]
 [-0.1  -0.24  1.   -0.26 -0.11]
 [ 0.03  0.01 -0.26  0.02  0.04]
 [ 0.02  0.02 -0.11  0.03  0.01]]
principal axes 1 and 4 (columns of E), first two rows of the 8 x 8 patch:
[[0.1  0.11 0.12 0.12 0.12 0.12 0.11 0.1 ]
 [0.11 0.12 0.13 0.13 0.13 0.13 0.12 0.11]] 
 [[ 0.16  0.16  0.12  0.04 -0.06 -0.16 -0.23 -0.25]
 [ 0.15  0.16  0.13  0.06 -0.03 -0.13 -0.2  -0.22]]
sklearn PCA(whiten=True): max |cov - I| = 0.0000

Neighbouring pixels are strongly correlated, and the eigenvalues fall by four orders of magnitude, from the smooth patterns that carry most of the patch variance down to a floor of about $10^{-4}$, the variance of the added sensor noise. Without a regularizer both transforms produce an exactly white covariance, but they do it differently. The PCA-whitened vectors have almost no resemblance to the patches (a mean cosine of 0.003 with the standardized patch), because their coordinates are amounts of global patterns: the first principal axis is a nearly uniform patch and the fourth a horizontal ramp across the whole patch, the global filters Bell and Sejnowski described. The ZCA matrix is symmetric, its outputs stay pixel-aligned (cosine 0.47), and the filter that produces one output pixel is a centre-surround kernel: a positive centre and negative neighbours, strongest along the row and column, that subtracts what the neighbours predict about the pixel. Adding $10^{-3}$ to the eigenvalues, ten times the noise floor, deliberately leaves the covariance non-white: directions whose variance is below the constant are no longer inflated, so the noise is not amplified, and the ZCA output resembles the image even more. scikit-learn’s whiten=True produces the PCA form, exactly white.

An outlier and a unit mistake

The third script renders 300 RGB images of a disk whose colour varies along one known direction (more red and less blue, or the reverse) under a small random brightness change and noise. It measures how well the first principal axis matches the true colour direction after adding overexposed frames, after storing the blue channel in 0–255 instead of 0–1, and after standardizing every pixel.

import numpy as np
np.set_printoptions(precision=3, suppress=True)

# 300 synthetic 16 x 16 RGB images (values in [0, 1]) of a disk whose colour varies
# along one known direction (more red, less blue), under a random brightness, plus noise.
rng = np.random.default_rng(2)
n, H, W = 300, 16, 16
yy, xx = np.mgrid[0:H, 0:W]
disk = ((yy - 8) ** 2 + (xx - 8) ** 2 < 5 ** 2)[..., None]

def render(n):
    t = rng.uniform(-1, 1, n)                                # the colour factor
    b = rng.uniform(-0.05, 0.05, n)                          # small brightness change
    col = np.array([0.5, 0.4, 0.5]) + t[:, None] * np.array([0.2, 0.0, -0.2])
    imgs = 0.3 + b[:, None, None, None] + disk * (col[:, None, None, :] - 0.3)
    return (imgs + rng.normal(0, 0.02, imgs.shape)).reshape(n, -1)

def first_axis(X):
    Xc = X - X.mean(axis=0)
    _, s, Vt = np.linalg.svd(Xc, full_matrices=False)
    return Vt[0], s ** 2 / (s @ s)                            # axis, explained-variance ratios

X = render(n)
v_clean, r_clean = first_axis(X)
true_dir = (disk * np.array([0.2, 0.0, -0.2])).ravel()
true_dir /= np.linalg.norm(true_dir)
print("clean:           |cos(axis 1, true colour direction)| = %.3f   explained ratio = %.3f"
      % (abs(v_clean @ true_dir), r_clean[0]))

# Outliers: overexposed frames in which every pixel saturated at 1.0.
for m in [1, 3]:
    Xo = np.vstack([X, np.ones((m, X.shape[1]))])
    v, r = first_axis(Xo)
    print(f"{m} overexposed frame(s): |cos(axis 1, true colour direction)| = {abs(v @ true_dir):.3f}"
          f"   explained ratio = {r[0]:.3f}")

# Units: the same clean images with the blue channel stored as 0-255 instead of 0-1.
Xu = X.reshape(n, H, W, 3).copy()
Xu[..., 2] *= 255
Xu = Xu.reshape(n, -1)
v, r = first_axis(Xu)
print("blue in 0-255:   |cos(axis 1, true colour direction)| = %.3f   blue channel's share of axis 1 = %.3f"
      % (abs(v @ true_dir), np.sum(v.reshape(H, W, 3)[..., 2] ** 2)))
# Standardizing every pixel (correlation-matrix PCA) removes the unit problem, since
# it gives the same Z for 0-1 and 0-255 data, but it also changes which directions win.
Z = (X - X.mean(0)) / X.std(0)
print("same Z for both units:", np.allclose(Z, (Xu - Xu.mean(0)) / Xu.std(0)))
v, r = first_axis(Z)
v = v * X.std(0)                                             # back to pixel units to compare
v /= np.linalg.norm(v)
flat = np.ones(X.shape[1]) / np.sqrt(X.shape[1])             # every pixel brighter by the same amount
print("standardized:    |cos(axis 1, true colour direction)| = %.3f   |cos(axis 1, uniform brightening)| = %.3f"
      % (abs(v @ true_dir), abs(v @ flat)))

Output:

clean:           |cos(axis 1, true colour direction)| = 0.999   explained ratio = 0.662
1 overexposed frame(s): |cos(axis 1, true colour direction)| = 0.957   explained ratio = 0.479
3 overexposed frame(s): |cos(axis 1, true colour direction)| = 0.014   explained ratio = 0.644
blue in 0-255:   |cos(axis 1, true colour direction)| = 0.703   blue channel's share of axis 1 = 1.000
same Z for both units: True
standardized:    |cos(axis 1, true colour direction)| = 0.092   |cos(axis 1, uniform brightening)| = 0.995

On clean data the first axis is the colour direction. One overexposed frame among 301 already tilts it, and three frames, 1% of the data, take it over completely: the first axis now points from the mean towards the white frames, nearly orthogonal to the colour direction. The explained-variance ratio gives no warning, since it is about as high as on clean data. The reason is that PCA maximizes a sum of squares, and each white frame lies far from the mean in every pixel at once.

Storing blue on a 255-times larger scale makes blue variance dominate the total, so the first axis lies entirely in the blue channel; its cosine with the true direction, $0.703 \approx 1/\sqrt{2}$, is just the blue half of a red–blue direction. Standardizing every pixel removes the unit problem, since the standardized data are identical for both encodings, but it raises the background pixels to the same weight as the disk; because brightness moves all 768 values together, the first standardized axis becomes uniform brightening, and the colour factor that dominated the raw data is no longer first. For pixel data in common units, rescaling the mis-scaled channel would have been the repair.

Properties and Behavior

Sign and basis ambiguity. An eigenvector is defined only up to sign: if $\mathbf{v}$ is a principal axis, so is $-\mathbf{v}$, with the scores negated. Different libraries, versions and solvers can return either, as the first example shows; scikit-learn 1.7.2 fixes the sign by making the entry of largest magnitude in each component positive (the svd_flip step in its source). When two eigenvalues are equal or nearly equal, as for the noise components in the first example, any orthonormal basis of their joint eigenspace is equally valid, and individual axes are unstable even though the subspace is well defined (Bishop, 2006, §12.1.2, notes that with equal eigenvalues any choice of principal direction gives the same error). Code that compares components across runs, or attaches meaning to the sign of a score, has to align them first; Point Cloud shows the same issue for normals estimated by PCA.

Uncorrelated, not independent. Scores along different axes are uncorrelated, which implies independence when the data are jointly Gaussian but not in general. Bell and Sejnowski (1997) contrasted PCA and ZCA, which only decorrelate, with independent component analysis, which seeks outputs that are independent, and found the independent components of natural images to be localized, oriented edge filters.

Sensitivity to outliers and to units. Because PCA minimizes squared distances, a few gross outliers can rotate the leading axes, as the overexposed frames did; detecting and removing such samples, or using a robust estimate of the covariance, is the remedy (our recommendation). The dependence on units is discussed under standardization.

Linear factors only. PCA finds a linear subspace. Factors that act non-linearly on pixels, such as translation, rotation, scale or viewpoint, spread over many components. Ke and Sukthankar (2004) describe PCA as very sensitive to registration error and apply it only to patches that SIFT’s detector has already aligned in position, scale and orientation; in their tests PCA-SIFT tolerated added position errors of 1 and 3 pixels but degraded when the orientation or scale was perturbed.

Unsupervised. PCA keeps directions of large variance whatever causes them. Bishop (2006, §12.1.3, Figure 12.7) shows a two-class example where the direction of maximum variance mixes the classes and Fisher’s linear discriminant, which uses the labels, separates them; Face Recognition describes the same contrast between eigenfaces and Fisherfaces.

Relation to linear autoencoders. A linear autoencoder trained with squared error spans the same principal subspace, though its units are mixtures of the components; that article derives and demonstrates the equivalence.

Limitations

The behaviours above are also its limitations: one global linear subspace, squared error, dependence on units and on alignment, and no use of labels. Two more matter in practice:

• Variance is not relevance. The leading components describe what varies most, which for images is often lighting, not the factor a task needs.
• Dense components. Every component mixes all the original variables, which can make it hard to interpret; sparse PCA, available in scikit-learn as SparsePCA, adds an $\ell_1$ penalty to obtain components with many zeros (scikit-learn 1.7.2 User Guide, §2.5.1.4).

Variants

• Kernel PCA: replaces inner products between samples by a nonlinear kernel and solves an eigenvalue problem for the $n \times n$ kernel matrix, which yields nonlinear components; Bishop (2006, §12.3) attributes it to Schölkopf, Smola and Müller (1998) and notes that it needs the eigenvectors of an $n \times n$ matrix, so large data sets need approximations.
• Factor analysis: replaces PCA’s isotropic noise $\sigma^2 I$ with a diagonal noise covariance, so that each variable has its own noise level. Tipping and Bishop (1999, §2–§3.3) explain that this makes factor analysis covariant under rescaling of individual variables, while PCA and probabilistic PCA are covariant under rotations of the data axes.
• Independent component analysis: seeks statistically independent rather than uncorrelated components (Bell and Sejnowski, 1997).

Uses in Computer Vision

Eigenfaces. Sirovich and Kirby (1987) characterized a picture of a face, to within an error bound, by a relatively low-dimensional vector (according to their abstract). Turk and Pentland (1991) built a near-real-time recognition system on the eigenvectors (principal components) of a set of face images, which they called eigenfaces and which need not correspond to features such as eyes or noses (their abstract); Face Recognition describes how recognition worked and what followed. With images of thousands of pixels and hundreds of training faces, the Gram-matrix route above keeps the computation small (our remark).

PCA-SIFT. Ke and Sukthankar (2004) kept SIFT’s keypoint detector and replaced its histogram descriptor with PCA of the normalized gradient patch. At each keypoint they extracted a $41 \times 41$ patch at the keypoint’s scale, rotated to its dominant orientation, concatenated the horizontal and vertical gradient maps into a vector of $2 \times 39 \times 39 = 3{,}042$ elements, and normalized it to unit length. The projection matrix was the leading eigenvectors of the covariance of 21,000 such patches collected from a separate, diverse set of images. Most of their results used 20 components instead of SIFT’s 128 values. They observed that the eigenvalues of gradient patches around keypoints decay much faster than those of randomly chosen patches, because the detector has already normalized position, scale and orientation, and that among the dimensions they compared, matching was best with 36 components and then declined: with all 3,042 it was only slightly better than with 12 (their Figure 7). Their explanation was that the first components encode what identifies the keypoint, and the later ones encode detail such as distortions from projective warps that does not help matching. They reported that the descriptor was more distinctive and more compact than SIFT’s and that matching was faster. Later evaluations, summarized in Feature Descriptors and the HOG paper, compared it with other descriptors.

Descriptor compression and whitening for retrieval. Jégou and Chum showed for bag-of-words and VLAD vectors that PCA helps retrieval through centring and whitening, not only by shortening the vectors, and that renormalizing the whitened vectors helps; the same chain was then applied to CNN descriptors. Embeddings covers that analysis and its results, and Image Retrieval the descriptors that use it, such as R-MAC. PCA also shortens descriptors before encoding, as in the Fisher-vector encoding of dense trajectories, and aligns randomized k-d trees with the principal axes of the data in approximate nearest-neighbour search.

Colour augmentation. AlexNet added to each training image random multiples of the principal components of the training set’s RGB pixel values, scaled by their eigenvalues, to imitate changes in the intensity and colour of illumination; the AlexNet paper and Data Augmentation articles give the details and its effect.

Geometry. The normal of a surface at a point of a point cloud is the eigenvector with the smallest eigenvalue of the covariance of its neighbours, used for example in point-to-plane iterative closest point.

History

Pearson (1901) posed the geometric problem for points whose coordinates all carry measurement error: fit the line or plane of closest fit by minimizing the sum of squared perpendicular distances, rather than the distances along one coordinate that regression minimizes when it treats one variable as known exactly. He showed that the best line passes through the centroid and runs along the major axis of the correlation ellipse, and solved the general case through the axes of an “ellipsoid of residuals” built from the standard deviations and correlations of the variables. Hotelling (1933), analysing correlated test scores, asked for a first component whose contributions to the variances of the variables have as great a total as possible, then a second, independent of the first, that takes as much of the residual variance as possible, and so on. He named this the method of principal components, identified its axes with the principal axes of the ellipsoids of uniform density of a normal population, and proposed an iterative computation: applying the correlation matrix repeatedly to a trial vector converges to the axis with the largest root, and subtracting each component’s contribution from the matrix lets the process find the next. He called it a vast saving of arithmetical effort over solving the characteristic equation; it is the power method with deflation. In vision, Sirovich and Kirby’s low-dimensional representation of face pictures (1987) and Turk and Pentland’s eigenfaces (1991) came first; whitened image patches for feature learning (1997), PCA-SIFT (2004) and the whitening of retrieval descriptors followed.