Curse of Dimensionality

The curse of dimensionality is the collection of ways in which problems become harder, often exponentially harder, as the number of variables that describe each example grows. Goodfellow, Bengio and Courville (2016, §5.11.1) single out the growth in the number of distinct configurations of the variables, and note that the curse arises in many places in computer science, especially in machine learning. For images the question is pressing: a $32 \times 32$ grey image is a point in a space of 1,024 dimensions, and a megapixel photograph a point in a space of millions. Methods that learn by looking at nearby examples, such as k-nearest neighbours, would need impossibly many images if those dimensions were all free to vary.

They are not. Images of real scenes occupy a tiny part of pixel space, and the number that matters for how much data a method needs is the intrinsic dimension of the data, the number of independent ways the images actually vary, rather than the ambient dimension, the length of the pixel vector. This article gathers the forms of the curse (combinatorial, geometric, statistical and computational), links the articles that treat particular forms in depth, describes how intrinsic dimension is estimated, and measures both ideas on synthetic images whose intrinsic dimension is known.

Definition

The curse of dimensionality is not one theorem but a family of related effects of high dimension $d$ :

Combinatorial: a grid with $v$ distinguishable values per variable has $v^d$ cells, so covering the space with examples needs a number of examples exponential in $d$ .
Geometric: almost all the volume of a high-dimensional ball lies near its surface, and a ball occupies a vanishing fraction of the cube around it.
Statistical: for a fixed number of examples, the neighbourhood that holds a given number of them grows with $d$ , so estimates that average over neighbours stop being local.
Distance concentration: under broad conditions the distances from a query to all data points become nearly equal.
Computational: exact search and sampling methods whose cost depends on partitioning or covering the space degrade with $d$ .

Which $d$ governs these effects depends on the data. For points that fill a $d$ -dimensional region uniformly it is the coordinate count; for data that lie on or near a lower-dimensional set inside a high-dimensional space, the statistical effects follow the intrinsic dimension of that set.

Intuition

Goodfellow and colleagues (§5.11.1, Figure 5.9) give the counting picture. To learn a function of one variable well enough to distinguish 10 regions, a learner needs examples in each of the 10 cells. With two variables and 10 values each there are 100 cells, with three 1,000, and with $d$ variables $O(v^d)$ . A learner that predicts by copying the answer from training examples in the same cell, or from the nearest one, is then left with mostly empty cells: the number of configurations is far larger than any training set.

The geometric forms say the same thing about volume. Take a unit ball and shrink it by 5% of its radius: in three dimensions the shrunken ball keeps 86% of the volume, in ten dimensions 60%, in a hundred dimensions 0.6%. Almost every point of a high-dimensional ball is near its boundary, so “the typical point” and “the centre” are far apart.

What rescues learning from images is that images are not spread through pixel space. Goodfellow and colleagues (§5.11.3, Figure 5.12) note that pixels drawn uniformly at random give static that essentially never resembles a face or any other structured scene, and that one can imagine the transformations that trace out a low-dimensional manifold of natural images: dimming or brightening the light, moving or rotating objects, changing their colours. The machine learning article treats this manifold hypothesis as one of the inductive biases vision relies on. If a set of images varies in only a few ways, a nearest-neighbour method behaves as if the data had those few dimensions, however many pixels each image has.

Formal Definition

Counting and sample size

For a function that must be resolved to $v$ values along each of $d$ input variables, a method that treats every cell independently needs $O(v^d)$ examples (Goodfellow et al., §5.11.1). The statistical version follows from local averaging. For $n$ points with density $f$ in $\mathbb{R}^m$ , the distance $T_k(x)$ from a point $x$ to its $k$ -th nearest neighbour satisfies approximately

\frac{k}{n} \approx f(x)\, V(m)\, T_k(x)^m, \qquad V(m) = \frac{\pi^{m/2}}{\Gamma(m/2 + 1)},

the fraction of the sample in a small ball being the density times the ball’s volume (Levina & Bickel, 2004, eq. 2). Solving for $T_k$ gives $T_k(x) \propto (k/n)^{1/m}$ : a fixed neighbour count is found within a radius that shrinks only as $n^{-1/m}$ . In one dimension, doubling the data halves the distance to the nearest neighbour; in $m$ dimensions it shrinks the distance by the factor $2^{-1/m}$ , and halving it needs $2^m$ times as much data. By our reasoning, an estimator that averages noise-free targets over the $k$ nearest neighbours has an error set by how far those neighbours are, to first order the gradient of the target times their mean offset, so its root-mean-square error also falls as $n^{-1/m}$ and its mean squared error as $n^{-2/m}$ .

The k-nearest neighbours article gives, from Hastie, Tibshirani and Friedman (§2.5), whom it reports tracing the phrase to Bellman (1961), the cube form that Murphy also presents (§16.1.2): a cube-shaped neighbourhood that captures a fraction $r$ of uniformly spread data has edge length $r^{1/p}$ in $p$ dimensions, so in ten dimensions it spans 80% of every coordinate’s range to hold 10% of the data. Murphy (§16.1.2) names two ways out: assume a form for the function (a parametric model), or use a metric that attends only to some of the dimensions.

Volume near the surface and the ball inside the cube

For any body $A$ in $\mathbb{R}^d$ , shrinking it by a factor $1 - \varepsilon$ scales its volume by $(1 - \varepsilon)^d \le e^{-\varepsilon d}$ , so at least a fraction $1 - e^{-\varepsilon d}$ of the volume lies in the shell of relative thickness $\varepsilon$ at the boundary (Blum, Hopcroft & Kannan, 2020, §2.3); for a ball, most of the volume is within $O(1/d)$ of the surface. The volume of the unit ball,

V(d) = \frac{\pi^{d/2}}{\Gamma(d/2 + 1)} = \frac{2\pi^{d/2}}{d\,\Gamma(d/2)},

tends to zero as $d$ grows, because $\Gamma(d/2)$ grows like a factorial while $\pi^{d/2}$ grows only exponentially (Blum et al., Lemma 2.6, which writes the second form). The cube $[-1, 1]^d$ around the ball has volume $2^d$ , so the ratio $V(d)/2^d$ collapses even faster: points drawn uniformly in the cube almost never fall in the ball it encloses.

Two further statements in Blum and colleagues’ Chapter 2 complete the picture. For points drawn uniformly from the unit ball, with probability $1 - O(1/n)$ all of them have length at least $1 - 2\ln n / d$ and every pair is nearly orthogonal, with $|x_i \cdot x_j| \le \sqrt{6 \ln n}/\sqrt{d - 1}$ (Theorem 2.8); and for Gaussian coordinates the distances between all pairs of points become essentially equal for large $d$ (§2.1–2.2). The latent space article measures the Gaussian version, the thin shell on which high-dimensional Gaussian codes live.

Intrinsic dimension and its maximum-likelihood estimate

Levina and Bickel (2004, §3) model the data as $X_i = g(Y_i)$ , where the $Y_i$ are drawn from a smooth density on $\mathbb{R}^m$ and $g$ is a smooth map into $\mathbb{R}^p$ with $m \le p$ ; $m$ is the intrinsic dimension. Treating the points near $x$ as a Poisson process whose rate grows as $t^{m-1}$ with the distance $t$ , they derive the maximum-likelihood estimate from the distances $T_1(x) \le \dots \le T_k(x)$ to the $k$ nearest neighbours (eq. 8):

\hat m_k(x) = \left[ \frac{1}{k - 1} \sum_{j=1}^{k-1} \log \frac{T_k(x)}{T_j(x)} \right]^{-1},

averaged over all $n$ points and then over a range $k = k_1, \dots, k_2$ (eq. 9); they fix $k_1 = 10$ and $k_2 = 20$ for all their results (§4). It follows from eq. (8) that the estimator uses only ratios of neighbour distances, so multiplying the data by a constant does not change it; its cost is that of finding $k_2$ neighbours for every point. They note that dividing by $k - 2$ instead of $k - 1$ makes it asymptotically unbiased: under the Poisson approximation the $k - 1$ values $m \log(T_k/T_j)$ behave as standard exponential order statistics, so their sum $U$ follows a Gamma $(k - 1, 1)$ distribution and $\mathbb{E}[U^{-1}] = 1/(k-2)$ (§3.1). The paper prints the factor as $m^{-1}$ ; by our derivation it must be $m$ , since $(T_j/T_k)^m$ are then uniform order statistics, and only with $m$ does $U$ have that distribution. It follows, by our arithmetic, that the $k - 1$ form has expectation $m(k-1)/(k-2)$ under that approximation, a factor of 1.125 at $k = 10$ and 1.056 at $k = 20$ , about 1.08 averaged over $k = 10, \dots, 20$ .

In simulations on $m$ -dimensional Gaussians the estimate showed a negative bias for large $k$ , which decreased with sample size and increased with dimension (Figure 1), and data near rather than on a manifold changed the estimate only when the spread off the manifold became very small (Figure 2a). They attribute the negative bias at high dimension to the large sample a small ball needs before it holds enough points, and to edge effects for data such as points in a cube (§5). On images they reported (Table 1) an estimate of 4.3 for 698 rendered face images of $64 \times 64$ pixels varying in illumination and two pose angles, and 3.1 for 481 frames of $480 \times 512$ pixels of a hand rotating along a one-dimensional path; they argue that several basis images are needed to represent the different poses, so an estimate above the number of physical parameters is reasonable.

Properties

The dimension that counts is intrinsic. The statistical forms of the curse follow the dimension of the set the data occupy, not the length of the vector. Muja and Lowe’s finding, summarized in the approximate nearest neighbor search article, that their approximate search lost its speedup over a linear scan beyond about 800 dimensions on uniform vectors but not on image patches of up to 4,096 dimensions, is one instance; the example below is another.
Distance concentration. Beyer and colleagues’ theorem, given in the k-nearest neighbours article, states when the farthest and nearest points of a data set become indistinguishable in relative terms; that article also measures how pixel noise erodes the contrast between them. The similarity and distance measures article discusses how the choice of $\ell_p$ norm interacts with concentration.
Hubness. In high dimensions some points appear among the $k$ nearest neighbours of many others, as the embeddings article describes. Radovanović, Nanopoulos and Ivanović (2010, §4.3) measured the skewness of these neighbour counts on 50 real data sets and found it correlated with the ambient dimension and more strongly with the intrinsic dimension given by the Levina–Bickel estimator (their Table 2).
Search and partitioning. Exact tree indexes visit a growing share of the data as the dimension grows; the approximate nearest neighbor search article gives Weber, Schek and Blott’s analysis, the k-d tree measurements and locality-sensitive hashing, two of whose cited papers carry the title “Towards Removing the Curse of Dimensionality”.
Sequential Monte Carlo. The particle filter article describes how the number of particles needed to avoid weight collapse grows exponentially with the variance of the observation log-likelihood, which typically grows with the state and measurement dimensions.
Kernels do not remove it. The support vector machine article reports Bishop’s warning that a kernel’s high-dimensional feature space is not a way around the curse, since the mapped features are constrained to a low-dimensional surface.
Random projection preserves distances in few dimensions. The Johnson–Lindenstrauss lemma, in the form Blum and colleagues prove (Theorem 2.11, with the proof they credit to Dasgupta and Gupta, §2.10), says that projecting $n$ points onto $k \ge 3 \ln n / (c \varepsilon^2)$ random Gaussian directions, for a constant $c$ they do not specify, keeps every pairwise distance within a factor $1 \pm \varepsilon$ of $\sqrt{k}$ times the original with probability at least $1 - 3/(2n)$ . The required $k$ depends on the number of points, not on the original dimension. scikit-learn 1.7.2’s johnson_lindenstrauss_min_dim uses a different statement of the bound, $k \ge 4 \ln n / (\varepsilon^2/2 - \varepsilon^3/3)$ for squared distances within $1 \pm \varepsilon$ ; its docstring gives 11,841 components for a million points at $\varepsilon = 0.1$ . The two bounds use $\varepsilon$ for different quantities and are not interchangeable.

Examples

Intrinsic dimension at a fixed number of pixels

The script below holds the ambient dimension fixed and varies the intrinsic dimension. Every image is $32 \times 32$ pixels, a point in $\mathbb{R}^{1024}$ . The image area is divided into a grid of 2 columns by 4 rows of $16 \times 8$ -pixel cells, and a set of intrinsic dimension $d \in \{1, 2, 4, 8\}$ draws Gaussian spots (standard deviation 1.5 pixels, peak 1) in the first $d$ cells; each spot’s parameter $u_i$ , uniform on $[0, 1]$ , slides it 6 pixels sideways inside its cell, and the other cells stay blank. Because the parameters move the spots rather than scale them, the images lie on a curved $d$ -dimensional manifold, not a flat subspace. Pixel coordinates are pixel centres $0, \dots, 31$ . The target to predict is $y = d^{-1/2} \sum_i \cos(\pi u_i)$ , whose variance is $1/2$ for every $d$ , so errors are reported as root-mean-square error divided by $\sqrt{1/2}$ (relative RMSE), where 1 is the error of always predicting the mean.

The script first checks that the spots of different cells barely overlap, that no spot is cut by the image border and that blank cells hold only the faint tails of neighbouring spots. It then runs three blocks:

(a) For seeds 0, 1 and 2, a training set of 4,000 images and a test set of 500, it measures for nested training sizes $n$ from 125 to 4,000 the mean distance from each test image to its nearest training image, and the error of k-nearest-neighbour regression with $k = 5$ . Both are computed on the pixels, which is the method under test, and on the parameter vectors $u$ , the exact geometry of the sample. It fits the slope of each quantity against $n$ on log–log axes and compares it with $-1/d$ , after checking its NumPy k-NN regression against scikit-learn’s KNeighborsRegressor.
(b) For $n = 2{,}000$ images per set and the same seeds, it computes the Levina–Bickel estimate with $k = 10, \dots, 20$ on the parameter vectors, on the noise-free pixels and on pixels with independent Gaussian noise of standard deviation 0.02, in the paper’s $k - 1$ form and, in brackets, the $k - 2$ form. Self-matches are excluded, and the NumPy neighbour distances are checked against scikit-learn’s NearestNeighbors(...).kneighbors(), which, called without an argument, does not count a point as its own neighbour (docstring in sklearn/neighbors/_base.py).
(c) It tabulates the fraction of the unit ball’s volume inside 95% of its radius, the outer 5% shell, and $\log_{10}$ of the ball-to-cube ratio $V(d)/2^d$ , computed in log space with gammaln so that $d = 1{,}024$ does not overflow, and checks the ball and shell fractions by Monte Carlo with $10^6$ points for $d \le 10$ , with a tolerance of four binomial standard errors.

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 10 seconds. Its output was identical over two runs.

import os
for v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"):
    os.environ[v] = "1"                       # one BLAS thread, set before NumPy is imported
import numpy as np
from scipy.special import gammaln
from sklearn.neighbors import NearestNeighbors, KNeighborsRegressor

# ---- Images with d known generating parameters, all in the same 32x32 = 1,024-pixel space ----
S, SIG = 32, 1.5
YY, XX = np.mgrid[0:S, 0:S].astype(float)

def spot_centres(u):
    """u: (n, d) in [0, 1]. Slot i sits in cell (column i % 2, row i // 2) of a 2 x 4 grid of
    16 x 8-pixel cells; its parameter slides the spot 6 px horizontally inside the cell."""
    i = np.arange(u.shape[1])
    return 16 * (i % 2) + 4.5 + 6 * u, np.broadcast_to(8 * (i // 2) + 3.5, u.shape)

def render(u):
    cx, cy = spot_centres(u)
    img = np.zeros((len(u), S, S))
    for i in range(u.shape[1]):
        img += np.exp(-((XX - cx[:, i, None, None]) ** 2 + (YY - cy[:, i, None, None]) ** 2) / (2 * SIG ** 2))
    return img.reshape(len(u), -1)

def target(u):                                  # smooth, Var = 1/2 for every d
    return np.cos(np.pi * u).sum(1) / np.sqrt(u.shape[1])

def sqd(A, B):                                  # squared Euclidean distances, Gram trick, float64
    return np.maximum((A * A).sum(1)[:, None] + (B * B).sum(1)[None, :] - 2 * A @ B.T, 0)

DIMS, SEEDS = (1, 2, 4, 8), (0, 1, 2)

# ---- Data check ----
g = np.linspace(0, 1, 13)
def one_spot(i, ui):                            # spot i alone at positions ui
    u = np.zeros((len(ui), i + 1)); u[:, i] = ui
    cx, cy = spot_centres(u)
    return np.exp(-((XX - cx[:, i, None, None]) ** 2 + (YY - cy[:, i, None, None]) ** 2) / (2 * SIG ** 2)).reshape(len(ui), -1)
S8 = [one_spot(i, g) for i in range(8)]
overlap = max(np.max((S8[i] / np.linalg.norm(S8[i], axis=1, keepdims=True)) @
                     (S8[j] / np.linalg.norm(S8[j], axis=1, keepdims=True)).T)
              for i in range(8) for j in range(8) if i != j)
mass = min(Si.sum(1).min() for Si in S8) / (2 * np.pi * SIG ** 2)
unused = 0.0
for d in DIMS[:-1]:
    img = render(np.random.default_rng(99).uniform(size=(200, d))).reshape(-1, S, S)
    for i in range(d, 8):
        r, c = i // 2, i % 2
        unused = max(unused, img[:, 8 * r:8 * r + 8, 16 * c:16 * c + 16].max())
print(f"data: images {S}x{S} = {S*S} px for every d; max spot overlap between slots {overlap:.4f} (<= 0.01); "
      f"min spot mass / 2 pi sigma^2 {mass:.4f} (>= 0.99); max pixel in unused cells {unused:.4f} (<= 0.02)")

# ---- (a) Nearest-neighbour distance and k-NN regression error against n ----
NS, NTEST, K = (125, 250, 500, 1000, 2000, 4000), 500, 5

def knn_predict(Dsq, ytr, k):
    idx = np.argpartition(Dsq, k, axis=1)[:, :k]
    return ytr[idx].mean(1)

# hand k-NN regression against scikit-learn (seed 0, d = 4, n = 1000)
rng = np.random.default_rng([0, 4, 0])
utr, ute = rng.uniform(size=(NS[-1], 4)), rng.uniform(size=(NTEST, 4))
Xtr, Xte = render(utr[:1000]), render(ute)
mine = knn_predict(sqd(Xte, Xtr), target(utr[:1000]), K)
sk = KNeighborsRegressor(n_neighbors=K).fit(Xtr, target(utr[:1000])).predict(Xte)
print(f"check: NumPy k-NN regression vs KNeighborsRegressor, max |difference| = {np.abs(mine - sk).max():.1e}")

print("\n(a) mean distance to the nearest training image, and k-NN (k=5) relative RMSE, against n; seed mean")
slopes, rmse_at_max = {}, {}
for d in DIMS:
    rows = {key: [] for key in ("nn_pix", "nn_par", "rm_pix", "rm_par")}
    for seed in SEEDS:
        rng = np.random.default_rng([seed, d, 0])
        utr, ute = rng.uniform(size=(NS[-1], d)), rng.uniform(size=(NTEST, d))
        Xtr, Xte, ytr, yte = render(utr), render(ute), target(utr), target(ute)
        Dp, Du = sqd(Xte, Xtr), sqd(ute, utr)
        r = {key: [] for key in rows}
        for n in NS:
            for sp, D in (("pix", Dp[:, :n]), ("par", Du[:, :n])):
                r["nn_" + sp].append(np.sqrt(D.min(1)).mean())
                r["rm_" + sp].append(np.sqrt(np.mean((knn_predict(D, ytr[:n], K) - yte) ** 2)) / np.sqrt(0.5))
        for key in rows:
            rows[key].append(r[key])
    rows = {key: np.array(v) for key, v in rows.items()}          # (seed, n)
    L = np.log(NS)
    slopes[d] = {key: np.array([np.polyfit(L, np.log(v[s]), 1)[0] for s in range(len(SEEDS))]) for key, v in rows.items()}
    rmse_at_max[d] = rows["rm_pix"][:, -1]
    print(f"d = {d}      n: " + "".join(f"{n:>8}" for n in NS))
    for key, name in (("nn_pix", "NN dist, pixels"), ("nn_par", "NN dist, params"),
                      ("rm_pix", "rel RMSE, pixels"), ("rm_par", "rel RMSE, params")):
        print(f"  {name:<17}" + "".join(f"{x:>8.4f}" for x in rows[key].mean(0)))

print("\nfitted log-log slopes against n: seed mean (min to max over seeds 0, 1, 2)")
print("  d    -1/d   NN dist, pixels          NN dist, params          rel RMSE, pixels         rel RMSE, params")
for d in DIMS:
    s = slopes[d]
    print(f"  {d}  {-1/d:7.4f}" + "".join(f"  {s[k].mean():7.4f} ({s[k].min():.4f} to {s[k].max():.4f})"
                                       for k in ("nn_pix", "nn_par", "rm_pix", "rm_par")))
print("\nH3 tolerances |mean param slope + 1/d| + 2 * seed range + 0.02, and pixel-slope deviation |s + 1/d|")
for d in DIMS:
    s = slopes[d]
    for a, b, h in (("nn_pix", "nn_par", "H3a"), ("rm_pix", "rm_par", "H3b")):
        tol = abs(s[b].mean() + 1 / d) + 2 * np.ptp(s[b]) + 0.02
        dev = abs(s[a].mean() + 1 / d)
        print(f"  d = {d} {h}: deviation {dev:.4f}, tolerance {tol:.4f} -> {'held' if dev <= tol else 'failed'}")
print("H2: relative RMSE (pixels) at n = 4000 by seed: " +
      "; ".join(f"d={d}: " + " ".join(f"{x:.4f}" for x in rmse_at_max[d]) for d in DIMS))
inc = all(np.all(rmse_at_max[a] < rmse_at_max[b]) for a, b in zip(DIMS, DIMS[1:]))
print(f"H2 (rises with d for every seed): {'held' if inc else 'failed'}")

# ---- (b) Levina-Bickel maximum-likelihood estimate of intrinsic dimension ----
NLB, K1, K2, NOISE = 2000, 10, 20, 0.02

def neighbours(X):                              # k2 nearest distances and indices, self excluded
    D = sqd(X, X); np.fill_diagonal(D, np.inf)
    idx = np.argpartition(D, K2, axis=1)[:, :K2]
    dd = np.take_along_axis(D, idx, 1); o = np.argsort(dd, axis=1)
    return np.sqrt(np.take_along_axis(dd, o, 1)), np.take_along_axis(idx, o, 1)

def levina_bickel(T, norm):
    """Levina & Bickel (2004) eq. (8)-(9); norm = 1 gives (k-1)/sum, norm = 2 gives (k-2)/sum."""
    est = [np.mean((k - norm) / np.log(T[:, k - 1:k] / T[:, :k - 1]).sum(1)) for k in range(K1, K2 + 1)]
    return np.mean(est)

print("\n(b) Levina-Bickel estimate, k = 10..20, n = 2000; seed mean (min to max), k-1 form [k-2 form]")
est = {}
mind = np.inf
for d in DIMS:
    for seed in SEEDS:
        rng = np.random.default_rng([seed, d, 1])
        u = rng.uniform(size=(NLB, d))
        X = render(u)
        Xn = X + rng.normal(0, NOISE, X.shape)
        for name, Z in (("params", u), ("pixels", X), ("noisy", Xn)):
            T, idx = neighbours(Z)
            mind = min(mind, T[:, 0].min())
            est[d, seed, name] = (levina_bickel(T, 1), levina_bickel(T, 2))
            if seed == 0 and name == "pixels":
                Ts, Is = NearestNeighbors(n_neighbors=K2).fit(Z).kneighbors()
                print(f"  check d = {d}: max |NumPy - sklearn distance| {np.abs(T - Ts).max():.1e}, "
                      f"index mismatches {int((idx != Is).sum())}, "
                      f"|estimate difference| {abs(levina_bickel(Ts, 1) - est[d, seed, name][0]):.1e}")
print(f"  smallest nearest-neighbour distance in any set: {mind:.3e} (> 1e-6)")
print("  d   parameters                 pixels, noise-free         pixels, noise 0.02")
for d in DIMS:
    line = f"  {d}"
    for name in ("params", "pixels", "noisy"):
        a = np.array([est[d, s, name][0] for s in SEEDS]); b = np.array([est[d, s, name][1] for s in SEEDS])
        line += f"   {a.mean():6.3f} ({a.min():.3f}-{a.max():.3f}) [{b.mean():.3f}]"
    print(line)
for d in DIMS[:3]:
    par = np.mean([est[d, s, "params"][0] for s in SEEDS]); tol = abs(par - d) + 0.1 * d
    dev = max(abs(est[d, s, "pixels"][0] - d) for s in SEEDS)
    print(f"H1 d = {d}: max over seeds |estimate - d| {dev:.3f}, tolerance {tol:.3f} -> {'held' if dev <= tol else 'failed'}")
bias = [np.mean([est[d, s, "pixels"][0] for s in SEEDS]) - d for d in DIMS]
dec = all(a > b for a, b in zip(bias, bias[1:]))
print("H1b signed bias (pixels, seed mean) " + " ".join(f"{b:+.3f}" for b in bias) +
      f" -> strictly decreasing: {'held' if dec else 'failed'}; negative at d = 8: {'held' if bias[-1] < 0 else 'failed'}")
gap = min(est[d, s, "noisy"][0] - est[d, s, "pixels"][0] for d in DIMS for s in SEEDS)
print(f"H4 noise raises the estimate for every d and seed: smallest rise {gap:.3f} -> {'held' if gap > 0 else 'failed'}")

# ---- (c) Volume in high dimensions: closed form, with a Monte Carlo check ----
print("\n(c)    d   inner 95% of radius   outer 5% shell   log10(ball / cube)")
for d in (1, 2, 3, 5, 10, 20, 50, 100, 1024):
    lr = (d / 2 * np.log(np.pi) - gammaln(d / 2 + 1) - d * np.log(2)) / np.log(10)
    print(f"  {d:>6}   {0.95 ** d:>18.3e}   {1 - 0.95 ** d:>14.6f}   {lr:>18.3f}")
rng = np.random.default_rng(0)
N = 1_000_000
print("  Monte Carlo, 10^6 points uniform in [-1, 1]^d: fraction in ball, and outer-shell fraction of those")
for d in (1, 2, 3, 5, 10):
    r = np.linalg.norm(rng.uniform(-1, 1, (N, d)), axis=1)
    p = np.exp(d / 2 * np.log(np.pi) - gammaln(d / 2 + 1) - d * np.log(2)); q = 1 - 0.95 ** d
    inb = r <= 1; ph, qh = inb.mean(), (r[inb] > 0.95).mean()
    tp, tq = 4 * np.sqrt(p * (1 - p) / N), 4 * np.sqrt(q * (1 - q) / inb.sum())
    print(f"  d = {d:>2}: ball {ph:.5f} vs {p:.5f} (|diff| {abs(ph-p):.5f}, tol {tp:.5f}) "
          f"shell {qh:.5f} vs {q:.5f} (|diff| {abs(qh-q):.5f}, tol {tq:.5f}) "
          f"-> {'ok' if abs(ph-p) <= tp and abs(qh-q) <= tq else 'FAIL'}")

Output:

data: images 32x32 = 1024 px for every d; max spot overlap between slots 0.0008 (<= 0.01); min spot mass / 2 pi sigma^2 0.9964 (>= 0.99); max pixel in unused cells 0.0111 (<= 0.02)
check: NumPy k-NN regression vs KNeighborsRegressor, max |difference| = 4.4e-16

(a) mean distance to the nearest training image, and k-NN (k=5) relative RMSE, against n; seed mean
d = 1      n:      125     250     500    1000    2000    4000
  NN dist, pixels    0.0298  0.0166  0.0076  0.0038  0.0019  0.0009
  NN dist, params    0.0040  0.0022  0.0010  0.0005  0.0003  0.0001
  rel RMSE, pixels   0.0210  0.0130  0.0051  0.0026  0.0013  0.0006
  rel RMSE, params   0.0210  0.0130  0.0051  0.0026  0.0013  0.0006
d = 2      n:      125     250     500    1000    2000    4000
  NN dist, pixels    0.3475  0.2443  0.1735  0.1220  0.0848  0.0600
  NN dist, params    0.0464  0.0325  0.0231  0.0162  0.0113  0.0080
  rel RMSE, pixels   0.0954  0.0605  0.0463  0.0306  0.0220  0.0160
  rel RMSE, params   0.0957  0.0604  0.0462  0.0306  0.0220  0.0160
d = 4      n:      125     250     500    1000    2000    4000
  NN dist, pixels    1.4548  1.2082  1.0211  0.8517  0.7128  0.5973
  NN dist, params    0.1987  0.1635  0.1375  0.1142  0.0953  0.0798
  rel RMSE, pixels   0.2721  0.1965  0.1568  0.1255  0.0986  0.0803
  rel RMSE, params   0.2703  0.1984  0.1573  0.1257  0.0987  0.0799
d = 8      n:      125     250     500    1000    2000    4000
  NN dist, pixels    3.5133  3.2138  2.9484  2.7140  2.4825  2.2620
  NN dist, params    0.5127  0.4601  0.4160  0.3791  0.3437  0.3111
  rel RMSE, pixels   0.4734  0.4041  0.3385  0.2954  0.2597  0.2303
  rel RMSE, params   0.4642  0.3977  0.3328  0.2940  0.2597  0.2297

fitted log-log slopes against n: seed mean (min to max over seeds 0, 1, 2)
  d    -1/d   NN dist, pixels          NN dist, params          rel RMSE, pixels         rel RMSE, params
  1  -1.0000  -1.0090 (-1.0403 to -0.9816)  -1.0091 (-1.0403 to -0.9816)  -1.0233 (-1.0792 to -0.9791)  -1.0233 (-1.0792 to -0.9791)
  2  -0.5000  -0.5074 (-0.5251 to -0.4884)  -0.5082 (-0.5259 to -0.4892)  -0.5092 (-0.5589 to -0.4814)  -0.5096 (-0.5591 to -0.4828)
  4  -0.2500  -0.2562 (-0.2582 to -0.2549)  -0.2625 (-0.2645 to -0.2610)  -0.3453 (-0.3701 to -0.3277)  -0.3464 (-0.3701 to -0.3305)
  8  -0.1250  -0.1261 (-0.1324 to -0.1225)  -0.1428 (-0.1507 to -0.1379)  -0.2085 (-0.2153 to -0.2015)  -0.2024 (-0.2168 to -0.1926)

H3 tolerances |mean param slope + 1/d| + 2 * seed range + 0.02, and pixel-slope deviation |s + 1/d|
  d = 1 H3a: deviation 0.0090, tolerance 0.1464 -> held
  d = 1 H3b: deviation 0.0233, tolerance 0.2434 -> held
  d = 2 H3a: deviation 0.0074, tolerance 0.1017 -> held
  d = 2 H3b: deviation 0.0092, tolerance 0.1822 -> held
  d = 4 H3a: deviation 0.0062, tolerance 0.0396 -> held
  d = 4 H3b: deviation 0.0953, tolerance 0.1955 -> held
  d = 8 H3a: deviation 0.0011, tolerance 0.0633 -> held
  d = 8 H3b: deviation 0.0835, tolerance 0.1457 -> held
H2: relative RMSE (pixels) at n = 4000 by seed: d=1: 0.0006 0.0007 0.0007; d=2: 0.0151 0.0171 0.0158; d=4: 0.0803 0.0816 0.0791; d=8: 0.2321 0.2317 0.2272
H2 (rises with d for every seed): held

(b) Levina-Bickel estimate, k = 10..20, n = 2000; seed mean (min to max), k-1 form [k-2 form]
  check d = 1: max |NumPy - sklearn distance| 2.1e-10, index mismatches 0, |estimate difference| 1.2e-09
  check d = 2: max |NumPy - sklearn distance| 6.8e-13, index mismatches 0, |estimate difference| 9.3e-15
  check d = 4: max |NumPy - sklearn distance| 2.9e-14, index mismatches 0, |estimate difference| 0.0e+00
  check d = 8: max |NumPy - sklearn distance| 1.0e-14, index mismatches 0, |estimate difference| 0.0e+00
  smallest nearest-neighbour distance in any set: 7.355e-08 (> 1e-6)
  d   parameters                 pixels, noise-free         pixels, noise 0.02
  1    1.077 (1.073-1.082) [0.995]    1.077 (1.073-1.082) [0.995]   103.780 (103.496-104.176) [95.829]
  2    2.128 (2.120-2.136) [1.967]    2.132 (2.125-2.140) [1.971]   27.126 (27.061-27.179) [24.984]
  4    4.001 (3.971-4.016) [3.697]    4.134 (4.105-4.152) [3.821]    6.423 (6.376-6.453) [5.932]
  8    7.309 (7.300-7.326) [6.754]    8.620 (8.616-8.626) [7.967]    9.286 (9.276-9.292) [8.582]
H1 d = 1: max over seeds |estimate - d| 0.082, tolerance 0.177 -> held
H1 d = 2: max over seeds |estimate - d| 0.140, tolerance 0.328 -> held
H1 d = 4: max over seeds |estimate - d| 0.152, tolerance 0.401 -> held
H1b signed bias (pixels, seed mean) +0.077 +0.132 +0.134 +0.620 -> strictly decreasing: failed; negative at d = 8: failed
H4 noise raises the estimate for every d and seed: smallest rise 0.660 -> held

(c)    d   inner 95% of radius   outer 5% shell   log10(ball / cube)
       1            9.500e-01         0.050000                0.000
       2            9.025e-01         0.097500               -0.105
       3            8.574e-01         0.142625               -0.281
       5            7.738e-01         0.226219               -0.784
      10            5.987e-01         0.401263               -2.604
      20            3.585e-01         0.641514               -7.609
      50            7.694e-02         0.923055              -27.813
     100            5.921e-03         0.994079              -69.729
    1024            1.545e-23         1.000000            -1220.255
  Monte Carlo, 10^6 points uniform in [-1, 1]^d: fraction in ball, and outer-shell fraction of those
  d =  1: ball 1.00000 vs 1.00000 (|diff| 0.00000, tol 0.00000) shell 0.05012 vs 0.05000 (|diff| 0.00012, tol 0.00087) -> ok
  d =  2: ball 0.78636 vs 0.78540 (|diff| 0.00096, tol 0.00164) shell 0.09780 vs 0.09750 (|diff| 0.00030, tol 0.00134) -> ok
  d =  3: ball 0.52438 vs 0.52360 (|diff| 0.00078, tol 0.00200) shell 0.14334 vs 0.14263 (|diff| 0.00072, tol 0.00193) -> ok
  d =  5: ball 0.16396 vs 0.16449 (|diff| 0.00053, tol 0.00148) shell 0.22710 vs 0.22622 (|diff| 0.00088, tol 0.00413) -> ok
  d = 10: ball 0.00248 vs 0.00249 (|diff| 0.00001, tol 0.00020) shell 0.39839 vs 0.40126 (|diff| 0.00287, tol 0.03933) -> ok

Checks on the data and the regressor. The data check held: the largest overlap between spots of different cells is 0.0008 against a bar of 0.01, the smallest spot mass is 0.9964 of the untruncated mass against a bar of 0.99, and the brightest pixel in a blank cell is 0.0111 against a bar of 0.02. The check of the NumPy k-NN regression against KNeighborsRegressor held, with a largest difference of 4.4e-16 against a bar of $10^{-10}$ .

Error and neighbour distance against data. At $n = 4{,}000$ the k-NN error rises from under 0.1% of the target’s spread for $d = 1$ to 23% for $d = 8$ , on the same number of pixels and the same number of training images; the hypothesis that the error at fixed $n$ rises with $d$ held for every seed. The fitted slopes of the nearest-neighbour distance follow $-1/d$ in pixel space as in parameter space, so a thirty-two-fold increase in data shrinks the typical neighbour distance about thirtyfold for one parameter and by about a third for eight. The pixel and parameter columns of the error are nearly the same at every $n$ : on this data, a nearest-neighbour regressor working on 1,024 pixels does as well as one given the true parameters, and its rate of improvement is set by $d$ , not by 1,024. The error slopes at $d = 4$ and $d = 8$ are steeper than $-1/d$ in both spaces, which our first-order reasoning above does not predict; since the parameter columns show the same steepening, it is a property of k-NN on this target at these sample sizes, not of the rendering. The hypotheses that the pixel-space distance and error slopes lie within tolerance of $-1/d$ held for every $d$ , but their tolerances, from 0.04 to 0.24, are wide, and they were set from the parameter-space slopes (see the note on process below), so these two checks are illustrations rather than tests.

Estimating the intrinsic dimension. The check of the NumPy neighbours against scikit-learn’s held for every $d$ : distances agreed to $2 \times 10^{-10}$ with no differing indices, and the estimates from the two sets of distances agreed to $1.2 \times 10^{-9}$ , within the bar of $10^{-8}$ . The check that no set contains a nearest-neighbour distance below $10^{-6}$ failed: the smallest distance printed is $7.4 \times 10^{-8}$ . In our reading this follows from the construction rather than from duplicates, since 2,000 independent uniform draws on an interval have smallest gaps of order $1/n^2 = 2.5 \times 10^{-7}$ of its length, so a $d = 1$ set can hold close but distinct pairs below $10^{-6}$ ; the output does not show which set held the minimum, and every printed estimate is finite. On the noise-free pixel sets the estimate stays within its tolerance of the true $d$ for $d = 1, 2, 4$ : that hypothesis held, with the same caveat that its tolerance came from the parameter-space estimates. The hypothesis that the bias, estimate minus $d$ , decreases with $d$ and is negative at $d = 8$ , the direction Levina and Bickel report, failed on both parts. In our reading, two things in the output account for it. The $k - 1$ form carries the upward factor of about 1.08 derived above, which at $d = 1$ and $d = 2$ outweighs any negative bias; the bracketed $k - 2$ column, without that factor, gives 0.995 and 1.97 on the parameters. And on the parameter vectors themselves the estimate at $d = 8$ is 7.31, a negative bias, while on the pixels it is 8.62: the pixel estimates exceed the parameter estimates at $d = 4$ and $d = 8$ . In our reading, the reason is that the rendering is not an isometry, so the ratios of neighbour distances that the estimator uses are not the same in pixel space as in parameter space once neighbours are a sizeable fraction of a spot’s width apart.

Pixel noise of standard deviation 0.02, 2% of a spot’s peak, raises the estimate for every $d$ and seed, to about 104 for $d = 1$ , and that hypothesis held. Independent noise adds about $2 \times 1024 \times 0.02^2 \approx 0.82$ to every squared pixel distance, so even two identical noise-free images end up about 0.9 apart, which by the printed distances in block (a) at $n = 2{,}000$ is larger than the typical nearest-neighbour spacing for $d = 1$ , 2 and 4; at the scale of 10 to 20 neighbours the data then look like a cloud in many dimensions. This is the scale dependence Levina and Bickel describe, a noisy line looking one- or two-dimensional depending on the neighbourhood size (§3), at its extreme: the smaller the true dimension, the more the noise dominates.

Volume. The table reproduces the geometric forms: the outer 5% of the radius holds 40% of a ten-dimensional ball’s volume and all but $10^{-23}$ of a 1,024-dimensional one, and the ball fills a fraction $10^{-2.6}$ of its cube at $d = 10$ and $10^{-1220}$ at $d = 1{,}024$ . The Monte Carlo check held for every $d$ it covers: both fractions agreed within four standard errors. For data that really did fill a 1,024-dimensional cube, no sample could reach the corners of a neighbourhood; the images above occupy a manifold of at most eight dimensions inside that cube, and the nearest-neighbour results follow its dimension.

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

The renderer uses spot positions as the parameters, so the manifold is curved.
The target was chosen so that its variance is the same for every $d$ .
The error is the RMSE rather than the mean squared error, so that its reference slope is $-1/d$ .
Distances and errors are printed in both pixel and parameter space.
The Levina–Bickel estimate is printed in the $k - 1$ and $k - 2$ forms and on the parameter vectors.
A check of the NumPy k-NN regression against scikit-learn was added.
Duplicates are excluded by construction and checked, not removed.
The sample sizes are capped at 4,000 and 2,000.
There is a single noise level.
The “outer shell” is taken as the outer 5% of the radius, with the dimensions and Monte Carlo sizes shown.
A data check was added.

The tolerances of the three hypotheses on the estimate and the slopes were written as formulas before any script ran, but their values were computed by a planning script, reach.py, which printed the parameter-space slopes and estimates of the same samples before the pilot. Those checks are therefore illustrations. The same planning output showed before the pilot that the exact parameter sample does not meet the bias hypothesis; that hypothesis was kept as written. After that planning script ran and before the pilot, the bound of the blank-cell data check was corrected from $10^{-3}$ , which the construction cannot meet, to 0.02. The pilot then ran the script shown; its smallest-distance check and its bias hypothesis failed and were reported to the editor, and the first run is the pilot script unchanged.

Common Misconceptions

“Images have a million dimensions, so nothing local can work on them.” The statistical cost follows the intrinsic dimension. In the example, 1,024 pixels and eight parameters behave like eight dimensions, and one parameter like one.
“More pixels means a worse curse.” Sampling the same scene more finely adds coordinates without adding directions of variation; the k-nearest neighbours article shows the contrast between near and far neighbours staying constant as noise-free gratings are sampled on finer grids. Noise is different: every pixel’s independent noise is a new direction, and in the example a little of it inflates the estimated dimension enormously.
“An intrinsic-dimension estimate is a property of the data alone.” It depends on the scale at which neighbourhoods are taken, on the sample size and on the estimator’s normalization: Levina and Bickel report a negative bias that grows with dimension, while in the example the $k - 1$ form overestimates at small $d$ and the curved rendering pushes the estimate above $d$ at $d = 8$ .
“High-dimensional space is roomy, so random vectors are far apart in every sense.” Random directions are nearly orthogonal and their distances nearly equal (Blum et al., Theorem 2.8, §2.1), which makes distances less informative, not more.
“Mapping to a higher-dimensional feature space escapes the curse.” See the warning reported in the support vector machine article: the mapped points still lie on a low-dimensional surface.

Where It Is Used

The curse shapes how vision systems are designed more than any single algorithm:

Nonparametric recognition and retrieval. k-nearest neighbours and image retrieval depend on neighbourhoods staying meaningful, and in practice run on compact representations rather than raw pixels.
Representations of lower dimension. Principal component analysis and learned embeddings describe images with far fewer numbers than pixels; Murphy (§16.2.2) motivates learned embeddings for distance-based methods as being less subject to the curse. Metric learning learns the distance itself.
Architectures with built-in structure. Locality and weight sharing in convolutional networks, discussed in the inductive bias article, restrict the functions a model can represent, so it needs fewer examples than a model that could represent any function of the pixels.
Search. Approximate nearest neighbor search exists because exact indexes degrade with dimension; random projections, the subject of the Johnson–Lindenstrauss lemma above, appear in its hashing schemes.
Generative models and their codes. The latent space article explains why interpolation in a Gaussian latent space must respect its thin shell.
Tracking and estimation. The particle filter article describes the exponential particle counts needed for high-dimensional states.
Choosing the size of a model of the data. Levina and Bickel present intrinsic-dimension estimates as a way to choose the embedding dimension for manifold methods, either directly or to narrow the range searched by cross-validation (§1).