Manifold Hypothesis & Curse of Dimensionality¶

Hands-on notebook covering two interconnected theory notes:

• 03-manifold-hypothesis.md — Manifold Hypothesis, Tangent Space, Local vs Global Structure
• 04-curse-of-dimensionality.md — Volume Concentration, Distance Concentration, Hubness, kNN Search Cost

Each section has a conceptual explanation, a runnable proof or demo, and at least one exercise to implement from scratch.

Conventions:

• X has shape (n_samples, n_features)
• All text is in English to match the project language rules for artifacts

---

Linked theory: research/03-manifold-hypothesis.md · research/04-curse-of-dimensionality.md · Chạy: uv run jupyter lab từ thư mục latent-anything-theory/ (cần các package trong pyproject.toml).

In [1]:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401 — registers 3-D projection
from collections import Counter

from sklearn.datasets import make_swiss_roll
from sklearn.decomposition import PCA
from sklearn.neighbors import NearestNeighbors
from sklearn.metrics import pairwise_distances

plt.style.use('seaborn-v0_8-whitegrid')
np.set_printoptions(precision=4, suppress=True)
rng = np.random.default_rng(42)

import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D # noqa: F401 — registers 3-D projection from collections import Counter from sklearn.datasets import make_swiss_roll from sklearn.decomposition import PCA from sklearn.neighbors import NearestNeighbors from sklearn.metrics import pairwise_distances plt.style.use('seaborn-v0_8-whitegrid') np.set_printoptions(precision=4, suppress=True) rng = np.random.default_rng(42)

---

Part 1 — Manifold Hypothesis¶

High-dimensional data actually lives on or near a low-dimensional manifold embedded in the ambient space.

Sections:

1. 1.1 Data on a low-dimensional submanifold — Swiss Roll demo
2. 1.2 Intrinsic dimension via Local PCA (exercise)
3. 1.3 Tangent space — local linear approximation
4. 1.4 Local structure vs global structure
5. 1.5 Euclidean distance vs geodesic distance

1.1 — Data on a Low-Dimensional Submanifold¶

The Swiss Roll is the canonical toy example:

• Ambient dimension: 3 — each point is an (x, y, z) triple
• Intrinsic dimension: 2 — every point is fully described by an angle t and a height y

The 3-D shape is just a flat 2-D sheet rolled up. Data are not spread uniformly through 3-D — they concentrate on a curved 2-D surface (the manifold).

The colour encodes the intrinsic parameter t (angle along the roll). Points with the same colour are at the same position along the intrinsic coordinate, regardless of where they appear in 3-D.

In [2]:

n_samples = 2000
X_roll, t_param = make_swiss_roll(n_samples=n_samples, noise=0.15, random_state=42)
# t_param: angle that parameterises the roll (one of the two intrinsic coordinates)

print(f'Ambient dimension   : {X_roll.shape[1]}D')
print(f'Intrinsic dimension : 2D  (angle t, height y)')
print(f'Data shape          : {X_roll.shape}')
print(f't_param range       : [{t_param.min():.2f}, {t_param.max():.2f}]')

n_samples = 2000 X_roll, t_param = make_swiss_roll(n_samples=n_samples, noise=0.15, random_state=42) # t_param: angle that parameterises the roll (one of the two intrinsic coordinates) print(f'Ambient dimension : {X_roll.shape[1]}D') print(f'Intrinsic dimension : 2D (angle t, height y)') print(f'Data shape : {X_roll.shape}') print(f't_param range : [{t_param.min():.2f}, {t_param.max():.2f}]')

Ambient dimension   : 3D
Intrinsic dimension : 2D  (angle t, height y)
Data shape          : (2000, 3)
t_param range       : [4.74, 14.13]

In [3]:

fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, projection='3d')
sc = ax.scatter(
    X_roll[:, 0], X_roll[:, 1], X_roll[:, 2],
    c=t_param, cmap='Spectral', s=12, alpha=0.7
)
plt.colorbar(sc, ax=ax, label='Intrinsic parameter t (angle)')
ax.set_title('Swiss Roll — 2D manifold embedded in 3D')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.tight_layout()
plt.show()

fig = plt.figure(figsize=(8, 6)) ax = fig.add_subplot(111, projection='3d') sc = ax.scatter( X_roll[:, 0], X_roll[:, 1], X_roll[:, 2], c=t_param, cmap='Spectral', s=12, alpha=0.7 ) plt.colorbar(sc, ax=ax, label='Intrinsic parameter t (angle)') ax.set_title('Swiss Roll — 2D manifold embedded in 3D') ax.set_xlabel('x') ax.set_ylabel('y') ax.set_zlabel('z') plt.tight_layout() plt.show()

1.2 — Intrinsic Dimension via Local PCA (Exercise 1)¶

Intrinsic dimension is the minimum number of free parameters needed to describe every point on the manifold.

Local PCA trick:

1. Take the $k$ nearest neighbours of point $p$.
2. Run PCA (via SVD) on the local neighbourhood.
3. Count how many principal components explain $\geq 95\%$ of local variance.
4. That count $\approx$ local intrinsic dimension.

The trick works because: if the neighbourhood is small enough, the manifold looks flat (linear) locally. A flat $d$-dimensional patch has exactly $d$ significant directions of variation.

In [4]:

def estimate_local_dimension(
    X: np.ndarray,
    point_idx: int,
    k: int = 30,
    variance_threshold: float = 0.95,
) -> int:
    """
    Estimate the local intrinsic dimension at X[point_idx] by running PCA
    on its k nearest neighbours.

    Returns the smallest d such that the top-d principal components explain
    >= variance_threshold of local variance.
    """
    # TODO:
    # 1. Compute Euclidean distances from X[point_idx] to all points.
    # 2. Select the k nearest neighbours (argsort, take first k).
    # 3. Extract the (k, n_features) neighbourhood matrix.
    # 4. Centre: subtract column-wise mean.
    # 5. SVD on the centred matrix — singular values s.
    # 6. Cumulative explained variance = cumsum(s^2) / sum(s^2).
    # 7. Return first index where cumulative variance >= threshold (1-based).
    # raise NotImplementedError('TODO: implement estimate_local_dimension')

    dists = np.linalg.norm(X - X[point_idx], axis=1)
    nbr_idx = np.argsort(dists)[:k]
    nbhd = X[nbr_idx]
    centred = nbhd - nbhd.mean(axis=0)
    _, s, _ = np.linalg.svd(centred, full_matrices=False)
    cumvar = np.cumsum(s ** 2) / (s ** 2).sum()
    return int(np.searchsorted(cumvar, variance_threshold) + 1)

def estimate_local_dimension( X: np.ndarray, point_idx: int, k: int = 30, variance_threshold: float = 0.95, ) -> int: """ Estimate the local intrinsic dimension at X[point_idx] by running PCA on its k nearest neighbours. Returns the smallest d such that the top-d principal components explain >= variance_threshold of local variance. """ # TODO: # 1. Compute Euclidean distances from X[point_idx] to all points. # 2. Select the k nearest neighbours (argsort, take first k). # 3. Extract the (k, n_features) neighbourhood matrix. # 4. Centre: subtract column-wise mean. # 5. SVD on the centred matrix — singular values s. # 6. Cumulative explained variance = cumsum(s^2) / sum(s^2). # 7. Return first index where cumulative variance >= threshold (1-based). # raise NotImplementedError('TODO: implement estimate_local_dimension') dists = np.linalg.norm(X - X[point_idx], axis=1) nbr_idx = np.argsort(dists)[:k] nbhd = X[nbr_idx] centred = nbhd - nbhd.mean(axis=0) _, s, _ = np.linalg.svd(centred, full_matrices=False) cumvar = np.cumsum(s ** 2) / (s ** 2).sum() return int(np.searchsorted(cumvar, variance_threshold) + 1)

In [5]:

probe_idx = rng.integers(0, len(X_roll), size=30)
local_dims = [estimate_local_dimension(X_roll, int(i), k=40) for i in probe_idx]

count_by_dim = Counter(local_dims)
top_dim, top_count = count_by_dim.most_common(1)[0]

print(f'Swiss Roll ambient dimension        : {X_roll.shape[1]}D')
print(f'Local intrinsic dim estimates       : {dict(sorted(count_by_dim.items()))}')
print(f'Most common estimate                : {top_dim}D  ({top_count}/{len(local_dims)} points)')
assert top_dim == 2, f'Expected intrinsic dim = 2, got {top_dim}'
print('Exercise 1.2 — looks good.')

probe_idx = rng.integers(0, len(X_roll), size=30) local_dims = [estimate_local_dimension(X_roll, int(i), k=40) for i in probe_idx] count_by_dim = Counter(local_dims) top_dim, top_count = count_by_dim.most_common(1)[0] print(f'Swiss Roll ambient dimension : {X_roll.shape[1]}D') print(f'Local intrinsic dim estimates : {dict(sorted(count_by_dim.items()))}') print(f'Most common estimate : {top_dim}D ({top_count}/{len(local_dims)} points)') assert top_dim == 2, f'Expected intrinsic dim = 2, got {top_dim}' print('Exercise 1.2 — looks good.')

Swiss Roll ambient dimension        : 3D
Local intrinsic dim estimates       : {2: 30}
Most common estimate                : 2D  (30/30 points)
Exercise 1.2 — looks good.

1.3 — Tangent Space $T_pM$¶

At any point $p$ on a smooth manifold $M$, the tangent space $T_pM$ is the best linear approximation of $M$ near $p$.

Formal definition:

$$T_pM = \bigl\{\, \dot{\gamma}(0) \;\big|\; \gamma : (-\varepsilon, \varepsilon) \to M,\; \gamma(0) = p,\; \gamma \text{ smooth} \,\bigr\}$$

It is the set of all velocity vectors of smooth curves on $M$ that pass through $p$.

Intuition: zoom in on $p$ until the curved manifold looks flat — the flat space you see is $T_pM$.

Why it matters for ML:

• All geometric operations (inner products, angles, distances) that only make sense in flat space are computed inside $T_pM$.
• The Riemannian metric $g_p$ is a bilinear form on $T_pM$ — not on the whole manifold.
• For a manifold embedded in $\mathbb{R}^n$, $T_pM$ is a literal linear subspace of $\mathbb{R}^n$.

Demo: The simplest non-trivial manifold is a circle $S^1$ (1-D manifold in 2-D space). Parametrisation: $\gamma(\theta) = (\cos\theta, \sin\theta)$. Tangent vector at $\theta$: $\dot{\gamma}(\theta) = (-\sin\theta, \cos\theta)$.

In [6]:

angles = np.linspace(0, 2 * np.pi, 300, endpoint=False)
circle = np.column_stack([np.cos(angles), np.sin(angles)])

theta = np.pi / 3                                              # probe point at 60 deg
p = np.array([np.cos(theta), np.sin(theta)])
tangent = np.array([-np.sin(theta), np.cos(theta)])            # unit tangent vector

fig, axes = plt.subplots(1, 2, figsize=(13, 5))

# --- global view ---
ax = axes[0]
ax.plot(circle[:, 0], circle[:, 1], 'b-', linewidth=2, alpha=0.5, label='Manifold M')
ax.scatter(*p, color='crimson', s=120, zorder=5, label=f'Point p (theta={theta:.2f} rad)')
t_global = np.linspace(-0.75, 0.75, 100)
tang_line = p[:, None] + tangent[:, None] * t_global
ax.plot(tang_line[0], tang_line[1], 'r-', linewidth=2.5, label='Tangent space T_pM')
ax.quiver(p[0], p[1], tangent[0] * 0.4, tangent[1] * 0.4,
          color='orange', scale=1, scale_units='xy', angles='xy',
          width=0.008, label='Tangent vector')
ax.set_xlim(-1.9, 1.9)
ax.set_ylim(-1.9, 1.9)
ax.set_aspect('equal')
ax.set_title('Global view — manifold is curved')
ax.legend(fontsize=9)

# --- local view: zoom in ---
ax = axes[1]
zoom = 0.28
local_mask = np.abs(angles - theta) < zoom
ax.plot(circle[local_mask, 0], circle[local_mask, 1], 'b-', linewidth=4, alpha=0.7,
        label='Local manifold patch')
ax.scatter(*p, color='crimson', s=120, zorder=5, label='Point p')
t_local = np.linspace(-0.28, 0.28, 80)
tang_local = p[:, None] + tangent[:, None] * t_local
ax.plot(tang_local[0], tang_local[1], 'r--', linewidth=2.5,
        label='T_pM (linear approx)')
ax.set_aspect('equal')
ax.set_title('Local view — manifold looks flat\n(T_pM is an excellent approximation here)')
ax.legend(fontsize=9)

plt.suptitle('Tangent Space on a Circle (1-D manifold in 2-D)', fontsize=13)
plt.tight_layout()
plt.show()

angles = np.linspace(0, 2 * np.pi, 300, endpoint=False) circle = np.column_stack([np.cos(angles), np.sin(angles)]) theta = np.pi / 3 # probe point at 60 deg p = np.array([np.cos(theta), np.sin(theta)]) tangent = np.array([-np.sin(theta), np.cos(theta)]) # unit tangent vector fig, axes = plt.subplots(1, 2, figsize=(13, 5)) # --- global view --- ax = axes[0] ax.plot(circle[:, 0], circle[:, 1], 'b-', linewidth=2, alpha=0.5, label='Manifold M') ax.scatter(*p, color='crimson', s=120, zorder=5, label=f'Point p (theta={theta:.2f} rad)') t_global = np.linspace(-0.75, 0.75, 100) tang_line = p[:, None] + tangent[:, None] * t_global ax.plot(tang_line[0], tang_line[1], 'r-', linewidth=2.5, label='Tangent space T_pM') ax.quiver(p[0], p[1], tangent[0] * 0.4, tangent[1] * 0.4, color='orange', scale=1, scale_units='xy', angles='xy', width=0.008, label='Tangent vector') ax.set_xlim(-1.9, 1.9) ax.set_ylim(-1.9, 1.9) ax.set_aspect('equal') ax.set_title('Global view — manifold is curved') ax.legend(fontsize=9) # --- local view: zoom in --- ax = axes[1] zoom = 0.28 local_mask = np.abs(angles - theta) < zoom ax.plot(circle[local_mask, 0], circle[local_mask, 1], 'b-', linewidth=4, alpha=0.7, label='Local manifold patch') ax.scatter(*p, color='crimson', s=120, zorder=5, label='Point p') t_local = np.linspace(-0.28, 0.28, 80) tang_local = p[:, None] + tangent[:, None] * t_local ax.plot(tang_local[0], tang_local[1], 'r--', linewidth=2.5, label='T_pM (linear approx)') ax.set_aspect('equal') ax.set_title('Local view — manifold looks flat\n(T_pM is an excellent approximation here)') ax.legend(fontsize=9) plt.suptitle('Tangent Space on a Circle (1-D manifold in 2-D)', fontsize=13) plt.tight_layout() plt.show()

1.4 — Local Structure vs Global Structure¶

	Local structure	Global structure
What	Tangent space at one point	How all tangent spaces connect
Shape	Always flat (Euclidean)	Curved, may have holes or handles
Tool	Riemannian metric on $T_pM$	Atlas of overlapping charts
Topology	Not captured	Fully captured

Charts and Atlas:

• A chart $\phi: U \to \mathbb{R}^d$ maps an open patch $U \subset M$ to flat coordinates.
• An atlas is a consistent collection of charts that covers all of $M$.
• Transition maps $\phi_j \circ \phi_i^{-1}$ must be smooth wherever patches overlap — this is what smooth manifold means.

Topological features visible only globally:

• Is $M$ closed (compact, no boundary) — like a sphere $S^2$?
• Does $M$ have handles (genus $g$) — like a torus $T^2$ with $g=1$?
• Is $M$ orientable — unlike a Möbius band or Klein bottle?

Why this matters for ML: t-SNE preserves local structure (neighbourhood layout) but may distort global topology (relative positions of clusters). UMAP uses a topological framework to preserve more global structure.

1.5 — Euclidean Distance vs Geodesic Distance¶

• Euclidean distance ($\ell_2$ in ambient space): straight line through the embedding space. Ignores manifold structure.
• Geodesic distance: length of the shortest path along the manifold.

On a curved manifold these can differ dramatically.

Swiss Roll example: points on adjacent coil layers may appear close in 3-D but are geodesically distant — to travel between them along the manifold you must traverse the entire spiral, not cut through the air.

Demo: find the pair (inner coil, outer coil) with the smallest 3-D Euclidean distance and compare it to the approximate geodesic distance (arc length along the spiral).

In [7]:

# Strategy: inner coil (small t) vs outer coil (large t); minimise Euclidean distance
idx_sorted_t = np.argsort(t_param)
inner_idx = idx_sorted_t[:60]    # smallest t = innermost coil layer
outer_idx = idx_sorted_t[-60:]   # largest t  = outermost coil layer

cross_dists = pairwise_distances(X_roll[inner_idx], X_roll[outer_idx])
i_best, j_best = np.unravel_index(cross_dists.argmin(), cross_dists.shape)
idx_a, idx_b = inner_idx[i_best], outer_idx[j_best]

euclidean = float(np.linalg.norm(X_roll[idx_a] - X_roll[idx_b]))

# Approximate geodesic: arc length = radius * delta_angle
r_avg = 0.5 * (
    np.sqrt(X_roll[idx_a, 0] ** 2 + X_roll[idx_a, 2] ** 2) +
    np.sqrt(X_roll[idx_b, 0] ** 2 + X_roll[idx_b, 2] ** 2)
)
geodesic_approx = float(abs(t_param[idx_b] - t_param[idx_a]) * r_avg)

print('Misleading pair — close in 3-D, far on the manifold:')
print(f'  Point A  t={t_param[idx_a]:.2f}  coords={X_roll[idx_a]}')
print(f'  Point B  t={t_param[idx_b]:.2f}  coords={X_roll[idx_b]}')
print()
print(f'  Euclidean distance (3-D)     : {euclidean:.3f}')
print(f'  Approx. geodesic distance    : {geodesic_approx:.3f}')
print(f'  Ratio geodesic / Euclidean   : {geodesic_approx / euclidean:.1f}x')
print()
print('  -> These points look close in 3-D but the path along the manifold is much longer.')

# Strategy: inner coil (small t) vs outer coil (large t); minimise Euclidean distance idx_sorted_t = np.argsort(t_param) inner_idx = idx_sorted_t[:60] # smallest t = innermost coil layer outer_idx = idx_sorted_t[-60:] # largest t = outermost coil layer cross_dists = pairwise_distances(X_roll[inner_idx], X_roll[outer_idx]) i_best, j_best = np.unravel_index(cross_dists.argmin(), cross_dists.shape) idx_a, idx_b = inner_idx[i_best], outer_idx[j_best] euclidean = float(np.linalg.norm(X_roll[idx_a] - X_roll[idx_b])) # Approximate geodesic: arc length = radius * delta_angle r_avg = 0.5 * ( np.sqrt(X_roll[idx_a, 0] ** 2 + X_roll[idx_a, 2] ** 2) + np.sqrt(X_roll[idx_b, 0] ** 2 + X_roll[idx_b, 2] ** 2) ) geodesic_approx = float(abs(t_param[idx_b] - t_param[idx_a]) * r_avg) print('Misleading pair — close in 3-D, far on the manifold:') print(f' Point A t={t_param[idx_a]:.2f} coords={X_roll[idx_a]}') print(f' Point B t={t_param[idx_b]:.2f} coords={X_roll[idx_b]}') print() print(f' Euclidean distance (3-D) : {euclidean:.3f}') print(f' Approx. geodesic distance : {geodesic_approx:.3f}') print(f' Ratio geodesic / Euclidean : {geodesic_approx / euclidean:.1f}x') print() print(' -> These points look close in 3-D but the path along the manifold is much longer.')

Misleading pair — close in 3-D, far on the manifold:
  Point A  t=4.99  coords=[ 1.5268  6.0607 -4.6736]
  Point B  t=13.85  coords=[ 3.9393  6.9152 13.14  ]

  Euclidean distance (3-D)     : 17.997
  Approx. geodesic distance    : 82.603
  Ratio geodesic / Euclidean   : 4.6x

  -> These points look close in 3-D but the path along the manifold is much longer.

In [8]:

fig = plt.figure(figsize=(9, 6))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(
    X_roll[:, 0], X_roll[:, 1], X_roll[:, 2],
    c=t_param, cmap='Spectral', s=8, alpha=0.3
)
# Highlight the two misleading points
for idx, label, color in [(idx_a, 'A (inner)', 'blue'), (idx_b, 'B (outer)', 'red')]:
    ax.scatter(*X_roll[idx], color=color, s=150, zorder=10, label=label)

# Draw a straight 3-D line (Euclidean path — cuts through space)
line_pts = np.linspace(X_roll[idx_a], X_roll[idx_b], 50)
ax.plot(line_pts[:, 0], line_pts[:, 1], line_pts[:, 2],
        'k--', linewidth=2, label=f'Euclidean ({euclidean:.1f})')

ax.set_title(f'Euclidean: {euclidean:.2f}   |   Geodesic approx: {geodesic_approx:.2f}')
ax.legend()
plt.tight_layout()
plt.show()

fig = plt.figure(figsize=(9, 6)) ax = fig.add_subplot(111, projection='3d') ax.scatter( X_roll[:, 0], X_roll[:, 1], X_roll[:, 2], c=t_param, cmap='Spectral', s=8, alpha=0.3 ) # Highlight the two misleading points for idx, label, color in [(idx_a, 'A (inner)', 'blue'), (idx_b, 'B (outer)', 'red')]: ax.scatter(*X_roll[idx], color=color, s=150, zorder=10, label=label) # Draw a straight 3-D line (Euclidean path — cuts through space) line_pts = np.linspace(X_roll[idx_a], X_roll[idx_b], 50) ax.plot(line_pts[:, 0], line_pts[:, 1], line_pts[:, 2], 'k--', linewidth=2, label=f'Euclidean ({euclidean:.1f})') ax.set_title(f'Euclidean: {euclidean:.2f} | Geodesic approx: {geodesic_approx:.2f}') ax.legend() plt.tight_layout() plt.show()

---

Part 2 — Curse of Dimensionality¶

In high dimensions, intuitions built in 2-D and 3-D break down catastrophically.

Sections:

1. 2.1 Volume Concentration — all volume migrates to the surface
2. 2.2 Distance Concentration — all pairwise distances become equal (exercise)
3. 2.3 Hubness — a few points become universal neighbours
4. 2.4 kNN Search Cost — finding neighbours requires scanning the whole space (exercise)
5. 2.5 Dimensionality reduction as the cure

2.1 — Volume Concentration¶

Consider the $d$-dimensional unit hypercube $[0,1]^d$.

The outer shell (all points within $\varepsilon$ of any face) has volume:

$$V_{\text{shell}} = 1 - (1 - 2\varepsilon)^d$$

Claim: for any fixed $\varepsilon > 0$, $V_{\text{shell}} \to 1$ as $d \to \infty$.

Proof: $(1 - 2\varepsilon)^d \leq e^{-2\varepsilon d} \to 0$. $\square$

Consequence: in high dimensions the interior of the cube is practically empty. Every data point lives near the boundary. Density estimation, interpolation, and integration become exponentially harder.

The same phenomenon holds for the hypersphere: volume concentrates in a thin equatorial band; the polar caps are exponentially small.

In [9]:

def hypercube_shell_fraction(d: int, epsilon: float) -> float:
    """Fraction of the unit d-cube volume in the outer shell of thickness epsilon."""
    inner_side = max(0.0, 1.0 - 2.0 * epsilon)
    return 1.0 - inner_side ** d


eps_vals = [0.01, 0.05, 0.10]
dims_table = [1, 2, 3, 5, 10, 20, 50, 100, 500, 1000]

header = f'{"Dim":>6}  |  ' + '  |  '.join(f'eps={e}' for e in eps_vals)
print(header)
print('-' * len(header))
for d in dims_table:
    row = '  |  '.join(f'{hypercube_shell_fraction(d, e):>9.6f}' for e in eps_vals)
    print(f'{d:>6}  |  {row}')

def hypercube_shell_fraction(d: int, epsilon: float) -> float: """Fraction of the unit d-cube volume in the outer shell of thickness epsilon.""" inner_side = max(0.0, 1.0 - 2.0 * epsilon) return 1.0 - inner_side ** d eps_vals = [0.01, 0.05, 0.10] dims_table = [1, 2, 3, 5, 10, 20, 50, 100, 500, 1000] header = f'{"Dim":>6} | ' + ' | '.join(f'eps={e}' for e in eps_vals) print(header) print('-' * len(header)) for d in dims_table: row = ' | '.join(f'{hypercube_shell_fraction(d, e):>9.6f}' for e in eps_vals) print(f'{d:>6} | {row}')

   Dim  |  eps=0.01  |  eps=0.05  |  eps=0.1
--------------------------------------------
     1  |   0.020000  |   0.100000  |   0.200000
     2  |   0.039600  |   0.190000  |   0.360000
     3  |   0.058808  |   0.271000  |   0.488000
     5  |   0.096079  |   0.409510  |   0.672320
    10  |   0.182927  |   0.651322  |   0.892626
    20  |   0.332392  |   0.878423  |   0.988471
    50  |   0.635830  |   0.994846  |   0.999986
   100  |   0.867380  |   0.999973  |   1.000000
   500  |   0.999959  |   1.000000  |   1.000000
  1000  |   1.000000  |   1.000000  |   1.000000

In [10]:

d_range = np.arange(1, 251)
colors_eps = ['#4e79a7', '#f28e2b', '#e15759']

fig, ax = plt.subplots(figsize=(9, 5))
for eps, color in zip(eps_vals, colors_eps):
    fracs = [hypercube_shell_fraction(d, eps) for d in d_range]
    ax.plot(d_range, fracs, color=color, linewidth=2, label=f'eps = {eps}')

ax.axhline(0.99, color='gray', linestyle='--', alpha=0.7, label='99% threshold')
ax.set_xlabel('Dimension d')
ax.set_ylabel('Fraction of volume in outer shell')
ax.set_title('Volume Concentration: all volume migrates to the surface as d grows')
ax.legend()
plt.tight_layout()
plt.show()

d_range = np.arange(1, 251) colors_eps = ['#4e79a7', '#f28e2b', '#e15759'] fig, ax = plt.subplots(figsize=(9, 5)) for eps, color in zip(eps_vals, colors_eps): fracs = [hypercube_shell_fraction(d, eps) for d in d_range] ax.plot(d_range, fracs, color=color, linewidth=2, label=f'eps = {eps}') ax.axhline(0.99, color='gray', linestyle='--', alpha=0.7, label='99% threshold') ax.set_xlabel('Dimension d') ax.set_ylabel('Fraction of volume in outer shell') ax.set_title('Volume Concentration: all volume migrates to the surface as d grows') ax.legend() plt.tight_layout() plt.show()

2.2 — Distance Concentration (Exercise 2)¶

In high dimensions, the relative contrast of pairwise distances collapses:

$$\text{RC}(X) = \frac{d_{\max} - d_{\min}}{d_{\min}} \xrightarrow{d \to \infty} 0$$

where $d_{\max}$ and $d_{\min}$ are the maximum and minimum pairwise Euclidean distances.

When $\text{RC} \to 0$, every point is approximately equidistant from every other — nearest-neighbour search loses meaning because there is no near anymore.

Formal result (Beyer et al., 1999): For i.i.d. samples under mild conditions:

$$\frac{d_{\max} - d_{\min}}{d_{\min}} \xrightarrow{p} 0 \quad \text{as } d \to \infty$$

The $\ell_p$ norm concentrates around its expected value $E[\|\mathbf{x}\|_p]$.

In [11]:

def relative_contrast(X: np.ndarray) -> float:
    """
    Compute the relative contrast of pairwise distances:

        (d_max - d_min) / d_min

    where d_max and d_min are the largest and smallest pairwise Euclidean
    distances among all n*(n-1)/2 unique point pairs.

    A value near 0 means all distances are indistinguishable.
    """
    # TODO:
    # 1. Compute the full (n, n) pairwise distance matrix D.
    #    Hint: sklearn pairwise_distances(X) already imported.
    # 2. Extract upper-triangle entries (k=1) to get unique pair distances.
    # 3. Return (max - min) / min.
    # raise NotImplementedError('TODO: implement relative_contrast')

    D = pairwise_distances(X)
    upper = D[np.triu_indices(len(X), k=1)]
    return float((upper.max() - upper.min()) / upper.min())

def relative_contrast(X: np.ndarray) -> float: """ Compute the relative contrast of pairwise distances: (d_max - d_min) / d_min where d_max and d_min are the largest and smallest pairwise Euclidean distances among all n*(n-1)/2 unique point pairs. A value near 0 means all distances are indistinguishable. """ # TODO: # 1. Compute the full (n, n) pairwise distance matrix D. # Hint: sklearn pairwise_distances(X) already imported. # 2. Extract upper-triangle entries (k=1) to get unique pair distances. # 3. Return (max - min) / min. # raise NotImplementedError('TODO: implement relative_contrast') D = pairwise_distances(X) upper = D[np.triu_indices(len(X), k=1)] return float((upper.max() - upper.min()) / upper.min())

In [12]:

X_2d   = rng.standard_normal((200, 2))
X_100d = rng.standard_normal((200, 100))

rc_2d   = relative_contrast(X_2d)
rc_100d = relative_contrast(X_100d)

print(f'Relative contrast in   2D : {rc_2d:.4f}')
print(f'Relative contrast in 100D : {rc_100d:.4f}')
assert rc_2d > rc_100d, 'High-d should have lower relative contrast'
print('Exercise 2.2 — looks good.')

X_2d = rng.standard_normal((200, 2)) X_100d = rng.standard_normal((200, 100)) rc_2d = relative_contrast(X_2d) rc_100d = relative_contrast(X_100d) print(f'Relative contrast in 2D : {rc_2d:.4f}') print(f'Relative contrast in 100D : {rc_100d:.4f}') assert rc_2d > rc_100d, 'High-d should have lower relative contrast' print('Exercise 2.2 — looks good.')

Relative contrast in   2D : 449.9325
Relative contrast in 100D : 0.7183
Exercise 2.2 — looks good.

In [13]:

n_pts_rc = 300
dims_rc  = [2, 5, 10, 20, 50, 100, 200, 500]
contrasts = []

for d in dims_rc:
    Xd = rng.standard_normal((n_pts_rc, d))
    contrasts.append(relative_contrast(Xd))

fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(dims_rc, contrasts, 'o-', color='steelblue', linewidth=2, markersize=7)
ax.set_xlabel('Dimension d')
ax.set_ylabel('Relative contrast  (d_max - d_min) / d_min')
ax.set_title('Distance Concentration: contrast collapses in high dimensions\n'
             '(-> 0 means all pairwise distances are equal; kNN loses meaning)')
ax.set_xscale('log')
plt.tight_layout()
plt.show()

n_pts_rc = 300 dims_rc = [2, 5, 10, 20, 50, 100, 200, 500] contrasts = [] for d in dims_rc: Xd = rng.standard_normal((n_pts_rc, d)) contrasts.append(relative_contrast(Xd)) fig, ax = plt.subplots(figsize=(9, 5)) ax.plot(dims_rc, contrasts, 'o-', color='steelblue', linewidth=2, markersize=7) ax.set_xlabel('Dimension d') ax.set_ylabel('Relative contrast (d_max - d_min) / d_min') ax.set_title('Distance Concentration: contrast collapses in high dimensions\n' '(-> 0 means all pairwise distances are equal; kNN loses meaning)') ax.set_xscale('log') plt.tight_layout() plt.show()

2.3 — Hubness¶

Define $N_k(p)$ = number of times point $p$ appears in the $k$-NN list of other points.

• Low dimensions: $N_k(p) \approx k$ for all $p$ (roughly equal).
• High dimensions: the $N_k$ distribution becomes heavily right-skewed:
  • A few hub points appear in almost everyone's neighbourhood.
  • Many points never appear as anyone's neighbour (isolated).

This is not a sampling artefact — it is a deterministic consequence of geometry.

Result (Radovanović et al., 2010): The skewness of $N_k$ grows monotonically with $d$ even for perfectly uniform data.

Impact: kNN classification and retrieval degrade; cluster centroids become pathological hubs; recommendation lists are dominated by a handful of items.

In [14]:

def knn_occurrences(X: np.ndarray, k: int = 10) -> np.ndarray:
    """Return how many times each point appears in others' k-NN lists."""
    nbrs = NearestNeighbors(n_neighbors=k + 1, algorithm='auto').fit(X)
    _, indices = nbrs.kneighbors(X)
    # indices[:, 0] is the point itself — drop it
    return np.bincount(indices[:, 1:].ravel(), minlength=len(X))


n_pts_hub = 600
k_hub     = 10
dims_hub  = [2, 10, 100, 1000]

fig, axes = plt.subplots(1, len(dims_hub), figsize=(16, 4), sharey=False)
for ax, d in zip(axes, dims_hub):
    Xd     = rng.standard_normal((n_pts_hub, d))
    counts = knn_occurrences(Xd, k=k_hub)
    ax.hist(counts, bins=30, color='#4e79a7', edgecolor='white', alpha=0.85)
    ax.axvline(k_hub, color='crimson', linestyle='--', linewidth=1.5,
               label=f'Expected = {k_hub}')
    ax.set_title(f'd = {d}', fontsize=12)
    ax.set_xlabel('Occurrences in kNN lists')
    ax.legend(fontsize=8)

axes[0].set_ylabel('Number of points')
fig.suptitle(
    'Hubness: N_k distribution becomes right-skewed in high dimensions\n'
    '(red line = expected count if all points were equally central)',
    y=1.03,
)
plt.tight_layout()
plt.show()

def knn_occurrences(X: np.ndarray, k: int = 10) -> np.ndarray: """Return how many times each point appears in others' k-NN lists.""" nbrs = NearestNeighbors(n_neighbors=k + 1, algorithm='auto').fit(X) _, indices = nbrs.kneighbors(X) # indices[:, 0] is the point itself — drop it return np.bincount(indices[:, 1:].ravel(), minlength=len(X)) n_pts_hub = 600 k_hub = 10 dims_hub = [2, 10, 100, 1000] fig, axes = plt.subplots(1, len(dims_hub), figsize=(16, 4), sharey=False) for ax, d in zip(axes, dims_hub): Xd = rng.standard_normal((n_pts_hub, d)) counts = knn_occurrences(Xd, k=k_hub) ax.hist(counts, bins=30, color='#4e79a7', edgecolor='white', alpha=0.85) ax.axvline(k_hub, color='crimson', linestyle='--', linewidth=1.5, label=f'Expected = {k_hub}') ax.set_title(f'd = {d}', fontsize=12) ax.set_xlabel('Occurrences in kNN lists') ax.legend(fontsize=8) axes[0].set_ylabel('Number of points') fig.suptitle( 'Hubness: N_k distribution becomes right-skewed in high dimensions\n' '(red line = expected count if all points were equally central)', y=1.03, ) plt.tight_layout() plt.show()

In [15]:

print(f'{"Dim":>6} | {"Max hub count":>14} | {"Isolated (0 occ)":>18} | {"Skewness":>10}')
print('-' * 58)
for d in [2, 10, 50, 100, 500, 1000]:
    Xd = rng.standard_normal((n_pts_hub, d))
    c  = knn_occurrences(Xd, k=k_hub).astype(float)
    isolated_pct = 100.0 * np.mean(c == 0)
    mu, sigma = c.mean(), c.std()
    skew = float(np.mean(((c - mu) / (sigma + 1e-12)) ** 3))
    print(f'{d:>6} | {int(c.max()):>14} | {isolated_pct:>17.1f}% | {skew:>10.3f}')

print(f'{"Dim":>6} | {"Max hub count":>14} | {"Isolated (0 occ)":>18} | {"Skewness":>10}') print('-' * 58) for d in [2, 10, 50, 100, 500, 1000]: Xd = rng.standard_normal((n_pts_hub, d)) c = knn_occurrences(Xd, k=k_hub).astype(float) isolated_pct = 100.0 * np.mean(c == 0) mu, sigma = c.mean(), c.std() skew = float(np.mean(((c - mu) / (sigma + 1e-12)) ** 3)) print(f'{d:>6} | {int(c.max()):>14} | {isolated_pct:>17.1f}% | {skew:>10.3f}')

   Dim |  Max hub count |   Isolated (0 occ) |   Skewness
----------------------------------------------------------
     2 |             18 |               0.0% |     -0.439
    10 |             37 |               2.8% |      0.941
    50 |            142 |              11.8% |      3.956
   100 |            196 |              14.5% |      4.818
   500 |            159 |              18.5% |      3.764
  1000 |            197 |              22.7% |      4.411

What you see¶

1. What the plot shows: Each histogram counts how often every point appears in other points' 10-NN lists, for dimensions $d = 2, 10, 100, 1000$. At $d=2$ the counts cluster tightly around the expected value of 10. As $d$ grows the distribution becomes increasingly right-skewed: a long tail of "hub" points that appear in many neighbourhoods, plus a pile-up at low counts for points that are almost never anyone's neighbour. The stats table confirms it — the maximum hub count and the skewness both climb with $d$, and the share of isolated (0-occurrence) points grows.

2. Why it looks this way: In high dimensions distances concentrate (previous experiment), so neighbour relations become unstable and a few points that sit slightly closer to the global mean get pulled into a disproportionate number of neighbour lists. This hubness is a deterministic geometric effect (Radovanović et al., 2010), not a sampling artefact — it shows up even for perfectly uniform or Gaussian data.

3. Key takeaway: kNN-based methods (classification, retrieval, clustering, recommendation) degrade in high dimensions not only because distances blur, but because neighbourhoods become lopsided — dominated by hubs while many points drop out entirely. This is a strong reason to reduce dimensionality before relying on neighbour structure.

2.4 — kNN Search Cost Explosion (Exercise 3)¶

To find $k$ nearest neighbours among $n$ points uniformly distributed in $[0,1]^d$, the side length of the hypercube window that must be searched is:

$$\ell(d) = \left(\frac{k}{n}\right)^{1/d}$$

Derivation: the window must cover fraction $k/n$ of total volume. $\ell^d = k/n \;\Rightarrow\; \ell = (k/n)^{1/d}$.

As $d \to \infty$:

$$\ell(d) = \left(\frac{k}{n}\right)^{1/d} = e^{\frac{\ln(k/n)}{d}} \xrightarrow{d \to \infty} e^0 = 1$$

The window must cover 100 % of the space — kNN degenerates to a full linear scan regardless of the data structure (tree, hash, quantised index, etc.).

In [16]:

def knn_search_side(d: int, k: int, n: int) -> float:
    """
    Return the side length of the hypercube window needed to expect k nearest
    neighbours among n uniform points in [0,1]^d.

    Derived from:  side^d = k/n
    """
    # TODO: implement the one-line formula above.
    # raise NotImplementedError('TODO: implement knn_search_side')
    return (k / n) ** (1.0 / d)


k_knn, n_knn = 10, 10_000

print(f'k = {k_knn} neighbours,  n = {n_knn} points')
print()
print(f'{"Dim":>5} | {"Side length":>12} | {"% of space":>12}')
print('-' * 38)
for d in [1, 2, 3, 5, 10, 20, 50, 100, 200, 500]:
    side = knn_search_side(d, k_knn, n_knn)
    print(f'{d:>5} | {side:>12.6f} | {side * 100:>11.2f}%')

def knn_search_side(d: int, k: int, n: int) -> float: """ Return the side length of the hypercube window needed to expect k nearest neighbours among n uniform points in [0,1]^d. Derived from: side^d = k/n """ # TODO: implement the one-line formula above. # raise NotImplementedError('TODO: implement knn_search_side') return (k / n) ** (1.0 / d) k_knn, n_knn = 10, 10_000 print(f'k = {k_knn} neighbours, n = {n_knn} points') print() print(f'{"Dim":>5} | {"Side length":>12} | {"% of space":>12}') print('-' * 38) for d in [1, 2, 3, 5, 10, 20, 50, 100, 200, 500]: side = knn_search_side(d, k_knn, n_knn) print(f'{d:>5} | {side:>12.6f} | {side * 100:>11.2f}%')

k = 10 neighbours,  n = 10000 points

  Dim |  Side length |   % of space
--------------------------------------
    1 |     0.001000 |        0.10%
    2 |     0.031623 |        3.16%
    3 |     0.100000 |       10.00%
    5 |     0.251189 |       25.12%
   10 |     0.501187 |       50.12%
   20 |     0.707946 |       70.79%
   50 |     0.870964 |       87.10%
  100 |     0.933254 |       93.33%
  200 |     0.966051 |       96.61%
  500 |     0.986279 |       98.63%

In [17]:

assert np.isclose(knn_search_side(1, 10, 10_000), 10 / 10_000)
assert np.isclose(knn_search_side(2, 10, 10_000), (10 / 10_000) ** 0.5)
assert knn_search_side(1000, 10, 10_000) > 0.99
print('Exercise 2.4 — looks good.')

assert np.isclose(knn_search_side(1, 10, 10_000), 10 / 10_000) assert np.isclose(knn_search_side(2, 10, 10_000), (10 / 10_000) ** 0.5) assert knn_search_side(1000, 10, 10_000) > 0.99 print('Exercise 2.4 — looks good.')

Exercise 2.4 — looks good.

In [18]:

dims_knn = np.arange(1, 201)
sides    = [knn_search_side(d, k_knn, n_knn) for d in dims_knn]

fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(dims_knn, [s * 100 for s in sides], color='firebrick', linewidth=2)
ax.fill_between(dims_knn, [s * 100 for s in sides], 0, alpha=0.10, color='firebrick')
ax.axhline(99, color='gray', linestyle='--', alpha=0.7, label='99% of space')
ax.set_xlabel('Dimension d')
ax.set_ylabel(f'% of space searched  (k={k_knn}, n={n_knn:,})')
ax.set_title('kNN Search Cost: neighbourhood covers the entire space in high dimensions')
ax.legend()
plt.tight_layout()
plt.show()

dims_knn = np.arange(1, 201) sides = [knn_search_side(d, k_knn, n_knn) for d in dims_knn] fig, ax = plt.subplots(figsize=(9, 5)) ax.plot(dims_knn, [s * 100 for s in sides], color='firebrick', linewidth=2) ax.fill_between(dims_knn, [s * 100 for s in sides], 0, alpha=0.10, color='firebrick') ax.axhline(99, color='gray', linestyle='--', alpha=0.7, label='99% of space') ax.set_xlabel('Dimension d') ax.set_ylabel(f'% of space searched (k={k_knn}, n={n_knn:,})') ax.set_title('kNN Search Cost: neighbourhood covers the entire space in high dimensions') ax.legend() plt.tight_layout() plt.show()

2.5 — Dimensionality Reduction as the Cure¶

All four curses share the same root cause: too many dimensions relative to the true structure of the data.

The fix: find and use the intrinsic low-dimensional representation.

Method	What it optimises	Linear?	Preserves
PCA	Variance	Yes	Global covariance structure
UMAP	Topological neighbourhood graph	No	Local + some global
t-SNE	KL divergence of neighbour distributions	No	Local clusters
Autoencoder	Reconstruction loss	No	Task-dependent

PCA is the simplest baseline. For the Swiss Roll (true intrinsic dim = 2), PCA partially unrolls it but cannot fully recover intrinsic coordinates because the manifold is non-linearly embedded. A nonlinear method (UMAP, Isomap) is needed to fully unroll the spiral.

The connection with Part 1: the manifold hypothesis justifies dimensionality reduction as a cure for the curse. If data truly lives on a $d$-dimensional manifold, projecting to $d$ intrinsic coordinates restores all geometric properties and makes kNN meaningful again.

In [19]:

pca_full = PCA(random_state=42).fit(X_roll)
explained   = pca_full.explained_variance_ratio_
cumulative  = np.cumsum(explained)

fig, axes = plt.subplots(1, 2, figsize=(13, 5))

# Scree plot
ax = axes[0]
ax.bar(range(1, len(explained) + 1), explained * 100,
       color='#4e79a7', alpha=0.8, label='Individual PC')
ax.step(range(1, len(cumulative) + 1), cumulative * 100,
        color='crimson', linewidth=2, where='mid', label='Cumulative')
ax.axhline(95, color='gray', linestyle='--', alpha=0.6, label='95% threshold')
ax.set_xlabel('Principal component')
ax.set_ylabel('Explained variance (%)')
ax.set_title('Swiss Roll — PCA Scree Plot')
ax.legend()

# 2-D projection
pca2  = PCA(n_components=2, random_state=42)
X_2d  = pca2.fit_transform(X_roll)

ax = axes[1]
sc = ax.scatter(X_2d[:, 0], X_2d[:, 1], c=t_param, cmap='Spectral', s=12, alpha=0.7)
plt.colorbar(sc, ax=ax, label='Intrinsic parameter t')
ax.set_xlabel('PC 1')
ax.set_ylabel('PC 2')
ax.set_title(
    f'PCA -> 2D  ({cumulative[1]*100:.1f}% variance)\n'
    '(linear: partially unrolls but spiral folds remain)'
)

plt.tight_layout()
plt.show()

print(f'PC1: {explained[0]*100:.1f}%   PC2: {explained[1]*100:.1f}%')
print(f'Total with 2 PCs: {cumulative[1]*100:.1f}%  (true intrinsic dim = 2)')

pca_full = PCA(random_state=42).fit(X_roll) explained = pca_full.explained_variance_ratio_ cumulative = np.cumsum(explained) fig, axes = plt.subplots(1, 2, figsize=(13, 5)) # Scree plot ax = axes[0] ax.bar(range(1, len(explained) + 1), explained * 100, color='#4e79a7', alpha=0.8, label='Individual PC') ax.step(range(1, len(cumulative) + 1), cumulative * 100, color='crimson', linewidth=2, where='mid', label='Cumulative') ax.axhline(95, color='gray', linestyle='--', alpha=0.6, label='95% threshold') ax.set_xlabel('Principal component') ax.set_ylabel('Explained variance (%)') ax.set_title('Swiss Roll — PCA Scree Plot') ax.legend() # 2-D projection pca2 = PCA(n_components=2, random_state=42) X_2d = pca2.fit_transform(X_roll) ax = axes[1] sc = ax.scatter(X_2d[:, 0], X_2d[:, 1], c=t_param, cmap='Spectral', s=12, alpha=0.7) plt.colorbar(sc, ax=ax, label='Intrinsic parameter t') ax.set_xlabel('PC 1') ax.set_ylabel('PC 2') ax.set_title( f'PCA -> 2D ({cumulative[1]*100:.1f}% variance)\n' '(linear: partially unrolls but spiral folds remain)' ) plt.tight_layout() plt.show() print(f'PC1: {explained[0]*100:.1f}% PC2: {explained[1]*100:.1f}%') print(f'Total with 2 PCs: {cumulative[1]*100:.1f}% (true intrinsic dim = 2)')

PC1: 40.2%   PC2: 31.3%
Total with 2 PCs: 71.6%  (true intrinsic dim = 2)

---

Summary¶

Part 1 — Manifold Hypothesis¶

Concept	Key takeaway
Manifold hypothesis	Real data concentrates on a low-dim submanifold — not filling the ambient space
Intrinsic vs ambient dim	True degrees of freedom $\ll$ embedding dimension; recoverable via local PCA
Tangent space $T_pM$	Best local linear approximation; host of the Riemannian metric
Local vs global structure	Local = flat (Euclidean patch); global = curved (topology matters)
Geodesic vs Euclidean	Euclidean ignores curvature — can be catastrophically misleading

Part 2 — Curse of Dimensionality¶

Phenomenon	Effect	Fix
Volume concentration	Interior is empty — data lives near the surface	Reduce dim first
Distance concentration	All pairwise distances equalise — kNN meaningless	Use intrinsic coordinates
Hubness	A few points dominate everyone's kNN lists	Reduce dim, then index
kNN search cost	Must scan entire space to find neighbours	Reduce dim first

The connection: the manifold hypothesis justifies dimensionality reduction as a cure for the curse. Projecting to intrinsic coordinates restores geometry and makes kNN meaningful.