Nonlinear dimension reduction

Math 123: Mathematical Aspects of Data Analysis – Fall 2026

Published

November 5, 2026

Where we are

Last time we finished normalized spectral clustering. We embedded each vertex \(i\) as the row \((\phi_1(i),\ldots,\phi_k(i))\), where \(\bfphi_\ell=\bfD^{-1/2}\bfv_\ell\) and \(\bfv_1,\ldots,\bfv_k\) are eigenvectors of the normalized Laplacian \(\bfL_{\mathrm{sym}}=\bfI-\bfD^{-1/2}\bfW\bfD^{-1/2}\) for its \(k\) smallest eigenvalues, and then rounded the rows to a partition with \(k\)-means.

Today we start a new topic: finding low-dimensional coordinates for data that lie on a curved surface, where PCA fails. We will see two methods, which make opposite choices about what to keep.

Laplacian eigenmaps (Belkin and Niyogi 2003) keep only local information: neighbours should stay close. The method uses the same eigenvectors as spectral clustering, kept as coordinates instead of rounded to a partition. We have already done the hard work, so we move quickly.
Isomap (Tenenbaum et al. 2000) keeps global information: every pairwise distance, measured along the data. It is built on classical multidimensional scaling, which recovers points from their pairwise distances and turns out to be PCA in disguise.

Throughout, \(G\) is a weighted graph with no isolated vertices, as in Lecture 16.

Why nonlinear?

PCA finds the best subspace, so it can only find structure that is linear. Many datasets have few degrees of freedom but are curved:

A swiss roll: a rectangle of points rolled up in \(\Rbb^3\). It is intrinsically two-dimensional, but any projection to a plane stacks the layers on top of each other.
Images of a moving stripe: a white vertical stripe on a black background, at \(n\) horizontal positions, wrapping around at the edge. Each image is a vector in \(\Rbb^{\text{pixels}}\), but the data have one degree of freedom and the shape of a circle. Two images with non-overlapping stripes are equally far apart, whatever their positions, so Euclidean distance only sees which images are nearby.
Faces under changing light: one face photographed with the light at different angles, two degrees of freedom (Section 6.1.3 of the textbook).

The working assumption, the manifold hypothesis, is that the data lie on or near a smooth \(m\)-dimensional surface in \(\Rbb^d\), with \(m\ll d\). On such a surface, Euclidean distance is trustworthy only between nearby points: two points on adjacent layers of the swiss roll are close in \(\Rbb^3\) but far apart along the roll. So we keep only local information, in a similarity graph, and ask for coordinates in which neighbours stay close.

See the swiss roll example for PCA failing, and the method below succeeding, on exactly this problem.

Laplacian eigenmaps

We want to place each vertex \(i\) at a point \(\bfy_i\in\Rbb^m\), with \(m\) small, so that vertices joined by heavy edges land close together. Collect the points as the rows of \(\bfY\in\Rbb^{n\times m}\) and measure the failure by the energy \[ E(\bfY)=\frac12\sum_{i,j=1}^nw_{ij}\|\bfy_i-\bfy_j\|_2^2=\operatorname{tr}(\bfY^\top\bfL\bfY), \] which is the Laplacian quadratic form applied to each column of \(\bfY\).

On its own, \(E\) is minimized by cheating: shrink everything to \(\bfY=\mathbf{0}\), collapse every point to one point, or repeat one good column \(m\) times. The same degree-weighted constraints as for the relaxed normalized cut rule these out, \[ \bfY^\top\bfD\mathbf{1}=\mathbf{0}, \qquad \bfY^\top\bfD\bfY=\bfI_m. \] The first removes translation, and the second fixes the scale of each coordinate and stops different coordinates from repeating each other.

Theorem 1 (The smoothest embedding) Let \(\bfL_{\mathrm{sym}}\) have eigenvalues \(\lambda_1\leq\cdots\leq\lambda_n\) with orthonormal eigenvectors \(\bfv_1,\ldots,\bfv_n\), and let \(\bfphi_\ell=\bfD^{-1/2}\bfv_\ell\). For every \(\bfY\in\Rbb^{n\times m}\) with \(\bfY^\top\bfD\mathbf{1}=\mathbf{0}\) and \(\bfY^\top\bfD\bfY=\bfI_m\), \[ \frac12\sum_{i,j=1}^nw_{ij}\|\bfy_i-\bfy_j\|_2^2\geq\lambda_2+\cdots+\lambda_{m+1}, \] with equality at \(\bfY=[\bfphi_2,\ldots,\bfphi_{m+1}]\).

Proof

Proof. This is the proof of the relaxed normalized cut, with \(m\) columns instead of one. Set \(\bfZ=\bfD^{1/2}\bfY\). Then \(\bfZ^\top\bfZ=\bfI_m\), \(\bfZ^\top\bfv_1=\mathbf{0}\) since \(\bfv_1=\bfD^{1/2}\mathbf{1}/\sqrt{\operatorname{vol}(V)}\), and \(\operatorname{tr}(\bfY^\top\bfL\bfY)=\operatorname{tr}(\bfZ^\top\bfL_{\mathrm{sym}}\bfZ)\). By Rayleigh–Ritz with the constraint \(\bfZ^\top\bfv_1=\mathbf{0}\), this is at least \(\lambda_2+\cdots+\lambda_{m+1}\), with equality at \(\bfZ=[\bfv_2,\ldots,\bfv_{m+1}]\).

Definition 1 (Laplacian eigenmap) The Laplacian eigenmap of \(G\) in \(m\) dimensions sends vertex \(i\) to \[ \bfy_i=\big(\phi_2(i),\ \ldots,\ \phi_{m+1}(i)\big)\in\Rbb^m, \] where \(\bfphi_\ell=\bfD^{-1/2}\bfv_\ell\) and \(\bfv_2,\ldots,\bfv_{m+1}\) are orthonormal eigenvectors of \(\bfL_{\mathrm{sym}}\) for its smallest eigenvalues after \(\lambda_1=0\).

A few things to notice.

Only neighbours enter. The energy involves \(w_{ij}\), which is negligible unless \(i\) and \(j\) are close in the data. Nothing asks distant points to be far apart; the constraint \(\bfY^\top\bfD\bfY=\bfI_m\) does that implicitly, by forbidding collapse.
Up to rotation. For orthogonal \(\bfQ\in\Rbb^{m\times m}\), the matrix \(\bfY\bfQ\) satisfies the same constraints and has the same energy, by the cyclic property of the trace. So the eigenmap is the smoothest embedding up to a rotation or reflection of \(\Rbb^m\).
Each coordinate is ranked by roughness. Coordinate \(\bfphi_\ell\) has energy \(\bfphi_\ell^\top\bfL\bfphi_\ell=\lambda_\ell\), so the first coordinates are the smoothest functions on the graph.

Spectral clustering is k-means on the eigenmap

Compare the spectral embedding with the Laplacian eigenmap (Definition 1). The eigenmap with \(m=k-1\) coordinates is the spectral embedding with the constant first coordinate \(\phi_1(i)=1/\sqrt{\operatorname{vol}(V)}\) removed, and removing it changes no distance between rows, and so nothing about \(k\)-means. So normalized spectral clustering is \(k\)-means on the Laplacian eigenmap with \(m=k-1\) coordinates, and the same eigenvectors serve two purposes:

Use	What we do with \(\bfphi_2,\ldots,\bfphi_{m+1}\)
Clustering	round to a partition with \(k\)-means, \(m=k-1\)
Dimension reduction	plot or analyze the coordinates directly

The proof of Theorem 1 also gives the \(k\)-cluster version of the relaxed normalized cut.

Proposition 1 (Relaxed normalized cut for \(k\) clusters) For every partition of \(V\) into non-empty sets \(C_1,\ldots,C_k\), \[ \sum_{\ell=1}^k\frac{\operatorname{cut}(C_\ell)}{\operatorname{vol}(C_\ell)}\geq\lambda_1+\cdots+\lambda_k, \] and the relaxed problem, over all \(\bfY\in\Rbb^{n\times k}\) with \(\bfY^\top\bfD\bfY=\bfI_k\), is solved by \(\bfY=[\bfphi_1,\ldots,\bfphi_k]\). For \(k=2\) the left side is \(\operatorname{Ncut}(C_1)\) and the right side is \(\lambda_2\).

Proof

Proof. Let \(\bfH\in\Rbb^{n\times k}\) have columns \(\bfh_\ell=\mathbf{1}_{C_\ell}/\sqrt{\operatorname{vol}(C_\ell)}\). The columns have disjoint supports, and \(\bfh_\ell^\top\bfD\bfh_\ell=\sum_{i\in C_\ell}\deg(i)/\operatorname{vol}(C_\ell)=1\), so \(\bfH^\top\bfD\bfH=\bfI_k\). By the Laplacian quadratic form, only edges leaving \(C_\ell\) contribute to \(\bfh_\ell^\top\bfL\bfh_\ell\), each with \(1/\operatorname{vol}(C_\ell)\), so \(\operatorname{tr}(\bfH^\top\bfL\bfH)=\sum_\ell\operatorname{cut}(C_\ell)/\operatorname{vol}(C_\ell)\). Now repeat the proof of Theorem 1 without the constraint \(\bfY^\top\bfD\mathbf{1}=\mathbf{0}\): with \(\bfZ=\bfD^{1/2}\bfY\), Rayleigh–Ritz gives \(\operatorname{tr}(\bfY^\top\bfL\bfY)\geq\lambda_1+\cdots+\lambda_k\) with equality at \(\bfZ=[\bfv_1,\ldots,\bfv_k]\). Apply this to \(\bfY=\bfH\).

Laplacian eigenmaps on data

For a point cloud \(\bfx_1,\ldots,\bfx_n\in\Rbb^d\), build a similarity graph and take its eigenmap.

The graph decides everything. Too many neighbours, or too large an \(\varepsilon\) in the Gaussian weights, and the swiss roll’s layers get connected across, which the embedding then faithfully reports. Too few, and the graph falls apart into pieces, and by the disconnected case the first eigenvectors just report the pieces. A sparse \(k\)-nearest-neighbour graph is the usual choice.
Noise: pairwise distances in high dimension are dominated by noise, so a common pipeline first reduces to a few dozen principal components and then builds the graph there.
Computation: with a sparse graph, only a few of the smallest eigenvectors are needed, and an iterative sparse eigensolver finds them without forming a dense \(n\times n\) decomposition.

Laplacian eigenmaps in NumPy

The eigenmap is the spectral clustering embedding without the constant first column.

import numpy as np


def normalized_spectral_embedding(W, k):
    d = W.sum(axis=1)
    L_sym = np.eye(len(d)) - W / np.sqrt(np.outer(d, d))
    lam, V = np.linalg.eigh(L_sym)  # increasing order
    return lam[:k], V[:, :k] / np.sqrt(d)[:, None]  # phi_l = D^{-1/2} v_l


lam, Phi = normalized_spectral_embedding(W, m + 1)
Y = Phi[:, 1:]  # rows are the points y_i in R^m

scikit-learn’s SpectralEmbedding does the same, with a nearest-neighbour graph by default.

Points from distances: classical MDS

The eigenmap only ever asks neighbours to stay close. The other natural idea is to ask for all distances to be right: pick a notion of distance \(\delta_{ij}\) between data points, then find \(\bfy_1,\ldots,\bfy_n\in\Rbb^m\) with \[ \|\bfy_i-\bfy_j\|_2\approx\delta_{ij} \qquad\text{for all }i,j. \] This is multidimensional scaling (MDS). Before choosing a clever \(\delta_{ij}\), we answer a cleaner question: if the \(\delta_{ij}\) really are Euclidean distances between some points, can we recover the points? Distances do not change under translation, rotation, or reflection, so the best we can hope for is the points up to these rigid motions.

The centering matrix

Write \[ \bfJ=\bfI_n-\frac1n\mathbf{1}\mathbf{1}^\top\in\Rbb^{n\times n} \] for the centering matrix. For a vector \(\bfa\in\Rbb^n\), \(\bfJ\bfa=\bfa-\bar a\mathbf{1}\) subtracts the mean \(\bar a\) of the entries from every entry. In other words, \(\bfJ\) is the orthogonal projection onto the vectors with mean zero, so \(\bfJ^\top=\bfJ\), \(\bfJ^2=\bfJ\), and \(\bfJ\mathbf{1}=\mathbf{0}\). Applied to a data matrix \(\bfA\in\Rbb^{n\times d}\) with rows \(\bfx_i^\top\), it centers every column, so \(\bfJ\bfA\) is the centered data matrix \(\bfX\).

Distances determine the Gram matrix

For points \(\bfx_1,\ldots,\bfx_n\in\Rbb^d\), let \(\bfDelta\in\Rbb^{n\times n}\) be the matrix of squared distances, \(\Delta_{ij}=\|\bfx_i-\bfx_j\|_2^2\).

Lemma 1 (Double centering) Let \(\bfX\in\Rbb^{n\times d}\) be the centered data matrix of \(\bfx_1,\ldots,\bfx_n\), and \(\bfDelta\) their matrix of squared distances. Then \[ \bfX\bfX^\top=-\frac12\bfJ\bfDelta\bfJ. \]

Proof

Proof. Let \(\bfA\) be the uncentered data matrix, so that \((\bfA\bfA^\top)_{ij}=\bfx_i^\top\bfx_j\), and let \(\bfg\in\Rbb^n\) have entries \(g_i=\|\bfx_i\|_2^2\). Expanding the square, \[ \Delta_{ij}=\|\bfx_i\|_2^2+\|\bfx_j\|_2^2-2\bfx_i^\top\bfx_j, \qquad\text{that is,}\qquad \bfDelta=\bfg\mathbf{1}^\top+\mathbf{1}\bfg^\top-2\bfA\bfA^\top. \] Since \(\bfJ\mathbf{1}=\mathbf{0}\), multiplying by \(\bfJ\) on both sides kills the two rank-one terms, and \[ -\frac12\bfJ\bfDelta\bfJ=\bfJ\bfA\bfA^\top\bfJ=(\bfJ\bfA)(\bfJ\bfA)^\top=\bfX\bfX^\top. \]

So the distances alone determine the Gram matrix \(\bfX\bfX^\top\) of the centered points, the matrix of all inner products \((\bfx_i-\bar\bfx)^\top(\bfx_j-\bar\bfx)\). For any list of dissimilarities \(\delta_{ij}\), symmetric with \(\delta_{ii}=0\), this suggests forming \[ \bfB=-\frac12\bfJ\bfDelta\bfJ, \qquad \Delta_{ij}=\delta_{ij}^2, \] and treating \(\bfB\) as the Gram matrix of the points we are looking for. The lemma says this is exactly right when the \(\delta_{ij}\) are Euclidean, and the converse holds too.

Proposition 2 (When are distances Euclidean?) Let \(\delta_{ij}=\delta_{ji}\geq0\) with \(\delta_{ii}=0\), and \(\bfB=-\frac12\bfJ\bfDelta\bfJ\) as above. There are points \(\bfy_1,\ldots,\bfy_n\in\Rbb^p\) with \(\|\bfy_i-\bfy_j\|_2=\delta_{ij}\) for all \(i,j\) if and only if \(\bfB\) is positive semidefinite with rank at most \(p\).

Proof

Proof. If such points exist, Lemma 1 gives \(\bfB=\bfY\bfY^\top\) for their centered data matrix \(\bfY\in\Rbb^{n\times p}\), which is PSD with rank at most \(p\).

Conversely, if \(\bfB\) is PSD with rank at most \(p\), write \(\bfB=\bfQ\Lambda\bfQ^\top\) by the spectral theorem, with only the first \(p\) eigenvalues possibly nonzero, and let \(\bfy_i^\top\) be the \(i\)-th row of \(\bfY=[\sqrt{\lambda_1}\bfq_1,\ldots,\sqrt{\lambda_p}\bfq_p]\), so \(\bfY\bfY^\top=\bfB\). Then \(\|\bfy_i-\bfy_j\|_2^2=B_{ii}+B_{jj}-2B_{ij}\). Writing out \(\bfB=-\frac12\big(\bfDelta-\frac1n\bfDelta\mathbf{1}\mathbf{1}^\top-\frac1n\mathbf{1}\mathbf{1}^\top\bfDelta+\frac1{n^2}\mathbf{1}\mathbf{1}^\top\bfDelta\mathbf{1}\mathbf{1}^\top\big)\), the middle two terms contribute \(c_i+c_j\) to \(B_{ij}\), where \(c_i=\frac1n\sum_l\Delta_{il}\) (using \(\bfDelta^\top=\bfDelta\)), and the last term contributes a constant; both cancel in \(B_{ii}+B_{jj}-2B_{ij}\). What is left is \(-\frac12(\Delta_{ii}+\Delta_{jj}-2\Delta_{ij})=\delta_{ij}^2\).

This is a classical result of Schoenberg (1935) and of Young and Householder (1938). It also gives a diagnostic: a negative eigenvalue of \(\bfB\) certifies that no configuration of points in any dimension has exactly the distances \(\delta_{ij}\).

The MDS embedding

When \(\bfB\) is not PSD, or has rank larger than \(m\), we settle for the configuration whose Gram matrix \(\bfY\bfY^\top\) is closest to \(\bfB\).

Theorem 2 (Classical MDS) Let \(\bfB=-\frac12\bfJ\bfDelta\bfJ\) have eigenvalues \(\lambda_1\geq\cdots\geq\lambda_n\) with orthonormal eigenvectors \(\bfq_1,\ldots,\bfq_n\), and suppose \(\lambda_m>0\). Then \[ \min_{\bfY\in\Rbb^{n\times m}}\|\bfB-\bfY\bfY^\top\|_F \] is attained at \(\bfY=[\sqrt{\lambda_1}\bfq_1,\ldots,\sqrt{\lambda_m}\bfq_m]=\bfQ_m\Lambda_m^{1/2}\).

Proof

Proof. The candidates \(\bfP=\bfY\bfY^\top\) are exactly the PSD matrices of rank at most \(m\). Split \(\bfB=\bfB_+-\bfB_-\) into its positive and negative parts, \[ \bfB_+=\sum_{\lambda_i>0}\lambda_i\bfq_i\bfq_i^\top, \qquad \bfB_-=\sum_{\lambda_i<0}|\lambda_i|\bfq_i\bfq_i^\top, \] both PSD. Expand the square, using \(\langle\bfM,\bfN\rangle=\operatorname{tr}(\bfM^\top\bfN)\) for the inner product behind the Frobenius norm: \[ \|\bfB-\bfP\|_F^2=\|\bfB_+-\bfP\|_F^2+\|\bfB_-\|_F^2+2\operatorname{tr}(\bfB_-\bfP)-2\operatorname{tr}(\bfB_+\bfB_-). \] The last term is zero, since \(\bfB_+\bfB_-\) only involves products \(\bfq_i^\top\bfq_j=0\) with \(i\neq j\). The third term is non-negative: \(\operatorname{tr}(\bfB_-\bfP)=\sum_{\lambda_i<0}|\lambda_i|\,\bfq_i^\top\bfP\bfq_i\geq0\) since \(\bfP\) is PSD. So \(\|\bfB-\bfP\|_F^2\geq\|\bfB_+-\bfP\|_F^2+\|\bfB_-\|_F^2\).

Since \(\bfB_+\) is PSD, its eigendecomposition is an SVD, so by Eckart–Young the first term is smallest at \(\bfP=\sum_{i=1}^m\lambda_i\bfq_i\bfq_i^\top\). This \(\bfP\) equals \(\bfY\bfY^\top\) for the claimed \(\bfY\), and has \(\operatorname{tr}(\bfB_-\bfP)=0\), so it attains the lower bound.

The \(\bfY\) produced is automatically centered: \(\bfB\mathbf{1}=\mathbf{0}\) because \(\bfJ\mathbf{1}=\mathbf{0}\), so every \(\bfq_i\) with \(\lambda_i\neq0\) is orthogonal to \(\mathbf{1}\). As with the eigenmap, \(\bfY\bfQ\) is just as good for any orthogonal \(\bfQ\in\Rbb^{m\times m}\), so the answer is only determined up to rotation and reflection, which is all we could ask for.

Algorithm 1 (Classical MDS) Input: dissimilarities \(\delta_{ij}\), symmetric with \(\delta_{ii}=0\), and a target dimension \(m\).

Form \(\bfDelta\) with \(\Delta_{ij}=\delta_{ij}^2\), and \(\bfB=-\frac12\bfJ\bfDelta\bfJ\).
Compute the \(m\) largest eigenvalues \(\lambda_1\geq\cdots\geq\lambda_m\) of \(\bfB\) and orthonormal eigenvectors \(\bfq_1,\ldots,\bfq_m\).
Return the rows of \(\bfY=[\sqrt{\lambda_1}\bfq_1,\ldots,\sqrt{\lambda_m}\bfq_m]\).

MDS on Euclidean distances is PCA

Corollary 1 (Classical MDS on Euclidean distances is PCA) If \(\delta_{ij}=\|\bfx_i-\bfx_j\|_2\) for points with centered data matrix \(\bfX=\bfU_r\Sigma_r\bfV_r^\top\), then classical MDS in \(m\leq r\) dimensions returns the PCA scores \[ \bfY=\bfU_m\Sigma_m=\bfX\bfV_m, \] up to the signs of the columns (and, if singular values tie, a rotation within each tied group).

Proof

Proof. By Lemma 1 and the SVD and eigendecomposition, \(\bfB=\bfX\bfX^\top=\bfU_r\Sigma_r^2\bfU_r^\top\). So \(\lambda_i=\sigma_i^2\) and \(\bfq_i=\pm\bfu_i\), and \(\sqrt{\lambda_i}\bfq_i=\pm\sigma_i\bfu_i=\pm\bfX\bfv_i\).

So PCA only ever needed the pairwise distances, never the coordinates themselves. That is what makes the next step possible: keep the PCA machinery, but feed it a better notion of distance.

Isomap

On the swiss roll, the Euclidean distance between points on adjacent layers is small, but the distance along the roll, the geodesic distance, is large. Unrolling the roll does not stretch it, so the geodesic distance between two points equals the ordinary distance between them in the flat rectangle. So if we knew the geodesic distances, Proposition 2 says they would be Euclidean distances in \(\Rbb^2\), and Theorem 2 would recover the flat rectangle exactly, up to a rigid motion.

We do not know the surface, but Lecture 13 already suggested a substitute: give each edge of a nearest-neighbour graph the length \(\|\bfx_i-\bfx_j\|_2\), and measure distance by shortest paths. Short edges are trustworthy, because nearby points are close along the surface too, and a path of short edges follows the surface instead of cutting through empty space.

Definition 2 (Graph distance) Let \(G\) be a connected graph on \(\bfx_1,\ldots,\bfx_n\) in which each edge \(\{i,j\}\) has length \(\|\bfx_i-\bfx_j\|_2\). The graph distance \(d_G(i,j)\) is the smallest total length of a path from \(i\) to \(j\) in \(G\).

Computing all the graph distances is a classical problem: Dijkstra’s algorithm, run once from every vertex, finds them for a \(k\)-nearest-neighbour graph in about \(O(n^2(k+\log n))\) operations.

Algorithm 2 (Isomap) Input: points \(\bfx_1,\ldots,\bfx_n\in\Rbb^d\), a number of neighbours \(k\), and a target dimension \(m\).

Build the \(k\)-nearest-neighbour graph, with an edge \(\{i,j\}\) whenever \(\bfx_j\) is among the \(k\) nearest points to \(\bfx_i\) or vice versa, of length \(\|\bfx_i-\bfx_j\|_2\). Check that it is connected.
Compute all graph distances \(\delta_{ij}=d_G(i,j)\).
Run classical MDS (alg. 1) on the \(\delta_{ij}\).

Tenenbaum, de Silva, Langford, and collaborators (Tenenbaum et al. 2000) show that, as the samples get denser and the neighbourhoods smaller, graph distances converge to geodesic distances, so Isomap recovers the flat coordinates of any surface that can be unrolled without stretching onto a convex region of \(\Rbb^m\). Each hypothesis points at a way the method breaks.

Shortcuts. If \(k\) is too large, or the noise too big, one edge can jump between adjacent layers of the roll. Every shortest path that can use it does, so a single bad edge can corrupt a large fraction of the distances at once.
Holes. If the flat region is not convex, say a rectangle with a hole in it, geodesics must bend around the hole. Then the geodesic distances are not Euclidean distances in the plane, \(\bfB\) has negative eigenvalues, and the embedding is distorted near the hole.
Surfaces that cannot be flattened. A sphere cannot be unrolled onto the plane without stretching, so no planar embedding preserves its geodesic distances, and MDS returns the best compromise.
Disconnected graphs. If \(k\) is too small, some graph distances are infinite, and MDS cannot run at all.
Cost. The matrix \(\bfB\) is dense and \(n\times n\), whereas the eigenmap only needs eigenvectors of a sparse Laplacian.

Classical MDS and Isomap in NumPy

import numpy as np
from scipy.sparse.csgraph import shortest_path
from sklearn.neighbors import kneighbors_graph


def classical_mds(delta, m):
    n = len(delta)
    J = np.eye(n) - np.ones((n, n)) / n
    B = -0.5 * J @ (delta**2) @ J
    lam, Q = np.linalg.eigh(B)  # increasing order
    lam, Q = lam[::-1][:m], Q[:, ::-1][:, :m]  # m largest
    return Q * np.sqrt(lam)  # Y = Q_m Lambda_m^{1/2}; assumes lam[m-1] > 0


def isomap(X, m, k=10):
    G = kneighbors_graph(X, k, mode="distance")  # edge lengths
    delta = shortest_path(G, method="D", directed=False)  # Dijkstra
    return classical_mds(delta, m)

scikit-learn’s Isomap does the same.

Local versus global

	Laplacian eigenmap	Isomap
Keeps	neighbours close	all geodesic distances
Matrix	sparse Laplacian \(\bfL_{\mathrm{sym}}\)	dense \(\bfB=-\frac12\bfJ\bfDelta\bfJ\)
Eigenvectors	smallest, after the constant one	largest
Scale of the output	fixed by \(\bfY^\top\bfD\bfY=\bfI_m\), so not meaningful	distances are meant to be read off

Both start from the same nearest-neighbour graph, and both end in an eigenvalue problem, but they ask opposite questions of it.

References

Belkin, Mikhail, and Partha Niyogi. 2003. “Laplacian Eigenmaps for Dimensionality Reduction and Data Representation.” Neural Computation 15 (6): 1373–96.

Maaten, Laurens van der, and Geoffrey Hinton. 2008. “Visualizing Data Using t-SNE.” Journal of Machine Learning Research 9: 2579–605.

McInnes, Leland, John Healy, and James Melville. 2018. UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction. https://arxiv.org/abs/1802.03426.

Tenenbaum, Joshua B., Vin de Silva, and John C. Langford. 2000. “A Global Geometric Framework for Nonlinear Dimensionality Reduction.” Science 290 (5500): 2319–23.