Manifold Learning and Nonlinear Dimensionality Reduction

A manifold is a topological space that locally resembles Euclidean space — near any point, it looks flat, even if globally it is curved or twisted, like the surface of a sphere or a rolled-up sheet of paper. The manifold hypothesis states that many real, high-dimensional datasets are not scattered arbitrarily through their ambient space but concentrate near such a low-dimensional manifold, with the ambient dimensions filled by redundant or noisy coordinates (pixels, sensor channels, one-hot codes). If this hypothesis holds, then the data has a much smaller intrinsic dimensionality $d$ than its ambient dimensionality $D$, and the goal of manifold learning is to recover a faithful $d$-dimensional coordinate chart. The canonical illustration is the Swiss roll: a two-dimensional sheet rolled into a spiral and embedded in three-dimensional space. Two points that are close in ambient (straight-line) Euclidean distance can be far apart along the sheet — they might sit on adjacent windings of the roll — while two points that are far apart in Euclidean distance can be close along the sheet's surface. Any method relying purely on ambient Euclidean distance, including PCA and classical Multidimensional Scaling, will therefore confuse "wound-up" and "unrolled" distances. Manifold learning methods explicitly correct for this by trusting Euclidean distance only locally, in small neighborhoods where the manifold really is approximately flat, and inferring global structure from how these local neighborhoods chain together. This yields the shared algorithmic template of the classical nonlinear embedding methods: (1) build a neighborhood graph from the data, typically via $k$-nearest-neighbors or an $\varepsilon$-ball; (2) define a notion of "structure to preserve" using this graph — either geodesic distances (Isomap) or local linear reconstruction weights (LLE); (3) solve a global optimization problem, always reducible to an eigenvalue decomposition, that finds low-dimensional coordinates best preserving that structure. Because these problems are convex quadratic forms constrained to orthonormal coordinates, they have a unique, globally optimal, closed-form solution — a major theoretical advantage over gradient-based neural embeddings, at the cost of being more brittle to the neighborhood graph's quality.

Step 1 — Build the neighborhood graph

Both algorithms begin identically. For each point $x_i$ in the dataset $\{x_1, \ldots, x_n\} \subset \mathbb{R}^D$, find its $k$ nearest neighbors under the Euclidean distance metric, and connect $x_i$ to each neighbor with an edge weighted by their Euclidean distance. The resulting weighted graph $G$ encodes only short-range distances, which are trustworthy because the manifold is locally flat (a first-order Taylor approximation of any smooth surface is a tangent plane).

Step 2a — Isomap: geodesics via shortest paths

Isomap's insight is that the geodesic distance — the distance measured while walking along the manifold's surface — is what should be preserved globally, not the straight-line ambient distance. It estimates this geodesic distance by the shortest path through the neighborhood graph, since a chain of short local hops approximates a walk along the curved surface. Formally, define

computed with Dijkstra's Algorithm from every node, or the all-pairs Floyd–Warshall algorithm. This produces an $n \times n$ matrix of estimated geodesic distances, which Isomap then feeds into classical multidimensional scaling (MDS) — the same eigenvalue machinery underlying PCA, applied to a distance matrix instead of a covariance matrix — to find coordinates in $\mathbb{R}^d$ whose pairwise Euclidean distances best match $d_G$. Classical MDS proceeds by double-centering the squared distance matrix to form the Gram-like matrix $B$, then eigendecomposing it:

where $D^{(2)}$ is the matrix of squared geodesic distances. The top $d$ eigenvectors of $B$, scaled by the square roots of their eigenvalues, are the output embedding coordinates.

Step 2b — LLE: local linear reconstruction

LLE takes a different but related view: instead of distances, it preserves the linear relationships among each point and its neighbors. Each point is assumed to be well-approximated by a linear (affine, weights summing to 1) combination of its $k$ nearest neighbors — true to first order on a smooth manifold. LLE first solves, for each $i$, a constrained least-squares problem for weights $W_{ij}$:

These weights are invariant to rotation, scaling, and translation of the neighborhood, so they capture purely the local geometric relationship, independent of ambient orientation. LLE then finds low-dimensional coordinates $y_i \in \mathbb{R}^d$ that are reconstructed by the same weights:

This is solved by finding the bottom $d+1$ eigenvectors of the sparse matrix $M = (I-W)^\top (I-W)$ and discarding the trivial constant eigenvector. ``` ALGORITHM Isomap(X, k, d) Input: data X = {x_1,...,x_n} in R^D, neighbors k, target dim d Output: embedding Y in R^{n x d} 1. Build k-NN graph G: for each x_i: N(i) = k nearest neighbors of x_i by Euclidean distance for each x_j in N(i): add edge (i,j) with weight ||x_i - x_j|| 2. Compute geodesic distances: D_G = all-pairs shortest paths on G # Dijkstra or Floyd-Warshall if G is disconnected: raise error or process largest component only 3. Classical MDS on D_G: S = elementwise square of D_G H = I - (1/n) ones(n,n) B = -0.5 H @ S @ H (lambda_1 >= ... >= lambda_n, v_1,...,v_n) = eigendecomposition(B) Y[:, a] = sqrt(lambda_a) * v_a for a = 1..d 4. return Y

ALGORITHM LLE(X, k, d) Input: data X = {x_1,...,x_n} in R^D, neighbors k, target dim d Output: embedding Y in R^{n x d} 1. for each x_i: N(i) = k nearest neighbors of x_i solve for w_i (indexed over N(i)) minimizing || x_i - sum_j w_ij * x_j ||^2 s.t. sum_j w_ij = 1

# closed form: build local Gram matrix C_jk = (x_i - x_j).(x_i - x_k)

# solve C w = 1 (ones vector), then normalize w to sum to 1 set W[i, N(i)] = w_i ; W[i, others] = 0 2. M = (I - W)^T @ (I - W) 3. Find eigenvectors of M with the d+1 SMALLEST eigenvalues discard the eigenvector with eigenvalue ~0 (the constant vector) Y = remaining d eigenvectors, stacked as columns 4. return Y ```

Isomap complexity. Building the $k$-NN graph over $n$ points in $D$ dimensions costs $O(n^2 D)$ with a naive scan, or $O(n \log n \cdot D)$ with tree-based structures in low-to-moderate $D$ (KD-trees degrade toward $O(nD)$ per query in high dimensions, so the naive bound is the honest worst case). Computing all-pairs shortest paths with Dijkstra's algorithm from each of $n$ sources, using a binary heap, costs $O(n(E + n)\log n)$ where $E = O(nk)$ edges, i.e. $O(n^2 k \log n)$ total (Floyd–Warshall is a simpler but slower $O(n^3)$ alternative). The final eigendecomposition of the $n \times n$ matrix $B$ for the top $d$ eigenvectors costs $O(d n^2)$ using iterative methods (Lanczos), or $O(n^3)$ for a full eigendecomposition. The dominant cost is the all-pairs shortest path step, $O(n^2 k \log n)$, making Isomap's overall complexity roughly $O(n^2 k \log n)$ for typical $k \ll n$.

LLE complexity. The $k$-NN search costs the same $O(n^2 D)$ (or better with trees). For each point, solving the local weight regression involves inverting or solving a $k \times k$ linear system, costing $O(k^3)$ per point and $O(nk^3)$ total — cheap since $k$ is a small constant (typically 5–50), independent of $n$ or $D$. The matrix $M = (I-W)^\top(I-W)$ is $n \times n$ but extremely sparse (each row has at most $k$ nonzeros before the product, $O(k^2)$ afterward), so its bottom $d+1$ eigenvectors can be found with sparse iterative eigensolvers (e.g., Lanczos/ARPACK) in roughly $O(dn k^2)$ to $O(d n^2)$ depending on conditioning, far cheaper in practice than Isomap's dense shortest-path step. LLE's dominant true bottleneck is usually the neighbor search, not the eigenproblem.

Why classical MDS is correct (Isomap's core guarantee). If the geodesic distance matrix $D_G$ were exactly the pairwise Euclidean distances of some true low-dimensional configuration $Y^\ast$, then double-centering recovers this configuration exactly, because for centered points, the squared distance matrix satisfies

and the double-centering operation $-\tfrac12 H D^{(2)} H$ exactly cancels the norm terms, leaving $B = Y^\ast (Y^\ast)^\top$, the Gram matrix. This is a rank-$d$ positive semidefinite matrix, and its eigendecomposition recovers $Y^\ast$ up to rotation — the same identity that underlies PCA acting on a covariance matrix. Isomap's correctness therefore rests on a separate geometric claim: that graph shortest-path distances converge to true geodesic distances as $n \to \infty$ and $k$ shrinks appropriately, which holds under sampling-density and curvature conditions (Bernstein, de Silva, Langford 2000 give explicit bounds). In finite samples this is an approximation, not an exact identity — the eigen-machinery is exact, but its input (the geodesic estimate) is not.

Why LLE's optimization is correct. The local weight-fitting step is an exact constrained least-squares problem with a closed-form solution (solve a small linear system, then renormalize), so no approximation error is introduced there beyond the linearity assumption itself. The embedding step is a Rayleigh–Ritz problem: minimizing $Y^\top M Y$ subject to $Y^\top Y = I$ (with the mean constraint removing the trivial null-space direction) is solved exactly by the eigenvectors of $M$ with smallest eigenvalues — this is the standard variational characterization of eigenvalues, the same principle used to justify PCA's optimality for reconstruction error. The one eigenvalue that is exactly zero corresponds to the all-ones (constant) eigenvector, which must be discarded because it embeds every point at the same location; keeping the next $d$ smallest nonzero eigenvectors gives the reconstruction-preserving coordinates with the least "distortion energy."

Tradeoffs, honestly stated. Both methods find the global optimum of their respective objectives — there is no local-minimum risk, unlike gradient descent in neural embeddings — but this global optimum can still be a poor answer if the input graph is wrong. Isomap fails catastrophically if the $k$-NN graph contains a "short-circuit" edge connecting two distant, unrelated regions of the manifold (e.g., through a hole, or across nearby-but-unconnected windings), because a single erroneous edge can corrupt many shortest-path distances. LLE is more locally robust to such edges (an error contaminates only nearby reconstruction weights) but is known to struggle when neighborhoods are sparse relative to curvature, or when $k$ exceeds $D$ (making the local reconstruction underdetermined, requiring regularization). Neither method extends naturally to new, unseen points without retraining (both are inherently transductive), and neither provides an inverse mapping back to the ambient space by default.

Consider $n = 6$ points sampled along a quarter-circle arc of radius 1 in $\mathbb{R}^2$ (a toy manifold curving in its ambient space), at angles $\theta = 0°, 18°, 36°, 54°, 72°, 90°$:

The true intrinsic coordinate is arc length: $s = \theta \cdot \pi/180$, ranging from $0$ to $1.571$. A straight-line (chord) distance from $A$ to $F$ is $\sqrt{2} \approx 1.414$, while the true geodesic (arc) distance is $\pi/2 \approx 1.571$ — a 10% underestimate that grows worse for larger arcs, exactly the distortion PCA/MDS on raw Euclidean distance would suffer.

Isomap step-by-step (using $k=2$ nearest neighbors, connecting each point to its immediate arc-neighbors): 1. Graph edges connect consecutive letters: $A$–$B$, $B$–$C$, $C$–$D$, $D$–$E$, $E$–$F$, each with Euclidean weight $\approx 0.314$ (chord length for $18°$ of arc). 2. Shortest path $A \to F$ sums five such hops: $d_G(A,F) \approx 5 \times 0.314 = 1.571$ — matching the true arc length almost exactly, because chords over small angles closely approximate arc length. Compare to the raw chord distance $1.414$: the graph geodesic has already corrected the curvature-induced underestimate. 3. Feeding the full $6\times 6$ geodesic distance matrix into classical MDS and keeping $d=1$ dimension recovers coordinates that are, up to a sign and additive constant, linear in $\theta$ — i.e., Isomap linearizes the arc into a 1-D coordinate, recovering exactly the arc-length parameterization that generated the data.

LLE step-by-step on the same data with $k=2$: for point $C$ (neighbors $B$ and $D$), solve for weights $w_B, w_D$ with $w_B + w_D = 1$ minimizing $\|x_C - w_B x_B - w_D x_D\|^2$. By the arc's local symmetry around $C$, the least-squares solution gives $w_B \approx w_D \approx 0.5$ (locally the arc looks like a straight segment through $C$, symmetric about it). Repeating for all interior points and assembling $M = (I-W)^\top(I-W)$, the second-smallest eigenvector of $M$ (the first is discarded as constant) yields embedding coordinates that again vary monotonically along the arc — LLE unrolls the same curve into an ordered 1-D line, using only local symmetry rather than any global shortest-path computation. Both methods agree on the qualitative answer for this simple example — recovering $\theta$ as the intrinsic coordinate — illustrating that they solve the same underlying problem (unroll a curved manifold) via distinct mechanisms (global geodesics vs. local linear reconstruction).

:::physicist[Sam T. Roweis and Lawrence K. Saul|Roweis 1972–2010|Saul active researcher|Machine learning, NYU and UCSD|Published Locally Linear Embedding in the same 2000 issue of Science as Isomap| Roweis and Saul proposed LLE as a strikingly different route to the same destination: rather than global shortest paths, they showed that preserving purely local linear reconstruction weights, combined with a sparse eigenproblem, recovers global nonlinear structure — a template later echoed in graph-based semi-supervised learning and spectral clustering.

Manifold learning underlies exploratory data analysis across the sciences: single-cell genomics pipelines use nonlinear embeddings to visualize cell-type continuums that PCA smears together; computer vision researchers use manifold structure to model pose and viewpoint manifolds in face and object recognition; robotics and motion-planning systems exploit low-dimensional configuration-space manifolds to plan efficiently in what is nominally a high-dimensional joint-angle space. The conceptual apparatus here — trust local distances, infer global structure — is the direct ancestor of the neighbor-graph constructions inside t-SNE and UMAP, both of which replace Isomap's rigid geodesic distances and LLE's linear reconstruction weights with softer, probabilistic neighborhood models, trading exact global optimality for scalability and better handling of noisy real-world neighborhoods. Spectral clustering, graph-based semi-supervised learning, and even certain graph neural network architectures inherit the same "build a graph, extract eigenvectors" template introduced by this generation of algorithms.

Ask the tutor: How does the choice of the neighborhood size $k$ concretely change the shape of an Isomap or LLE embedding on a dataset with varying sampling density?