Dimensionality Reduction

Dimensionality reduction replaces each example, an image of $D$ pixels or a descriptor of $D$ numbers, by a shorter vector of $k \ll D$ numbers that keeps what a later step needs. There are two ways to do it. Feature selection keeps $k$ of the original variables and discards the rest; feature extraction (Guyon and Elisseeff call it feature construction) computes $k$ new variables from all of them, as principal component analysis does with linear combinations. The reduced vectors are used to compress and store data, to make a learner faster or more accurate, to measure distances that mean something, and to draw a data set in two or three dimensions. This article lays out the taxonomy of methods, covers feature selection and the classical manifold methods (multidimensional scaling, Isomap and locally linear embedding), explains how the target dimension is chosen, and runs an example in which the directions of largest variance contain none of the signal a classifier needs.

Definition

Given examples $\mathbf{x}_1, \dots, \mathbf{x}_N \in \mathbb{R}^D$ , dimensionality reduction produces $\mathbf{z}_1, \dots, \mathbf{z}_N \in \mathbb{R}^k$ with $k < D$ . Murphy (2022, ch. 20) distinguishes two forms of the output:

a parametric map $\mathbf{z} = f(\mathbf{x}; \theta)$ that can be applied to any input, including images that were not in the training set;
a nonparametric mapping, an embedding that gives a vector $\mathbf{z}_n$ for each training example and for no other point.

He notes that the second form is used mostly for visualization, while the first can also serve as a preprocessing step for a later learner, for example a linear classifier trained on $\mathbf{z}$ instead of $\mathbf{x}$ .

Guyon and Elisseeff (2003, footnote 1) call the raw inputs variables and quantities constructed from them features, and use the two words interchangeably when the distinction does not matter. In images the variables are pixels; a selection method returns a set of pixel indices, an extraction method a matrix or a function.

Intuition

A $16 \times 16$ grey image is a point in a 256-dimensional space, but the images a camera produces of one scene type do not fill that space. The pixels are tied together by the objects, the lighting and the optics, so the set of images that actually occurs has far fewer degrees of freedom than the pixel count. Murphy (2022, §20.4.2) calls the assumption that natural high-dimensional data lie near a low-dimensional manifold the manifold hypothesis, and illustrates it with 360 rotated copies of one digit image that PCA projects onto a circle. Curse of Dimensionality explains why the dimension that matters for learning is this intrinsic one, not the number of pixels.

What “keeps what is needed” means depends on the purpose. Reconstruction asks that $\mathbf{x}$ can be rebuilt from $\mathbf{z}$ ; distance preservation asks that $\lVert \mathbf{z}_i - \mathbf{z}_j \rVert$ track $\lVert \mathbf{x}_i - \mathbf{x}_j \rVert$ or a distance along the manifold; prediction asks that a label can still be read from $\mathbf{z}$ . Guyon and Elisseeff (2003, §5) separate the first goal, an unsupervised problem close to data compression, from the third, a supervised one. The goals can disagree: the directions along which images vary most can be lighting and exposure, which a recognition task wants to ignore. The example below builds such a case on purpose.

Formal Definition

Selection and extraction

Feature selection chooses an index set $S \subset \{1, \dots, D\}$ with $|S| = k$ and keeps $\mathbf{z} = \mathbf{x}_S$ . Since the number of subsets grows combinatorially, selection is a search problem; Guyon and Elisseeff (2003, §4.1) note that finding the best subset is NP-hard, citing Amaldi and Kann (1998). Linear feature extraction chooses a $D \times k$ matrix $W$ and sets $\mathbf{z} = W^\top(\mathbf{x} - \bar{\mathbf{x}})$ ; selection is the special case in which the columns of $W$ are coordinate axes. Nonlinear extraction replaces $W^\top$ by a nonlinear function, either learned as a network (autoencoder) or defined implicitly by an embedding of the training set.

Reconstruction: PCA

Murphy (2022, §20.1) defines PCA by the reconstruction error of a linear encoder and decoder,

\mathcal{L}(W) = \frac{1}{N} \sum_{n=1}^{N} \lVert \mathbf{x}_n - W W^\top \mathbf{x}_n \rVert_2^2 ,

for centred data and $W^\top W = I$ , minimized by the $k$ leading eigenvectors of the sample covariance. Principal Component Analysis derives this, its equivalence with maximum variance, whitening and probabilistic PCA; a linear autoencoder trained with squared error spans the same subspace, as Autoencoder shows.

Distances: classical multidimensional scaling

Multidimensional scaling (MDS) looks for points $\mathbf{z}_i$ whose pairwise distances match given dissimilarities $d_{ij}$ , which need not come from vectors at all (Murphy, 2022, §20.4.4). Classical MDS works with the centred Gram matrix $\tilde{K}_{ij} = \langle \mathbf{x}_i - \bar{\mathbf{x}}, \mathbf{x}_j - \bar{\mathbf{x}} \rangle$ and minimizes the strain $\lVert \tilde{K} - \tilde{Z}\tilde{Z}^\top \rVert_F^2$ (Murphy’s eq. 20.105). When only Euclidean distances are known, $\tilde{K}$ is recovered by double centring the squared distances (eq. 20.109):

\tilde{K}_{ij} = -\frac{1}{2}\Big( d_{ij}^2 - \frac{1}{N}\sum_{l=1}^{N} d_{il}^2 - \frac{1}{N}\sum_{l=1}^{N} d_{jl}^2 + \frac{1}{N^2}\sum_{l=1}^{N}\sum_{m=1}^{N} d_{lm}^2 \Big),

and the embedding is $\tilde{Z} = U_k S_k^{1/2}$ from the top $k$ eigenvectors and eigenvalues of $\tilde{K}$ . Murphy shows that on Euclidean distances this embedding equals the PCA scores. Metric MDS generalizes to any dissimilarity by minimizing the stress $\sqrt{\sum_{i<j} (d_{ij} - \hat{d}_{ij})^2 / \sum_{ij} d_{ij}^2}$ , with $\hat{d}_{ij} = \lVert \mathbf{z}_i - \mathbf{z}_j \rVert$ (eq. 20.111); non-metric MDS matches only the rank order of the dissimilarities, alternating isotonic regression and gradient steps; the Sammon mapping divides each squared error by $d_{ij}$ so that small distances count more, which Murphy describes as non-convex and as putting arguably too much weight on getting very small distances exactly right (§20.4.4.2–§20.4.4.4).

Neighbourhood graphs: Isomap, LLE and Laplacian eigenmaps

Straight-line distances cut across a curved manifold: two points on neighbouring turns of a rolled sheet are close in the ambient space and far apart along the sheet. Isomap (Tenenbaum, de Silva and Langford, 2000, as Murphy, 2022, §20.4.5 describes it) builds the $K$ -nearest-neighbour graph of the data, replaces each distance by the shortest path through that graph (computed with Dijkstra’s algorithm), and applies classical MDS to the resulting geodesic distances. Murphy reports that noise can create false edges in the graph, “short circuits” that significantly distort the embedding, and that choosing a very small neighbourhood does not solve this, since it can fragment the manifold into many disconnected regions.

Locally linear embedding (LLE; Roweis and Saul, 2000, as Murphy, 2022, §20.4.8 describes it) assumes that the manifold is locally linear. It first writes each point as an affine combination of its $K$ nearest neighbours,

\hat{W} = \arg\min_{W} \sum_{i=1}^{N} \Big\lVert \mathbf{x}_i - \sum_{j} w_{ij}\mathbf{x}_j \Big\rVert^2 \quad \text{subject to } w_{ij} = 0 \text{ unless } \mathbf{x}_j \text{ is a neighbour of } \mathbf{x}_i, \quad \sum_j w_{ij} = 1,

and then finds low-dimensional points that the same weights reconstruct, $\min_Z \sum_i \lVert \mathbf{z}_i - \sum_j \hat{w}_{ij}\mathbf{z}_j \rVert^2$ , whose solution is given by the eigenvectors of $(I - \hat{W})^\top (I - \hat{W})$ with the smallest nonzero eigenvalues (eqs. 20.117–20.120). The sum-to-one constraint rules out the trivial solution $W = 0$ . Because each point uses only its neighbours, the eigenproblem is sparse; Murphy found LLE’s results on his running example less good than Isomap’s but notes that it tends to be somewhat less sensitive to short-circuiting.

Laplacian eigenmaps (spectral embedding) minimize $\sum_{(i,j)} W_{ij}\lVert \mathbf{z}_i - \mathbf{z}_j \rVert^2$ over a weighted neighbour graph under a normalization constraint, which leads to eigenvectors of the graph Laplacian (Murphy, 2022, §20.4.9); Spectral Clustering covers the Laplacian and its eigenvectors. Kernel PCA, which Murphy groups with these methods, is covered in Kernel Methods, and t-SNE matches neighbour probabilities rather than distances. Murphy also describes maximum variance unfolding, which maximizes the spread of the embedding subject to keeping the distances between graph neighbours, posed as a semidefinite program (§20.4.7).

Properties

A taxonomy

The methods differ along several independent axes. The table places the common ones; each linked article owns its method.

Method	Selection or extraction	Linear	Uses labels	Out-of-sample map	Geometry kept
Univariate filters (F-test, mutual information)	selection	yes (coordinate axes)	yes	yes	single variables
Wrappers and embedded selection (RFE, sequential selection, lasso)	selection	yes	yes	yes	what the predictor uses
PCA	extraction	yes	no	yes	global variance
Linear Discriminant Analysis	extraction	yes	yes	yes	class separation, at most $c - 1$ directions
Random projection (Curse of Dimensionality)	extraction	yes	no	yes	pairwise distances, approximately
Classical MDS	extraction	yes on Euclidean distances (equals PCA)	no	none in scikit-learn (its `MDS` runs SMACOF and has no `transform`)	global distances
Kernel PCA	extraction	no	no	yes	variance in a feature space
Isomap	extraction	no	no	approximate	geodesic distances
LLE, Laplacian eigenmaps	extraction	no	no	LLE yes, spectral embedding no	local neighbourhoods
t-SNE	extraction	no	no	no	neighbour probabilities
Autoencoder	extraction	no	no	yes	reconstruction

The out-of-sample column follows scikit-learn 1.7.2’s source: TSNE and SpectralEmbedding have no transform method; Isomap.transform links new points into the training set’s neighbour graph and projects the resulting geodesic distances onto the training embedding; LocallyLinearEmbedding has a transform; and MDS is the iterative SMACOF algorithm, metric or non-metric, with only fit_transform. Murphy (2022, §20.4.3, Table 20.1) adds a computational axis: kernel PCA and Isomap solve dense eigenproblems, LLE and Laplacian eigenmaps sparse ones, while t-SNE and autoencoders have non-convex objectives. In Linear Discriminant Analysis the number of directions is bounded by the number of classes: scikit-learn’s n_components must be at most $\min(c - 1, D)$ , one direction for two classes.

Feature selection: filters, wrappers and embedded methods

Guyon and Elisseeff (2003, §4) divide subset selection methods into three families. They call variable ranking a filter following the classification of Kohavi and John (1997) (§2.1), and credit Kohavi and John with popularizing wrappers (§4.1):

Filters score variables, or subsets of variables, as a preprocessing step that does not depend on the predictor. The simplest filter ranks single variables by a score $S(i)$ and keeps nested subsets of the top-ranked ones (§2.1). Scores include the squared correlation with the target, which for two classes is closely related to Fisher’s criterion and the t-test, single-variable classifiers, and the mutual information between the variable and the target (§2.2–§2.4). The authors credit ranking with computational scalability (it needs $D$ scores and a sort) and statistical robustness: it introduces bias but may have considerably less variance than subset search (§2.1).
Wrappers treat the learning machine as a black box and score candidate subsets by its validated performance (§4.1). Exhaustive search quickly becomes intractable. The authors single out greedy forward selection and backward elimination: both produce nested subsets, and they write that greedy search strategies seem to be computationally advantageous and robust against overfitting.
Embedded methods select variables while the model is trained (§4.1–§4.3). Their examples are decision trees such as CART, which choose variables as they split, backward elimination by the smallest weight magnitude in linear SVMs, and objectives with a sparsity penalty, from the $\ell_1$ norm to approximations of the $\ell_0$ count. Linear Regression shows the lasso selecting pixels, and in the Viola–Jones detector each round of boosting selects one Haar feature from a large pool.

scikit-learn 1.7.2 provides each family (User Guide §1.13): VarianceThreshold removes low-variance features without looking at labels; SelectKBest and SelectPercentile keep the best-scoring features under f_classif (the ANOVA F-value), chi2 (for non-negative features) or mutual_info_classif; RFE and RFECV eliminate features recursively by a fitted model’s coefficients or importances; SelectFromModel keeps features above an importance threshold, for example the nonzero weights of an $\ell_1$ -penalized model; and SequentialFeatureSelector runs forward or backward selection scored by cross-validation. The user guide notes that F-tests measure only linear dependency, while mutual information captures any dependency but, being nonparametric, needs more samples. In the 1.7.2 source, mutual_info_classif uses a $k$ -nearest-neighbour estimator with n_neighbors=3, rescales continuous features to unit variance (scale(..., with_mean=False)) and adds noise of $10^{-10}$ times each rescaled feature’s mean absolute value (at least 1), so it takes a random_state; its scores are in nats. Every selection step is part of the learning procedure and has to be refit inside each fold of cross-validation; Training, Validation, and Test Sets shows how badly a selection made before splitting can mislead.

Variables that only help together

Guyon and Elisseeff (2003, §3) give small examples that limit what single-variable ranking can see. Two variables with identical class-conditional distributions are not redundant when their noise is independent: their average separates the classes better by a factor $\sqrt{2}$ (their Figure 1). Strong correlation between two variables does not imply that one is redundant (Figure 2). A variable that is “useless by itself” can improve separation when used with another, and in an XOR-like arrangement two variables that are each useless separate the classes together (Figure 3). Their general conclusion is that methods that score variables individually cannot tell which combination performs best, which is the argument for wrappers and embedded methods. A footnote to their check list cautions that it is heuristic and that the only recommendation almost surely valid is “to try the simplest things first” (footnote 2).

Feature extraction can be supervised or not. Guyon and Elisseeff (§5) give two reasons to construct features without labels even for a supervised problem: some problems come with more unlabelled than labelled data, and unsupervised feature selection is less prone to overfitting. The cost, which the example below shows, is that an unsupervised method has no way to know which variation matters.

Choosing the dimension

No single rule exists, and the right $k$ depends on the purpose.

Variance and reconstruction. For PCA, the cumulative explained variance, the scree plot and its elbow, Murphy’s profile likelihood and probabilistic PCA’s likelihood are covered under Parameters in Principal Component Analysis, together with Murphy’s (2022, §20.1.4.1) point that held-out reconstruction error keeps falling as $k$ grows and so cannot choose $k$ .
Intrinsic dimension. Estimators such as Levina and Bickel’s measure how many dimensions the data actually occupy; Curse of Dimensionality explains the estimate and shows how pixel noise inflates it.
Distance distortion. For random projection, the Johnson–Lindenstrauss lemma gives a dimension that keeps all pairwise distances within a chosen factor; scikit-learn’s GaussianRandomProjection draws its matrix entries from $\mathcal{N}(0, 1/k)$ and can set $k$ from that bound (n_components='auto', controlled by eps), which its docstring warns can be very conservative.
The downstream task. When the reduced vectors feed a predictor, $k$ is a hyperparameter: Guyon and Elisseeff’s check list (§1) suggests training a sequence of predictors on increasing subsets of the ranked features and asking whether a smaller subset matches or improves performance, with the subset size chosen by model selection (§6). The validation rules of Cross-Validation apply.

Examples

Selection against extraction when the class signal has low variance

The example asks what each kind of reduction keeps when the image variation that is largest has nothing to do with the class. It uses synthetic $16 \times 16$ images (pixel row $r$ and column $c$ from 0 to 15, $D = 256$ ), with every generating parameter known:

\mathbf{x} = 0.5 + b + g\,x_n(c) + \text{blob} - y\,a\,\text{notch} + \boldsymbol{\varepsilon}, \qquad x_n(c) = \frac{c - 7.5}{7.5}.

Blob: fixed in every image, $0.4\exp(-((r-7.5)^2 + (c-7.5)^2)/(2 \cdot 3^2))$ .
Nuisance: a brightness offset $b \sim \mathcal{N}(0, 0.2^2)$ and a left–right gradient $g \sim \mathcal{N}(0, 0.2^2)$ , drawn independently for each image. They carry no class information and contribute almost all of the pixel variance.
Class: images with $y = 1$ have a $2 \times 2$ notch of depth $a = 0.1$ at rows 6–7, columns 9–10, on the blob’s right flank. The classes are exactly balanced.
Noise: independent $\mathcal{N}(0, 0.05^2)$ per pixel.

Because both classes are Gaussian with the same covariance $\Sigma = 0.04\,\mathbf{1}\mathbf{1}^\top + 0.04\,\mathbf{x}_n\mathbf{x}_n^\top + 0.0025\,I$ and means that differ by $\boldsymbol{\delta} = -a\,\text{notch}$ , the Bayes-optimal rule is linear and its error is exactly $\Phi(-\Delta/2)$ with $\Delta^2 = \boldsymbol{\delta}^\top\Sigma^{-1}\boldsymbol{\delta}$ . The same formula, applied to the projections $W^\top\mathbf{x}$ with $W^\top\Sigma W$ and $W^\top\boldsymbol{\delta}$ , gives the oracle error of any linear representation: the lowest error any classifier could reach from that representation. Comparing a learned error with its oracle error separates what the reduction lost from what the classifier failed to learn.

For each seed (0 to 4), the script draws 2,000 training and 5,000 test images, reduces them to $k \in \{1, 2, 4, 8, 16\}$ dimensions fitted on the training images, and trains the same classifier, standardization followed by LogisticRegression(C=1.0, max_iter=5000), on each representation:

F filter: the $k$ pixels with the largest f_classif score.
MI filter: the $k$ pixels with the largest mutual_info_classif score (n_neighbors=3, random_state = seed).
PCA: the first $k$ of 16 components (svd_solver="full").
Random projection: GaussianRandomProjection(n_components=k, random_state=seed).
LDA: LinearDiscriminantAnalysis(n_components=1), its single direction, printed on a separate row.

The measures are the test error (mean, minimum and maximum over the five seeds), the oracle error of the learned subspace, and the fraction of total pixel variance kept, $\operatorname{tr}(Q^\top S Q)/\operatorname{tr}(S)$ , where $S$ is the covariance of the test images and $Q$ an orthonormal basis of the representation’s span (selected pixel axes, PCA components, random directions, the LDA vector). Per seed it prints the error of the exact Bayes rule on the test images, each filter’s top four pixels, the notch pixels’ F and MI scores against the largest score of any other pixel, and the comparison for H3.

Six checks are printed: C1, the data match their description (exact class counts, the standard deviations of $b$ and $g$ within 0.01 of 0.2, the residual after removing the described terms equal to the drawn noise, with standard deviation within $5 \times 10^{-4}$ of 0.05 and mean within $5 \times 10^{-4}$ of 0); C2, f_classif equals a one-way ANOVA F computed by hand, largest relative difference at most $10^{-8}$ ; C3, PCA’s explained variances equal the top eigenvalues of the training covariance from np.linalg.eigh, to $10^{-8}$ relative; C4, the random-projection entries times $\sqrt{k}$ have mean within 0.05 of 0 and standard deviation within 0.05 of 1; C5, the label-free check that the task is learnable but not trivial: the exact Bayes error is below 0.10 and the Bayes rule’s test error is within three binomial standard errors of it; C6, no learned error is below its oracle error by more than two binomial standard errors.

The hypotheses, each decided on every seed, are:

H1a: PCA keeps the largest fraction of the test variance at every $k$ .
H1b: at $k = 1$ and 2, PCA’s error exceeds both filters’ errors by more than a tolerance.
H1c: at $k = 1$ and 2, no other method’s error exceeds PCA’s by more than the tolerance. Together, H1b and H1c say that PCA has the highest error at small $k$ .
H2: each filter’s top four pixels are exactly the four notch pixels.
H3: LDA’s single direction is as good as any $k$ -dimensional representation: its error is at most the best error of that seed plus the tolerance.

The tolerance is two standard errors of the difference of two test errors, $2\sqrt{(e_1(1-e_1) + e_2(1-e_2))/5000}$ . The script uses NumPy 2.3.5, SciPy 1.16.3 and scikit-learn 1.7.2 with one BLAS thread and runs in about 5 seconds; two runs gave identical output.

import os
for v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"):
    os.environ[v] = "1"                       # one BLAS thread, set before NumPy
import numpy as np
from scipy.stats import norm
from sklearn.feature_selection import f_classif, mutual_info_classif
from sklearn.decomposition import PCA
from sklearn.random_projection import GaussianRandomProjection
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline

SEEDS = range(5)                              # first-run seeds
N_TR, N_TE, KS = 2000, 5000, (1, 2, 4, 8, 16)

# ---- generator: 16 x 16 images, a fixed blob, a strong nuisance, a faint notch
H = W = 16; D = H * W
rr, cc = np.mgrid[0:H, 0:W]
XN = ((cc - 7.5) / 7.5).ravel()                              # ramp from -1 to 1
BLOB = 0.4 * np.exp(-((rr - 7.5)**2 + (cc - 7.5)**2) / (2 * 3.0**2)).ravel()
NOTCH = np.zeros((H, W)); NOTCH[6:8, 9:11] = 1; NOTCH = NOTCH.ravel()
A, S_B, S_G, SIGMA = 0.1, 0.2, 0.2, 0.05
NOTCH_SET = set(np.flatnonzero(NOTCH))

def make(n, rng):
    y = rng.permutation(np.repeat([0, 1], n // 2))
    b = rng.normal(0, S_B, n)                                # brightness offset
    g = rng.normal(0, S_G, n)                                # left-right gradient
    eps = rng.normal(0, SIGMA, (n, D))
    X = 0.5 + b[:, None] + g[:, None] * XN + BLOB - A * y[:, None] * NOTCH + eps
    return X, y, b, g, eps

# exact class-conditional covariance and mean difference, known from the generator
SIG_W = S_B**2 * np.ones((D, D)) + S_G**2 * np.outer(XN, XN) + SIGMA**2 * np.eye(D)
DELTA = -A * NOTCH

def oracle_err(Wm):
    """Exact error of the Bayes rule that sees only the projections Wm^T x."""
    d = Wm.T @ DELTA
    return norm.cdf(-np.sqrt(d @ np.linalg.solve(Wm.T @ SIG_W @ Wm, d)) / 2)

def kept(Wm, S):
    """Fraction of total test-pixel variance kept by projecting onto span(Wm)."""
    Q, _ = np.linalg.qr(Wm)
    return np.trace(Q.T @ S @ Q) / np.trace(S)

def tol(e1, e2):                              # two standard errors of a difference
    return 2 * np.sqrt((e1 * (1 - e1) + e2 * (1 - e2)) / N_TE)

def logreg():
    return make_pipeline(StandardScaler(), LogisticRegression(C=1.0, max_iter=5000))

BAYES = oracle_err(np.eye(D))
print(f"exact Bayes error {BAYES:.5f}; seeds {list(SEEDS)}; "
      f"{N_TR} training and {N_TE} test images per seed")
methods = ("F-filter", "MI-filter", "PCA", "RandProj")
err = {m: {k: [] for k in KS} for m in methods}
orc = {m: {k: [] for k in KS} for m in methods}
var = {m: {k: [] for k in KS} for m in methods}
lda_e, lda_o, lda_v = [], [], []
checks = {f"C{i}": True for i in range(1, 7)}
h = {"H1a": [], "H1b": [], "H1c": [], "H2": [], "H3": []}
top4 = {"F-filter": [], "MI-filter": []}

for seed in SEEDS:
    rng = np.random.default_rng(seed)
    Xtr, ytr, btr, gtr, etr = make(N_TR, rng)
    Xte, yte, bte, gte, ete = make(N_TE, rng)
    # C1: the data match the generator's description
    b_all, g_all = np.r_[btr, bte], np.r_[gtr, gte]
    e_all = np.r_[etr.ravel(), ete.ravel()]
    resid = np.r_[(Xtr - 0.5 - btr[:, None] - gtr[:, None] * XN - BLOB
                   + A * ytr[:, None] * NOTCH).ravel(),
                  (Xte - 0.5 - bte[:, None] - gte[:, None] * XN - BLOB
                   + A * yte[:, None] * NOTCH).ravel()]
    c1 = (np.bincount(ytr).tolist() == [N_TR // 2] * 2
          and np.bincount(yte).tolist() == [N_TE // 2] * 2
          and abs(b_all.std() - S_B) <= 0.01 and abs(g_all.std() - S_G) <= 0.01
          and abs(resid.std() - SIGMA) <= 5e-4 and abs(resid.mean()) <= 5e-4
          and np.allclose(resid, e_all, atol=1e-12))
    checks["C1"] &= c1
    S_te = np.cov(Xte, rowvar=False)

    # filters: scores on the training images
    F, _ = f_classif(Xtr, ytr)
    m0, m1 = Xtr[ytr == 0].mean(0), Xtr[ytr == 1].mean(0)
    n0, n1 = (ytr == 0).sum(), (ytr == 1).sum()
    ssw = ((Xtr[ytr == 0] - m0)**2).sum(0) + ((Xtr[ytr == 1] - m1)**2).sum(0)
    mt = Xtr.mean(0)
    ssb = n0 * (m0 - mt)**2 + n1 * (m1 - mt)**2
    F_hand = ssb / (ssw / (N_TR - 2))
    checks["C2"] &= np.max(np.abs(F_hand - F) / F) <= 1e-8
    MI = mutual_info_classif(Xtr, ytr, n_neighbors=3, random_state=seed)
    order = {"F-filter": np.argsort(-F, kind="stable"),
             "MI-filter": np.argsort(-MI, kind="stable")}
    for m in order:
        top4[m].append(set(order[m][:4].tolist()) == NOTCH_SET)

    pca = PCA(n_components=16, svd_solver="full").fit(Xtr)
    ev = np.linalg.eigh(np.cov(Xtr, rowvar=False))[0][::-1][:16]
    checks["C3"] &= np.max(np.abs(pca.explained_variance_ - ev) / ev) <= 1e-8

    rp_entries = []
    for k in KS:
        rp = GaussianRandomProjection(n_components=k, random_state=seed).fit(Xtr)
        rp_entries.append(rp.components_.ravel() * np.sqrt(k))
        reps = {"F-filter": np.eye(D)[:, order["F-filter"][:k]],
                "MI-filter": np.eye(D)[:, order["MI-filter"][:k]],
                "PCA": pca.components_[:k].T,
                "RandProj": rp.components_.T}
        for m, Wm in reps.items():
            Ztr, Zte = (Xtr - pca.mean_) @ Wm, (Xte - pca.mean_) @ Wm
            e = 1 - logreg().fit(Ztr, ytr).score(Zte, yte)
            o = oracle_err(Wm)
            err[m][k].append(e); orc[m][k].append(o); var[m][k].append(kept(Wm, S_te))
            checks["C6"] &= e >= o - 2 * np.sqrt(o * (1 - o) / N_TE)
    rp_entries = np.concatenate(rp_entries)
    checks["C4"] &= abs(rp_entries.std() - 1) <= 0.05 and abs(rp_entries.mean()) <= 0.05

    lda = LinearDiscriminantAnalysis(n_components=1).fit(Xtr, ytr)
    w = lda.scalings_[:, :1]
    assert np.allclose(lda.transform(Xte), (Xte - lda.xbar_) @ w)
    e = 1 - logreg().fit(lda.transform(Xtr), ytr).score(lda.transform(Xte), yte)
    lda_e.append(e); lda_o.append(oracle_err(w)); lda_v.append(kept(w, S_te))
    checks["C6"] &= e >= lda_o[-1] - 2 * np.sqrt(lda_o[-1] * (1 - lda_o[-1]) / N_TE)

    # C5: the Bayes rule from the known generator, applied to the test images
    wb = np.linalg.solve(SIG_W, DELTA)
    mu0 = 0.5 + BLOB
    s = (Xte - mu0 - DELTA / 2) @ wb
    e_bayes = np.mean((s > 0).astype(int) != yte)
    checks["C5"] &= abs(e_bayes - BAYES) <= 3 * np.sqrt(BAYES * (1 - BAYES) / N_TE)
    checks["C5"] &= BAYES < 0.10
    print(f"seed {seed}: Bayes rule test error {e_bayes:.4f}; "
          f"top-4 F pixels {sorted(divmod(int(i), W) for i in order['F-filter'][:4])}; "
          f"top-4 MI pixels {sorted(divmod(int(i), W) for i in order['MI-filter'][:4])}")
    print(f"   notch F {np.round(F[NOTCH == 1], 2).tolist()}, best other F "
          f"{np.max(F[NOTCH == 0]):.2f}; notch MI (nats) {np.round(MI[NOTCH == 1], 4).tolist()}, "
          f"best other MI {np.max(MI[NOTCH == 0]):.4f}")

    # hypotheses, decided on this seed
    i = len(lda_e) - 1
    h["H1a"].append(all(var["PCA"][k][i] > max(var[m][k][i] for m in methods if m != "PCA")
                        and (k != 1 or var["PCA"][1][i] > lda_v[i]) for k in KS))
    h["H1b"].append(all(err["PCA"][k][i] - err[m][k][i] > tol(err["PCA"][k][i], err[m][k][i])
                        for k in (1, 2) for m in ("F-filter", "MI-filter")))
    h["H1c"].append(all(err[m][k][i] - err["PCA"][k][i] <= tol(err["PCA"][k][i], err[m][k][i])
                        for k in (1, 2) for m in ("F-filter", "MI-filter", "RandProj")))
    h["H2"].append(top4["F-filter"][i] and top4["MI-filter"][i])
    best = min(err[m][k][i] for m in methods for k in KS)
    h["H3"].append(lda_e[i] <= best + tol(best, lda_e[i]))
    print(f"   best k-dim error {best:.4f}, LDA error {lda_e[i]:.4f}, "
          f"tol {tol(best, lda_e[i]):.4f}")

def mmm(a):
    a = np.asarray(a)
    return f"{a.mean():.4f} [{a.min():.4f}, {a.max():.4f}]"

print("\ntest error of logistic regression: mean [min, max] over seeds;"
      " oracle = exact error of the best rule in the same subspace (mean)")
print(f"{'method':10s} {'k':>3s}  {'test error':>24s}  {'oracle':>7s}  {'variance kept':>13s}")
for m in methods:
    for k in KS:
        print(f"{m:10s} {k:3d}  {mmm(err[m][k]):>24s}  {np.mean(orc[m][k]):7.4f}  "
              f"{np.mean(var[m][k]):13.6f}")
print(f"{'LDA':10s} {1:3d}  {mmm(lda_e):>24s}  {np.mean(lda_o):7.4f}  {np.mean(lda_v):13.6f}")

print("\nchecks:", " ".join(f"{c} {'held' if v else 'FAILED'}" for c, v in checks.items()))
for name, v in h.items():
    print(f"{name}: {'held' if all(v) else 'failed'} overall; "
          f"held on {sum(v)} of {len(v)} seeds, failed on {len(v) - sum(v)}")

exact Bayes error 0.02378; seeds [0, 1, 2, 3, 4]; 2000 training and 5000 test images per seed
seed 0: Bayes rule test error 0.0238; top-4 F pixels [(6, 9), (6, 10), (7, 9), (7, 10)]; top-4 MI pixels [(0, 11), (6, 10), (9, 14), (10, 0)]
   notch F [89.84, 87.22, 94.38, 83.27], best other F 2.46; notch MI (nats) [0.0222, 0.0335, 0.0263, 0.0237], best other MI 0.0337
   best k-dim error 0.0248, LDA error 0.0336, tol 0.0067
seed 1: Bayes rule test error 0.0240; top-4 F pixels [(6, 9), (6, 10), (7, 9), (7, 10)]; top-4 MI pixels [(0, 14), (6, 9), (6, 10), (7, 10)]
   notch F [138.44, 120.07, 122.25, 118.37], best other F 2.61; notch MI (nats) [0.0505, 0.0402, 0.0212, 0.0421], best other MI 0.0380
   best k-dim error 0.0276, LDA error 0.0340, tol 0.0069
seed 2: Bayes rule test error 0.0278; top-4 F pixels [(6, 9), (6, 10), (7, 9), (7, 10)]; top-4 MI pixels [(6, 9), (6, 10), (7, 9), (12, 11)]
   notch F [128.08, 121.35, 122.44, 114.9], best other F 1.81; notch MI (nats) [0.0546, 0.0471, 0.0379, 0.0178], best other MI 0.0358
   best k-dim error 0.0292, LDA error 0.0404, tol 0.0073
seed 3: Bayes rule test error 0.0208; top-4 F pixels [(6, 9), (6, 10), (7, 9), (7, 10)]; top-4 MI pixels [(5, 12), (6, 0), (7, 9), (7, 10)]
   notch F [99.76, 86.63, 93.34, 91.94], best other F 1.65; notch MI (nats) [0.019, 0.0139, 0.0454, 0.0492], best other MI 0.0263
   best k-dim error 0.0232, LDA error 0.0358, tol 0.0068
seed 4: Bayes rule test error 0.0250; top-4 F pixels [(6, 9), (6, 10), (7, 9), (7, 10)]; top-4 MI pixels [(6, 9), (6, 10), (7, 9), (7, 10)]
   notch F [136.49, 136.42, 124.91, 138.71], best other F 3.73; notch MI (nats) [0.0258, 0.0248, 0.0433, 0.0391], best other MI 0.0244
   best k-dim error 0.0276, LDA error 0.0378, tol 0.0071

test error of logistic regression: mean [min, max] over seeds; oracle = exact error of the best rule in the same subspace (mean)
method       k                test error   oracle  variance kept
F-filter     1   0.4062 [0.4028, 0.4086]   0.4065       0.003196
F-filter     2   0.4043 [0.4006, 0.4058]   0.4048       0.006339
F-filter     4   0.4056 [0.4046, 0.4080]   0.4046       0.013012
F-filter     8   0.2238 [0.1742, 0.2676]   0.2244       0.029297
F-filter    16   0.1487 [0.1042, 0.1776]   0.1477       0.063578
MI-filter    1   0.4262 [0.4020, 0.5000]   0.4253       0.003697
MI-filter    2   0.4046 [0.3992, 0.4120]   0.4043       0.007015
MI-filter    4   0.2590 [0.1790, 0.4050]   0.2592       0.014874
MI-filter    8   0.1008 [0.0804, 0.1358]   0.1015       0.029511
MI-filter   16   0.0528 [0.0460, 0.0688]   0.0501       0.059099
PCA          1   0.4966 [0.4866, 0.5020]   0.4984       0.690477
PCA          2   0.4929 [0.4846, 0.5040]   0.4981       0.956318
PCA          4   0.0265 [0.0232, 0.0292]   0.0261       0.957298
PCA          8   0.0272 [0.0236, 0.0306]   0.0260       0.957976
PCA         16   0.0278 [0.0242, 0.0306]   0.0259       0.959347
RandProj     1   0.4838 [0.4524, 0.4938]   0.4850       0.001891
RandProj     2   0.4657 [0.4142, 0.4948]   0.4664       0.005689
RandProj     4   0.4328 [0.4100, 0.4656]   0.4297       0.014072
RandProj     8   0.3822 [0.3514, 0.3984]   0.3820       0.026165
RandProj    16   0.3321 [0.3116, 0.3588]   0.3238       0.052076
LDA          1   0.0363 [0.0336, 0.0404]   0.0340       0.000736

checks: C1 held C2 FAILED C3 held C4 held C5 held C6 held
H1a: held overall; held on 5 of 5 seeds, failed on 0
H1b: failed overall; held on 4 of 5 seeds, failed on 1
H1c: held overall; held on 5 of 5 seeds, failed on 0
H2: failed overall; held on 1 of 5 seeds, failed on 4
H3: failed overall; held on 1 of 5 seeds, failed on 4

Checks. C1, C3, C4 and C5 held. The exact Bayes error is 0.02378, and the exact rule reached 0.0208 to 0.0278 on the five test sets, so the class can be learned from the pixels but not perfectly. C2 failed: the script prints only “C2 FAILED”, not the largest relative difference it compared with $10^{-8}$ , so this article does not claim that f_classif matches the hand-computed ANOVA F. The filter rankings that follow come from f_classif itself. C6 held: no learned error fell below its subspace’s oracle error by more than its bar.

Variance is not the signal. H1a held on all five seeds. PCA kept 0.690 of the test variance with one component and 0.956 with two. On average no other representation kept more than 0.064 at any $k$ , and LDA’s direction kept 0.000736. Yet PCA’s one- and two-component errors (0.4966 and 0.4929) are at chance, with oracle errors of 0.4984 and 0.4981: no classifier could do better from those components. The brightness offset and the gradient vary along two image directions, the constant image and the ramp. In our reading, those two directions are what the first two components follow. The printed oracle error then falls to 0.0261 at $k = 4$ , close to the exact Bayes error, and PCA’s learned errors at $k = 4$ , 8 and 16 (0.0265 to 0.0278) are the lowest of any reduction. Once the nuisance has been absorbed, a few more components contain nearly all the class information. H1b failed: it held on 4 of 5 seeds and failed on 1. The output does not name that seed, but the printed ranges show one failure: the MI filter’s one-pixel error reached 0.5000 on some seed, and PCA’s one-component error is at most 0.5020 on any seed, so there PCA’s margin is at most 0.0020, below the tolerance of about 0.020 that the formula gives at errors near 0.5. With one failing seed, that is the one (in the pilot, 0.5098 against at most 0.5150). H1c held on all five seeds: random projection’s errors at $k = 1$ and 2 (0.4838 and 0.4657 on average) stay below or within the tolerance of PCA’s.

What the filters select. For the F filter, the four notch pixels scored 83.27 to 138.71 against at most 3.73 for any other pixel, and they were its top four on every seed. For the MI filter, the notch pixels scored 0.0139 to 0.0546 nats against 0.0244 to 0.0380 for the best other pixel, and its top four were the notch only on seed 4. H2 therefore failed: it held on 1 of 5 seeds and failed on 4. Even the F filter’s perfect choice gives an error of about 0.405 at $k = 1$ , 2 and 4, matching oracle errors of 0.4046 to 0.4065. By construction, the four notch pixels hold the class signal but also the image’s brightness and gradient, and without other pixels nothing in the representation measures those. The F filter’s errors fall only when it adds pixels outside the notch (0.2238 at $k = 8$ , 0.1487 at $k = 16$ ), and so do its oracle errors. The MI filter’s errors are lower from $k = 4$ on: 0.2590, 0.1008 and 0.0528 against the F filter’s 0.4056, 0.2238 and 0.1487. This is reported as printed; the output does not show why the MI filter’s selections serve better. In general terms, Guyon and Elisseeff’s point is that a variable useless by itself can be useful with others (§3.3); the example does not test whether that is the reason here. Random projection keeps the least variance of the $k$ -dimensional reductions at every $k$ except 4, and its error falls slowly, from 0.4838 at $k = 1$ to 0.3321 at $k = 16$ , following its oracle error.

One supervised direction. LDA’s single direction reached an average error of 0.0363, and its oracle error is 0.0340. The direction estimated from 2,000 images in 256 dimensions is itself worse than the exact discriminant, apart from any loss in the classifier. H3 failed: it held on 1 of 5 seeds (seed 1, 0.0340 against 0.0276 + 0.0069) and failed on 4; on the other seeds LDA’s error exceeded the best $k$ -dimensional error, always a PCA representation’s, by more than the tolerance, for instance 0.0404 against 0.0292 + 0.0073 on seed 2. The best error on every seed (at most 0.0292) belongs to PCA, since every other $k$ -dimensional reduction has a minimum of at least 0.0460. Linear Discriminant Analysis discusses how covariance estimates behave with few samples per dimension and how shrinkage helps.

Note on process. The example departs from the batch plan’s outline in eight ways, all fixed before the pilot:

The generator’s values (amplitudes, noise, the notch on the blob’s flank, the fixed blob), the sizes and the seeds were chosen for the example; the plan gave none.
The oracle error of each learned subspace, and check C6, were added.
Checks C1–C5 were added; C5 makes the plan’s label-free Bayes check concrete.
The plan’s hypothesis “PCA keeps the most variance at every k but has the highest error at small k” was split into H1a–H1c, with “small $k$ ” = 1 and 2 and “highest” meaning within the tolerance. “The filters select the notch pixels first” (H2) means that each filter’s top four pixels are exactly the notch pixels, and H3’s “as good as” means within the tolerance of the best error.
The variance kept is measured on the test images by orthogonal projection onto each representation’s span.
The classifier is standardized logistic regression with $C = 1$ .
The pilot ran seeds 100–104; the first run, shown, ran seeds 0–4.
The XOR variant of Guyon and Elisseeff’s Figure 3 was not run, as the plan said.

Every check and hypothesis here is an illustration, not a test. Every bar and tolerance (C1–C6, the tolerance formula, “small $k$ ” = 1 and 2, and the split of H1) was written after a planning script, reach.py, had printed the exact quantities the hypotheses rest on:

the exact Bayes error, 0.02378;
the oracle errors of one notch pixel, (6, 9), and of the four notch pixels, 0.4059 and 0.4046;
the oracle errors of the population principal subspaces, 0.498 at $k = 1$ and 2 and 0.02378 for $k \ge 4$ ;
the oracle errors of the pilot’s random projections, 0.475 to 0.498 at $k = 1$ and 0.460 to 0.488 at $k = 2$ ;
the population variance kept by the first one and two principal components, 0.694 and 0.956, and the largest kept by any 16 pixels, 0.0894;
the approximate expected F score of a notch pixel, 108 to 114, against 21.40 for the 99.9% quantile of the largest of 252 null scores treated as independent;
the exact mutual information between notch pixel (6, 9) and the class, 0.02757 nats.

H1c’s tolerance form was chosen knowing that random projection’s oracle errors at small $k$ are close to PCA’s. The pilot ran seeds 100–104, and its outcomes differ from those shown in four places:

C6 failed in the pilot and held in the first run.
H2 held on none of the pilot’s seeds, against 1 of 5 here; it failed overall in both.
At $k = 2$ the MI filter’s mean error was below the F filter’s in the pilot (0.3242 against 0.4022) but not here (0.4046 against 0.4043). At $k \ge 4$ it was below in both.
Random projection kept the least mean variance of the $k$ -dimensional reductions at every $k$ except 16 in the pilot (the MI filter kept less there), and at every $k$ except 4 here.

C2 failed in both runs. The other checks held in both, and H1a, H1b, H1c and H3 held on the same numbers of seeds in both. The pilot’s hypothesis lines printed the number of seeds on which a hypothesis held under a misleading label (“failed on 4 of 5 seeds” meant held on 4). With the editor’s approval the label was changed, before the first run, to the “held on n of 5 seeds, failed on m” form shown. Nothing else changed between the pilot and the first run except the seeds.

Common Misconceptions

“The components that keep the most variance keep the most information.” They keep the most information about the pixels in the squared-error sense, which is what PCA optimizes. For a task, the relevant variation can lie in directions with tiny variance, as in the example, where two components keep 96% of the variance and none of the class signal.
“Ranking variables one at a time finds the best subset.” Guyon and Elisseeff’s examples show variables that help only together, and in the example the notch pixels alone cannot discount the brightness that every pixel shares.
“Selection and extraction are interchangeable.” Selection keeps measurable, interpretable variables and can lower measurement and storage costs, which Guyon and Elisseeff (§5) give as the reason to prefer it when those costs matter. Extraction can combine information spread over many variables, but every new variable depends on all the inputs.
“A 2-D embedding is a reduced representation for any image.” t-SNE, spectral embedding and scikit-learn’s MDS give coordinates only for the training points; Isomap’s out-of-sample map is an approximation through the training graph. A pipeline that must reduce new images needs a parametric method or such an extension.
“Selecting features before cross-validation is harmless.” The selection has seen the held-out labels; it must be refit within each training fold, as Cross-Validation explains.
“Classical MDS is a different algorithm from PCA.” On Euclidean distances the two give the same embedding (Murphy, 2022, §20.4.4.1); MDS becomes different when it is given other dissimilarities or a stress objective.

Where It Is Used

Descriptors and retrieval. PCA compresses local and global descriptors before matching, as in PCA-SIFT and, with whitening, VLAD, and product quantization compresses vectors further for approximate nearest neighbor search; random projections appear in its hashing schemes.
Faces and appearance models. Eigenfaces represent a face by PCA coefficients, and Fisherfaces replace them by LDA directions, as Face Recognition describes.
Feature selection in detectors. The Viola–Jones detector selects a few Haar features out of a very large pool by boosting; Random Forest measures pixel importance, which can be used to rank variables.
Visualization of representations. t-SNE and related embeddings are used to inspect what a network’s features separate; Embeddings and Representation Learning discuss what a good representation needs beyond being short.
Learning on fewer dimensions. A reduction step can precede a classical learner such as k-nearest neighbors or a support vector machine, as Murphy (2022, ch. 20) describes for a linear classifier; Curse of Dimensionality explains why the dimension that matters to such learners is the intrinsic one.