Diffusion Maps: A Spectral Approach

• Diffusion maps are a spectral manifold-learning method that leverages random-walk operators to capture intrinsic geometry and stochastic dynamics.
• They compute leading eigenpairs to approximate diffusion distances, preserving connectivity through multiple short paths in the data.
• The technique is versatile, extending to applications like dimensionality reduction, state estimation, and PDE discretization in high-dimensional settings.

Diffusion maps are a spectral manifold-learning method that builds coordinates adapted to the intrinsic geometry and stochastic dynamics of high-dimensional data. The construction starts from a kernel on sampled data, defines a diffusion or random-walk operator on the resulting graph, and uses its leading eigenpairs to represent smooth, slowly varying modes of the underlying space. Euclidean distances in the embedding approximate diffusion distances, so the method captures connectivity through many short paths rather than only ambient Euclidean proximity. At the same time, diffusion maps are best understood as a data-driven spectral representation of intrinsic geometry rather than as a complete charting method: the appropriate low-dimensional parametrization can lie in the span of multiple diffusion coordinates rather than coincide with the leading coordinates themselves (Shnitzer et al., 2017, Candanedo et al., 30 Mar 2026).

1. Core construction and diffusion geometry

The standard setting assumes data points $x_1,\dots,x_n$ sampled from, or observed through a smooth function of, a low-dimensional Riemannian manifold $M$ embedded in a high-dimensional ambient space. A common affinity is the Gaussian kernel

$$k_\varepsilon(x,y)=\exp\!\Big(-\frac{\|x-y\|^2}{4\varepsilon}\Big),$$

together with the empirical density estimate

$q(x_i)=\sum_{j=1}^n k_\varepsilon(x_i,x_j).$

To mitigate non-uniform sampling, diffusion maps use the Coifman–Lafon reweighting

$$k_\varepsilon^{(\alpha)}(x_i,x_j)=\frac{k_\varepsilon(x_i,x_j)}{q(x_i)^\alpha q(x_j)^\alpha},$$

followed by row-normalization,

$$p(x_i,x_j)=\frac{k_\varepsilon^{(\alpha)}(x_i,x_j)}{\sum_{\ell=1}^n k_\varepsilon^{(\alpha)}(x_i,x_\ell)}.$$

Different conventions for the kernel width appear in the literature, but the role of $\varepsilon$ is consistently to set locality and spectral resolution (Shnitzer et al., 2017, Candanedo et al., 30 Mar 2026, Beier et al., 28 Jan 2026).

The central spectral problem is

$$P\phi_i=\lambda_i\phi_i,\qquad 1=\lambda_0\ge \lambda_1\ge \lambda_2\ge \cdots .$$

The leading nontrivial eigenvectors capture large-scale geometry and slow dynamics. For diffusion time $t$, the diffusion distance is

$$D_t(x,y)^2=\sum_{i\ge 1}\lambda_i^{2t}\big(\phi_i(x)-\phi_i(y)\big)^2,$$

and the diffusion map embedding is

$M$0

When sufficiently many modes are retained, Euclidean distance in $M$1-space equals diffusion distance; truncating to the leading $M$2 coordinates yields a low-dimensional representation that preserves diffusion distances at scale $M$3 (Shnitzer et al., 2017, Candanedo et al., 30 Mar 2026).

Many implementations diagonalize the similar symmetric matrix

$M$4

rather than the non-symmetric Markov matrix $M$5. This preserves the spectrum while simplifying numerical linear algebra. The same spectral machinery underlies closely related methods such as normalized graph Laplacians and Laplacian Eigenmaps, but diffusion maps emphasize a Markov-semigroup interpretation and a multi-scale notion of geometry through the diffusion time $M$6 (Candanedo et al., 30 Mar 2026, Shan et al., 2022).

2. Continuum limits and stochastic interpretations

As $M$7 and $M$8, with appropriate normalization, diffusion maps approximate differential operators on the underlying manifold. With $M$9, the generator converges to the Laplace–Beltrami operator. With other choices of $$k_\varepsilon(x,y)=\exp\!\Big(-\frac{\|x-y\|^2}{4\varepsilon}\Big),$$0, the limit becomes a weighted Laplacian or a drift–diffusion operator tied to the sampling distribution. A standard form is

$$k_\varepsilon(x,y)=\exp\!\Big(-\frac{\|x-y\|^2}{4\varepsilon}\Big),$$1

with stationary density $$k_\varepsilon(x,y)=\exp\!\Big(-\frac{\|x-y\|^2}{4\varepsilon}\Big),$$2. In overdamped Langevin settings, $$k_\varepsilon(x,y)=\exp\!\Big(-\frac{\|x-y\|^2}{4\varepsilon}\Big),$$3 yields the Fokker–Planck-compatible generator, whereas $$k_\varepsilon(x,y)=\exp\!\Big(-\frac{\|x-y\|^2}{4\varepsilon}\Big),$$4 preserves maximal density influence (Shnitzer et al., 2017, Trstanova et al., 2019, Beier et al., 28 Jan 2026).

This operator viewpoint explains why diffusion maps recover slow variables. For gradient flows with isotropic diffusion,

$$k_\varepsilon(x,y)=\exp\!\Big(-\frac{\|x-y\|^2}{4\varepsilon}\Big),$$5

the dominant diffusion coordinates empirically approximate eigenfunctions of the generator and capture metastable sets and transitions. From a Koopman perspective, if $$k_\varepsilon(x,y)=\exp\!\Big(-\frac{\|x-y\|^2}{4\varepsilon}\Big),$$6, then

$$k_\varepsilon(x,y)=\exp\!\Big(-\frac{\|x-y\|^2}{4\varepsilon}\Big),$$7

so the expected evolution of each coordinate is linear in time after spectral lifting (Shnitzer et al., 2017).

The same interpretation is frequently phrased in terms of the heat semigroup. On a manifold, $$k_\varepsilon(x,y)=\exp\!\Big(-\frac{\|x-y\|^2}{4\varepsilon}\Big),$$8 satisfies $$k_\varepsilon(x,y)=\exp\!\Big(-\frac{\|x-y\|^2}{4\varepsilon}\Big),$$9, and discrete diffusion operators are constructed to approximate this semigroup. This perspective motivates both diffusion distance and recent parameter-tuning criteria based directly on semigroup consistency rather than visualization heuristics (Shan et al., 2022).

3. Spectral basis versus charting

Diffusion maps are routinely used for nonlinear dimensionality reduction, but recent work argues that this description is incomplete. The eigenfunctions $q(x_i)=\sum_{j=1}^n k_\varepsilon(x_i,x_j).$0 form a basis of smooth modes adapted to intrinsic geometry; they are not, in general, the coordinates of a chart. Coordinate functions on a manifold are generally not eigenfunctions of the Laplace–Beltrami operator, so a meaningful parametrization typically requires a linear combination of several diffusion modes. This distinction is crucial in practice (Candanedo et al., 30 Mar 2026).

A Swiss-roll study with oracle sheet coordinates makes this point explicit. Isomap most efficiently recovers the low-dimensional chart at small embedding dimension, UMAP provides an intermediate tradeoff, and diffusion maps become accurate only after combining multiple diffusion modes. The correct chart lies in the span of diffusion coordinates, but standard diffusion maps do not by themselves identify the appropriate combination. On that example, two chart coordinates require tens to hundreds of diffusion modes before the oracle readout becomes near-exact (Candanedo et al., 30 Mar 2026).

Practice-oriented analysis further shows that the first components are not necessarily the most relevant ones. Scaling, normalization, quasi-discrete variables, and redundancy can all alter which diffusion coordinates dominate. In a Swiss-roll example, higher components may be polynomial functions of the first component when one direction dominates in scale; a later component can carry the missing complementary direction. To address this, the paper introduces Neural Reconstruction Error, which scores subsets and orderings of diffusion components by training an approximate inverse map and measuring reconstruction fidelity (Beier et al., 28 Jan 2026).

A common misconception is therefore that a spectral gap or the first $q(x_i)=\sum_{j=1}^n k_\varepsilon(x_i,x_j).$1 nontrivial eigenvectors automatically reveal a $q(x_i)=\sum_{j=1}^n k_\varepsilon(x_i,x_j).$2-dimensional chart. The available evidence instead supports a narrower claim: diffusion maps furnish a manifold-adapted spectral basis and a geometry-aware distance, while chart recovery remains an additional inverse problem (Candanedo et al., 30 Mar 2026, Beier et al., 28 Jan 2026).

4. Algorithms, parameter tuning, and scalability

A standard computational pipeline consists of constructing a kernel, normalizing it, building either a dense or sparse graph, computing leading eigenpairs, choosing a diffusion time $q(x_i)=\sum_{j=1}^n k_\varepsilon(x_i,x_j).$3, and truncating to an embedding dimension $q(x_i)=\sum_{j=1}^n k_\varepsilon(x_i,x_j).$4. Practical performance depends strongly on $q(x_i)=\sum_{j=1}^n k_\varepsilon(x_i,x_j).$5, $q(x_i)=\sum_{j=1}^n k_\varepsilon(x_i,x_j).$6, $q(x_i)=\sum_{j=1}^n k_\varepsilon(x_i,x_j).$7, the graph sparsity pattern, and the spectral truncation rule. Too small an $q(x_i)=\sum_{j=1}^n k_\varepsilon(x_i,x_j).$8 produces noisy or disconnected graphs; too large an $q(x_i)=\sum_{j=1}^n k_\varepsilon(x_i,x_j).$9 oversmooths the geometry. In sparse implementations, $$k_\varepsilon^{(\alpha)}(x_i,x_j)=\frac{k_\varepsilon(x_i,x_j)}{q(x_i)^\alpha q(x_j)^\alpha},$$0-nearest-neighbor graphs and sparse kernels reduce memory and computation, while spectral gaps and cross-validated downstream performance are common ways to choose $$k_\varepsilon^{(\alpha)}(x_i,x_j)=\frac{k_\varepsilon(x_i,x_j)}{q(x_i)^\alpha q(x_j)^\alpha},$$1 (Shnitzer et al., 2017, Beier et al., 28 Jan 2026).

The diffusion time $$k_\varepsilon^{(\alpha)}(x_i,x_j)=\frac{k_\varepsilon(x_i,x_j)}{q(x_i)^\alpha q(x_j)^\alpha},$$2 controls scale: larger $$k_\varepsilon^{(\alpha)}(x_i,x_j)=\frac{k_\varepsilon(x_i,x_j)}{q(x_i)^\alpha q(x_j)^\alpha},$$3 downweights higher-frequency modes through $$k_\varepsilon^{(\alpha)}(x_i,x_j)=\frac{k_\varepsilon(x_i,x_j)}{q(x_i)^\alpha q(x_j)^\alpha},$$4. A principled recent proposal selects $$k_\varepsilon^{(\alpha)}(x_i,x_j)=\frac{k_\varepsilon(x_i,x_j)}{q(x_i)^\alpha q(x_j)^\alpha},$$5 by enforcing the semigroup property. Using a symmetric positive semidefinite operator $$k_\varepsilon^{(\alpha)}(x_i,x_j)=\frac{k_\varepsilon(x_i,x_j)}{q(x_i)^\alpha q(x_j)^\alpha},$$6, the semigroup error

$$k_\varepsilon^{(\alpha)}(x_i,x_j)=\frac{k_\varepsilon(x_i,x_j)}{q(x_i)^\alpha q(x_j)^\alpha},$$7

is minimized over a useful short-to-moderate time range. The rationale is that the correct $$k_\varepsilon^{(\alpha)}(x_i,x_j)=\frac{k_\varepsilon(x_i,x_j)}{q(x_i)^\alpha q(x_j)^\alpha},$$8 is the scale at which the discrete operator behaves most like a manifold diffusion semigroup, while avoiding the large-$$k_\varepsilon^{(\alpha)}(x_i,x_j)=\frac{k_\varepsilon(x_i,x_j)}{q(x_i)^\alpha q(x_j)^\alpha},$$9 over-smoothed regime (Shan et al., 2022).

Out-of-sample extension is typically handled by Nyström interpolation. For a new point $$p(x_i,x_j)=\frac{k_\varepsilon^{(\alpha)}(x_i,x_j)}{\sum_{\ell=1}^n k_\varepsilon^{(\alpha)}(x_i,x_\ell)}.$$0,

$$p(x_i,x_j)=\frac{k_\varepsilon^{(\alpha)}(x_i,x_j)}{\sum_{\ell=1}^n k_\varepsilon^{(\alpha)}(x_i,x_\ell)}.$$1

which enables online embedding without recomputing the full eigendecomposition. This basic idea underlies several acceleration schemes. Nyström-accelerated diffusion maps achieved speedups of roughly two to four times when approximating dominant diffusion components, while Landmark Diffusion Maps reduce out-of-sample complexity from $$p(x_i,x_j)=\frac{k_\varepsilon^{(\alpha)}(x_i,x_j)}{\sum_{\ell=1}^n k_\varepsilon^{(\alpha)}(x_i,x_\ell)}.$$2 to $$p(x_i,x_j)=\frac{k_\varepsilon^{(\alpha)}(x_i,x_j)}{\sum_{\ell=1}^n k_\varepsilon^{(\alpha)}(x_i,x_\ell)}.$$3, $$p(x_i,x_j)=\frac{k_\varepsilon^{(\alpha)}(x_i,x_j)}{\sum_{\ell=1}^n k_\varepsilon^{(\alpha)}(x_i,x_\ell)}.$$4, and report up to 50-fold speedups in molecular examples with less than $$p(x_i,x_j)=\frac{k_\varepsilon^{(\alpha)}(x_i,x_j)}{\sum_{\ell=1}^n k_\varepsilon^{(\alpha)}(x_i,x_\ell)}.$$5 errors in manifold reconstruction fidelity (Shnitzer et al., 2017, Erichson et al., 2018, Long et al., 2017).

Compression can also occur at the level of data regions rather than landmark points. “Compressed Diffusion” computes diffusion geometry from a compressed process between data regions using an inverse-density Gaussian-correlation kernel, then interpolates pointwise embeddings from region-level coordinates. This suggests a broader view in which diffusion geometry can be approximated by operator compression, landmarking, or region aggregation, provided the coarse representation preserves local diffusion structure (Gigante et al., 2019).

5. Generalizations of the basic framework

Several extensions replace the classical Laplace–Beltrami or sampling-density-driven limit by more targeted operators. One line of work introduces Target Measure Diffusion Maps and Local Kernel Diffusion Maps, which replace the drift induced by the sampling density with the gradient of an arbitrary density of interest or approximate the forward and backward generators of general non-degenerate Itô diffusions with prescribed drift and diffusion coefficients. The local kernel

$$p(x_i,x_j)=\frac{k_\varepsilon^{(\alpha)}(x_i,x_j)}{\sum_{\ell=1}^n k_\varepsilon^{(\alpha)}(x_i,x_\ell)}.$$6

encodes general drift $$p(x_i,x_j)=\frac{k_\varepsilon^{(\alpha)}(x_i,x_j)}{\sum_{\ell=1}^n k_\varepsilon^{(\alpha)}(x_i,x_\ell)}.$$7 and diffusion tensor $$p(x_i,x_j)=\frac{k_\varepsilon^{(\alpha)}(x_i,x_j)}{\sum_{\ell=1}^n k_\varepsilon^{(\alpha)}(x_i,x_\ell)}.$$8, and the associated normalization removes bias from arbitrary sampling densities (Banisch et al., 2017).

Another line generalizes diffusion maps from manifolds to fibre bundles and tangent bundles. Horizontal Diffusion Maps model a dataset with pairwise structural correspondences as a fibre bundle equipped with a connection, and their asymptotics reveal horizontal and vertical differential operators associated with sub-Riemannian geometry. Hypoelliptic Diffusion Maps on tangent bundles approximate operators of the form $$p(x_i,x_j)=\frac{k_\varepsilon^{(\alpha)}(x_i,x_j)}{\sum_{\ell=1}^n k_\varepsilon^{(\alpha)}(x_i,x_\ell)}.$$9, interpolating between horizontal transport-dominated and fibre-averaged regimes through the bandwidth ratio $\varepsilon$0 (Gao, 2016, Gao, 2015).

Parameterized or changing datasets motivate still another extension. For a family of kernels $\varepsilon$1 indexed by parameters, “dynamic diffusion distance” compares $\varepsilon$2 and $\varepsilon$3 through the $\varepsilon$4 distance between their parameter-specific diffusion kernels. The framework defines change-of-basis operators $\varepsilon$5 that transport diffusion embeddings into a common space, and a global diffusion distance that measures operator-level change across the parameter space (Coifman et al., 2012).

Diffusion maps have also been coupled to state-space modeling. The Diffusion Maps Kalman Filter constructs “virtual state coordinates” from measurements, uses a modified Mahalanobis distance to approximate latent geodesic distances, and exploits the approximately linear drift of diffusion coordinates to build a linear Kalman filter in the learned coordinate system. The result is a non-parametric state-estimation framework for high-dimensional nonlinear stochastic systems with gradient-flow structure (Shnitzer et al., 2017).

On manifolds with boundary, a weak variational re-analysis shows that diffusion maps approximate Dirichlet energies and naturally recover Neumann eigenvalue problems. By estimating the distance-to-boundary function directly from the point cloud, the framework can impose Dirichlet, Neumann, and mixed boundary conditions for elliptic and parabolic PDEs without constructing a mesh or triangulation (Vaughn et al., 2019).

6. Applications across scientific domains

In molecular systems, diffusion maps approximate the generator of Langevin dynamics from simulation data and identify slowly evolving principal modes. Global analyses use leading eigenfunctions to detect metastable sets and compute committor functions; local analyses within metastable regions use quasi-stationary distributions to characterize intra-basin dynamics. This local–global perspective was demonstrated on alanine dipeptide and deca-alanine, and combined with biasing methods to accelerate sampling of the stationary Boltzmann–Gibbs distribution (Trstanova et al., 2019).

In stochastic state estimation, diffusion maps have been used to linearize nonlinear latent dynamics in a data-driven coordinate system. The Diffusion Maps Kalman Filter was evaluated on nonlinear object tracking with Gaussian and Poisson observations and on rodent hippocampal neural recordings. In the neural setting, applying the method directly yielded a representation of animal position from CA1/CA3 or CA1/MEC activity, and the resulting coordinates correlated strongly with two-dimensional position while outperforming raw diffusion maps and PCA in trajectory reconstruction (Shnitzer et al., 2017).

In wireless localization, diffusion coordinates derived separately from RSSI-based distances and from blueprint distances define a shared intrinsic space in which equipment and candidate locations can be matched. The method solves the resulting assignment problem with the Hungarian algorithm. In a realistic factory-floor layout with $\varepsilon$6 positions, it achieved perfect or near-perfect matching at $\varepsilon$7, and the same framework remained effective across grids, biaxial layouts, long-thin strips, and three-dimensional configurations (Ghafourian et al., 2019).

In neuroimaging, diffusion maps were used to cluster fMRI ICA spatial maps in a small-$\varepsilon$8, large-$\varepsilon$9 regime. In the reported experiment, $$P\phi_i=\lambda_i\phi_i,\qquad 1=\lambda_0\ge \lambda_1\ge \lambda_2\ge \cdots .$$0 maps with $$P\phi_i=\lambda_i\phi_i,\qquad 1=\lambda_0\ge \lambda_1\ge \lambda_2\ge \cdots .$$1 voxels were embedded by a Gaussian-kernel diffusion map before spectral clustering. The resulting two-cluster partition matched agglomerative hierarchical clustering and $$P\phi_i=\lambda_i\phi_i,\qquad 1=\lambda_0\ge \lambda_1\ge \lambda_2\ge \cdots .$$2-means, while producing more compact clusters when needed and offering an interpretable low-dimensional visualization (Sipola et al., 2013).

More recent applications emphasize operator learning rather than embedding alone. Diffusion map particle systems use spectral approximations of Langevin generators together with Laplacian-adjusted Wasserstein gradient descent for generative modeling, with no offline training and minimal tuning; the method can outperform competing approaches on data sets of moderate dimension (Li et al., 2023). “Learning functions through Diffusion Maps” uses diffusion geometry to extend real-valued functions from sampled manifolds to the ambient space, introduces a low-rank SVD strategy for distance compression, and reports that the resulting method outperforms classical feedforward neural networks and interpolation methods in sparse CT reconstruction in both accuracy and efficiency (Gomez, 3 Sep 2025).

Taken together, these results support a stable characterization. Diffusion maps are not merely a visualization heuristic or a fixed dimensionality-reduction recipe. They are a spectral framework for constructing geometry-aware coordinates, approximating diffusion generators, and organizing high-dimensional data through intrinsic connectivity. Their distinctive strength lies in the operator they approximate and in the family of tasks that become accessible once that operator is available: embedding, filtering, clustering, state estimation, PDE discretization, metastability analysis, and function extension.