Key words. Stochastic differential equations with switching and jumps, hybrid systems, randomized Hamiltonian Monte Carlo, discrete Poisson equation, computing ergodic limits.

AMS subject classifications. 65C30, 60H35, 60H10, 37H10

1 Introduction↩︎

Sampling techniques are widely used to evaluate risk and measure unknowns, finding parameters within statistical frameworks, compute thermodynamic variables in physics simulations and designing new global optimization methods. On the other hand, regime-switching jump-diffusion stochastic differential equations (SDEs) find a large number of applications in biology, insurance, control etc. (for more in-depth exposition see [1], [2]). The underlying process in such hybrid systems has two components \((Y(t), \mathfrak{s}(t))\) where \(Y\) represents the jump-diffusion and \(\mathfrak{s}\) describes switching mechanism. This paper utilizes the framework of switching stochastic dynamics for the purpose of sampling from mixture distributions as in [3]. Mixture distributions play an important role in statistical modeling and Bayesian inference [4], [5].

Arguably, discrete-time Hamiltonian Monte Carlo [6] and its variants [7], [8] are among the most widely used sampling methods. In this paper, we take a continuous-time perspective on Hamiltonian Monte Carlo. We consider hybrid-Hamiltonian dynamics which includes a switching mechanism and is interlaced with randomized refreshments. We write the resulting regime-switching jump dynamics as SDEs driven by Poisson random measure with added regime-switching mechanism in the drift term so that the target measure is a mixture Boltzmann-Gibbs distribution. Then the contributions of this paper are as follows:

• We propose new numerical integrators based on splitting technique to discretize switching Hamiltonian-type SDEs where continuous-time finite state space conditionally Markov chain is simulated via uniformization technique or Gillespie’s algorithm. The proposed schemes improve the error incurred in estimating ergodic averages to \(\mathcal{O}(h^2)\) from \(\mathcal{O}(h)\) bias resulting from switching Langevin algorithm proposed in [3] (here \(h\) represents the discretization step).

• We present a new approach based on discrete Poisson equation associated with numerical scheme itself for obtaining error bounds in integral probability metric between desired measure \(\mu\) and the invariant measure \(\mu_h\) of Markov chain governed by the rules of proposed numerical schemes, i.e. \[\begin{align} \left|\int_{\mathbb{R}^d} \varphi \mathop{}\!\mathrm{d}\mu - \int_{\mathbb{R}^d} \varphi \mathop{}\!\mathrm{d}\mu_h \right| \leq Ch^2, \end{align}\] where \(\varphi\) belongs to a class of test functions allowing for bounded measurable functions, Lipschitz functions etc., and \(C\) is independent of \(h\) and is dependent on the solution of a discrete Poisson equation.

The other contributions which are of independent interest and serve as the building blocks of the above main results are as follows. We first show the geometric ergodicity of continuous-time dynamics. Writing the dynamics as SDEs allows us to use tools available in stochastic numerics literature to develop numerical integrators for randomized Hamiltonian Monte Carlo dynamics. A lot of attention has been given in designing and analyzing numerical integrators for kinetic Langevin dynamics, for example [9]–[18] but the same can not be said about Hamiltonian dynamics interlaced with Poisson refreshments for which we propose splitting schemes and also show their uniform in \(h\) geometric ergodicity. We now position our work within the literature on numerical integration of switching SDEs. Existing convergence results for numerical schemes for regime-switching processes cover several different settings. For switching diffusions without jumps, [19] proved weak convergence of a numerical scheme without rate, [20] established an optimal 1/2 rate in \(L^1\) and [21] investigated continuous time Markov chain approximation scheme where the generator itself has to be discretized. For regime-switching jump diffusions, convergence of the Euler method is investigated in [22] in \(L^2\) with decreasing step sizes. Here, we take splitting approach for numerical approximation which allows to use uniformization technique or Gillespie’ methods to simulate continuous time chain governing switching mechanism.

We denote with \((a\cdot b)\) scalar product between two vectors \(a, b \in \mathbb{R}^d\) and with \(|a|\) norm of \(a \in \mathbb{R}^d\). We denote with \(\mathcal{N}(0, I_d)\) \(d\)-dimensional standard Gaussian distribution and with \(\mathcal{U}(a,b)\) the uniform distribution on \((a,b)\) with \(a< b \in \mathbb{R}\). We represent \(I_d\) and \(\mathbf{0}_d\) the identity matrix and zero matrix, respectively. For a matrix \(P \in \mathbb{R}^{n \times n}\), we denote its matrix norm as \(\|P\|\).

2 Switching Hamiltonian Monte Carlo↩︎

Andersen’s thermostat [23] models a molecular system using Hamiltonian dynamics with random Poisson jump updates in the momentum variable of a randomly chosen particle. These random refreshments of randomly selected particles imitate collisions. In the same spirit, the authors of [24] introduced randomized Hamiltonian Monte Carlo (RHMC), where Hamiltonian dynamics is interlaced with Poisson refreshments, to sample from the Boltzmann-Gibbs measure defined as \[\begin{align} \mu(\text{d}x, \text{d}v) \propto e^{-U(x) - |v|^2/2} \mathop{}\!\mathrm{d}x \mathop{}\!\mathrm{d}v . \end{align}\] In RHMC, starting at time \(t_0\), a random time \(\delta t\) is drawn from an exponential distribution with mean \(1/\lambda\), and the system evolves according to standard Hamiltonian dynamics until \(t_1 = t_0 + \delta t\). At time \(t_1\), the position remains unchanged, but the velocity is partially refreshed by replacing it with \[V(t_1) = \cos(\phi)V(t_1^-) + \sin(\phi)Z,\] where \(\phi \in (0, \pi/2]\) determines the amount of randomness added and \(Z \sim \mathcal{N}(0, I_d)\). We can write the RHMC dynamics in terms of SDEs governed by Poisson random measure as follows: \[\label{eqn95Poishmc} \begin{align} \mathop{}\!\mathrm{d}X(t) & = V(t) \mathop{}\!\mathrm{d}t, \\ \mathop{}\!\mathrm{d}V(t) & = -\nabla U(X(t)) \mathop{}\!\mathrm{d}t + \int_{\mathbb{R}^d}(- V(t^{-}) + \cos(\phi)V(t^{-}) + \sin(\phi) z) N(\mathop{}\!\mathrm{d}t ,\mathop{}\!\mathrm{d}z), \end{align}\tag{1}\] where \(N(\mathop{}\!\mathrm{d}t, \mathop{}\!\mathrm{d}z)\) denotes the Poisson random measure with intensity \(\lambda \nu( \mathop{}\!\mathrm{d}z)\mathop{}\!\mathrm{d}t\), \(\nu( \mathop{}\!\mathrm{d}z) \propto e^{-|z|^2/2} \mathop{}\!\mathrm{d}z\). Under the assumptions of local one-sided Lipschitz continuity and super-linear growth of \(-\nabla U\), there exists a unique strong solution of 1 (see [25]). As can be verified \(\rho\) satisfies the corresponding Fokker-Planck partial integro-differential equation (PIDE):

\[\begin{align} &- (v \cdot \nabla_x \rho) + (\nabla U(x) \cdot \nabla_v \rho) + \frac{\lambda}{(2 \pi)^{d/2} \sin^d (\phi)} \int_{\mathbb{R}^d} \rho(x, z) e^{-\frac{|v - \cos(\phi) z|^2}{2 \sin^2(\phi)}}\mathop{}\!\mathrm{d}z - \lambda\rho(x,v)= 0. \end{align}\]

Applications requiring sampling from mixture distributions and the computation of their ergodic averages arise in several statistical models and hybrid-control systems. For this purpose, we present below a switching Hamiltonian Monte Carlo method. Let \(\mathcal{M}=\{1,\ldots ,m\}\) denote the finite set of regimes, let \(U(x; m) > 0\), \(x \in \mathbb{R}^d\), \(m \in \mathcal{M}\), be sufficiently smooth functions and consider a given finite mixture distribution with sufficiently smooth density \[\rho(x) \propto \sum_{i=1}^{m} \rho(x; i),\;\;\text{ where} \;\;\rho(x; i) = w_i \exp(-U(x; i)),\] \(w_i > 0\) are weights. With a slight abuse of notation which should not lead to any confusion, we introduce the mixture Boltzmann-Gibbs density as \[\begin{align} \rho(x,v) = \frac{1}{\mathcal{Z}} \sum_{i=1}^{m} \rho(x, v; i),\;\; \text{where }\;\; \rho(x,v;i) = w_i \exp(-U(x; i) - |v|^2/2), \end{align}\] and \(\mathcal{Z}\) is the normalization constant. Let us now introduce the transition-rate matrix governing the regime-switching mechanism. \[\label{swi95hmc95Qmatrix} {Q}(x)=\left[ \begin{array}{cccc} -\mathfrak{q}_{1}(x) & \mathfrak{q}_{12}(x) & \cdots & \mathfrak{q}_{1m}(x) \\ \mathfrak{q}_{21}(x) & - \mathfrak{q}_{2}(x) & \cdots & \mathfrak{q}_{2m}(x) \\ \vdots & \vdots & \cdots & \vdots \\ \mathfrak{q}_{m1}(x) & \mathfrak{q}_{m2}(x) & \cdots & - \mathfrak{q}_{m}(x)\end{array}\right] ,\quad x\in \mathbb{R}^d,\tag{2}\] where \(\mathfrak{q}_{i}(x)\) and \(\mathfrak{q}_{ij}(x)\) are continuous non-negative functions bounded in \(\mathbb{R}^{d}\) and the following relation holds: \[\sum_{j\neq i}\mathfrak{q}_{ij}(x)=\mathfrak{q}_{i}(x),\;\;i\in \mathcal{M}. \label{hmc95switching95eqn95q}\tag{3}\] Therefore, the regime-switching jump SDEs governing switching Hamiltonian Monte Carlo dynamics are \[\label{eqn95switching95Poishmc} \begin{align} \mathop{}\!\mathrm{d}X(t) & = V(t) \mathop{}\!\mathrm{d}t \\ \mathop{}\!\mathrm{d}V(t) & = -\nabla U(X(t); \mathfrak{s}(t)) \mathop{}\!\mathrm{d}t + \int_{\mathbb{R}^d}(- V(t^{-}) + \cos(\phi)V(t^{-}) + \sin(\phi) z) N(\mathop{}\!\mathrm{d}t ,\mathop{}\!\mathrm{d}z), \end{align}\tag{4}\] where \(N(\mathop{}\!\mathrm{d}t, \mathop{}\!\mathrm{d}z)\) is as in 1 and \(\mathfrak{s}(t)\) is a continuous-time finite state conditionally Markov process with right-continuous sample paths taking values in \(\mathcal{M}\) and \[\begin{align} \mathbb{P}(\mathfrak{s}(t+\delta ) &= j\mid\mathfrak{s}(t)= i ,X(t)=x)=\mathfrak{q}_{ij}(x)\delta +o(\delta ),\;i\neq j, \label{eq:sde2} \\ \mathbb{P}(\mathfrak{s}(t+\delta ) &= i\mid \mathfrak{s}(t)=i,X(t)=x)=1-\mathfrak{q}_{i}(x)\delta +o(\delta ) \notag \end{align}\tag{5}\] as \(\delta \downarrow 0.\) The Poisson random measure \(N(\mathop{}\!\mathrm{d}t, \mathop{}\!\mathrm{d}z)\) and the continuous-time conditional Markov chain governing the transitions of \(\mathfrak{s}(t)\) are independent. Here, the corresponding Fokker-Planck equation with adjoint operator \(\mathcal{A}^*\) is \[\begin{align} \mathcal{A}^* \rho(x,v; \mathfrak{s}) = & -(v \cdot\nabla_x \rho(x,v;\mathfrak{s})) + (\nabla U(x; \mathfrak{s})\cdot\nabla_v \rho(x,v; \mathfrak{s}) ) \\ & + \frac{\lambda}{(2 \pi)^{d/2} \sin^d (\phi)} \int_{\mathbb{R}^d} \rho(x, z;\mathfrak{s}) e^{-\frac{|v - \cos(\phi) z|^2}{2 \sin^2(\phi)}}\mathop{}\!\mathrm{d}z - \lambda\rho(x,v;\mathfrak{s}) \\ & - \mathfrak{q}_\mathfrak{s}(x) \rho(x,v; \mathfrak{s}) + \sum_{j \neq \mathfrak{s}} \mathfrak{q}_{j\mathfrak{s}}(x) \rho(x, v; j). \end{align}\]

To ensure that 3 leaves the desired mixture distribution invariant, i.e. \(\mathcal{A}^*\rho = 0\) we require the following relation to hold: \[-\mathfrak{q}_\mathfrak{s}(x) \rho(x, v; \mathfrak{s}) + \sum_{j \neq \mathfrak{s}} \mathfrak{q}_{j\mathfrak{s}}(x) \rho(x ,v; j) = 0, \quad \mathfrak{s}\in \mathcal{M},\] or equivalently \(\sum_{j \neq \mathfrak{s}} [\mathfrak{q}_{j\mathfrak{s}}(x) \rho(x, v; j) - \mathfrak{q}_{\mathfrak{s}j}(x) \rho(x, v; \mathfrak{s})] = 0\), \(\mathfrak{s}\in \mathcal{M}\). One way it is satisfied if for \(j \neq \mathfrak{s}\) \[\frac{\mathfrak{q}_{j\mathfrak{s}}(x)}{\mathfrak{q}_{\mathfrak{s}j}(x)} = \frac{\rho(x , v; \mathfrak{s})}{\rho(x,v; j)}. \label{swi95eqn952469}\tag{6}\]

2.0.0.1 Geometric ergodicity.

Under Assumption 1, the regime-switching dynamics with jumps 4 is well-posed, i.e. there exists a unique strong solution [26], [27]. To establish geometric ergodicity via Meyn-Tweedie approach, we need two ingredients: the Lyapunov drift condition and the minorization condition. The Lyapunov drift condition follows from [24] under Assumption 2. However, to prove minorization condition, diffusion matrix is explicitly used in [26] (also see [27], [28]). Likewise ellipticity or local ellipticity of diffusion matrix is used in [29]. Therefore, we can not rely on these quoted results for geometric ergodicity of 4 . However, 4 , for a fixed switch \(\mathfrak{s}(t) =k\) for all \(t\), satisfies minorization condition [24]. We will use this fact to show minorization condition for 4 and thus establishing its geometric ergodicity.

Let us consider 4 again. We fix the switches \(\mathfrak{s}(t) = k\) and \(\mathfrak{s}(t) = l\) for all \(t \geq 0\) and denote the fixed switch (i.e., without any regime-switching mechanism) processes with \((X^{(k)} , V^{(k)})\) and \((X^{(l)}, V^{(l)})\), respectively. Thanks to [24], there exist times \(t_1 > 0\) and \(t_3 > 0\) such that the unswitched processes \((X^{(k)} , V^{(k)})\) and \((X^{(l)}, V^{(l)})\) satisfy minorization on a compact set \(\mathcal{C}\). We take \(t_2 >0\) to be an arbitrary intermediate transition time whose choice we will make later. Let \(t^* = t_1 + t_2 + t_3\). We apply the Chapman-Kolmogorov equation twice, expanding over all possible intermediate states \((y, m)\) at time \(t_1\), and \((z, n)\) at time \(t_1 + t_2\), to get \[\begin{align} P_{t^*}&\big((y_0, k), A \times \{l\}\big) \\ &= \sum_{i \in \mathcal{M}} \sum_{n \in \mathcal{M}} \int_{\mathbb{R}^{2d}} \int_{\mathbb{R}^{2d}} P_{t_1}\big((y_0, k), \mathrm{d}y \times \{i\}\big) P_{t_2}\big((y, i), \mathrm{d}z \times \{n\}\big) P_{t_3}\big((z, n), A \times \{l\}\big). \end{align}\] Restricting the first sum to \(i = k\), the second sum to \(n = l\), and the spatial integral over \(z\) to the compact set \(\mathcal{C}\) yields \[\begin{align} \label{eq:CK95bound} P_{t^*}&\big((y_0, k), A \times \{l\}\big) \nonumber\\ & \geq \int_{\mathbb{R}^{2d}} \int_{\mathcal{C}} P_{t_1}\big((y_0, k), \mathrm{d}y \times \{k\}\big) P_{t_2}\big((y, k), \mathrm{d}z \times \{l\}\big) P_{t_3}\big((z, l), A \times \{l\}\big). \end{align}\tag{7}\] We can get a lower bound on first and third terms of the integrand by conditioning on the event of no switching. For the third term, starting at \(z \in \mathcal{C}\) and remaining in \(l\) for duration \(t_3\), we apply the unswitched minorization condition to arrive at \[\label{eq:bound3} P_{t_3}\big((z, l), A \times \{l\}\big) \geq e^{-\bar{q}_{\max} t_3} P^{(l)}_{t_3}(z, A) \geq e^{-\bar{q}_{\max} t_3} \epsilon_{l} \nu_{l}(A),\tag{8}\] where \(\bar{q}_{\max} = \sup_{x \in \mathbb{R}^d}\max_{i \in \mathcal{M}}\mathfrak{q}_i(x)\). Similarly, for the first term, conditioning on starting at \(y_0 \in \mathcal{C}\) and remaining in \(k\) for duration \(t_1\) yields \[\label{eq:bound1} P_{t_1}\big((y_0, k), \mathrm{d}y \times \{k\}\big) \ge e^{-\bar{q}_{\max} t_1} P^{(k)}_{t_1}(y_0, \mathrm{d}y) \geq e^{-\bar{q}_{\max} t_1} \epsilon_{k} \nu_{k}(\mathrm{d}y).\tag{9}\] Using 8 and 9 in 7 gives \[\begin{align} P_{t^*}&\big((y_0, k), A \times \{l\}\big) \geq e^{-\bar{q}_{\max}(t_1+t_3)} \epsilon_k \epsilon_{l} \nu_{l}(A) \int_{\mathbb{R}^{2d}} P_{t_2}\big((y, k), \mathcal{C} \times \{l\}\big) \nu_k(\mathrm{d}y). \end{align}\] Let \(B_k\subset\operatorname{int}\mathcal{C}\), i.e. interior \(\mathcal{C}\). From [24], we know that \(\nu_k\) is a probability measure on \(\mathcal{C}\) with \(\nu_k(B_k)>0\). By continuity of the Hamiltonian flow, for \(t_2>0\) sufficiently small, the continuous component starting from any \(y\in B_k\) remains in \(\mathcal{C}\) on \([0,t_2]\) conditional on the event that there is no velocity-refreshment jump, which has probability \(e^{-\lambda t_2}>0\). If \(l=k\), the event of no switching of regime and no velocity-refreshment jump gives that \(P_{t_2}\big((y,k),\mathcal{C}\times\{k\}\big)>0\), \(y\in B_k\). If \(l\neq k\), then using continuity and strict positivity of \(\mathfrak{q}_{kl}\), the event that exactly one switch \(k\to l\) occurs during \([0,t_2]\), no further switch occurs, and no velocity-refreshment jump occurs has strictly positive probability. On this event the continuous component remains in \(\mathcal{C}\) and the terminal regime is \(l\). Hence, \(P_{t_2}\big((y,k),\mathcal{C}\times\{l\}\big)>0\), \(y\in B_k\). Consequently, \[\begin{align} \int_{\mathbb{R}^{2d}} P_{t_2}\big((y,k),\mathcal{C}\times\{l\}\big)\nu_k(\mathrm dy) \geq \int_{B_k} P_{t_2}\big((y,k),\mathcal{C}\times\{l\}\big)\nu_k(\mathop{}\!\mathrm{d}y)>0. \end{align}\] Therefore, \[\begin{align} P_{t^*}&\big((y_0, k), A \times \{l\}\big) \geq \epsilon \tilde{\nu}(A\times\{l\}), \end{align}\] where \(\tilde{\nu}(\mathop{}\!\mathrm{d}z, \{j\}) = \nu_l(\mathop{}\!\mathrm{d}z) I_{j = l}\) and \(\epsilon = \min_{k,l \in \mathcal{M}}e^{-\bar{q}_{\max}(t_1+t_3)} \epsilon_k \epsilon_{l} \int_{B_k} P_{t_2}\big((y, k), \mathcal{C} \times \{l\}\big) \nu_k(\mathrm{d}y).\)

3 Numerical integrators↩︎

In Subsections 3.1 and 3.2, we present numerical integrators for randomized HMC dynamics 1 and regime-switching HMC dynamics 4 , respectively.

3.1 [E, B] schemes↩︎

Splitting techniques are a standard way to construct higher-order weak or strong approximations of stochastic dynamics. The idea is to decompose the SDEs into components that are either exactly solvable or computationally straightforward to approximate. This approach has successfully been utilized for Hamiltonian ODEs and Langevin SDEs [12], [14], [15]. In this subsection, we consider the SDEs 1 and partition the dynamics in the following manner: \[\begin{pmatrix} \mathop{}\!\mathrm{d}X(t) \\ \mathop{}\!\mathrm{d}V(t) \end{pmatrix} = \underbrace{\begin{pmatrix} 0 \\ -\nabla U(X(t))\mathop{}\!\mathrm{d}t \end{pmatrix}}_{\mathcal{B}} + \underbrace{\begin{pmatrix} V(t) \mathop{}\!\mathrm{d}t \\ \displaystyle\int_{\mathbb{R}^d} (-V(t^-) + \cos(\phi) V(t^{-}) + \sin( \phi) z) N(\mathop{}\!\mathrm{d}t, \mathop{}\!\mathrm{d}z) \end{pmatrix}}_{\mathcal{E}},\] where each of these parts can be solved exactly.

Fix a final time \(T\). Consider a uniform partition of \([0,T]\) with discretization parameter \(h \in (0,1)\) such that the number of steps is \(\lfloor T/h\rfloor\). Let \(N_t\) denote the Poisson process with intensity \(\lambda\). Let \(N_h\) be the total number of jumps in the interval \([t_k, t_k + h)\) and conditioned on \(N_h = n\), let \(\tau_1, \tau_2, \dots, \tau_{n}\) be the exact times when jumps occur which are distributed as the order statistics simulated from \(\mathcal{U}(t_k, t_k+h)\). We define the boundaries as \(\tau_0 = t_k\) and \(\tau_{n+1} = t_k +h\). Let \(Z_i \sim \mathcal{N}(0, I_d)\), \(i = 1,\dots, n\). On any given interval \([\tau_i, \tau_{i+1})\), \(i =0,\dots, n\), the velocity in \(\mathcal{E}\) is constant and given by \[\zeta_i = v \cos^{i}(\phi) + \sin(\phi)\sum_{j =1}^{i}\cos^{i-j}(\phi)Z_j.\] This also means that integral \(\int_{0}^{h}V(s) \mathop{}\!\mathrm{d}s\) is simply a summation over the inter-jump intervals. Having all the ingredients, we can write the exact solution of \(\mathcal{B}\) and \(\mathcal{E}\) as follows: \[\begin{align} \mathscr{S}_h^{\mathrm{B}}(x,v) &= (x, v - h\nabla U(x)), \label{B95update95swi95HMC} \\ \mathscr{S}_h^{\mathrm{E}}(x,v) &= \left(x + \textstyle\sum\limits_{i=0}^{N_h}\zeta_{i} (\tau_{i+1} - \tau_{i}), \zeta_{N_h}\right), \end{align}\tag{10}\] and, for brevity, we use same notation to denote the conditional update \[\begin{align} \mathscr{S}_h^{\mathrm{E}}(x,v; N_h = n) &= \left(x + \textstyle\sum\limits_{i=0}^{n}\zeta_{i} (\tau_{i+1} - \tau_{i}), \zeta_{n}\right). \end{align}\] Since our main interest is in sampling, the update rule for \(\mathcal{E}\) step can be simplified via the Cholesky decomposition. Specifically, the same joint distribution of the update in \((x,v)\) can be preserved by generating two Gaussian random vectors \(\xi_1, \xi_2 \sim \mathcal{N}(0, I_{d})\), providing the following representation: \[\begin{align} \label{E95update95swi95HMC} \mathscr{S}_h^{\mathrm{E}} &= \left(x + m_x + \sin(\phi)\bigg( \frac{\sigma_{xv}}{\sqrt{\sigma_{vv}}} \xi_1 + \sqrt{\sigma_{xx} - \frac{\sigma_{xv}^2}{\sigma_{vv}}} \xi_2\bigg) , v\cos^{n}(\phi) + \sqrt{1 - \cos^{2n}(\phi)}\xi_1\right), \end{align}\tag{11}\] where \[\begin{align} m_x = v\sum_{i=0}^{n}\cos^{i}(\phi)(\tau_{i+1} - \tau_{i}), \quad \sigma_{vv} = \frac{1 - \cos^{2n}(\phi)}{\sin^2(\phi)}, \quad \sigma_{xx} = \sum_{l=1}^{n}a_l^2, \quad \sigma_{xv} =\sum_{l=1}^{n} a_{l}\cos^{n-l}(\phi), \end{align}\] with \(a_l = \sum_{i = l}^{n}(\tau_{i+ 1} - \tau_{i})\cos^{i-l}(\phi)\). Therefore, \(a_{l} = \tau_{l+1}- \tau_{l} + \cos(\phi) a_{l+1}\) with \(a_n = t_k + h - \tau_{n}\). We assign \(\sigma_{vv} = 1\) and \(\sigma_{xx} = \sigma_{xv} = 0\) if \(n = 0\). \[\begin{array}{cc} \text{[EBE]} & \text{[BEB]} \\ [6pt] \begin{align} (x_{1/2}, \hat{v}_{1/2}) &= \mathscr{S}_{h/2}^{\mathrm{E}}(x, v), \\ (x_{1/2},v_{1/2}) &= \mathscr{S}_{h}^{\mathrm{B}}(x_{1/2}, \hat{v}_{1/2}), \\ (x_{1}, v_{1}) &= \mathscr{S}_{h/2}^{\mathrm{E}}(x_{1/2}, v_{1/2}). \end{align} & \begin{align} (x,\hat{v}_{1/2}) &= \mathscr{S}_{h/2}^{\mathrm{B}}(x, v), \\ (x_{1}, v_{1/2}) &= \mathscr{S}_{h}^{\mathrm{E}}(x, \hat{v}_{1/2}), \\ (x_1,v_{1}) &= \mathscr{S}_{h/2}^{\mathrm{B}}(x_{1}, v_{1/2}). \end{align} \end{array}\]

Instead of decomposing into deterministic flow and jump component as is done for PDMPs [30], in the proposed schemes above we split the deterministic drift into two subflows and include one of them in the jump component leading to \(\mathrm{[BEB]}\) and \(\mathrm{[EBE]}\) integrators. This construction is motivated by viewing RHMC as a stochastic differential equation with jumps 1 . In contrast, the papers devoted to numerical integrators underlying Hamiltonian Monte Carlo or its variants focus mainly on the discretization of Hamiltonian dynamics over a fixed duration as the refreshments arrive at fixed times (see for example [7], [31]–[33]). The presented integrators have clear resemblance to \(\mathrm{[UBU]}\) integrator for discretizing kinetic Langevin SDEs proposed in [15] and investigated in [16], [34]. This means that several new numerical integrators can be designed for 1 as has been done for kinetic Langevin dynamics. Such a study will be of interest considering path-wise properties of 1 are different from kinetic Langevin SDEs.

3.2 [E, B, S] schemes↩︎

Here, we include the switching mechanism in the [E, B] schemes introduced in the previous subsection. During the S-step, the position \(x\) is kept fixed. This results in a continuous-time Markov chain (CTMC) \(\mathfrak{s}(t)\) with finite state space whose generator matrix is \(Q\) (from 2 ). Let \(\Pi(t)\) denote the row vector of probabilities of the chain at time \(t\). Then, its evolution is described by \[\begin{align} \frac{\mathop{}\!\mathrm{d}}{\mathop{}\!\mathrm{d}t}\Pi(t) = \Pi(t) Q(x), \qquad \Pi(t) = \Pi(0) \exp({t Q(x)}), \quad t \in [0,h], \label{swi95hmc9534650} \end{align}\tag{12}\] where \(\exp(\cdot)\) represents the matrix exponential. Consequently, the semigroup generated by \(\mathrm{S}\) acts on functions \(g = g(x,v,m)\) as \[(P_t^\mathrm{S} g)(x,v,m) = \sum_{j \in \mathcal{M}} (\exp({tQ(x)}))_{mj} \, g(x,v,j).\] To incorporate the switching mechanism in the [E, B] scheme, we need a random sample from \(\Pi(h)\) i.e. a sample from the probability distribution of the continuous time Markov chain \(\mathfrak{s}\) (note that \(x\) is fixed) at time \(h\). The practicability of the above solution is limited by the required computational cost and memory usage if the size of matrix \(Q\) is large. There have been several strategies proposed in the literature to simulate CTMC : (i) Gillespie method [35]–[37], commonly known as stochastic simulation algorithm (SSA) or direct jump method, and (ii) Uniformization method [38]. We can also use a higher order approximation which induces \(\mathcal{O}(h^\alpha)\), \(\alpha \geq 3\) error in each step and hence maintaining second order of convergence of resulting [E, B, S] schemes. Since we are interested in a random sample at fixed time \(h\), we will first discuss the uniformization technique to simulate S step. This is exact in distribution like the Gillespie method. The idea behind the algorithm is to transform the CTMC, where transition to next state depends on current state, to a discrete time Markov chain via random time change. Here, for the sake of completeness, we provide a brief derivation. Let \(P = I + \frac{1}{\eta}Q\) where \(\eta \geq \max_{i} \mathfrak{q}_i\). This implies \[\begin{align} \Pi(t) &= \Pi(0) \exp((\eta P - \eta I )t) = \Pi(0) \exp(-\eta I t) \exp(\eta P t) \\ & = e^{-\eta t}\Pi(0) \exp(\eta P t) = e^{-\eta t}\sum_{n=0}^{\infty} \frac{(\eta t)^{n}}{n !}\Pi(0)P^n. \end{align}\] The uniformization allows us to simulate a continuous time Markov chain with generator matrix \(Q\) by simulating discrete time Markov chain with transition matrix \(P\). If \(\mathfrak{s}\) is our original continuous time Markov chain and we denote the discrete time Markov chain with transition matrix \(P\) as \(\mathfrak{m}\), then \(\mathfrak{s}(t) = \mathfrak{m}_{N^{\eta}(t)}\) where \(N^{\eta}(t)\) is a Poisson clock with intensity \(\eta\). Let \(P_{ij}\) denote the \(ij\)-th entry in matrix \(P\). Therefore, given \(x\), the update in switching Markov chain is given by \[\mathscr{S}_h^{\mathrm{S}}(\mathbf{s}) = \mathfrak{m}_{n}, \text{ where } \begin{cases} \text{simulate } n := N^{\eta} \sim \text{Poisson}(\eta h), \\ \mathfrak{m}_i \sim P_{\mathfrak{m}_{i-1} \cdot} \text{ for } i = 1, \dots, n, & P = I + \frac{1}{\eta}Q,\; \mathfrak{m}_0 = \mathbf{s}, \end{cases}\] where \(\mathfrak{m}_i \sim P_{\mathfrak{m}_{i-1} \cdot}\) means \[\mathfrak{m}_i = \min \left\{ j \in \mathcal{M} \;\Bigg| \;\sum_{l=1}^{j} P_{\mathfrak{m}_{i-1} l} \geq u_i \right\}, \text{ with } u_i \sim \mathcal{U}(0, 1).\] Having introduced the uniformization technique above, we may now combine it with the \([\text{E,\;B}]\) schemes in order to define the class of \([\text{E,\;B,\;S}]\) schemes for switching Hamiltonian Monte Carlo dynamics 4 . In this setting, there can be six possible symmetric concatenations, namely [SEBES], [SBEBS], [BESEB], [BSESB], [ESBSE] and [EBSBE]. Except [SBEBS] and [EBSBE] all other schemes require \(M+1\) gradient evaluations for \(M\) steps. We write below two possible concatenations [SEBES] and [BESEB]. \[\label{swi95hmc95sebes} \begin{array}{cc} \text{[SEBES]} & \text{[BESEB]} \\ [6pt] \begin{align} \mathfrak{s}_1 &= \mathscr{S}_{h/2}^{\mathrm{S}}(\mathfrak{s}_0\; ; \; x), \\ (x_{1/2}, \hat{v}_{1/2}) &= \mathscr{S}_{h/2}^{\mathrm{E}}(x, v), \\ (x_{1/2}, v_{1/2}) &= \mathscr{S}_{h}^{\mathrm{B}}(x_{1/2}, \hat{v}_{1/2}\; ;\; \mathfrak{s}_1), \\ (x_{1}, v_{1}) &= \mathscr{S}_{h/2}^{\mathrm{E}}(x_{1/2}, v_{1/2}),\\ \mathfrak{s}_2 &= \mathscr{S}_{h/2}^{\mathrm{S}}(\mathfrak{s}_1\; ; \; x_1). \end{align} & \begin{align} (x, \hat{v}_{1/2}) &= \mathscr{S}_{h/2}^{\mathrm{B}}(x, v\;;\; \mathfrak{s}_0), \\ (x_{1/2}, v_{1/2}) &= \mathscr{S}_{h/2}^{\mathrm{E}}(x, \hat{v}_{1/2}), \\ \mathfrak{s}_1 &= \mathscr{S}_{h}^{\mathrm{S}}(\mathfrak{s}_0\; ; \; x_{1/2}), \\ (x_1, \hat{v}_{1}) &= \mathscr{S}_{h/2}^{\mathrm{E}}(x_{1/2}, v_{1/2}),\\ (x_1, v_1) &= \mathscr{S}_{h/2}^{\mathrm{B}}(x_1, \hat{v_1}\; ; \; \mathfrak{s}_1). \end{align} \end{array}\tag{13}\]

If \(\eta\) is large, simulation based on uniformization may have large number of virtual jumps. Instead of uniformizing the chain, we may simulate its jump times via Gillespie’s construction. Recall that \(\mathfrak{q}_\mathfrak{s}= \sum_{j\neq \mathfrak{s}}\mathfrak{q}_{\mathfrak{s}j}\) denotes the total outgoing rate from state \(\mathfrak{s}\). Starting from \(\mathfrak{m}_0 = \mathfrak{s}\), \(\hat{t}_0 = t_k\), define recursively \[\begin{align} \label{swi95hmc9534613} \tau_{1} = \mathcal{E}(\mathfrak{q}_{\mathfrak{m}_{0}}),\quad \hat{t}_{1} := \hat{t}_0 + \tau_{1}. \end{align}\tag{14}\] If \(\hat{t}_1 > t_k + h\) then the chain does not jump during the interval \([t_k,t_k +h]\), and the S step is simply \(\mathscr{S}_h = \mathfrak{m}_{0} = \mathfrak{s}\). Else the next regime is sampled as \[\begin{align} \label{swi95hmc9534614} \mathbb{P}(\mathfrak{m}_{1} = j\; | \; \mathfrak{m}_{0} = \mathfrak{s}, x) = \frac{\mathfrak{q}_{\mathfrak{s}j}(x)}{\mathfrak{q}_\mathfrak{s}(x)}\qquad j \neq \mathfrak{s}. \end{align}\tag{15}\] Equivalently, with \(u_1\sim\mathcal{U}(0,1)\), \[\begin{align} \label{swi95hmc9534615} \mathfrak{m}_{1} = \min\left\{ j\in\mathcal{M}\setminus\{\mathfrak{s}\} \;\Bigg|\; \sum_{\ell\in\mathcal{M}\setminus\{\mathfrak{s}\},\,\ell\leq j} \frac{\mathfrak{q}_{\mathfrak{s}\ell}(x)}{\mathfrak{q}_{\mathfrak{s}}(x)} \ge u_1 \right\}. \end{align}\tag{16}\] The procedure 14 16 is repeated for \(i =1,\dots, n\) where \(n\) is the largest integer such that \(\hat{t}_{n}\le t_k + h\), i.e. until the next proposed jump time exceeds \(t_k + h\). Therefore, the Gillespie update can be written as \[\begin{align} \mathscr{S}_h^{\mathrm{S}}(\mathbf{s}) = \mathfrak{m}_{n}. \end{align}\] We can simply replace the S step in the [E, B, S] schemes with Gillespie’s construction.

One simplification is possible if we take the switching rates of the form \[\begin{align} \mathfrak{q}_{ij} = c p_j, \quad i \neq j, \end{align}\] where \(c\) is some positive constant, \(p_j \geq 0\) and \(\sum_{j=1}^m p_j = 1\). As is clear, we have \[\begin{align} \mathfrak{q}_{ii} = -\sum_{j \neq i}\mathfrak{q}_{ij} = - c(1 -p_i). \end{align}\] Then the generator can be written as \(Q(x)=c\bigl(\mathbf{1}_m p(x)^{\top}-I_m\bigr)\) where \(p(x)=(p_1(x),\ldots,p_m(x))^\top\) and \(\mathbf{1}_m \in \mathbb{R}^{m}\) with each component \(1\). Since \(\mathbf{1}p(x)^\top\) is an idempotent matrix, the transition matrix is explicit, i.e. \[\begin{align} \exp(hQ(x)) = e^{-ch}I_m+\bigl(1-e^{-ch}\bigr)\mathbf{1}_m p(x)^\top . \end{align}\] In component-wise form, we have \[\begin{align} \mathbb{P}(\mathfrak{s}_1=j\mid \mathfrak{s}_0=i,x) = e^{-ch}I({\{i=j\}}) + \bigl(1-e^{-ch}\bigr)p_j(x). \end{align}\] Thus the S step can be sampled exactly without uniformization or Gillespie method.

3.2.0.1 Comparison with switching Langevin algorithm from [3].

To sample from mixture distributions, the following stochastic differential equation with state-dependent switching is considered in [3]: \[\begin{align} &\mathop{}\!\mathrm{d}X(t) = -\frac{1}{2}\nabla U(X(t); \mathfrak{s}(t))\mathop{}\!\mathrm{d}t + \mathop{}\!\mathrm{d}B(t), \;\;\mathbb{P}(\mathfrak{s}(t+\delta)=j \mid \mathfrak{s}(t)=i) = \mathfrak{q}_{ij}(X(t))\delta + o(\delta), \quad j \neq i, \end{align}\] where \(\mathfrak{s}(t)\) is the conditionally Markov switching state and \(B(t)\) is a standard Brownian motion, and the switching rates \(\mathfrak{q}_{ij}(x)\) are from 2 . Switching Langevin algorithm proposed in [3] is based on an explicit Euler scheme with constant time step \(h\): \[\begin{align} \label{swi95lang95algo951} X_{k+1} &= X_k - \frac{h}{2} \nabla U(X_k; \mathfrak{s}_k) + \sqrt{h} \xi_{k+1}, \quad \xi_{k+1} \sim \mathcal{N}(0, I_d), \end{align}\tag{17}\] where state-dependent regime-switching chain is approximated as \[\begin{align} \label{swi95lang95algo952} \mathfrak{s}_{k+1} &= \begin{cases} j \neq \mathfrak{s}_k & \text{with prob } \mathfrak{q}_{\mathfrak{s}_k j}(x_k)h, \\ \mathfrak{s}_k & \text{with prob } 1 - \mathfrak{q}_{\mathfrak{s}_k}(x_k)h. \end{cases} \end{align}\tag{18}\] In comparison with 18 , we use uniformization technique or Gillespie’s algorithm to simulate regime changes. Moreover, thanks to splitting technique relying on exact simulation of switching chain, error incurred in \(\mathrm{[E, B, S]}\) scheme is \(\mathcal{O}(h^2)\) (with one force evaluation in most of the concatenations) as compared to \(\mathcal{O}(h)\) in 17 18 .

4 Main Results↩︎

Our first main result establishes geometric ergodicity of numerical scheme \(\mathrm{[SEBES]}\) proposed in previous section for 4 . We ensure that there exists a constant \(h_0 \in (0,1)\) such that the geometric ergodicity result is uniform in \(h \in (0,h_0)\). This is along the lines of [39], [40] in contrast to geometric ergodicity obtained in [12], [41] for a fixed discretization parameter.

We will use the following Lyapunov function in our analysis [24], [41]: \[\begin{align} \label{eqn95swi95hmc95lyapunov} W (x,v) = 1 + \mathcal{V}(x,v) = 1 + a|x|^2 + b|v|^2 + c( x \cdot v) \text{ with } ab > c^2/4. \end{align}\tag{19}\] The further conditions needed on \(a, b\) and \(c\) are discussed in the proof of Lemma 2.

The theorem is proved once we prove the following two lemmas.

The next result is establishing the order of convergence of the proposed schemes when employed to compute ergodic averages. To this end, we need some notation. For functions on \(\mathbb{R}^{2d}\times\mathcal{M}\), define \[\begin{align} L^\infty_\mathcal{V} := \left\{ \varphi : \sup_{(x,v,i)} \frac{|\varphi(x,v,i)|}{1+\mathcal{V}(x,v)} < \infty \right\},\; L^{1,\mathcal{V}} := \left\{ g : \sum_{i\in\mathcal{M}} \int_{\mathbb{R}^{2d}} |g(x,v,i)|(1+\mathcal{V}(x,v))\,\mathop{}\!\mathrm{d}x\,\mathop{}\!\mathrm{d}v < \infty \right\}. \end{align}\] We denote the associated norms by \(|\varphi|_{L^\infty_\mathcal{V}}\) and \(|g|_{L^{1,\mathcal{V}}}\), respectively. Let \(f\in C_c^\infty(\mathbb{R}^{2d} \times \mathcal{M})\). For fixed \((x,v,i)\), the operators \(\mathrm{E}\) and \(\mathrm{B}\), corresponding to \(\mathscr{S}^{\mathrm{E}}\) and \(\mathscr{S}^{\mathrm{B}}\) steps, act on the \((x,v)\) variables only as \[\begin{align} (\mathrm{B}f)(x,v,i) &= -(\nabla U(x;i)\cdot \nabla_v f(x,v,i)),\\ (\mathrm{E}f)(x,v,i) &= (v\cdot \nabla_x f(x,v,i)) +\lambda\big(\mathcal{R} f(x,v,i)-f(x,v,i)\big), \end{align}\] where the refreshment operator \(\mathcal{R}\) is \((\mathcal{R} f)(x,v,i) = \int_{\mathbb{R}^d} f\big(x,\,v\cos\phi+\sin\phi\,z,\,i\big) G(z)\mathop{}\!\mathrm{d}z\) with \(G\) being the standard Gaussian density on \(\mathbb{R}^d\). Likewise, the operator corresponding to switching chain is \[(\mathrm{S}f)(x,v,i) = -\mathfrak{q}_i(x)f(x,v,i)+\sum_{j\neq i} \mathfrak{q}_{ij}(x) f(x,v,j).\] Since \(\mathrm{E}\) contains a nonlocal refreshment term, we interpret \([\mathrm{E},\mathrm{B}]\) as the commutator of linear operators rather than as a Lie bracket of vector fields. Namely, on the common test space \(C_c^\infty(\mathbb{R}^{2d})\), \[[\mathrm{E},\mathrm{B}]f := \mathrm{E}(\mathrm{B}f)-\mathrm{B}(\mathrm{E}f).\] More generally, all nested brackets appearing below are understood in this operator-theoretic sense.

We can write the corresponding Markov semigroups as \[\begin{align} (P^{\mathrm{E}}_tf)(x,v) = \sum_{n=0}^\infty e^{-\lambda t} \frac{(\lambda t)^n}{n!} \mathbb{E} \big[ f\big( \mathscr{S}_t^{\mathrm{E},(n)}(x, v) \big) \mid N_t = n \big], \quad (P^{\mathrm B}_t f)(x,v) = f\big(x, v-t\nabla U(x)\big), \end{align}\] and \[(P^{\mathrm{S}}_t f)(x,v,i) = \sum_{j\in\mathcal{M}} (\exp{(tQ(x))})_{ij} f(x,v,j).\]

To quantify the discretization bias in estimating \(\mu(\varphi)\), where \(\varphi\in L_\mathcal{V}^\infty\), we use the following discrete Poisson equation associated with the geometrically ergodic Markov chain on \(\mathbb{R}^{2d}\times\mathcal{M}\) whose one-step transition operator we denote with \(\mathcal{P}_h\): \[\begin{align} \label{dis95Poisson95equation95SEBES} \Big(\frac{\mathcal{P}_h-I}{h}\Big)\Psi = \varphi-\mu_h(\varphi). \end{align}\tag{20}\] This equation admits a solution [42], i.e. \(\Psi = -h\sum_{k=0}^\infty (\mathcal{P}_h)^k\big(\varphi-\mu_h(\varphi)\big)\), and due to uniform in \(h\) geometric ergodicity result (i.e. Theorem 1), we also have the following bound: \[\begin{align} \label{swi95hmc95poiss95bound} |\Psi(x,v,i)|\le \frac{hK}{1- e^{-\gamma h}}\,|\varphi|_{L_\mathcal{V}^\infty}\,(1+\mathcal{V}(x,v)), \end{align}\tag{21}\] where \(K\) and \(\gamma\) are independent of \(h\). We impose the following assumption.

We say that \(\mathrm{E}^*\) and \(\mathrm{B}^*\) are the formal adjoints of \(\mathrm{E}\) and \(\mathrm{B}\) with respect to the Lebesgue pairing on \(\mathbb{R}^{2d}\times\mathcal{M}\) if, for every \(f\in C_c^\infty(\mathbb{R}^{2d}\times\mathcal{M})\) and every smooth function \(g\) for which the pairings below are absolutely integrable and integration by parts produces no boundary terms, we have \[\begin{align} \sum_{i\in\mathcal{M}}\int_{\mathbb{R}^{2d}} (\mathrm{E}f)(x,v,i)\,g(x,v,i)\,\mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}v &= \sum_{i\in\mathcal{M}}\int_{\mathbb{R}^{2d}} f(x,v,i)\,(\mathrm{E}^* g)(x,v,i)\,\mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}v, \\ \sum_{i\in\mathcal{M}}\int_{\mathbb{R}^{2d}} (\mathrm{B}f)(x,v,i)\,g(x,v,i)\,\mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}v &= \sum_{i\in\mathcal{M}}\int_{\mathbb{R}^{2d}} f(x,v,i)\,(\mathrm{B}^* g)(x,v,i)\,\mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}v. \end{align}\] In particular, this applies to the functions \(g\) arising later in the proof of Theorem 2, such as \(\rho\), \(\mathrm{E}^*\rho\), \(\mathrm{B}^*\rho\), and their finite compositions, since by Assumption 3 these functions are smooth and have sufficiently fast decay.

For two probability measures \(\nu_1,\nu_2\) on \(\mathbb{R}^{2d}\times\mathcal{M}\), define the weighted total variation distance associated with \(\mathcal{V}\) by \[\begin{align} \|\nu_1-\nu_2\|_{\mathcal{V}} := \sup_{\|\varphi\|_{L^\infty_{\mathcal{V}}}\leq 1} \Big| \int \varphi\,\mathop{}\!\mathrm{d}\nu_1-\int \varphi\,\mathop{}\!\mathrm{d}\nu_2 \Big|. \end{align}\] We also define the usual total variation distance and the bounded-Lipschitz distance by \[\begin{align} \|\nu_1-\nu_2\|_{\mathrm{TV}} := \sup_{\|\varphi\|_\infty\leq 1} \big| \textstyle \int \varphi\,\mathop{}\!\mathrm{d}\nu_1-\int \varphi\,\mathop{}\!\mathrm{d}\nu_2 \big|\; \text{ and }\; d_{\mathrm{BL}}(\nu_1,\nu_2) := \sup\limits_{\|\varphi\|_\infty\leq 1,\; \mathrm{Lip}(\varphi)\leq 1} \left| \int \varphi\,\mathop{}\!\mathrm{d}\nu_1-\int \varphi\,\mathop{}\!\mathrm{d}\nu_2 \right|, \end{align}\] respectively.

The proof of above corollary follows from Theorem 2 since we have uniform in \(h\) bound on the solution of discrete Poisson equation.

One approach to obtain weak error expansions/bounds for ergodic estimates is to use the solution to a parabolic PDE: \[\frac{\partial u}{\partial t} = \mathcal{L}u, \quad u(0, z) = \varphi(z),\] where \(\mathcal{L}\) is the operator associated with the underlying Markov process. This technique, based on backward Kolmogorov equation, has been employed in [3], [13], [43]–[45] for hypoelliptic and elliptic stochastic differential equations assuming \(U \in C^{\infty}\). This approach requires that the derivatives of the semigroup corresponding to SDE decay exponentially in time. There is one other approach based on Poisson partial differential equation. It has been employed in [46]–[49]. To get weak error bounds using this methodology, we require Schauder estimates on the solution of Poisson equation.

It is clear that both techniques, which have been used for weak error analysis of numerical schemes to compute ergodic averages, require certain types of bounds on the solutions of parabolic or elliptic PDEs. The paper [3] investigating switching Langevin algorithm conveniently relies on [44] for such bounds. In our case, we have a first order partial integro-differential equation (PIDE). We would need such bounds if we were to follow the above mentioned approaches for error analysis. Since we do not have elliptic or hypoelliptic diffusion in the dynamics, we are not aware of any such results on the derivatives of the solution of corresponding first order PIDE. Therefore, we devise an approach based on discrete Poisson equation associated with the numerical scheme itself to obtain the error bounds. Discrete Poisson equation has been used as a tool for studying additive functionals, especially when proving limit theorems (see [50], [51] and [52]), and for convergence of stochastic approximation and reinforcement learning algorithms [53]–[55]. One important consequence of the presented approach is that we get convergence in weighted total variation norm and \(1\)-Wasserstein distance (see Corollary 1) which is not straightforward in PDEs based approaches.

For the discussion that follows let \(\mathcal{L}\) denote the generator of kinetic Langevin dynamics. The error analysis for the numerical schemes of kinetic Langevin equation in [40] also uses the operator \((\mathcal{P}_{h}- I)/h\) where \(\mathcal{P}_h\) is the one-step Markov operator of the numerical scheme. However, the analysis still relies on two Poisson PDEs : one involving generator \(\mathcal{L}\) and other with its \(L^2\) adjoint \(\mathcal{L}^*\). In [56], discrete Poisson equation with operator \((e^{\mathcal{L}h} - I)/h\) has been employed to obtain non-asymptotic bounds for numerical schemes approximating kinetic Langevin equation but assuming strongly convex \(U\). In contrast, we use discrete Poisson equation associated with numerical scheme itself i.e. we have \(\mathcal{P}_h^{\mathrm{EBE}} = e^{h\mathrm{E}/2}e^{h\mathrm{B}}e^{h\mathrm{E}/2}\) and similarly for \(\mathcal{P}_h^{\mathrm{SEBES}}\) to analyze the error. We leave it for future work to obtain non-asymptotic bounds on error incurred in invariant measure approximation using numerical integrators for randomized HMC 1 , switching HMC 4 or kinetic Langevin equation with strongly convex potential \(U\) (or \(U_i\) in case with switching) using the discrete Poisson equation based approach proposed in this paper.

5 Proofs↩︎

5.1 Proof of Theorem 1↩︎

By applying Harris theorem [59] on fixed time skeleton kernel \(\mathcal{Q} = (\mathcal{P}_h)^{K_\varepsilon}\), we obtain the desired result.

5.2 Proof of Theorem 2↩︎

Consider a measurable function \(\varphi\) defined on \(\mathbb{R}^{2d} \times \mathcal{M}\) and bounded by \(1 +\mathcal{V}\). Integrating 20 with respect to \(\rho\) yields \[\begin{align} \sum_{i\in\mathcal{M}} \int_{\mathbb{R}^{2d}} \left(\frac{\mathcal{P}_h-I}{h}\right)\Psi(x, v,i) \rho(x,v,i) \mathop{}\!\mathrm{d}x \mathop{}\!\mathrm{d}v = \mu(\varphi) - \mu_h(\varphi). \end{align}\] Because \(\mathcal{P}_h\) is a positive operator and the Poisson solution satisfies \[\begin{align} |\Psi(x,v,i)|\leq |\Psi|_{L^\infty_\mathcal{V}}(1+\mathcal{V}(x,v)), \end{align}\] we may use the Lyapunov drift condition \(\mathcal{P}_h(1+\mathcal{V})\le c_1 \mathcal{V}+c_2\) to bound the pairing against the invariant density on \(\mathbb{R}^{2d}\times\mathcal{M}\): \[\begin{align} \sum_{i\in\mathcal{M}}\int_{\mathbb{R}^{2d}} |\mathcal{P}_h\Psi(x,v,i)|\,\rho(x,v,i)\mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}v \nonumber &\leq |\Psi|_{L^\infty_\mathcal{V}} \sum_{i\in\mathcal{M}}\int_{\mathbb{R}^{2d}} \big(\mathcal{P}_h(1+\mathcal{V})\big)(x,v,i)\,\rho(x,v,i)\mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}v \nonumber \\ &\leq |\Psi|_{L^\infty_\mathcal{V}} \sum_{i\in\mathcal{M}}\int_{\mathbb{R}^{2d}} \big(c_1\mathcal{V}(x,v)+c_2\big)\,\rho(x,v,i)\mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}v <\infty, \label{swi95eqn95new9554655} \end{align}\tag{37}\] where the finiteness follows from the exponential decay of \(\rho\). The adjoint \(\mathcal{P}^*_{h}\) is defined via duality \[\begin{align} \sum_{i\in\mathcal{M}}\int_{\mathbb{R}^{2d}} f(x,v,i) \mathcal{P}_h^*\rho(x,v,i) \mathop{}\!\mathrm{d}x \mathop{}\!\mathrm{d}v = \sum_{i\in\mathcal{M}}\int_{\mathbb{R}^{2d}} \mathcal{P}_h f(x,v,i) \rho(x,v,i) \mathop{}\!\mathrm{d}x \mathop{}\!\mathrm{d}v, \end{align}\] where \(f\) belongs to a class of bounded and measurable functions. Since each substep of \(\mathrm{SEBES}\) maps absolutely continuous measures to absolutely continuous measures, \(\mathcal{P}^*\rho\) can be identified as density. Because of 37 , Fubini’s theorem can be used to extend duality for \(f = \Psi\), and therefore, we get \(\mathcal{P}^{*}_h \rho \in L^{1, \mathcal{V}}\) and \[\begin{align} \label{swi95eq95adjoint} \sum_{i\in\mathcal{M}} \int_{\mathbb{R}^{2d}} \Psi(x,v,i) \left(\frac{\mathcal{P}_h^*-I}{h}\right)\rho(x,v,i) \mathop{}\!\mathrm{d}x \mathop{}\!\mathrm{d}v =\mu(\varphi) - \mu_h(\varphi). \end{align}\tag{38}\] For the \(\mathrm{B}\)-step, Taylor’s theorem along the explicit flow yields \[P^{\mathrm{B}}_sf = f+s\mathrm{B}f+\frac{s^2}{2}\mathrm{B}^2f+\frac{s^3}{6}\mathrm{B}^3f +\int_0^s\frac{(s-r)^3}{3!}P_r^\mathrm{B}\mathrm{B}^4f \mathop{}\!\mathrm{d}r.\] The function \(\mathrm{E}^j f\), \(0<j \leq 4\) is smooth and has at most polynomial growth. It is not difficult to show the moment bounds for \(\mathrm{SEBES}\) under Assumptions 1-2 following the similar set of arguments as used in establishing Lyapunov drift condition and therefore, we have \[\sup_{i\in \mathcal{M}}\mathbb{E}\int_0^h |\mathrm{E}^{j+1}f(X_s,V_s, i)|\mathop{}\!\mathrm{d}s<\infty, \quad 0<j\leq 4,\] therefore, we get \[\begin{align} P^{\mathrm{E}}_sf = f+ s \mathrm{E}f+\frac{s^2}{2}\mathrm{E}^2f+\frac{s^3}{6}\mathrm{E}^3f + \int_0^{s}\frac{(s-r)^3}{3!} P_r^\mathrm{E} \mathrm{E}^4 f\mathop{}\!\mathrm{d}r,\quad 0\leq s \leq h. \end{align}\] Denote \[\begin{align} \mathrm{R}_\mathrm{E} f = \int_0^{h/2}\frac{(h/2-s)^3}{3!}\,P_s^\mathrm{E} \mathrm{E}^4 f\,\mathop{}\!\mathrm{d}s, \qquad \mathrm{R}_\mathrm{B} f =\int_0^{h}\frac{(h-s)^3}{3!}\,P_s^\mathrm{B} \mathrm{B}^4 f\,\mathop{}\!\mathrm{d}s. \end{align}\] For a given \((x,v) \in \mathbb{R}^d\), we treat \(i \to f(x,v, i)\) as \(\mathbb{R}^m\) valued function. We apply Taylor’s theorem to \(P_s^{\mathrm{S}} := e^{s Q}\) to get \[P_{s}^\mathrm{S}f = f+ s\mathrm{S}f+\frac{s^2}{2}\mathrm{S}^2 f+\frac{s^3}{6}\mathrm{S}^3 f+\mathrm{R}_\mathrm{S}f,\] with \(\mathrm{S}:= Q(x)\) and \(\mathrm{R}_\mathrm{S}f = \int_0^{s}\frac{(s-r)^3}{3!}P_r^\mathrm{S}\mathrm{S}^4 f \mathop{}\!\mathrm{d}r\). Consider the symmetric splitting \[\mathcal{P}_h := P_{h/2}^\mathrm{S}P_{h/2}^\mathrm{E}P_h^\mathrm{B}P_{h/2}^\mathrm{E}P_{h/2}^\mathrm{S}.\] Denoting \(P_{s}^\mathrm{F}= \mathrm{T}_\mathrm{F}+ \mathrm{R}_\mathrm{F}\) where \(\mathrm{T}_{\mathrm{F}} := \mathrm{I} + s\mathrm{F} + \frac{s^2}{2}\mathrm{F}^2 + \frac{s^3}{6}\mathrm{F}^3\) with \(\mathrm{F} \in \{\mathrm{E}, \mathrm{B}, \mathrm{S}\}\) and \(s = h\) if \(\mathrm{F} = \mathrm{B}\) and \(s = h/2\) otherwise. Then \[\begin{align} \mathcal{P}_h = (\mathrm{T}_\mathrm{S}+ \mathrm{R}_\mathrm{S})(\mathrm{T}_\mathrm{E}+ \mathrm{R}_\mathrm{E})(\mathrm{T}_\mathrm{B}+ \mathrm{R}_\mathrm{B})(\mathrm{T}_\mathrm{E}+ \mathrm{R}_\mathrm{E})(\mathrm{T}_\mathrm{S}+ \mathrm{R}_\mathrm{S}). \end{align}\] Expanding the product and collecting the terms with no remainder factor, we obtain \[\mathrm{T}_\mathrm{S}\mathrm{T}_\mathrm{E}\mathrm{T}_\mathrm{B}\mathrm{T}_\mathrm{E}\mathrm{T}_\mathrm{S} = I+h\mathcal{A} +\frac{h^2}{2}\mathcal{A}^2 +h^3\Bigl(\frac{1}{6}\mathcal{A}^3+\mathcal{S}\Bigr)+ \mathrm{R}_1,\] where \(\mathcal{A}=\mathrm{E}+\mathrm{B}+ \mathrm{S}\) and \(\mathrm{R}_1\) contains all terms of order at least \(h^4\) arising from \(\mathrm{T}_\mathrm{S}\mathrm{T}_\mathrm{E}\mathrm{T}_\mathrm{B}\mathrm{T}_\mathrm{E}\mathrm{T}_\mathrm{S}\). All the terms containing at least one of \(\mathrm{R}_\mathrm{S},\mathrm{R}_\mathrm{E},\mathrm{R}_\mathrm{B}\) are grouped into \(\mathrm{R}_2\). Consequently, we have the following expansion for \(f \in C_c^{\infty}(\mathbb{R}^{2d}\times \mathcal{M})\): \[\begin{align} \label{eq:SEBES-expansion} \mathcal{P}_h f = f + h\mathcal{A} f + \frac{h^2}{2}\mathcal{A}^2 f +h^3\Big(\frac{1}{6}\mathcal{A}^{3}+ \mathcal{S}\Big)f +\mathrm{R}_1 f+ \mathrm{R}_2 f. \end{align}\tag{39}\] where \(\mathcal{S} = \frac{1}{2}[\mathrm{B},[\mathrm{B},\mathrm{E}]] -\frac{1}{24}[\mathrm{E},[\mathrm{E},\mathrm{B}]] +\frac{1}{12}[\mathrm{B}+\mathrm{E},[\mathrm{B}+\mathrm{E},\mathrm{S}]] -\frac{1}{24}[\mathrm{S},[\mathrm{S},\mathrm{B}+\mathrm{E}]]\). We now pair the above expansion against the Gibbs density \(\rho\) in order to derive the corresponding adjoint expansion. Since \(f\in C_c^\infty(\mathbb{R}^{2d}\times\mathcal{M})\), all terms produced by the operators under consideration are smooth and have at most polynomial growth. Together with the fact that \(\rho\) is smooth and its derivatives decay sufficiently fast, this implies that all duality pairings and remainder integrals below are absolutely integrable. Consequently, Fubini’s theorem may be applied to interchange the order of integration, and integration by parts yields no boundary contributions (due to sufficiently fast decay of \(\rho\) and its derivatives, see Assumption 3). The weighted spaces \(L^{1,\mathcal{V}}\) and \(L_\mathcal{V}^\infty\) are used only for the subsequent remainder estimates and final duality argument. With these preliminaries in place, the adjoints can be obtained. For terms like \(\mathrm{T}_\mathrm{B}\) it is straightforward to obtain its adjoint which we write below \[\begin{align} \mathrm{T}_\mathrm{B}^* = I+h\mathrm{B}^*+\frac{h^2}{2}(\mathrm{B}^*)^2+\frac{h^3}{6}(\mathrm{B}^*)^3. \end{align}\] Likewise, we can obtain adjoints for \(\mathrm{T}_\mathrm{E}\), \(\mathrm{T}_\mathrm{S}\) and \(\mathrm{R}_1\) via repeated application of integration by parts and Fubini’s theorem. Likewise, for reminder terms, for example: \[\begin{align} \mathrm{R}_\mathrm{S}^*\rho = \int_0^{h/2}\frac{(h/2-s)^3}{3!}(\mathrm{S}^*)^4 (P_s^\mathrm{S})^*\rho \mathop{}\!\mathrm{d}s. \end{align}\] It is more involved to get adjoint of one of the terms appearing in \(\mathrm{R}_2\), for example \[\begin{align} (\mathrm{R}_\mathrm{S}\mathrm{R}_\mathrm{E}\mathrm{R}_\mathrm{B}\mathrm{T}_\mathrm{E}\mathrm{T}_\mathrm{S})f = \int_0^{h/2}\int_0^{h/2}\int_0^h q(s,r,u) P_s^\mathrm{S}\mathrm{S}^4 P_r^\mathrm{E}\mathrm{E}^4 P_u^\mathrm{B}\mathrm{B}^4 \mathrm{T}_\mathrm{E}\mathrm{T}_\mathrm{S}f \mathop{}\!\mathrm{d}u \mathop{}\!\mathrm{d}r \mathop{}\!\mathrm{d}s, \end{align}\] where \(q(s , r, u) =(h/2-s)^3(h/2-r)^3(h-u)^3/(3!)^3\). Hence, by reversing the operator order, its adjoint is \[\begin{align} (\mathrm{R}_\mathrm{S}\mathrm{R}_\mathrm{E}\mathrm{R}_\mathrm{B}\mathrm{T}_\mathrm{E}\mathrm{T}_\mathrm{S})^*\rho &= \int_0^{h/2}\int_0^h\int_0^{h/2} q(u,r,s)\mathrm{T}_\mathrm{S}^* \mathrm{T}_\mathrm{E}^* (\mathrm{B}^*)^4 (P_r^\mathrm{B})^* (\mathrm{E}^*)^4 (P_s^\mathrm{E})^* (\mathrm{S}^*)^4 (P_u^\mathrm{S})^*\rho \mathop{}\!\mathrm{d}s\mathop{}\!\mathrm{d}r\mathop{}\!\mathrm{d}u. \end{align}\] All remaining adjoint terms are obtained in the same manner. Therefore, pairing 39 against \(\rho\), we obtain \[\begin{align} \sum_{i\in\mathcal{M}}&\int_{\mathbb{R}^{2d}} (\mathcal{P}_h f)(x,v,i)\,\rho(x,v,i)\mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}v = \sum_{i\in\mathcal{M}}\int_{\mathbb{R}^{2d}} f(x,v,i)\rho(x,v,i)\mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}v \nonumber\\ & + h\sum_{i\in\mathcal{M}}\int_{\mathbb{R}^{2d}} f_i\mathcal{A}^*\rho_i\mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}v +\frac{h^2}{2}\sum_{i\in\mathcal{M}}\int_{\mathbb{R}^{2d}} f_i(\mathcal{A}^*)^2\rho_i\mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}v +h^3\sum_{i\in\mathcal{M}}\int_{\mathbb{R}^{2d}} f_i\Big(\frac{1}{6}(\mathcal{A}^*)^3+\mathcal{S}^*\Big)\rho_i\mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}v \nonumber\\ & +\sum_{i\in\mathcal{M}}\int_{\mathbb{R}^{2d}} f_i\mathrm{R}_1^*\rho_i\mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}v +\sum_{i\in\mathcal{M}}\int_{\mathbb{R}^{2d}} f_i\mathrm{R}_2^*\rho_i\mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}v, \label{swi95eq95SEBES95adjoint95pairing} \end{align}\tag{40}\] where \(\mathcal{A}^*=\mathrm{S}^*+\mathrm{E}^*+\mathrm{B}^*\), \(\mathcal{S}^* = \frac{1}{12}[\mathrm{B}^*,[\mathrm{B}^*,\mathrm{E}^*]] -\frac{1}{24}[\mathrm{E}^*,[\mathrm{E}^*,\mathrm{B}^*]] +\frac{1}{12}[\mathrm{B}^*+\mathrm{E}^*,[\mathrm{B}^*+\mathrm{E}^*,\mathrm{S}^*]] -\frac{1}{24}[\mathrm{S}^*,[\mathrm{S}^*,\mathrm{B}^*+\mathrm{E}^*]]\) and, for convenience, we have denoted \(f_i := f(x,v, i)\) and \(\rho_i := \rho(x,v,i)\). Thanks to the invariance of Gibbs measure for exact dynamics \((\mathcal{A}^{*})^j \rho = 0\) under condition 6 , 40 becomes \[\begin{align} \sum_{i\in\mathcal{M}}&\int_{\mathbb{R}^{2d}} (\mathcal{P}_h f)(x,v,i)\,\rho(x,v,i)\mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}v = \sum_{i\in\mathcal{M}}\int_{\mathbb{R}^{2d}} f(x,v, i)\rho(x,v, i)\mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}v \nonumber\\ & +h^3\sum_{i\in\mathcal{M}}\int_{\mathbb{R}^{2d}} f_i\mathcal{S}^*\rho_i\mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}v +\sum_{i\in\mathcal{M}}\int_{\mathbb{R}^{2d}} f_i\mathrm{R}_1^*\rho_i\mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}v +\sum_{i\in\mathcal{M}}\int_{\mathbb{R}^{2d}} f_i\mathrm{R}_2^*\rho_i\mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}v. \label{swi95eqn9554662} \end{align}\tag{41}\] Since the density \(\rho\) has sufficiently fast decaying derivatives, applying any finite combination of the adjoint operators \(\mathrm{S}^*,\mathrm{E}^*,\mathrm{B}^*\) produces functions which retain this decay. Since \(\mathcal{V}\) is quadratic, these functions belong to \(L^{1,\mathcal{V}}\). By passing the weighted norm inside the exact time-integral representations of the remainder terms \(\mathrm{R}_i^* \rho\), the bound below follows: \[|\mathrm{R}_i^*\rho|_{L^{1,\mathcal{V}}}\le Ch^4, \qquad i=1,2.\] And, similarly \(\mathcal{S}^*\rho\in L^{1,\mathcal{V}}\). To construct a pairing with the Poisson solution \(\Psi\), we approximate \(\Psi\) with compactly supported functions. To this end, let \(\chi_n\in C_c^\infty(\mathbb{R}^{2d})\) satisfy \(0\le \chi_n\le 1\) and \(\chi_n\to 1\) point-wise as \(n \to \infty\), and define \[g_n(x,v,i):=\chi_n(x,v) \Psi(x,v,i).\] As our Lyapunov function is locally integrable and since \(|\Psi|\le C(1+\mathcal{V})\) for all \(i\in \mathcal{M}\), we have \(g_n\in L^1(\mathbb{R}^{2d}\times \mathcal{M})\). Let \(\eta\in C_c^\infty(\mathbb{R}^{2d})\) be a standard mollifier with unit mass, and define \(\eta_\delta(z):=\delta^{-2d}\eta(z/\delta)\). For each \(n\), choose \(\delta_n>0\) such that for every \(i \in \mathcal{M}\) (see [60]), we have \[\label{swi95eqn9554664} |\eta_{\delta_n}*g_n-g_n|_{L^1(\mathbb{R}^{2d})}<2^{-n},\tag{42}\] Define \(f_n(x,v,i):=(\eta_{\delta_n}*g_n(\cdot, \cdot, i))(x,v)\in C_c^\infty(\mathbb{R}^{2d}\times \mathcal{M}).\) Then, due to Tonelli’s theorem, 42 and finiteness of \(\mathcal{M}\), we get \[\begin{align} \int_{\mathbb{R}^{2d}}\sum_{n=1}^\infty \sum_{i \in \mathcal{M}} |f_n-g_n| \mathop{}\!\mathrm{d}y = \sum_{n=1}^{\infty}\sum_{i \in \mathcal{M}} |f_n - g_n|_{L^1(\mathbb{R}^{2d})} < \infty, \end{align}\] hence, \(f_n(x,v,i)-g_n(x,v,i)\to0\) for a.e. \((x,v) \in \mathbb{R}^{2d}\) and for every \(i \in \mathcal{M}\). Since \(g_n=\chi_n\Psi\to\Psi\) pointwise, we have established \[\label{eq95SEBES95fn95to95Psi} f_n(x,v,i)\to \Psi(x,v,i) \qquad\text{for a.e. } (x,v) \in \mathbb{R}^{2d} \text{ and every } i\in \mathcal{M}.\tag{43}\] We next show that \((f_n)\) is uniformly bounded by \(1+\mathcal{V}\). Since \(\Psi\in L^\infty_\mathcal{V}\), \(|\Psi(y,i)|\leq |\Psi|_{L^\infty_\mathcal{V}}(1+\mathcal{V}(y))\). Hence, for \(z=(x,v)\), \[\begin{align} |f_n(z, i)| &= \left| \int_{\mathbb{R}^{2d}}\eta_{\delta_n}(z-y)\chi_n(y)\Psi(y, i)\mathop{}\!\mathrm{d}y \right| \leq |\Psi|_{L^\infty_\mathcal{V}} \int_{\mathbb{R}^{2d}}\eta_{\delta_n}(z-y)(1+\mathcal{V}(y))\mathop{}\!\mathrm{d}y. \end{align}\] Since \(\eta_{\delta_n}\) is compactly supported and \(\delta_n\leq 1\), whenever \(\eta_{\delta_n}(z-y)\neq 0\) one has \(|z-y|\le1\). For quadratic \(\mathcal{V}\), it is not difficult to show that there exists \(C>0\) such that \(1+\mathcal{V}(y)\le C(1+\mathcal{V}(z))\) whenever \(|z-y|\leq 1\). Therefore, \(|f_n(z,i)|\le C |\Psi|_{L^\infty_\mathcal{V}}(1+\mathcal{V}(z))\). Thus for every \(\varrho\in L^{1,\mathcal{V}}\), \[|f_n\varrho| \leq C |\Psi|_{L^\infty_\mathcal{V}}(1+\mathcal{V})|\varrho| \in L^1(\mathbb{R}^{2d}),\] Therefore, using 43 and dominated convergence theorem, we obtain for every \(\varrho \in L^{1, \mathcal{V}}\) \[\label{eq:SEBES-pairing-limit} \sum_{i \in \mathcal{M}}\int_{\mathbb{R}^{2d}} f_n (x,v,i) \varrho(x,v,i)\mathop{}\!\mathrm{d}x \mathop{}\!\mathrm{d}v \longrightarrow \sum_{i \in \mathcal{M}} \int_{\mathbb{R}^{2d}} \Psi(x,v,i) \varrho(x,v,i) \mathop{}\!\mathrm{d}x \mathop{}\!\mathrm{d}v,\tag{44}\] as \(n \to \infty\). As we already discussed that functions \(\varrho=\frac{P_h^*-I}{h}\rho_i\), \(\varrho=\mathcal{S}^*\rho_i\) or \(\varrho=R_j^*\rho_i\), \(j =1, 2\) belong to \(L^{1,\mathcal{V}}\), therefore applying 44 , we may pass to the limit in the expansion 41 and conclude that \[\begin{align} \sum_{i \in \mathcal{M}}\int_{\mathbb{R}^{2d}} \Psi(x,v,i) \left(\frac{\mathcal{P}_h^*-I}{h}\right)\rho(x,v,i)\mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}v &= h^2\sum_{i \in \mathcal{M}}\int_{\mathbb{R}^{2d}} \Psi(x,v,i)\mathcal{S}^*\rho(x,v,i)\mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}v \nonumber \\ & +\frac{1}{h}\sum_{i \in \mathcal{M}} \sum_{j=1,2}\int_{\mathbb{R}^{2d}} \Psi(x,v,i) \mathrm{R}_j^*\rho(x,v,i)\mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}v. \label{eq95SEBES95strong95pairing} \end{align}\tag{45}\] Finally, since \(|R_j^*\rho|_{L^{1,\mathcal{V}}}\le Ch^4\), \(j=1,2\) and \(\Psi\in L^\infty_\mathcal{V}\), then using 38 , we get \[\mu(\varphi) - \mu_h(\varphi) = h^2\sum_{i\in \mathcal{M}}\int_{\mathbb{R}^{2d}} \Psi(x,v,i) \mathcal{S}^*\rho (x,v,i)\mathop{}\!\mathrm{d}x\mathop{}\!\mathrm{d}v + \mathcal{O}(h^3).\] This completes the proof.

6 Numerical illustration↩︎

In this section, we present a numerical experiment validating the second order convergence of SEBES scheme. With \(\hat{\varphi}_{N}\), we denote the Monte Carlo approximation of \(\mathbb{E}\varphi(X_{N})\), i.e. \(\hat{\varphi}_{N} = \frac{1}{M}\sum_{i=1}^{M} \varphi(X^{(i)}_{N})\) where \(X^{(i)}_{N}\) represents independent samples of the numerical solution. We borrow this experiment from [3] to compare with switching Langevin algorithm. We consider a mixture of two-dimensional Gaussian distributions: \[\begin{align} \rho(x) = \frac{1}{\mathcal{Z}} \left[\alpha_1 \exp\left( -\frac{1}{2}(x-m_1)^{T}\Sigma_1^{-1}(x-m_1) \right) + \alpha_2 \exp\left( -\frac{1}{2}(x-m_2)^{T}\Sigma_2^{-1}(x-m_2) \right) \right], \quad x \in \mathbb{R}^2 , \end{align}\] where \(\mathcal{Z}\) denotes the normalizing constant, \(\alpha_1 = 0.7\), \(\alpha_2 = 0.5\), \(m_1 = (1, 1)\), \(m_2 = (-2, -1)\), and \[\begin{align} \Sigma_1 &= \begin{pmatrix} 2.0 & 0.1 \\ 0.1 & 0.5 \end{pmatrix}, & \Sigma_2 &= \begin{pmatrix} 1.0 & -0.1 \\ -0.1 & 1.0 \end{pmatrix}. \end{align}\] In this experiment, we take \(T = 200\) for both SEBES and switching Langevin algorithm. For SEBES, we choose \(\lambda = 1\) and \(\phi = \pi/4\).

7 Concluding remarks↩︎

In this paper, we have considered splitting discretizations of continuous time Hamiltonian dynamics interlaced with Poisson refreshment and interacting with a regime-switching chain. We have proved second-order convergence in computing ergodic averages with respect to mixture probability distributions. Since we do not have access to derivative bounds on the solution of Kolmogorov equations or Poisson equations which are first order partial integral differential equations in our case, we have developed a new approach based on discrete Poisson equation linked with numerical integrators. This methodology allows for larger class of test functions and therefore also providing convergence in total variation and \(1\)-Wasserstein distances.

This work opens up several interesting avenues of future research:

• An important line of future work is to utilize the discrete Poisson equation based error analysis, presented in this paper, to obtain non-asymptotic bounds for numerical integrators for Langevin dynamics, randomized HMC and switching HMC under the assumption of strong convexity on potential function \(U\) (\(U_i\) in mixture distribution setting).

• Another significant step for the switching HMC is to better understand the behavior of switching mechanism, perhaps via coupling techniques, so that an optimal choice of switching rates can be made.

• There have been a lot of focus in developing and analyzing the numerical integrators for kinetic Langevin dynamics. Therefore, development of numerical integrators for randomized HMC is also an interesting topic of future work considering both are stochastic generalizations of deterministic Hamiltonian dynamics with significant differences in path-wise properties but serving the same purpose.

References↩︎