Hybrid stochastic dynamics, adaptive diffusion, Langevin dynamics, sampling algorithm, Fokker–Planck equation, Gibbs distribution

1 Introduction↩︎

Sampling from Gibbs distributions is a fundamental task in scientific and statistical computing. Let \(\Omega\subseteq\bbR^d\) (\(d\ge 1\)) be the given state space. Given a potential energy landscape \(F:\Omega\to \bbR\) and a parameter \(\varepsilon>0\), the objective is to sample from the Gibbs distribution \[\label{EQ:pi} \pi(\bx)=\frac{1}{Z} e^{-F(\bx)/\eps}, \qquad Z:=\int_\Omega e^{-F(\bx)/\eps}\,\,\mathrm{d}\bx.\tag{1}\] A standard approach is to simulate the overdamped Langevin dynamics \[\label{EQ:Langevin} dX_t=-\nabla F(X_t)\,dt+\sqrt{2\eps}\,dW_t,\tag{2}\] where \(W_t\) is a standard \(d\)-dimensional Brownian motion. The invariant measure of 2 is precisely \(\pi\). This dynamics forms the basis of many modern sampling algorithms; see, for instance, [1]–[8] and references therein.

Despite its simplicity and theoretical appeal, overdamped Langevin dynamics often suffers from severe metastability in multi-well energy landscapes. That is, when the potential \(F\) contains multiple local minima separated by high barriers, the process may remain trapped in one metastable basin for exponentially long times before transitioning to another. This phenomenon leads to slow mixing and poor sampling efficiency, especially in low-temperature (i.e., small \(\eps\)) regimes. To address this difficulty, various accelerated sampling strategies have been proposed [9]–[11], including underdamped dynamics [12]–[20], preconditioned Langevin dynamics [21]–[25], adaptive diffusion methods [26], [27], tempering, replica exchange and annealing-based methods [28]–[34], and many more [35]–[42]. A common principle underlying many successful approaches is to modify the local geometry of the dynamics in order to facilitate barrier crossing while preserving the target Gibbs law 1 .

In this work, we propose a framework for hybrid dynamics sampling, in which two distinct diffusion dynamics are run in different regions of the state space. More precisely, the state space is partitioned into two subregions, for example, according to a level set of the potential, and different stochastic dynamics are prescribed in the two regions. The coefficients are chosen so that the resulting hybrid process still preserves the original Gibbs distribution globally. The motivation is to retain the favorable properties of overdamped Langevin dynamics in regions where the potential already provides efficient confinement, while replacing the dynamics in metastable regions by other dynamics whose effective geometry is better suited for rapid exploration.

Our primary example is a hybridization between the standard overdamped Langevin dynamics and the derivative-free adaptive diffusion introduced in [26], [43], although the framework developed here is more general. The earlier adaptive-variance sampler showed that, on bounded domains, a carefully chosen state-dependent diffusion coefficient can flatten the effective geometry inside nonconvex potential wells and thereby accelerate metastable transitions. However, using this derivative-free dynamics globally is not always appropriate: on unbounded domains, removing the Langevin drift everywhere may destroy the confining mechanism at infinity, and the resulting process need not converge to the target Gibbs distribution. The hybrid construction is also motivated by a simple numerical observation (see Section 5.1). If the derivative-free diffusion coefficient is used on the whole space, then the coefficient grows exponentially in high-energy regions, and the resulting numerical scheme can lose confinement immediately.

The hybrid dynamics proposed here localizes this acceleration mechanism. Inside a prescribed interior region, the drift is removed, and the diffusion coefficient is modified so that the Gibbs measure remains invariant; in particular, the weighted conductance density \(a(\bx)\pi(\bx)\) becomes constant there. Outside this region, the dynamics revert to overdamped Langevin dynamics, preserving confinement and convergence.

The main part of this work provides a rigorous convergence theory for the hybrid dynamics. Because the coefficients are discontinuous across the switching interface, the corresponding Fokker–Planck equation naturally takes the form of a transmission problem. We introduce a smooth regularization of the coefficients, establish uniform entropy dissipation estimates independent of the regularization parameter, and prove compactness and convergence of weak solutions to the limiting hybrid Fokker–Planck system. This yields exponential convergence of the hybrid dynamics toward the Gibbs measure in relative entropy under a logarithmic Sobolev inequality assumption.

We also analyze the metastability properties of the hybrid dynamics in a radial double-well landscape, as an example of situations where the hybrid dynamics achieve the better mixing speed of the two separate dynamics. Using the radial Poisson problem and Laplace asymptotics, we compare the mean transition times between metastable wells for the hybrid dynamics and for overdamped Langevin dynamics. Our analysis shows that, while the overdamped Langevin transition time is governed by the classical saddle barrier height, the hybrid dynamics replaces this barrier by the switching level of the interface. Consequently, if the switching interface is chosen below the saddle level, precisely when \(F_r-F_g<H-F_u\) (with \(F_g\), \(F_u\), \(H\), and \(F_r\) the energies of the lower well, upper well, saddle, and switching interface, respectively), the hybrid dynamics achieves an exponentially faster transition rate than the overdamped Langevin dynamics.

The remainder of the paper is organized as follows. We introduce the hybrid dynamics and its regularization in Section 2. We then analyze the regularized dynamics uniformly in \(\delta\) and then pass to the sharp-interface limit in Section 3. The key estimate is the entropy dissipation identity, which gives compactness and transfers the exponential convergence rate to the limiting hybrid dynamics. In Section 4, we compare the metastable transition behavior of the hybrid and overdamped Langevin dynamics in a radial double-well potential to show that the hybrid dynamics accelerates metastable transition in this setting. Numerical simulations are presented in Section 5 to verify that the hybrid dynamics indeed sample the target distribution. Concluding remarks are offered in Section 6.

2 Hybrid dynamics and diffuse-interface regularization↩︎

We introduce the sharp-interface hybrid sampler and then define a smooth diffuse-interface approximation. The regularization is chosen so that the Gibbs density remains invariant for every \(\delta>0\), while the coefficients remain uniformly controlled as \(\delta\to0\).

2.1 Sharp-interface hybrid dynamics↩︎

In our hybrid approach, we divide the state space \(\Omega\) into two regions. Two different dynamics are run in the regions. We simultaneously consider two cases:

• \({\rm (Case~I)}\)

  : The state space \(\Omega\) is a bounded domain, with boundary \(\partial\Omega\in \cC^2\).

• \({\rm (Case~II)}\)

  : The state space is \(\Omega=\bbR^d\).

We assume that the potential function \(F\) has the following properties.

• \(F\) is sufficiently smooth, has a bounded gradient, and is bounded from below. More precisely: \[\label{EQ:F32Assumptions} F\in\cC^2(\Omega),\qquad \nabla F\in L^\infty(\Omega), \qquad and\qquad \inf_\Omega F >-\infty\,.\tag{3}\]

Note that the property \(\nabla F\in L^\infty(\Omega)\) is necessary as it does not directly follow from the first property in \({\rm (Case~II)}\). We assume throughout that the normalization constant \[Z=\int_\Omega e^{-F(\bx)/\eps}\,\,\mathrm{d}\bx<\infty.\]

In this work, we use two related ways of partitioning the domain. The main analysis in [SEC:Hybrid-Regularized,SEC:all_properties] treats the case where the switching interface is an entire level set of the potential. For a given \(F_0\in\mathbb{R}\), define \[\Gamma:=\{\bx\in\Omega:F(\bx)=F_0\}.\] In the full-level-set setting, we assume that \(\Gamma\) is sufficiently smooth, say \(\Gamma\in\mathcal{C}^1\), and separates \(\Omega\) into two regions \[\Omega_- := \{\bx\in\Omega:F(\bx)<F_0\}, \qquad \Omega_+ := \{\bx\in\Omega:F(\bx)>F_0\}.\] The whole level set \(\Gamma\) is then used as the interface between the two dynamics.

The exit-time analysis in 4 treats a radially symmetric, possibly nonconvex case, where the relevant level set may have several connected components, and the switching interface is chosen to be the outermost one, namely the component adjacent to the exterior region. Until 4, we work exclusively in the full-level-set setting.

Our goal is to be able to run two different sampling schemes in \(\Omega_-\) and \(\Omega_+\) while preserving the target distribution in the whole space \(\Omega\). More precisely, we run the overdamped Langevin dynamics 2 in \(\Omega_+\), and a derivative-free adaptive diffusion dynamics in \(\Omega_-\).

For a given \(\eps\), we define the functions \(a:\Omega\to \bbR\) and \(b:\Omega\to \bbR\) as \[\label{EQ:Coefficients32True} a(\bx)= \begin{cases} a_-:=\eps e^{(F(\bx)-F_0)/\eps} & \bx\in \Omega_-\\ a_+:=\eps, & \bx\in \Omega\setminus \Omega_- \end{cases}\,, \quad b(\bx)= \begin{cases} b_-:=0, & \bx\in \Omega_-\\ b_+:=-\nabla F(\bx), & \bx\in \Omega\setminus\Omega_- \end{cases}\,.\tag{4}\] Note that we scale the diffusion coefficient inside \(\Omega_-\) in a way such that \(a(\bx)\) is continuous across \(\Gamma\). The continuity of \(a\) is an essential property to enforce; see Remark [REM:Continuity32of32a].

The hybrid dynamics we are interested in are then \[\label{EQ:SDE32HB} \,\mathrm{d}X_t^\pm=b_\pm(X_t^\pm)\,\mathrm{d}t+\sqrt{2a_\pm(X_t^\pm)}\,\,\mathrm{d}W_t,\quad in\;\Omega_\pm\,.\tag{5}\] In \({\rm (Case~I)}\), that is, when \(\Omega\) is bounded, the process in \(\Omega_+\) is reflected at \(\partial\Omega\) to enforce the no-flux condition at \(\partial\Omega\).

The sampling dynamics in the inner region \(\Omega_-\) is purely diffusive in nature, in the sense that the drift term is turned off (since \(b_-=0\)). To keep the stationary distribution Gibbsian, the diffusion coefficient is chosen to have the particular form of \(a_-\). This leads to the fact that the weighted conductance density \(a_-(\bx)\pi(\bx)\) is constant: \(a_-(\bx)\pi(\bx) = \eps Z^{-1}e^{-F_0/\eps}\). Thus, although the invariant measure remains Gibbsian, the effective conductivity governing the associated Dirichlet form becomes spatially uniform in \(\Omega_-\); see more detailed discussions in 4 and  [26]. Therefore, this dynamics removes the exponential weighting induced by the potential inside the \(\Omega_-\) and replaces it with a flat diffusion geometry. Hence, for a potential whose landscape below the level set \(\Gamma\) is complex, the hybrid dynamics can potentially achieve better overall performance than running the overdamped Langevin dynamics in the whole state space \(\Omega\); see numerical examples in Section 5.

The hybrid process \(X_t\) is a Markov process on \(\Omega\). Its infinitesimal generator is \[\label{EQ:Generator} \cL f(\bx) =a(\bx) \Delta f + b(\bx) \cdot \nabla f(\bx)\,.\tag{6}\] The domain of \(\cL\), that is, which functions \(f\) count as admissible test functions, requires that \(f\) and \(a \bn_{\Gamma_-}\cdot\nabla f\) be continuous across \(\Gamma\). Therefore, the adjoint Fokker–Planck operator \(\cL^*\) acts on densities \(\rho\) with the dual transmission conditions. That is, the Fokker–Planck system for the hybrid dynamics is \[\label{EQ:FP32HB} \begin{array}{rcll} \partial_t\rho_\pm &=& -\nabla\cdot(b_\pm \rho_\pm)+\Delta(a_\pm \rho_\pm), & \text{in }\Omega_\pm\\[1ex] (a_-\rho_-)|_{\Gamma_-} &=& (a_+\rho_+)|_{\Gamma_+}, & \text{on}\;\Gamma\\ \bn_\Gamma\cdot J_-(\rho_-)|_{\Gamma_-} &=& \bn_\Gamma\cdot J_+(\rho_+)|_{\Gamma_+}, &\text{on}\;\Gamma \end{array}\tag{7}\] where the flux \[\label{EQ:Flux} J_\pm(\rho_\pm):=b_\pm \rho_\pm-\nabla(a_\pm\rho_\pm)\,,\tag{8}\] and \(\bn_\Gamma\) is the unit normal to \(\Gamma\) (pointing, say, from \(\Omega_-\) to \(\Omega_+\)). Since the coefficients \((a,b)\) are discontinuous across \(\Gamma\), the limiting Fokker–Planck equation takes the form of a transmission problem. The matching conditions in 7 express conservation of the weighted conductance density and the probability flux across the switching interface. These are the natural transmission conditions for a conservative diffusion process in a heterogeneous environment. With the continuity of \(a(\bx)\) across the interface \(\Gamma\), the first interface condition implies the continuity of the density across the interface: \(\rho|_{\Gamma_-}=\rho|_{\Gamma_+}\).

When \(\Omega\) is bounded, the no-flux boundary condition on the distribution level is described as \[\label{EQ:FP32BC} \bn_{\partial\Omega}\cdot J_+(\rho_+)|_{\partial\Omega}=0\,.\tag{9}\]

We intentionally choose the interface to switch the dynamics as a level set of the potential rather than an arbitrary hypersurface. While this choice appears to be quite restrictive, it has some advantages. First, it ensures that the matching condition for continuity of \(a\) can be enforced by a single scalar constant \(\gamma:=\eps e^{-F_0/\eps}\), since \(F\) is constant on \(\Gamma\). Second, from the sampling perspective, the value of \(F\) provides a natural energetic criterion for deciding where to modify the dynamics: the hybridization is activated in regions where the potential is below a prescribed threshold and metastability is expected to dominate.

2.2 Diffuse-interface regularization↩︎

The hybrid sampling dynamics we proposed has the drift coefficient that is discontinuous across the switching interface \(\Gamma\), and the associated Fokker–Planck equation 7 therefore takes the form of a transmission problem. This sharp-interface system is not very friendly for analyzing the performance of the system.

Rather than working directly with the sharp-interface system, we approximate the interface by a thin transition layer of width \(\cO(\delta)\), for some small \(\delta>0\), across which the coefficients vary smoothly. The regularized dynamics can therefore be viewed as a diffuse-interface approximation of the hybrid process. As \(\delta\to 0\), the transition layer collapses and the sharp-interface hybrid dynamics is formally recovered. To be precise, let \(\mathrm{dist}(\bx,\Gamma)\) denote signed distance (negative inside \(\Omega_-\)). For a given \(\delta>0\), we denote \[\chi_\delta(\bx)=\chi(\mathrm{dist}(\bx, \Gamma)/\delta)\] with a fixed smooth profile \(\chi:\bbR\to[0,1]\) satisfying \(\chi\equiv 1\) on \((-\infty,-1]\) and \(\chi\equiv 0\) on \([1,\infty)\). See 1 for an illustration of a candidate \(\chi\).

We then introduce the smoothed coefficients as \[\label{EQ:Coefficients32Regularized} \begin{array}{rcl} a_\delta(\bx) &=& \chi_\delta(\bx)\,\eps e^{(F(\bx)-F_0)/\eps}+(1-\chi_\delta(\bx))\,\eps\,,\\[1ex] b_\delta(\bx) &=& \nabla a_\delta(\bx)-\dfrac{a_\delta(\bx)}{\eps}\,\nabla F(\bx)\,. \end{array}\tag{10}\] Note that instead of directly regularizing both \(a\) and \(b\), our strategy here is to regularize the diffusion coefficient to \(a_\delta\) and then select \(b_\delta\) through the Gibbs-preserving relation. This is of key importance as it allows the dynamics to generate the same Gibbs distribution in the whole state space \(\Omega\) for every \(\delta>0\); see more discussions at the end of this section.

Formally, as \(\delta\to0\), \[a_\delta(\bx)\to a(\bx),\qquad b_\delta(\bx)\to b(\bx)\] for every \(\bx\in\Omega\setminus\Gamma\), and hence almost everywhere in \(\Omega\).

2.3 Basic properties of the regularized dynamics↩︎

We first record the basic bounds on the regularized coefficients. These bounds are uniform in the regularization parameter and are the main reason the diffuse-interface approximation can be used to pass to the sharp-interface limit.

The regularized dynamics is therefore \[\label{EQ:SDE32Regularized} \,\mathrm{d}X_t=b_\delta(X_t)\,\mathrm{d}t+\sqrt{2a_\delta(X_t)}\,\,\mathrm{d}W_t,\tag{11}\] and the associated Fokker–Planck equation is \[\label{EQ:FP32Regularized} \begin{array}{rcll} \partial_t\rho_\delta &=& -\nabla\cdot(b_\delta\rho_\delta)+\Delta(a_\delta\rho_\delta), & \text{in }\Omega. \end{array}\tag{12}\] We also define the regularized flux \[J_\delta(\rho_\delta):=b_\delta\rho_\delta-\nabla(a_\delta\rho_\delta)\,.\] The regularized coefficients \((a_\delta, b_\delta)\) are smooth and uniformly elliptic. Hence, the regularized Fokker–Planck equation becomes a classical uniformly parabolic drift-diffusion equation, allowing the use of standard entropy methods, compactness arguments, and parabolic regularity theory to analyze the system. We can then pass to the limit \(\delta\to 0\) to understand the original sharp-interface transmission problem 7 .

When \(\Omega\) is bounded, that is, in \({\rm (Case~I)}\), the no-flux boundary condition 9 is imposed for the regularized system: \[\label{EQ:FP32Regularized32BC} J_\delta(\rho_\delta)\cdot \bn_{|\partial\Omega} = 0 \qquad \text{on }\partial\Omega.\tag{13}\]

A key observation is that the Gibbs density \(\pi\) defined in 1 is a stationary solution for the regularized system, independent of \(\delta\), since: \[J_\delta(\pi)=b_\delta\pi-\nabla(a_\delta\pi)\equiv 0,\quad \forall \delta>0.\] This happens because we choose a special way to regularize \(a\) and adjust \(b\) accordingly in 10 . A naive independent regularization of \(a\) and \(b\) will change the Gibbs distribution inside the thin layer around the interface \(\Gamma\).

From a computational perspective, the regularized dynamics may also be viewed as a practical implementation of the hybrid sampler. Indeed, any numerical discretization of a sharp switching interface introduces an effective smoothing at the mesh scale, so the parameter \(\delta\) can also be interpreted as an interface-resolution parameter.

The goal of the next section is to establish that the hybrid dynamics preserves the target Gibbs distribution and converges to it at an exponential rate, uniformly with respect to the regularization parameter \(\delta\).

The main difficulty is that the limiting hybrid system 7 is a transmission problem with discontinuous coefficients across the interface \(\Gamma\), for which standard parabolic theory does not directly apply. Our strategy is therefore to analyze the regularized system 12 and then pass to the limit \(\delta\to 0\).

For each \(\delta>0\), the regularized coefficients \((a_\delta,b_\delta)\) are smooth and uniformly elliptic (Lemma 1), so the regularized Fokker–Planck equation 12 is a classical uniformly parabolic drift-diffusion equation. It admits smooth classical solutions with smooth initial data [44]. In the rest of this work, we assume that the initial data we consider, \(\rho_{0}\), satisfies \[\label{EQ:IC32ASS} \begin{array}{lllll} \hyperlink{case1}{{\rm (Case~I)}}: & \rho_0\in C^\infty(\overline{\Omega}), & \rho_0>0, & \int_\Omega\rho_0\,\,\mathrm{d}\bx=1, & \mathrm{KL}(\rho_0\,\|\,\pi)<+\infty\\[1ex] \hyperlink{case2}{{\rm (Case~II)}}: & \rho_0\in\cS(\bbR^d),& \rho_0>0,& \int_{\bbR^d}\rho_0\,\,\mathrm{d}\bx=1,& \mathrm{KL}(\rho_0\,\|\,\pi)<+\infty \end{array}\tag{14}\] where \(\cS(\bbR^d)\) is the standard Schwartz space of rapidly decreasing functions on \(\bbR^d\).

We recall this classical result here for convenience.

3 Entropy estimates, sharp-interface limit, and exponential convergence↩︎

In this section, we analyze the regularized dynamics uniformly in \(\delta\) and then pass to the sharp-interface limit. The key estimate is the entropy dissipation identity, which gives compactness and transfers the exponential convergence rate to the limiting hybrid dynamics.

3.1 Uniform entropy dissipation↩︎

For convenience, we first rewrite 12 into the standard divergence (Onsager) form \[\label{EQ:FP32Onsager32Form} \partial_t\rho_\delta = \nabla\cdot\left(a_\delta(\bx)\,\rho_\delta\,\nabla\log\frac{\rho_\delta}{\pi}\right) \qquad \text{in }\Omega,\tag{15}\] using the identity \(\nabla(a_\delta\rho) - b_\delta\rho = a_\delta\rho\,\nabla\log\!\left(\frac{\rho}{\pi}\right)\). The boundary condition 13 , required when \(\Omega\) is bounded, that is, in \({\rm (Case~I)}\), then takes the form \[\label{EQ:FP32Onsager32Form32BC} \left(a_\delta\rho_\delta\,\nabla\log(\rho_\delta/\pi)\right)\cdot \bn_{|\partial\Omega}=0, \qquad on\;\partial\Omega\,.\tag{16}\] We define the standard relative entropy energy \[\label{EQ:KL} \cE(\rho):=\mathrm{KL}(\rho\|\pi)=\int_\Omega \rho\log\frac{\rho}{\pi}\,\,\mathrm{d}\bx.\tag{17}\] Then the PDE 15 can be interpreted as a gradient flow of the entropy energy \(\cE(\rho_\delta)\) in a Wasserstein-type metric \(\mathsf W_\delta\) with position-dependent mobility \(m_\delta(\bx,\rho):=a_\delta \rho\) [43]. That is, \[\label{eq:GF95metric} \partial_t\rho_\delta = \nabla_{\mathsf W_\delta}\cE(\rho_\delta) \quad\Longleftrightarrow\quad \partial_t\rho_\delta=\nabla\cdot \!\left(a_\delta\rho_\delta\nabla\frac{\delta \cE}{\delta \rho_\delta}\right) = \nabla\cdot \!\left(a_\delta\rho_\delta\,\nabla \log\frac{\rho_\delta}{\pi}\right)\,.\tag{18}\] Here \(\mathsf W_\delta\) is the metric whose tangent vectors \(\sigma=-\nabla\cdot(a_\delta\rho\nabla\phi)\) carry the weighted kinetic energy \(\int_\Omega a_\delta\rho|\nabla\phi|^2\,\,\mathrm{d}\bx\); since \(\delta\cE/\delta\rho=\log(\rho/\pi)+1\), the steepest descent of \(\cE\) is precisely the Onsager form 15 . When \(a_\delta\) is constant, \(\mathsf W_\delta\) reduces, up to a time rescaling, to the usual quadratic Wasserstein metric [47], [48].

In the rest of the work, we assume that the potential \(F\) is such that the stationary distribution \(\pi\) in 1 satisfies the following properties:

• (Log-Sobolev Inequality (LSI)) For all sufficiently smooth \(\rho\) with \(\dint_\Omega\rho\, \,\mathrm{d}\bx=1\) and \(\rho\ll \pi\), we have \[\label{EQ:LSI} \mathrm{KL}(\rho\|\pi)\le \frac{1}{2\lambda_{\mathrm{LSI}}}\int_\Omega \rho\left | \nabla\log\frac{\rho}{\pi}\right |^2 d \bx,\tag{19}\] for some LSI constant \(\lambda_{\mathrm{LSI}}>0\).

• (Exponential Moment Condition) In \({\rm (Case~II)}\), we have that \[\label{EQ:Moment32Cond} \dint_{\Omega} e^{\alpha|\bx|^2}\pi(\bx)\,\mathrm{d}\bx<\infty,\quad for some\;\alpha>0\,.\tag{20}\]

Log-Sobolev inequalities for \(\pi\), such as 19 , are available under various scenarios. On bounded domains, LSI for \(\pi\) is often available under mild regularity, though the LSI constants may depend on \(\Omega\), \(F\), and \(\eps\) [49]–[51].

The exponential moment condition 20 ensures that when \(\mathrm{KL}(\rho\|\pi)\) is finite, \(\rho\) has finite second moment. This is classical [48], [52]. We reproduce it in the following lemma.

We define the entropy dissipation functional (or the weighted Fisher information) \[\mathcal{I}_{a_\delta}(\rho_\delta\|\pi):=\int_\Omega a_\delta(\bx)\,\rho_\delta(\bx)\left|\nabla\log\frac{\rho_\delta}{\pi}\right|^2\,\mathrm{d}\bx\,.\] The following result is standard; see, for instance,  [49], [53]–[56]. We recall it for convenience.

The key point of Lemma 3 is that, thanks to the uniform ellipticity \(a_\delta\ge a_{\min}>0\), the dissipation \(\cI_{a_\delta}\) controls the (weighted) Fisher information uniformly in \(\delta\), so the LSI yields a \(\delta\)-independent entropy decay rate. The two resulting uniform bounds, on the Fisher information and on its time integral ?? , furnish, respectively, spatial and (weak) time regularity of \(\rho_\delta\). In 3.2 we turn these into compactness of \(\{\rho_\delta\}\), from which we extract a convergent subsequence and pass to the limit \(\delta\to0\).

3.2 Passage to sharp-interface limit↩︎

The goal of this section is to pass to the limit \(\delta \to 0\) in the regularized Fokker–Planck equation and recover a weak solution of the hybrid transmission problem introduced in Section 2. We follow the standard strategy of using the compactness provided by the entropy dissipation estimates to extract a convergent subsequence and then identifying the limit by passing to the weak formulation [57].

First, we observe that the assumption about the initial condition in 14 immediately yields the following bounds, which we will use throughout this section.

3.2.1 Compactness of solution sequences↩︎

The uniform entropy dissipation we see gives two key bounds on the interval \([0,T]\). The first is that \(\sqrt{\rho_\delta}\) is bounded in \(L^2((0,T);H^1_{\mathrm{loc}})\), and the other is that \(\partial_t\rho_\delta\) is bounded in \(L^2((0,T);H^{-1}_{\mathrm{loc}})\). We establish these in the following two lemmas. The strong local compactness of \(\rho_\delta\) for positive times will allow us to take \(\delta\to 0\) using the lower-semicontinuity of \(\mathrm{KL}\) divergence.

We record a \(\delta\)-uniform local sup bound on \(\rho_\delta\), used repeatedly below. The point is that this bound depends only on the ellipticity ratio and the drift bound, hence stays uniform as \(\delta\to0\), even though the classical (Schauder) regularity constants degenerate since \(\nabla^2 a_\delta=\cO(\delta^{-2})\).

We first establish the following \(H^1\) control of \(\sqrt{\rho_\delta}\) for \(t>0\).

We now sharpen Lemma 6 by establishing a uniform \(H^1\) bound on the product \(a_\delta \rho_\delta\). This bound will be the key ingredient in deducing the continuity of \(\rho\) across the interface \(\Gamma\) in the sharp-interface limit; see Theorem 3.

The next step is to show the following time-derivative bound in \(H^{-1}_{\mathrm{loc}}\).

3.2.2 The sharp-interface limit↩︎

We now show that \(\rho_\delta(t)\) converges to \(\rho(t)\), in an appropriate sense, for some \(\rho(t)\) that solves 7 in the distributional sense.

It turns out that the limit of \(\rho_\delta\) as \(\delta\to 0\) is indeed a distributional solution to 7 .

It is important to realize that the transmission conditions across \(\Gamma\) in 7 are encoded in the weak formulation. Whenever the limiting density has enough regularity for traces to be defined, this weak formulation recovers the classical transmission conditions across \(\Gamma\). Moreover, we can use Lemma 7 to deduce a corresponding \(H^1\) regularity statement for the limit \(\rho\), and consequently the continuity of \(\rho\) across the interface \(\Gamma\).

In fact, the continuity statement ?? holds for almost every \(t \in (0,T)\). After invoking the weak \(L^1\)-continuity result of Theorem 4 in the next subsection, one may upgrade this to: for every \(t \in [0,T]\), the trace identity ?? holds in \(H^{1/2}(\Gamma)\). Note also that the conclusion ?? concerns the trace of \(\rho\) on \(\Gamma\), not the trace of its gradient. The normal derivative \(\bn_\Gamma \cdot \nabla \rho\) is generally not continuous across \(\Gamma\): indeed, by the flux-matching condition \([\bn_\Gamma \cdot J(\rho)]_\Gamma = 0\) encoded in 7 , together with the discontinuity of \(b\) across \(\Gamma\), the normal flux \(\nabla(a\rho) - b\rho\) must adjust its normal-derivative contribution to compensate for the jump in \(b\).

3.3 Exponential rate for hybrid dynamics↩︎

We now use the uniform entropy dissipation estimates and the result in the previous subsection to show that the exponential \(\mathrm{KL}\) rate does not deteriorate as \(\delta\to 0\), that is, the hybrid dynamics converges to the Gibbs distribution.

We first prove the following weak \(L^1\)-continuity result for the limit \(\rho\), the weak solution to 7 .

We are finally ready to show that the hybrid dynamics converges to the equilibrium exponentially fast.

The rate from Theorem 5 is \(2a_{\min}\lambda_{\text{LSI}}\) where \(a_{\min} = \eps \min\{1, e^{(F_{\min} - F_0)/\eps}\}\). If the interior well is deeper than the switching level (that is, \(F_{\min}<F_0\), which is the interesting case, as otherwise the hybridization region does not contain the wells), \(a_{\min}\) is exponentially small in \(1/\eps\), and the rate from Theorem 5 is worse than the analogous bound for pure Langevin. Thus, Theorem 5 should be understood as a qualitative “rate does not deteriorate under regularization" result rather than a quantitative improvement over pure Langevin. The genuine quantitative advantage of the hybrid construction is captured at the level of mean exit times and is established in Section 4 next.

4 Mean exit time comparison↩︎

The analysis in the previous sections shows that the hybrid dynamics preserves the Gibbs distribution and converges to equilibrium at an exponential rate comparable to that of the regularized system. The entropy-based convergence results, however, do not fully capture the mechanism responsible for slow mixing in multimodal landscapes, that is, the presence of metastable states separated by energy barriers.

We now complement the entropy-based analysis by studying the metastability properties of the hybrid dynamics. In particular, we investigate how the hybrid construction modifies the transition mechanism between metastable wells and whether this can potentially lead to an improvement in transition times. We will focus on the limit object whose existence was established in Theorem 2. Therefore, we will use the explicit discontinuous coefficients without further justification.

We now introduce the second partitioning convention. In the preceding sections, the switching interface was the full level set \[\Gamma_{F_0}:=\{\bx\in\Omega:F(\bx)=F_0\},\] which separates \(\Omega\) into two regions. In the metastability setting considered below, the potential is radially symmetric and may be nonconvex, so \(\Gamma_{F_0}\) may have several connected components. We therefore take the switching interface to be only the outermost connected component, denoted by \(\Gamma_{F_0}^{\mathrm{out}}\), namely the component adjacent to the exterior region, and set \(\Gamma:=\Gamma_{F_0}^{\mathrm{out}}\). The regions \(\Omega_-\) and \(\Omega_+\) are then defined as the interior and exterior regions separated by this component: \(\Omega_-\) is the region enclosed by \(\Gamma\), while \(\Omega_+\) is the exterior region. The hybrid coefficients are defined by the same sharp-interface construction as before, but with this outermost component serving as the interface. Since \(F\) is radially symmetric in this setting, the partition can equivalently be described in terms of the radial variable \(r=\|\bx\|\).

4.1 Radially symmetric setting↩︎

We will compare the hybrid dynamics with the standard overdamped Langevin dynamics in a radially symmetric landscape. We consider the case when the state space \(\Omega=\bbR^d\).

Let \(r>0\) be given, and \(U(s)\in \cC^2([0,\infty))\) be a potential that has a double-well structure on \([0,r)\) and is strictly convex on \([r,\infty)\); see Figure 2. We construct the potential \(F\) from \(U\) as \[F(\bx)=U(|\bx|)\,.\] We assume that \(U'(0)=0\) and that \(U\) satisfies the usual even radial compatibility conditions at \(s=0\), ensuring that \(F(\bx)=U(|\bx|)\) is \(\cC^2\) on \(\mathbb{R}^d\). The switching surface is the sphere \[\Gamma_r:=\{\bx\in\mathbb{R}^d:|\bx|=r\},\] which separates the inner and outer regions \[\Omega_-:=\{\bx\in\mathbb{R}^d:|\bx|<r\}, \qquad \Omega_+:=\{\bx\in\mathbb{R}^d:|\bx|>r\}.\] The hybrid dynamics 5 coincides with overdamped Langevin dynamics in \(\Omega_+\) and replaces the inner dynamics by a drift-free diffusion with appropriately chosen diffusion coefficient: \[a(\bx) = \begin{cases} a_-:=\eps e^{(F(\bx)-F_r)/\eps}, & \bx \in \Omega_-,\\ a_+:=\eps, & \bx \in \Omega_+, \end{cases} \qquad b(\bx) = \begin{cases} b_-:=0, & \bx \in \Omega_-,\\ b_+:=-\nabla F(\bx), & \bx \in \Omega_+, \end{cases}\] where \[F_r:=U(r),\] so that \(F(\bx)=F_r\) for every \(\bx\in\Gamma_r\).

With the potential \(U\) selected here, the overdamped Langevin dynamics is left unchanged outside the sphere \(\Gamma_r\), where the potential is confining and already drives the process back toward the interior. The main question is whether the switching of the dynamics to the adaptive diffusion in \(\Omega_-\) improves the mechanism responsible for transitions and, if so, what it does to the overall convergence to equilibrium.

4.2 The conductance density \(a\pi\)↩︎

For a reversible diffusion with generator 6 and invariant density \(\mu(\,\mathrm{d}\bx)=\pi(\bx)\,\,\mathrm{d}\bx\), the associated Dirichlet form is \[\mathcal{E}(f,f) = \int_{\mathbb{R}^d}a(\bx)|\nabla f(\bx)|^2\,d\mu(\bx) = \int_{\mathbb{R}^d}a(\bx)\pi(\bx)|\nabla f(\bx)|^2\,\,\mathrm{d}\bx.\] Thus the relevant coefficient in all variational problems is the product \[w(\bx):=a(\bx)\pi(\bx),\] which we view as the conductance density. Regions where \(w\) is large make variations of the committor energetically expensive, whereas regions where \(w\) is small allow the committor to change at low Dirichlet cost. Thus, small conductance regions form bottlenecks for transitions. This point is worth emphasizing: from the perspective of capacities, committors, and relaxation estimates, the process does not see \(a\) and \(\pi\) separately. What enters the variational formulas is their product \(a\pi\). Therefore, the right quantity to compare between the two dynamics is not the diffusion coefficient \(a\) alone but the conductance density \(a\pi\).

For the hybrid dynamics, the matching condition gives \[a(\bx)\pi(\bx) = \frac{\eps}{Z}e^{-F_r/\eps} \qquad\text{for }\bx\in\Omega_-.\] In particular, the conductance density is constant in the entire inner region. This simple identity is the main structural reason that the hybrid dynamics behaves differently from overdamped Langevin. The dynamics still preserves the same invariant law, but the geometry relevant to transitions is no longer governed by the wells and saddle of \(F\) in the same way as in overdamped Langevin dynamics.

4.3 Mean exit time comparison↩︎

We next compare the transition from the upper well to the global minimum. The comparison is intended to isolate the main advantage of the hybrid construction. In a double-well landscape, the difficult event is the passage from the upper local minimum to the deeper well. For overdamped Langevin, this passage is controlled by the saddle height. For the hybrid dynamics, the expectation is that the saddle no longer appears in the leading exponential barrier.

We consider a radially symmetric potential with both local minima and global minima. Assume that \(U\in \cC^2([0,\infty))\), \(s^{d-1}e^{-U(s)/\eps}\in L^1((0,\infty))\) for all sufficiently small \(\eps>0\), and that there exist radii \[0<r_g<r_s<r_u<r\] such that:

1. \(U'(r_g)=U'(r_s)=U'(r_u)=0\) and these are the only critical points of \(U\) in \((0,r]\);

2. the critical points are nondegenerate, namely, \(U''(r_g)>0\), \(U''(r_u)>0\) and \(U''(r_s)<0\);

3. \(U'(s)>0\) on \((r_g,r_s)\cup(r_u,\infty)\) and \(U'(s)<0\) on \((0,r_g)\cup(r_s,r_u)\);

4. \(U(s)\ge c_0s^\alpha-C_0\) for all sufficiently large \(s\), for some \(c_0,\alpha>0\) and \(C_0\in\mathbb{R}\).

We set \[F_g:=U(r_g),\qquad F_u:=U(r_u),\qquad H:=U(r_s),\qquad F_r:=U(r).\] and assume \(F_g < F_u < H\). See Figure 2 for an illustration of the function \(U\).

The transition of interest starts from the upper well shell \[\{\bx\in \mathbb{R}^d: |\bx|=r_u\}\] and ends on the lower well shell \[\{\bx\in \mathbb{R}^d: |\bx|=r_g\}.\]

5 Numerical simulations↩︎

We now present numerical simulations illustrating two features of the proposed hybrid dynamics. First, we verify empirically that the hybrid and regularized hybrid dynamics sample the desired Gibbs distribution. Second, we illustrate the accelerated redistribution of probability mass across metastable regions predicted by the mean-exit-time analysis of Section 4.

We are interested in comparing results from the overdamped Langevin dynamics (OL), the hybrid dynamics 5 (HD), and the regularized hybrid dynamics with regularization parameter \(\delta\) 10 (\(\delta\)-HD). We discretize the continuous systems using the standard Euler-Maruyama scheme. In particular, with the time stepsize \(\eta\), the overdamped Langevin becomes \[X_{n+1}=X_n-\eta \nabla F(X_n)+\sqrt{2\eta\eps}\,\xi_n, \qquad \xi_n\sim \cN(\bzero,\bI).\] For the \(\delta\)-regularized hybrid dynamics, we have \[X_{n+1}=X_n + \eta b_\delta(X_n) +\sqrt{2\eta a_\delta(X_n)}\,\xi_n\,.\] The limiting hybrid dynamics gives the piecewise Euler-Maruyama update \[X_{n+1}= \begin{cases} X_n+\sqrt{2\eta\,\eps\exp\left((F(X_n)-F_0)/\eps\right)}\,\xi_n, & X_n \in\Omega_-,\\[1ex] X_n-\eta \nabla F(X_n)+\sqrt{2\eta\eps}\,\xi_n, & X_n \in\Omega_+, \end{cases}\] where \(F_0\) is the chosen value that determines the interface \(\Gamma\).

5.1 1D sampling examples↩︎

We first show some one-dimensional simulations.

5.1.0.1 Example I (localized initialization).

In the first numerical example, we consider the one-dimensional double-well potential (see Figure 3): \[F(x)=\frac{1}{10}\bigl(x^2-c^2\bigr)^2+0.1, \qquad c=1.5.\] The corresponding Gibbs equilibrium, for a given temperature parameter \(\eps\), is \[\pi_G(x)=Z^{-1}\exp\{-F(x)/\eps\}\,.\] We take a small \(\eps=0.08\) to make the sampling challenging. The interior region for the derivative-free dynamics is taken to be \(\Omega_-:=\{x\in\bbR \mid |x|<2\}\). The switching interface consists of the two outer points \(x=\pm L\) with \(L=2\). These points lie on the level \(F_0=F(\pm L)=0.40625\). Thus this example uses the outermost-component convention: the hybridized region is the interval enclosed by the outermost level-set points.

For the coefficient regularization scheme 10 , we take \[S(t)=10t^3-15t^4+6t^5,\qquad 0\le t\le 1,\] and define the \(C^2\) cutoff \[\chi_\delta(x)= \begin{cases} 1, & |x|\le L-\delta,\\ 1-S\!\left(\dfrac{|x|-(L-\delta)}{\delta}\right), & L-\delta<|x|<L,\\ 0, & |x|\ge L. \end{cases}\] We test \(\delta\in\{0.02,0.05,0.10,0.20\}\). The limiting HD method corresponds to the discontinuous-interface limit \(\delta=0\).

The simulations are “cold-started": they are initialized from a very concentrated distribution \(\cN(3,0.01^2)\). All methods use the same time step \(\eta=10^{-3}\), are initialized from \(N(3,0.01^2)\), and are evolved up to \(T=2000\). For each method, we use \(4500\) trajectories. The reported KL and \(\chi^2\) divergences are computed from normalized histogram bin masses against the normalized grid approximation of \(\pi_G\), using histogram spacing \(\Delta x=0.05\). The convergence curves use pooled histograms over windows of \(50\) saved frames, corresponding to a physical window length of \(50\) time units.

We show in Figure 4 the empirical densities given by the different schemes as well as the corresponding Gibbs density. At the final time of the simulation, the HD and \(\delta\)-HD produced the densities that are closer to the Gibbs than the overdamped Langevin. Figure 5 shows the KL convergence history. The regularized HD curves are nearly indistinguishable from the limiting HD curve at this time step and sample size. After a short initial transient, the HD curves overtake OL around \(t=83\) and remain below OL for the rest of the simulation. In particular, all HD variants reach \(\mathrm{KL}\le 10^{-2}\) by time approximately \(810\), while OL does not reach this accuracy by the final time \(T=2000\). At the final time, OL has KL divergence \(3.53\times 10^{-2}\), whereas the HD variants have final KL divergences between \(3.88\times 10^{-4}\) and \(4.51\times 10^{-4}\).

The comparison across \(\delta\), in both Figure 5 and Table 1, indicates that the observed acceleration persists under \(C^2\) regularization of the HD interface. The regularized dynamics replace the abrupt switch by a thin transition layer in which both \(a_\delta\) and the correction drift \(b_\delta\) vary smoothly, but the observed convergence is essentially unchanged for the tested widths. This indicates that the main effect comes from flattening the weighted conductance \(a(\bx)\pi(\bx)\) inside the hybridized region. The drift-free interior dynamics removes the usual energetic bottleneck in the Dirichlet form, while the exterior OL dynamics preserves confinement and pulls trajectories back toward the target region. The small differences among the \(\delta\)-curves are comparable to the Monte Carlo resolution of the pooled histogram estimator. In this example, regularization improves the smoothness of the numerical dynamics without degrading the sampling advantage of HD over OL.

We also tested the fully derivative-free dynamics obtained by applying the interior HD diffusion coefficient on the entire real line, \[X_n+\sqrt{2\eta\,\eps\exp\left((F(X_n)-F_0)/\eps\right)}\,\xi_n.\] This experiment is not included as a density plot because it is numerically unstable at the same step size and initialization used above. Starting from \(N(3,0.01^2)\), all \(4500\) trajectories leave the plotting window \([-4,4]\) after the a few time steps. Thus the pure derivative-free dynamics does not provide a meaningful sampling baseline for this setup. The observation instead reinforces the role of the hybrid switch: large diffusion is useful inside the low-energy region, but outside this region the OL part is needed to provide confinement.

5.1.0.2 Example II (uniform initialization).

In the second numerical example, we consider the asymmetric potential given by (see Figure 6): \[F(x)=3(x^2-1)^2+0.3(x+1)^2-0.4\,.\] We choose the interior region for the derivative-free dynamics as \(\Omega_-:=(-1.2753, 1)\). This means that \(F_0=0.8\) for the potential. Note that \(\Omega_-\) is not symmetric with respect to \(x=0\) due to the non-symmetry of \(F\).

We start the simulations from a uniform density \(\cU([-2.5, 2.5])\). We compare here only the hybrid dynamics and the overdamped Langevin dynamics. Figure 7 shows the comparison of the empirical snapshots of the densities \(\wh\rho^N(t,\cdot)\) at \(t=1, 2, 10\) and \(30\) from \(N=2\times 10^3\) trajectories. Both samplers reproduce the two-peaked structure of \(\pi\) fairly well by the time \(t=30\). However, we see a clear advantage of the hybrid dynamics in this case from the evolution of the snapshots and the \(\mathrm{KL}\) divergence \(\mathrm{KL}(\wh\rho^N(t,\cdot)\|\pi)\) in Figure 8. The hybrid dynamics quickly redistribute the extra mass from the smaller peak at \(x=1\) to the larger peak at \(x=-1\), resulting in faster convergence than the overdamped Langevin dynamics.

5.2 2D sampling examples↩︎

We next consider a two-dimensional example designed to test whether HD can rapidly redistribute probability mass among several metastable radial wells.

5.2.0.1 Example III (radially symmetric potential).

Let \(r=|\bx|\) and \(q=r^2=|\bx|^2\). We define \[\theta(q)=6\pi(q-0.04),\] and \[G(q) =0.90\left(1-\cos\theta(q)\right) +0.45\left(1-\cos(2\theta(q))\right) +0.35\sin^2\theta(q) +0.18\sin(3\theta(q)).\] With \(\delta=0.08\), set \[\psi_\delta(s)= \begin{cases} 0, & s\le 0,\\[0.5ex] \delta^2\left(t^5-3t^4+3t^3\right), \qquad t=s/\delta, & 0<s<\delta,\\[0.5ex] s^2, & s\ge \delta. \end{cases}\] The radial profile is \[V_{\rm rad}(r)=G(r^2)+80\,\psi_\delta(r-1),\] and the two-dimensional potential is \[V(\bx)=V_{\rm rad}(|\bx|) =G(|\bx|^2)+80\,\psi_\delta(|\bx|-1).\] Thus \(V\in \cC^2(\mathbb{R}^2)\), has several oscillatory radial wells in the unit disk, and grows quadratically for \(|\bx|\ge 1+\delta\). The target Gibbs distribution is \[\pi_\eps(\bx)=Z_\eps^{-1}\exp(-V(\bx)/\eps), \qquad \eps=0.35 .\] Figure 9 shows the radial profile, the HD level-set region, and the Gibbs density.

Here the HD interior diffusion coefficient is \(D(\bx)=\eps\exp((V(\bx)-F_0)/\eps)\), so that the hybrid update uses increment \(\sqrt{2\eta D(X_n)}\xi_n\) on \(\{V<F_0\}\) and the OL update on \(\{V\ge F_0\}\). In all experiments below \(F_0=2\), so \(D(\bx)=\eps\) continuously on the interface \(V(\bx)=F_0\). Both OL and HD methods use step size \(\eta=10^{-3}\), final time \(T=200\), and \(16000\) particles per run. Reported \(\mathrm{KL}\) values are computed from normalized histogram cell masses against the normalized grid approximation of \(\pi_\eps\). We use six independent runs and also report the KL of the pooled histogram.

We test two initial distributions: \[\mu_0^{\mathrm{well}}:\quad X_0=0.18(\cos\Theta,\sin\Theta)+0.02Z, \qquad \Theta\sim\mathrm{Unif}(0,2\pi),\quad Z\sim N(0,I_2),\] and \[\mu_0^{\mathrm{disk}}:\quad X_0\sim \mathrm{Unif}\{\bx\in\mathbb{R}^2:\;|\bx|\le 2\}.\] The first initial distribution is concentrated near the innermost well, while the second one places mass well outside the target region.

The results are shown in Figure 10 and summarized in Table 2. For the near-well initial distribution \(\mu_0^{\mathrm{well}}\), HD reduces the pooled KL divergence from order one to \(7.37\times10^{-2}\) by final time, while OL remains near \(1.21\). For the broad uniform-disk initialization \(\mu_0^{\mathrm{disk}}\), HD reaches essentially the same accuracy: the best recorded pooled KL is \(7.21\times10^{-2}\) at \(t=60\), and the final pooled KL is \(7.44\times10^{-2}\). In contrast, OL again remains near \(1.21\). The final density plots show that OL quickly contracts mass toward the inner well and ignore the outer target rings, whereas HD equilibrates the radial shell masses much more effectively.

5.2.0.2 Example IV (anisotropic potential).

We next deform the radial example into a genuinely anisotropic two-dimensional example by replacing the Euclidean radius with an elliptic radius. This tests whether HD still rapidly redistributes probability mass when the metastable wells are organized along elliptic, rather than circular, level sets. Define \[s(\bx)=\sqrt{\left(\frac{x_1}{1.45}\right)^2+ \left(\frac{x_2}{0.75}\right)^2}, \qquad \bx=(x_1,x_2)\in\mathbb{R}^2,\] and set \(g(\bx)=s(\bx)^2\). Using the same oscillatory profile \(G\) and the same \(C^2\) confining regularization \(\psi_\delta\) from Example III, with \(\delta=0.08\), we define \[V(\bx)=G(g(\bx))+80\,\psi_\delta(s(\bx)-1).\] Thus \(V\in\mathcal{C}^2(\mathbb{R}^2)\), has several oscillatory wells along elliptic shells, and grows quadratically as \(s(\bx)\to\infty\). The target Gibbs distribution is \[\pi_\eps(\bx)=Z_\eps^{-1}\exp(-V(\bx)/\eps), \qquad \eps=0.35 .\] Figure 11 shows the elliptic-radius profile, the HD level-set region, and the Gibbs density.

We use the same OL and HD discretizations as in Section 5.2. In all experiments below \(F_0=2\), so \[D(\bx)=\eps\exp((V(\bx)-F_0)/\eps)\] satisfies \(D(\bx)=\eps\) continuously on the interface \(V(\bx)=F_0\). Both methods use step size \(\eta=10^{-3}\), final time \(T=200\), and \(16000\) particles per run. \(\mathrm{KL}\) values are computed from normalized histogram cell masses against the normalized grid approximation of \(\pi_\eps\). We use six independent runs and report the \(\mathrm{KL}\) of the pooled histogram.

We test two initial distributions: \[\mu_0^{\mathrm{ell}}:\quad X_0=\bigl(1.45\cdot 0.18\cos\Theta,\;0.75\cdot 0.18\sin\Theta\bigr)+0.02Z, \qquad \Theta\sim\mathrm{Unif}(0,2\pi),\quad Z\sim N(0,I_2),\] and \[\mu_0^{\mathrm{disk}}:\quad X_0\sim \mathrm{Unif}\{\bx\in\mathbb{R}^2:\;|\bx|\le 2\}.\] The first initial distribution is concentrated near the innermost elliptic well, while the second one places mass well outside the target region.

The results are shown in Figure 12 and summarized in Table 3. For the near-well-elliptic initialization \(\mu_0^{\mathrm{ell}}\), HD reduces the pooled KL divergence to \(1.09\times 10^{-1}\) by the final time, while OL remains near \(1.27\). For the broad uniform-disk initialization \(\mu_0^{\mathrm{disk}}\), HD again achieves essentially the same accuracy, with a final pooled KL of \(1.09\times 10^{-1}\), whereas OL remains near \(1.27\). The final density plots show the same behaviors of the two dynamics as in the radial example: OL moves mass toward the innermost well and fails to put enough mass on the outer elliptic shells. At the same time, HD equilibrates mass across all shell structures.

6 Concluding remarks↩︎

This work introduces a hybrid dynamics framework for sampling from Gibbs distributions by combining two distinct stochastic dynamics in different regions of the state space. The construction is designed so that the resulting process preserves the target Gibbs distribution while modifying the dynamics’ effective geometry in regions where metastability is expected to dominate.

We developed a rigorous analysis of the hybrid dynamics through a regularization approach. We showed that the hybrid dynamics converges exponentially fast to the Gibbs distribution. We also analyzed the metastability properties of the hybrid dynamics in a radially symmetric landscape. We showed that the transition mechanism between metastable states is fundamentally altered. While overdamped Langevin dynamics is governed by the saddle barrier, the hybrid dynamics replaces it with the interface’s switching level. As a consequence, an appropriate choice of the switching surface leads to exponentially faster transition times, demonstrating a clear advantage of the hybrid construction in multimodal landscapes.

A central question for the future is how to choose the switching interface in a systematic and possibly adaptive manner. In complex high-dimensional landscapes, it would be desirable to identify regions where the hybridization yields the greatest improvement, potentially using data-driven or learning-based approaches. Also, the present framework combines overdamped Langevin dynamics with a derivative-free diffusion. It would be of interest to extend the hybrid approach to incorporate other dynamics, such as underdamped Langevin processes, preconditioned diffusions, or nonreversible dynamics, and to understand how these choices affect both convergence rates and metastability.

Data availability statements↩︎

Data sets generated during the current study are available from the corresponding author on reasonable request.

Declarations↩︎

The authors declare no competing interests. This work is partially supported by the National Science Foundation through grants DMS-2208504 (BE), DMS-1937254 (KR), DMS-2309802 (KR), and DMS-2409855 (YY). YY was also partially supported by Office of Naval Research through grant N00014-24-1-2088.

References↩︎

---

1. Department of Mathematics and the Oden Institute, The University of Texas at Austin, Austin, TX 78731; engquist@oden.utexas.edu↩︎

2. Department of Applied Physics and Applied Mathematics, Columbia University, New York, NY 10027; kr2002@columbia.edu↩︎

3. Department of Mathematics, Cornell University, Ithaca, NY 14850; yunan.yang@cornell.edu↩︎

4. This work is partially supported by the National Science Foundation through grants DMS-2208504 (BE), DMS-1937254 (KR), DMS-2309802 (KR), and DMS-2409855 (YY). YY was also partially supported by Office of Naval Research through grant N00014-24-1-2088.↩︎