Manifold Learning: MDS, Isomap, LLE, and Diffusion Maps

In 2005 the US House of Representatives held 671 roll calls among 435 members. Diaconis and colleagues dropped the first two — a call of the House by States and the election of the Speaker — and kept the 401 Representatives who voted on at least 90% of the rest: 220 Republicans, 180 Democrats and 1 Independent, Rep. Bernie Sanders. Classical MDS on that 401 x 669 matrix returned a three-dimensional embedding. It separated into what they call “twin horseshoes”, one per party.

The curve is the kind of picture that gets a story attached to it: a political spectrum bending back on itself, moderates meeting at the far end. The story is unnecessary, and in 2008 they proved it. Their abstract: “We show that, in general, a latent ordering of the data gives rise to these patterns when one only has local information. That is, when only the interpoint distances for nearby points are known accurately.” The horseshoe is a property of the distance model, not of the legislators.

Meilă and Zhang state the reading rule that follows: “when an embedding looks like a 'horseshoe', this may not represent a property of the data, but an artifact signalling that one of the data dimensions is collapsed, or poorly reflected in the embedding (Diaconis et al., 2008)”.

Hand the same matrix to the other families and each answers a different geometric question. Classical MDS fits coordinates to a supplied dissimilarity matrix. Isomap follows shortest paths through a neighbourhood graph. LLE preserves the weights that reconstruct each point from its neighbours. Diffusion maps preserve connectivity under repeated random-walk transitions. All four can return a plausible two-dimensional picture of the same data. Calling them collectively “nonlinear PCA” hides the assumption that decides which picture means anything.

Manifold methods differ less by their final plot than by the geometry optimized before plotting. Each objective also comes with a setting that has to be chosen before any picture exists.

Classical MDS finds coordinates whose Euclidean distances approximate a supplied dissimilarity matrix. The dissimilarity matrix is the input, which is why the 401 x 669 roll-call matrix and its twin horseshoes are an MDS result rather than a fact about the House. Isomap approximates geodesic distance through a neighbourhood graph and then embeds those graph distances; its dial is the neighbourhood, K = 6 for the 698 faces, epsilon = 5 versus 4.3–4.6 for the noisy Swiss roll. Locally linear embedding preserves the coefficients that reconstruct each point from nearby points; its dials are K and the regularizer that keeps those coefficients unique. Diffusion maps preserve connectivity under repeated random-walk transitions; their dials are the kernel width, sigma = 2 in the Goldberg experiments, and the normalization exponent alpha, which decides which operator is being approximated at all.

A method can satisfy its own objective while distorting a relation another method values. In each family the cost has been measured rather than merely predicted.

MDS fits distances from an explicit dissimilarity matrix. That makes global stress visible and admits non-Euclidean input in some variants, but it discovers no geodesics on its own. On the 2005 House matrix it produced the twin horseshoes that were later proved to come from having only local distances accurate.

Isomap targets global manifold geometry through shortest paths, and pays for it twice. De Silva and Tenenbaum named the two bottlenecks when they introduced Landmark Isomap at NIPS in 2002: the all-pairs shortest-path matrix, and an MDS eigenvalue calculation on a full N x N matrix. “In contrast, the eigenvalue computations in LLE and Laplacian Eigenmaps are sparse (hence considerably cheaper).” Their fix replaces the N x N distance matrix with an n x N matrix of distances to n landmark points. Meilă and Zhang put the asymptotics on it: “Computation complexity of Isomap is O(n3), with a most computational burden for computing all pairs of shortest path distance. Space complexity is O(n2). Since Isomap works with dense matrices, this space complexity cannot be improved.” The shipped implementation agrees, giving “The overall complexity of Isomap is O[D log(k) N log(N)] + O[N^2(k + log(N))] + O[d N^2]” against “The overall complexity of standard LLE is O[D log(k) N log(N)] + O[D N k^3] + O[d N^2]” (scikit-learn, “Manifold learning”). And Isomap remains the family that a single false edge can reorganize.

LLE focuses on local patch geometry and avoids explicit global distance fitting. Its local step, though, is ill-posed rather than merely delicate. Once the neighbourhood size K exceeds the input dimensionality D, each point is reconstructed perfectly from its neighbours and the weights stop being unique. The local Gram matrix then has to be conditioned by adding a multiple of the identity. Saul and Roweis used a regularizer of Δ = 0.1 for the synthetic manifolds of their Figure 1, where K = 8 exceeded D = 3, and reported in 2003 that “The results are stable over a wide range of values but do break down as K becomes too small or large.” The same conditioning is a tuning knob in the library: “reg float, default=1e-3 — Regularization constant, multiplies the trace of the local covariance matrix of the distances”, with n_neighbors defaulting to 5. A default of 5 neighbours and 1e-3 is a modelling decision someone already made for you.

Diffusion maps can suppress short-scale noise and offer a scale parameter through diffusion time. They are also not one construction. In PNAS in 2005, Coifman, Lafon and five colleagues showed that they are a one-parameter family indexed by a normalization exponent: “One can easily see that the interesting cases are (i) α = 0, corresponding to the classical normalized graph Laplacian; (ii) α = 1, yielding the Laplace–Beltrami operator; and (iii) α = 1/2 yielding the backward Fokker–Planck operator.” At alpha = 1 the sampling density plays no role. Meilă and Zhang record the consequence: “The seminal work Coifman & Lafon (2006), which introduced renormalization, showed that in this case, the eigenvectors of Lnorm, Lrw are biased by the sampling density and that renormalization removes this bias.” That is why every honest report of a diffusion-map coordinate states the exponent, as the grid experiments above do: “We ran the DFM algorithm on the same data sets. We used the normalization constant α = 1 and the kernel width σ = 2; the results for σ = 1,4, and 8 are similar.” Whether a coordinate reports geometry or reports where you happened to sample is decided by a number chosen before the run.

Applying one of these usefully means naming both the hypothesized manifold and the evidence that would challenge it. The published failures name that evidence for you.

The workflow separates geometric assumptions from visual preference. Every step below has a number attached to it somewhere in this lesson.

1. Define the preserved relation. Choose global dissimilarity, graph geodesic, local reconstruction or diffusion connectivity before looking at any plot. The 401 x 669 roll-call horseshoe is what happens when the relation is chosen implicitly.

2. Build competing neighbourhoods. Vary the metric, the scaling, the neighbour count, the radius and the approximate-search settings. K = 6 recovered the faces; K = 8 with D = 3 forced Saul and Roweis to add Δ = 0.1; epsilon moved from 5 to the 4.3–4.6 band between a broken and a working Swiss roll.

3. Inspect graph pathology. Measure components, degree distribution, hubs, shortcuts and isolated observations. A single short-circuit edge can alter many entries of the geodesic distance matrix at once.

4. Quantify preservation, and quantify cost. Use stress, reconstruction error, neighbour retention and domain-specific continuity checks. Budget for O(n^3) time and O(n^2) space if the method is Isomap, or move to landmarks — an n x N matrix instead of N x N.

5. Challenge the interpretation. Repeat across seeds, samples, batches and plausible preprocessing choices, and test the aspect-ratio failure specifically. An embedding that looks one-dimensional may be a two-dimensional manifold the normalized-output family refuses to report.

MDS, Isomap, LLE and diffusion maps do not offer interchangeable routes to the same hidden plane. Each constructs a different contract between the original space and the embedding, and each contract has a published price. A horseshoe that a theorem attributes to the method rather than to the House of Representatives. A two-dimensional grid collapsed to a curve once its aspect ratio passed 2, with no sample size that repairs it. An embedding reorganized by one false edge under 2% noise. Reconstruction weights that stop being unique the moment K exceeds D. A dense N x N eigenproblem whose O(n^2) space cannot be improved. A normalization exponent that decides whether a coordinate describes geometry or sampling density.

Select the method from that contract, not from the most attractive scatterplot. When no contract matches the task, a simpler representation may be more honest.