Source-linked AI summary
Nonlocal Hamiltonian Dynamics on Sparse L\'evy Graphs: Spectral Analysis and Multimodal Sampling
Miaolei Zheng, Ting Gao, Jinqiao Duan
TL;DR
Multimodal samplers can struggle to balance probability mass across isolated modes. This paper develops sparse, density-level nonlocal Hamiltonian transport with direct long-range exchange, and reports improved mode coverage and occupancy balance on tested benchmarks compared with several baselines.
Problem
When modes are separated by low-density regions, local MCMC proposals can linger around one mode or miss modes with significant probability mass under finite computational budgets.
Source: §1 · p. 2Method
The method evolves density and edge momentum through damped nonlocal Hamiltonian dynamics on a sparse graph with local and L´evy-type long-range edges.
Source: §5 · p. 22Results
On tested multimodal benchmarks, the method achieves better mode coverage and more balanced occupancy than MCMC, MALA, and underdamped Langevin dynamics under the reported computational budgets.
Source: §5 · pp. 22–23Takeaways & Limitations
The sparse graph construction reduces per-step complexity to O(N(klocal + m)) while preserving the heavy-tailed interaction structure used for nonlocal mode exploration.
Source: §5 · p. 22The fixed spatial discretization is subject to the curse of dimensionality, and the L´evy kernel parameters and damping coefficient require problem-specific tuning.
Source: §5 · p. 23
Abstract
from arXiv · showhide
We develop a sparse graph method for transporting probability mass toward multimodal target distributions through damped nonlocal Hamiltonian dynamics. The formulation combines logarithmic-mean mobility with symmetric L\'evy-type interaction weights, coupling the evolving density to an edge momentum field. A graph constructed from nearest-neighbor connections and sampled long-range edges provides direct mass exchange between spatially separated regions. Once the graph is constructed, the density evolution is deterministic, and each update costs linear in the number of nodes and the long-range sampling budget. Linearization around the target distribution yields a damped oscillator governed by a weighted graph Laplacian. Its spectrum characterizes the interaction between nonlocal connectivity and inertia, with the L\'evy exponent alpha tuning the nonlocal connectivity: the spectral gap determines the optimal asymptotic damping, while the largest eigenvalue governs the time-step stability. Experiments on synthetic multimodal distributions demonstrate improved mode balance and more stable mode coverage relative to first-order and MCMC baselines. The resulting framework provides a sparse implementation of nonlocal inertial density transport for sampling problems with low-dimensional spatial structure.
1 Introduction
The paper proposes sparse, density-level nonlocal Hamiltonian transport to address mode exploration in multimodal targets. It reports improved mode balance against tested baselines and characterizes the role of graph spectra in relaxation.
1 Introduction
- Local MCMC proposals [4] and gradient-based samplers can struggle to move between modes separated by low-density regions; SVGD interactions [9] are limited by bandwidth.Graph-based acceleration [11] is also constrained when distant modes are weakly connected.Source: §1 · p. 2
- L´evy-type jumps permit direct communication across distance, motivating their combination with accelerated graph-based transport in this framework.Prior sampling approaches have combined L´evy jumps with Hamiltonian dynamics and heavy-tailed noise.Source: §1 · p. 2
- The proposed method improves mode coverage and balances mode occupancy on tested multimodal benchmarks relative to first-order, MCMC, MALA, and ULD baselines.The method evolves an explicit density on a fixed spatial grid and is evaluated in low dimensions.Source: §1 · p. 2
- Sparse graphs combine nearest-neighbor connections with sampled L´evy-type long-range edges, reducing per-step cost from O(N2) to O(N(klocal + m)).The stated cost applies once the graph is constructed.Source: §1 · p. 2
- The paper characterizes linearized relaxation through the weighted graph spectrum, deriving optimal damping from the spectral gap and distinguishing relaxation from discretization errors.Source: §1 · p. 2
2 Nonlocal Wasserstein Geometry and Hamiltonian Dynamics
The section develops damped Hamiltonian transport in nonlocal Wasserstein geometry, then expresses the density–momentum dynamics on a weighted graph. The graph formulation couples mass-conserving density evolution to damped edge dynamics.
2.1 Preliminaries
- The preliminaries define the target through KL divergence and introduce a symmetric jump kernel, logarithmic-mean edge density, and nonlocal continuity equation.Source: §2.1 · p. 3
2.2 Hamiltonian formulation and damped dynamics
- The Hamiltonian formulation combines kinetic and potential energy; linear damping is added because the undamped system conserves energy and does not relax to equilibrium.In the overdamped limit, the system reduces to first-order nonlocal Wasserstein gradient flow.Source: §2.2 · pp. 4–5
2.3 Graph discretization
- The graph discretization uses node masses and symmetric edge weights, with edge potential differences defining a coupled density–velocity system.Its antisymmetric edge flux drives mass-conserving density evolution, while the second equation governs edge dynamics under damping and nonlocal convection.Source: §2.3 · p. 5
3 Convergence analysis
The convergence analysis establishes local behavior around equilibrium, temporal accuracy and stability conditions, and probabilistic control of sparse-graph approximation error. It gives a high-probability finite-horizon error bound combining sparse forcing, temporal discretization, and continuous-time relaxation, with a nonlocal Wasserstein bound under an additional neighborhood condition.
3.1 Continuous-time mixing
- The energy–dissipation identity makes the total Hamiltonian non-increasing and motivates a discrete nonlocal Fisher information, but does not quantify the convergence rate.The analysis therefore turns to linearization near equilibrium to characterize convergence rates.Source: §3.1 · p. 6
3.1.1 Linearization around equilibrium
- Near equilibrium, the self-adjoint graph Laplacian’s spectral gap controls local convergence; optimal damping is γ*=2√λ1, with a second-order advantage when 0<λ1<1.Under Assumptions 3.1 and 3.2, solutions starting sufficiently near (π, 0) exist globally, remain in a local region, and converge exponentially.Source: §3.1.1 · pp. 7–8
- With zero initial momentum, the local exponential convergence result also gives a convergence rate in the nonlocal Wasserstein distance WNL.Its energy decay rate is λc=2r1(γ)>0.Source: §3.1.1 · p. 9
3.2 Temporal discretization error
- The semi-implicit Euler scheme has O(∆t) global ℓ2 error on finite intervals, while linear-mode stability must use the graph’s actual λmax; mesh scaling alone gives no graph-specific prefactor.The λmax stability condition applies to the unregularized linear update.Source: §3.2 · p. 10
3.3 Sparse graph approximation error
- Under the graph-structure and small-sampling assumptions, Proposition 3.1 bounds sparse-versus-full operator error with high probability; Figure 3.1 illustrates the predicted m^-1/2 estimator scaling.The illustrated full-graph errors also include spatial discretization and domain-truncation effects.Source: §3.3 · p. 12, Figure 3.1 · p. 12
- With probability at least 1 − δp, the sparse forcing error is bounded along the trajectory under uniform density and velocity bounds and an admissible Bernstein threshold.The stated forcing-error constant is independent of m and K under the uniform state bounds.Source: §3.3 · pp. 12–13
- The finite-horizon total error combines sparse-graph forcing, temporal discretization, and exponentially decaying relaxation, with constants independent of m, Δt, and k for fixed T.The bound holds with probability at least 1 − δp under the stated trajectory and bounded-region assumptions.Source: §3.3 · p. 13
- Under Theorem 3.2’s assumptions and a neighborhood condition, the total error bound transfers to the nonlocal Wasserstein metric; at γ = 2√λ1, r1(γ) = √λ1.WNL and λ1 refer to the same weighted graph.Source: §3.3 · p. 13
4 Numerical Experiments
The experiments test how fractional order and damping affect sampling, then compare the method with baseline samplers on synthetic and real-world targets. Results show broad mode coverage across synthetic benchmarks, while real-data performance varies with α and the target geometry.
4.1 Experimental Setup
- Experiments compare MCMC, MALA, and ULD on three synthetic distributions, MNIST, and alanine dipeptide, with shared grid and graph settings.The common setup uses a 50×50 grid, 10 local neighbors, and three long-range edges per node at α = 1.0.Source: §4 · p. 13, §4.1 · p. 14
4.2 Effect of Nonlocality
- Smaller α produces faster, longer-range transport while preserving complete mode coverage; larger α yields more localized transport and slower convergence in the banana experiment.The long-jump ratio falls from 0.145 at α = 0.1 to 0.016 at α = 1.5, and transport length decreases from 1.118 to 0.451.Source: §4.2 · pp. 14–15
4.3 First- and Second-Order Dynamics
- For the banana target, the fastest second-order decay occurs near γ∗ = 2√λ1, while excessive damping produces a strongly overdamped, first-order-like regime.The overdamped limit γ → ∞ reduces the second-order dynamics to the first-order flow.Source: §4.3 · p. 16
4.3.1 Spectral and Damping Analysis
- In this discretization, smaller α is associated with a larger spectral gap; first-order relaxation is earlier when λ1 > 1, while second-order relaxation is earlier when λ1 < 1.The critically damped second-order and first-order amplitude decay rates are √λ1 and λ1, respectively; the reported relationship is specific to this discretization and normalization.Source: §4.3.1 · pp. 17–18
4.4 Comparison with Other Samplers
- Across three synthetic targets, the method achieves better mode coverage than tested baselines, including all four mixture modes with balanced weights.On the single ring it covers the manifold, and on the banana it recovers curved tails; the four-Gaussian mixture shows the largest improvement.Source: Figure 4.6 · p. 18, Figure 4.7 · p. 18, Figure 4.8 · p. 18
- Within the equal 40-s budget, graph setup weakens initial wall-clock performance, but its cost is amortized and final KL and mode-weight-error trends are favorable among compared methods.The timing averages five seeds and includes graph construction for the proposed method and warm-up for particle baselines.Source: §4.4 · p. 19, Figure 4.9 · p. 19
4.4.1 Mode Exploration and Collapse Prevention
- The second-order flow maintains nearly balanced occupancy on the four-mode mixture and two-well banana, unlike baselines that show imbalance, collapse, or remain in one well.The first-order nonlocal flow visits both banana wells but with imbalanced occupancy; particle baselines remain in one well over the final 5,000 samples.Source: §4.4.1 · p. 20, Figure 4.10 · p. 20
4.5 Experiments on Real-World Data
- On MNIST’s 2D PCA latent space, smaller α gives broader mode coverage; on alanine, α = 0.5 captures both major states and has the lowest reported final KL among tested values.The real-data experiments are restricted to two-dimensional latent spaces because the grid discretization scales as N = O(h^-d), making it suited primarily to low-dimensional structure or meaningful reaction coordinates.Source: §4.5 · pp. 21–22, Figure 4.13 · p. 22
5 Conclusion
The paper presents a density-level nonlocal Wasserstein Hamiltonian flow on sparse graphs, combining local and Lévy-type long-range edges. It reports improved mode coverage and balance on tested benchmarks, with spectral behavior shaped by a trade-off between nonlocality and inertia.
5 Conclusion
- The nonlocal Wasserstein Hamiltonian flow evolves density directly and enables long-range mode communication without traversing low-density barriers.It combines a graph-discretized second-order system with Lévy jump kernels.Source: §5 · p. 22
- The experiments report improved mode coverage and more balanced mode occupancy on tested multimodal benchmarks.Source: §5 · p. 22
- Stronger nonlocality enlarges the spectral gap, speeding first-order flow but weakening the second-order advantage; for λ1 < 1, second-order flow is asymptotically faster.Source: §5 · p. 23
- The fixed spatial discretization faces the curse of dimensionality, while Lévy-kernel parameters and damping require problem-specific tuning.Future work includes high-dimensional extensions, adaptive graph construction, and optimal damping beyond the linearized regime.Source: §5 · p. 23
A Variational Derivation
The variational derivation computes first variations of the Hamiltonian with respect to momentum and density. Symmetry of the interaction kernel simplifies the momentum variation and the density-dependent mobility terms.
A Variational Derivation
- The Hamiltonian depends on Φ through its quadratic term, and the derivation takes its first variation against a test function.Source: A Variational Derivation · p. 23
- The density variation accounts for the logarithmic-mean mobility’s dependence on both density arguments.The displayed variation separates contributions at x and y.Source: A Variational Derivation · p. 24
- Symmetry of J and the identity relating the two partial derivatives of θ simplify the second term in the density variation.Source: A Variational Derivation · p. 24
B Proofs of theoretical results
The appendix gathers technical proofs for the paper’s convergence and approximation results. It covers linearization, local convergence, temporal discretization, and sparse graph approximation.
B Proofs of theoretical results
- The appendix’s proofs support the convergence analysis, including linearization, local convergence, temporal discretization error, and sparse graph approximation error.Source: B Proofs of theoretical results · p. 24
B.1 Local convergence and spectral analysis
The analysis linearizes the dynamics around equilibrium and reduces density perturbations to damped modes of a positive semidefinite graph operator. Under a spectral gap and sufficiently small initial perturbations, it establishes local exponential convergence and identifies optimal damping.
B.1 Local convergence and spectral analysis
- Linearization around equilibrium reduces the density perturbation to a damped second-order equation with stiffness operator A = BC.The expansion uses log(1 + φi) = φi + O(φ2), and the nonlinear remainder has zero value and derivative at equilibrium.Source: §B.1 · pp. 25–26
- The operator A is self-adjoint and positive semidefinite, linking its spectrum to the local dissipation rate of relative entropy.Source: §B.1 · p. 25
- Each eigenmode follows a damped oscillator, and the slowest mode λ1 determines the global linearized decay rate.Source: §B.1 · p. 26
- The linearized decay rate is maximized at γ*=2√λ1, where the rate equals √λ1.Source: §B.1 · p. 26
- For sufficiently small initial perturbations under the stated assumptions, the nonlinear solution exists globally and remains in the local neighborhood.The bootstrap estimate yields decay at rate r1(γ), with a 1 + t factor at critical damping.Source: §B.1 · p. 27
- At critical damping, the mixing-time bound scales as √λ1 Tmix ≤ log(1/ε) + O(log log(1/ε)).Source: §B.1 · p. 27
- Near equilibrium, the nonlocal Wasserstein metric is locally equivalent to a normalized mass-error norm, supporting convergence estimates in that metric.Source: §B.1 · pp. 27–28
B.2 Temporal discretization error (Lemma 3.3)
The temporal error analysis assumes bounded density and velocity and local Lipschitz control near equilibrium. Under these conditions, the time-discrete approximation has a first-order error over the stated time horizon.
B.2 Temporal discretization error (Lemma 3.3)
- The local error analysis assumes Lipschitz bounds for nonlinear terms near equilibrium (π, 0).The density and velocity are bounded along the trajectory, supporting the stated nonlinear-term estimates.Source: §B.2 · pp. 28–29
- The error recursion combines local truncation error with coupled density and velocity estimates, absorbing a velocity term when γ_k Δt ≤ 1/2.The proof then applies discrete Gronwall to obtain the global bound.Source: §B.2 · pp. 29–30
- The global error is bounded by C_temp Δt e^(C t_k), with C_temp independent of Δt and k.The result follows by unrolling the recursion from zero initial error and applying discrete Gronwall.Source: §B.2 · p. 30
B.3 Sparse graph error
The sparse-graph analysis establishes unbiased, symmetric operator weights and bounds the resulting sampling error. It then controls the accumulated sparse-graph forcing error along the full-graph trajectory and combines it with temporal and spatial errors.
B.3 Sparse graph error
- The graph construction combines deterministic KNN edges with sampled long-range edges, excluding duplicate and already-connected selections.When k_local + m ≪ N, the interaction evaluation avoids the O(N^2) cost of a fully connected graph.Source: §B.3 · p. 30
- The sparse operator is unbiased and symmetric, and its nodewise approximation error has mean zero.The variance bound is O(1/m) under the stated sampling assumptions.Source: §B.3 · pp. 30–31
- The uniform sparse-operator error bound scales with the sampling budget m and the log of the node count and confidence level.The concentration argument chooses a threshold proportional to the square root of log(4N/δ)/(ĉ_1 m).Source: §B.3 · p. 32
- The sparse-graph forcing error is controlled across all nodes and time levels by concentration bounds and a union bound, yielding the stated trajectory estimate.The full trajectory is independent of the sparse-graph draws; the bound accounts for N nodes and K + 1 time levels.Source: §B.3 · p. 33
- The full error bound combines sparse-graph forcing, temporal discretization, and the spatial error bound.The sparse-graph contribution is controlled using a Lipschitz constant independent of m and the sampling realization.Source: §B.3 · p. 33
C Numerical Implementation Details
The numerical experiments use a two-dimensional computational grid and a sparse graph with KNN and sampled long-range edges. The section also specifies benchmark timing protocols and the measured components of algorithmic time.
C Numerical Implementation Details
- Long-range edges are sampled with probabilities proportional to r^-(2+α), while local KNN edges retain their base weights.Repeated selections are merged; the stated unbiasedness and sampling-error estimates use the exact inclusion weights from Section 3.3.Source: C Numerical Implementation Details · p. 34
- The equal-wall-clock comparison uses a 40-second budget per method and seed on one CPU, with five independent seeds.The proposed method’s measured time includes target setup, graph construction, graph-variable initialization, and density evolution.Source: C Numerical Implementation Details · p. 34
- For the proposed method, measured algorithmic time includes target setup, sparse-graph construction, graph-variable initialization, and density evolution; experiments use a 50 × 50 grid, α = 1, KNN query size 10, m = 3, and Δt = 5 × 10^-4.The damping is γ* = 2√λ1, with λ1 computed separately for each sparse-graph realization.Source: C Numerical Implementation Details · p. 34
- Particle baselines use a 4-second warm-up and 36-second production phase, with adaptation restricted to warm-up.Reported algorithmic times exclude KDE construction, KL evaluation, mode analysis, plotting, and file I/O.Source: C Numerical Implementation Details · p. 34