Double Diffusion Maps: Sequential Manifold Learning

• Double Diffusion Maps is a two-stage manifold learning framework that first extracts intrinsic coordinates from ambient data and then refines them via latent space harmonization.
• It constructs a spectral basis that facilitates robust out-of-sample extension, lifting, and generative sampling, overcoming limitations of single-stage diffusion maps and POD.
• Applications span scientific computation, inverse regression, and generative modeling, achieving high accuracy in interpolating and reconstructing high-dimensional yet low-dimensional data.

Double Diffusion Maps is a manifold learning framework comprising two sequential applications of diffusion maps—one in the ambient space to uncover intrinsic coordinates, and a second in the discovered latent space—to construct a functionally rich basis for modeling, interpolation, lifting, and generative sampling on high-dimensional yet intrinsically low-dimensional data. This approach is designed to address challenges in reduced-order modeling, generative learning, forward and inverse regression, and multimodal/multiview correspondence, particularly when classical single-stage manifold learning methods or standard reduced-order models (e.g., POD) are insufficient for accurate out-of-sample extension and functional completeness.

1. Mathematical Foundations of Double Diffusion Maps

Double Diffusion Maps is rooted in the spectral theory of Markov kernels and graph Laplacians. The methodology comprises the following stages:

1. First Diffusion Map (Discovery Step):
   • Given $N$ samples $X=\{x_1, \ldots, x_N\}\subset\mathbb{R}^d$, construct a kernel $K(x_i, x_j) = \exp(-\|x_i-x_j\|^2/2\epsilon)$.
   • Apply density correction: $K̃ = P^{-\alpha} K P^{-\alpha}$ with $P_{ii} = \sum_j K_{ij}$ and exponent $\alpha\in[0,1]$ (commonly $\alpha=1$).
   • Row-normalize to a Markov matrix $A_{ij} = K̃(x_i,x_j) / \sum_\ell K̃(x_i,x_\ell)$.
   • Eigendecomposition: $A \phi_\ell = \lambda_\ell \phi_\ell$ with $|\lambda_0| ≥ |\lambda_1| ≥ \cdots$.
   • The diffusion map at time $t$ is $$\Psi_t(x_i) = (\lambda_1^t\phi_1(x_i), \ldots, \lambda_k^t\phi_k(x_i))$$.
   • Identify the $k$ non-harmonic eigenvectors (commonly via eigengaps or local linear regression residuals) to define the latent coordinates $$z_i = (\phi_1(x_i),\ldots,\phi_k(x_i)) \in \mathbb{R}^k$$.
2. Second Diffusion Map (Latent Space Harmonization):
   • Take $Z = \{z_i\} \subset \mathbb{R}^k$ from the first stage and define a new kernel $K^*(z_i,z_j)=\exp(-\|z_i-z_j\|^2/(2\epsilon_2))$.
   • (Optionally) Normalize and row-normalize to obtain a Markov matrix $A^*$.
   • Eigendecompose $A^* \psi_j = \mu_j \psi_j$. The $d$ leading modes $\psi_j$ are called latent harmonics; they form a data-driven basis for smooth functions on the latent manifold.

This two-stage construction yields a spectral basis that is well-suited for lifting (interpolation from latent to ambient space), restriction (projection of functions into latent coordinates), and enables principled out-of-sample extension via the Nyström or geometric harmonic formulae. All spectral computations depend critically on the choice of kernel bandwidths $\epsilon, \epsilon_2$, and the number of non-harmonic modes $k,d$.

2. Lifting and Out-of-Sample Extension via Latent Harmonics

A key innovation in Double Diffusion Maps is the use of the second-stage eigenbasis for global, out-of-sample lifting of latent representations to arbitrary functions defined on the original data set (including full ambient reconstructions):

• Nyström extension: For a new ambient point $x_\text{new}$, eigenfunctions are extended by computing appropriate kernel-weighted averages of their values in training data, normalized by the corresponding eigenvalue.
• Latent harmonic lifting: For any new latent coordinate $z_\text{new}\in\mathbb{R}^k$, one projects the ambient coordinate function $f(\cdot)$ onto the latent harmonic basis:

$c_j = \sum_{i=1}^N f(x_i) \psi_j(z_i)$

$$R(f)(z_\text{new}) = \sum_{j=0}^{d-1} c_j \Psi_j(z_\text{new}), \quad \Psi_j(z_\text{new}) = \mu_j^{-1} \sum_{i=1}^N K^*(z_\text{new}, z_i) \psi_j(z_i)$$

• The above can be applied componentwise to reconstruct $x^\text{new}$, or to reconstruct any other function of interest in the ambient domain.

The use of latent harmonics as a global interpolation basis ensures completeness and smoothness properties that are superior to one-stage (single DM) counterparts, as the latter are often contaminated by eigenvector multiplicities and functional harmonics.

3. Simulation and Generative Modeling Methodologies

Double Diffusion Maps enables multiple strategies for integrating reduced dynamical models or constructing generative samples on manifolds:

Methodology	Key Steps	Use Cases
Back-and-Forth (BF)	Alternate lifting to ambient, evaluating dynamics, and pushing derivative using chain rule in latent space.	ODE/PDE simulation with online lifting
Grid Tabulation (GT)	Precompute chain-ruled dynamics $\dot\varphi$ on a latent grid, interpolate on-the-fly as needed.	Fast integration; amortized cost
Tabulation with Latent Harmonics Interpolation (TaLHI)	Interpolate $\dot\varphi$ over entire manifold using latent harmonics; integrate $d\varphi/dt$ only in latent space.	Global surrogates; accuracy critical

In generative modeling (Giovanis et al., 5 Mar 2025, Giovanis et al., 2 Jun 2025), after learning the latent manifold, either a score-based diffusion model or an Itô SDE is used to sample from the learned latent density. The sampled latent points are then lifted to ambient space via the latent harmonics extension, enabling high-fidelity sampling consistent with the true (possibly unknown) manifold geometry.

4. Applications: Scientific Computation, Regression, and Generative Sampling

Several prominent applications of Double Diffusion Maps illustrate its versatility:

• Scientific computation in latent space: For both dissipative PDE systems (e.g., Chafee-Infante equation) and stiff ODEs (e.g., hydrogen combustion), reduced trajectory integration in latent space using BF, GT, or TaLHI yields errors $<1\%$ in latent coordinates and mean-square errors $O(10^{-6})$ upon lifting to full state (Evangelou et al., 2022).
• Parameter-to-state and inverse regression: In chemical vapor deposition (Cu–CVD) reactor modeling, forward (param$\rightarrow$state), inverse (state$\rightarrow$param), and partial-to-partial regression tasks are performed with mean errors below $0.5\%$ (Koronaki et al., 2023).
• Generative manifold sampling: On S-curve and Hermite polynomial function manifolds, generative pipelines adopting Double Diffusion Maps demonstrate recovery of both density and geometric modes with $L_1$ and $L_2$ errors below $1\%$ and $2\%$, respectively, even under moderate noise (Giovanis et al., 5 Mar 2025, Giovanis et al., 2 Jun 2025).
• Mitigating overfitting in small-data regimes: The two-stage construction enriches the available basis in low-sample/high-dimensional scenarios, preventing the loss of functional completeness when the number of samples approaches the latent dimension.

5. Algorithmic Workflow and Computational Considerations

Essential computational steps and considerations for Double Diffusion Maps include:

• Kernel bandwidth selection: Typically via the median-pairwise-distance heuristic; both $\epsilon$ and $\epsilon_{2}$ must be tuned separately for ambient and latent stages.
• Embedding dimension ($k$, $d$): Selected by spectral gap or using local linear regression residual analysis to isolate non-harmonic modes.
• Spectral truncation threshold ($\delta$): Mode retention is based on eigenvalue magnitude relative to the leading mode (e.g., $\delta = 10^{-3}$); ensures numerical stability in GH extension.
• Computational complexity: Each kernel/spectral construction is $O(N^3)$ in naive dense implementations, but $O(N \log N)$ or $O(NdL)$ with kernel/sparse approximations and low intrinsic manifold dimensions.
• Partial data embedding: To reconstruct or predict from partial observations, at least $2L+1$ generic samples are required by Whitney's embedding theorem for invertibility (Koronaki et al., 2023).
• Robustness and stability: The method is robust against noise due to the smoothing properties of kernel construction and the averaging in the Nyström/GH lifting. However, poorly chosen bandwidths or overly aggressive spectral truncation can result in numerical instabilities.

6. Comparisons with Alternative Methods

Double Diffusion Maps provides several advantages over conventional dimensionality reduction and reduced-order modeling techniques:

• Gappy-POD: Requires substantially more modes to match geometric fidelity on curved manifolds. Masking (partial observation) breaks orthogonality and creates ill-conditioning, whereas the nonlinear coordinates of Double Diffusion Maps offer insensitivity to sensor placement, provided the unique embedding criterion is met (Koronaki et al., 2023).
• Single-stage Diffusion Maps: Standard diffusion-map bases are often incomplete for function reconstruction due to harmonic contamination; Double Diffusion Maps' second stage restores completeness, crucial for generative tasks and smooth nonlinear interpolation (Giovanis et al., 2 Jun 2025).
• Multimodal/Multi-view Synchronization: Through joint approximate diagonalization, Double Diffusion Maps can yield shared low-dimensional bases across distinct data modalities, improving clustering and retrieval on real and synthetic datasets (Eynard et al., 2012).

7. Limitations, Hyperparameter Selection, and Practical Guidelines

• Limitations: The two-stage spectral procedure adds computational cost, particularly in high $N$ regimes and for high-dimensional embeddings. Hyperparameters ($\epsilon,\epsilon_{2},\delta$) require careful tuning; GH lifting may be sensitive to spectral truncation, and out-of-sample extension can be unreliable if the second kernel is not chosen appropriately (Giovanis et al., 2 Jun 2025).
• Hyperparameter guidelines:
  • $\epsilon, \epsilon_2$: median-of-distances or cross-validation heuristics.
  • $k, d$: select based on spectral gap and local linear regression ($r_k$ residual metric).
  • $\delta$: $10^{-3}$ recommended to capture $>99\%$ of kernel energy while suppressing noise amplification.
  • For generative pipelines, bandwidth $h$ (for KDE) is often set via Silverman’s rule.
• Implementation tips: Data should be normalized prior to kernel computation, k-nearest-neighbor sparsification is effective for scaling, and only as many harmonics as needed for function approximation accuracy should be retained.

In summary, Double Diffusion Maps delivers a robust, versatile, and mathematically grounded framework for nonlinear dimension reduction, regression, simulation, and generative modeling of high-dimensional data on manifolds, supporting both scientific computation and modern machine learning pipelines across a range of application domains.