Piecewise-Deterministic Markov Chain Monte Carlo

Introduction

Markov chain Monte Carlo (MCMC) methods are the tools of choice to sample non-standard probability distributions. In high-dimensional scenarios, the celebrated Metropolis–Hastings algorithm performs usually poorly and alternative algorithms are required. Two of the most popular alternatives are slice sampling and Hamiltonian Monte Carlo (HMC) methods which have had much empirical success over recent years. More recently, continuous-time non-reversible MCMC algorithms based on Piecewise-Deterministic Markov Processes (PDMP) schemes have also appeared in the literature in applied probability , automatic control , physics , statistics and machine learning . In physics, these schemes have become quickly popular as they provide state-of-the-art performance when applied to the simulation of large scale physical models. They also show promise for statistics applications, in particular for high dimensional sparse graphical models and big data .

However, the PDMP-based schemes currently available suffer from shortcomings which limit both their applicability and performance. To ensure invariance with respect to the target distribution, one needs to be able to simulate these continuous-time processes exactly. In practice, this restricts severely the deterministic dynamics one can use: all the existing algorithms use a simple linear dynamics that does not exploit the geometry of the target. Moreover, exact simulation of the event times is problem specific and may be impossible in certain scenarios. This prevents the development of a generic software implementation of these techniques.

In this paper, we address these limitations by developing novel continuous-time and discrete-time Piecewise-Deterministic Markov Chain Monte Carlo (PD-MCMC) techniques which bring together HMC, PDMP and generalized Metropolis–Hastings.

First, we show that it is possible to develop continuous-time PD-MCMC algorithms relying on Hamiltonian dynamics. In this context, exact simulation of the resulting PDMP remains possible for an important class of target distributions. The resulting algorithms provide an alternative to elliptical slice sampling-type algorithms . We also exploit a generalized version of Metropolis–Hastings algorithm (see, e.g., ) satisfying a skewed detailed balance condition to derive novel schemes.

Second, we introduce novel discrete-time PD-MCMC algorithms. These non-reversible algorithms can be thought of as a discretized version of continuous-time PD-MCMC but preserve the target distribution as invariant distribution for all discretization steps. These schemes are not only able to exploit complex dynamics, such as approximate Hamiltonian dynamics arising from symplectic integrators, but it is also always possible to simulate the event times. Moreover some versions of these discrete-time algorithms do not even require being able to compute the gradient of the log-target. These methods enjoy the same attractive features as their continuous-time counterparts: they can leverage any representation of the target as a product of non-negative factors. Additionally they can use unbiased estimators of the log-target distribution and its gradient and still provide algorithms with the correct invariant distribution.

The rest of the paper is organised as follows. In Section 2 we review continuous-time PDMPs, provide sufficient conditions to ensure invariance of a PDMP with respect to a given target distribution, discuss existing PD-MCMC algorithms and finally introduce novel algorithms relying on Hamiltonian dynamics. In Section 3, we introduce the class of discrete-time PDMP and provide sufficient conditions to ensure invariance of a PDMP with respect to a given target distribution which parallel the ones obtained in the continuous-time scenarios. We review existing and describe novel discrete-time PD-MCMC algorithms. Section 4 is dedicated to the efficient implementation of discrete-time algorithms using subsampling and prefetching ideas while Section 5 proposes discrete-time algorithms to handle scenarios where the target is intractable but its logarithm and the logarithm of its gradient can be estimated unbiasedly. Empirical performance of some of these schemes are reviewed in Section 6. Appendix A contains all the proofs of validity of the proposed algorithms while weak convergence of a specific discrete-time scheme to a PDMP is proven in Appendix B.

Continuous-Time PDMP and PD-MCMC

an Ordinary Differential Equation (ODE) with differentiable drift $\phi:\mathcal{Z}\rightarrow\mathcal{Z}$, i.e.,

satisfying the semi-group property $\Phi_{s}\circ\Phi_{t}=\Phi_{s+t}$ and such that $t\mapsto\Phi_{t}\left(z\right)$ is càdlàg,

a Markov transition kernel $Q$ from $\mathcal{Z}$ to $\mathcal{Z}$ where the state at event time $t$ is given by $z_{t}\sim Q\left(z_{t^{-}},\cdot\right)$, $z_{t^{-}}$ being the state of the process just before the event.

Algorithm 1 describes how to simulate the path of a PDMP.

To be able to exactly simulate a PDMP, we thus need to be able to simulate from the distribution (3) and compute the flow (4). Finally we also need to be able to simulate from the transition kernel $Q$. In important scenarios, exact simulation of the event times can be performed using inversion of the integrated rate function as in or using adaptive thinning procedures as in .

We now introduce the generator associated with the PDMP. For functions in the domain of the generator, it is defined by

Under suitable regularity conditions [15, Theorem 26.14], it can be shown that this generator is given by

where $\left\langle a,b\right\rangle$ denotes the scalar product between vectors $a,b$ and $|a|^{2}=\left\langle a,a\right\rangle$. The first term on the right hand side of (6) arises from the deterministic dynamics while the second term corresponds to the jump component of the process.

2 From PDMP to PD-MCMC

Assume we are interested in sampling from a given target probability distribution on the Borel space $\left(\mathcal{Z},\mathcal{B}\left(\mathcal{Z}\right)\right)$. If we want to use a PDMP mechanism to sample this target distribution, this PDMP needs at least to admit this distribution as invariant distribution. We provide here sufficient conditions to ensure this is satisfied. If additionally the PDMP is ergodic, this will allow us to estimate consistently expectations with respect to the invariant distribution.

From now onward, the target distribution will be assumed to have a strictly positive density $\rho\left(z\right)$ with respect to the Lebesgue measure $\text{d}z$ where

Invariance with respect to $\rho$ will be satisfied if

for all functions $f$ in the domain of the generator [15, Proposition 34.7]. From (6), this means that we need

However, using integration by parts, we obtain

where $\nabla\cdot\phi\left(z\right):=\sum_{i=1}^{n}\partial_{i}\phi_{i}\left(z\right)$ is the divergence of the vector field $\phi$. Hence, a sufficient condition to ensure invariance of a PDMP with respect to $\rho$ is to have

We provide here useful sufficient conditions on $\phi,$ $\lambda$, and $Q$ to ensure $\rho$-invariance of the associated PDMP, without making any structural assumptions on these objects.

Based on these assumptions, straightforward calculations show that the following result holds.

Assume (A(A1)). Then the PDMP admits $\rho$ as invariant distribution.

2.2 Sufficient conditions for local methods

Assume that $H\left(z\right)$ can be decomposed as follows

where potentially each $H_{i}\left(z\right)$ only depends on a subset of the components of $z$. In this context, like in standard MCMC, we might be interested in using a transition kernel which is a mixture of $n$ kernels performing local updates. This can be achieved in the PDMP framework by introducing an event rate of the form

where $Q_{i}$ are Markov transition kernels. Let us write $\left[n\right]:=\left\{1,2,...,n\right\}$. To simulate the event times of the resulting PDMP, one can associate a clock to each index $i\in[n]$ and use a priority queue . When it is possible to bound $\{\lambda_{i};i\in\left[n\right]\}$ locally in time, more elaborate thinning strategies have been developed in [10, Section 3.3.2] and .

Based on these structural assumptions on $\lambda$ and $Q$, we can provide useful sufficient “local” conditions on $\phi,$ $\{\lambda_{i}:i\in\left[n\right]\}$ and $\{Q_{i}:i\in[n]\}$ to ensure that invariance of the associated PDMP with respect to $\rho$ is satisfied.

Conditions on $\phi,$ $\{\lambda_{i}:i\in\left[n\right]\}$ and $\{Q_{i}:i\in[n]\}$

There exists a $\rho$-preserving mapping $\mathcal{S}:\mathcal{Z}\to\mathcal{Z}$.

The event rates $\{\lambda_{i}:i\in\left[n\right]\}$ satisfy

For all $i\in\left[n\right]$, the transition kernel $Q_{i}$ satisfies

If the functions $\{H_{i}:i\in\left[n\right]\}$ are differentiable then Assumption A(A2).2 is satisfied for a divergence-free vector field, i.e. $\nabla\cdot\phi=0$, if for all $i\in\left[n\right]$

Assume (A(A2)). Then the PDMP admits $\rho$ as invariant distribution.

2.3 Sufficient conditions for doubly stochastic methods

In this context, we consider an event rate of the form

where $Q_{\omega}$ is a Markov transition kernel from $\mathcal{Z}$ to $\mathcal{Z}$. In Section 2.2.2, (11), (12) and (13) simply correspond to (17), (18) and (19) if we select $\mu$ as the measure such that $\mu\left(\left\{i\right\}\right)=1$ for all $i\in\Omega=[n]$. The sufficient conditions of the previous section can be directly generalized.

Conditions on $\phi,$ $\{\lambda_{\omega}:\omega\in\Omega\}$ and $\{Q_{\omega}:\omega\in\Omega\}$

There exists a $\rho$-preserving mapping $\mathcal{S}:\mathcal{Z}\to\mathcal{Z}$.

The event rates $\{\lambda_{\omega}:\omega\in\Omega\}$ satisfy

For all $\omega\in\Omega$, the transition kernel $Q_{\omega}$ satisfies

If $\mu$ is a probability measure and the derivative $\nabla H_{\omega}\left(z\right)$ is well-defined for almost all $\omega\in\Omega$ then under weak regularity conditions, it follows from (17) that $\nabla H_{\omega}\left(z\right)$ is an unbiased estimate of $\nabla H\left(z\right)$ when $\omega\sim\mu$ and Assumption A(A3).2 will be satisfied for a divergence-free field if

We will refer to this class of PD-MCMC as “doubly stochastic” in reference to doubly-stochastic Poisson processes.

Assume (A(A3)). Then the PDMP admits $\rho$ as invariant distribution.

3 Existing PD-MCMC algorithms

so the resulting flow is analytically tractable and given by

In this case, we have $\nabla\cdot\phi=0$. Additionally, all these algorithms rely on $\mathcal{S}(x,v)=(x,-v)$ which can be viewed as a time reversal, so (9) becomes

These algorithms differ in the way the event rate and the transition kernels are specified. We just give a few examples here and refer the reader to the list of references for other examples.

This algorithm proposed in exploits any additive decomposition of the potential $U$, i.e.

For $\lambda_{\textrm{ref}}>0$, it uses the event rate

The BPS algorithm has been further extended to the scenario where one has access to an unbiased estimate of $\nabla U$; see and [20, Section 4.4.2]. The validity of this algorithm can be established as an application of the results of Section 2.2.3. We are not aware of any implementation of this algorithm in scenarios where $\mu$ is not an atomic measure with finite support, in which case the algorithm is the local BPS.

3.2 Zig-Zag sampler

This algorithm proposed in uses for $\psi$ the uniform distribution on $\left\{-1,1\right\}^{d}$ In this scenario, $\rho\left(\text{d}z\right)$ does not admit a density with respect to Lebesgue measure but the results discussed previously can be directly extended to this scenario. . It relies on the following event rates

while the transition kernel is selected as

It is also possible to further exploit any additive decomposition of $U\left(x\right)$ within this framework and this has been used to develop an efficient sampling algorithm for big data . Again, it is easy to show that Assumption A(A2) is satisfied.

3.3 BPS sampler with randomized bounces

Alternatives to bounces of the form (28) have been proposed where one uses

and $\psi\left(v\right)=g\left(|v|\right)$. In this case, Assumption A(A1).3 is verified if

Here $\psi$ will be the standard multivariate normal distribution. We consider the scenario where $\lambda\left(x,v\right)=\left\langle\nabla U\left(x\right),v\right\rangle_{+}$ as in the global BPS. To present the various methods proposed in the literature, a decomposition of the velocity similar to that adopted in is useful:

where $n_{\perp}$ and $n_{\parallel}$ are unit norm vectors such that

All the randomized bounce procedures return a vector $v^{\prime}$

where $\langle n^{\prime}_{\perp},n_{\parallel}\rangle=0$. With this notation, we obtain $\lambda\left(x,-v^{\prime}\right)=a_{\parallel}^{\prime}{}_{+}|\nabla U\left(x\right)|.$

Let $\chi\left(k\right)$ and $\chi^{2}\left(k\right)$ be the $\chi$ and $\chi^{2}$ distributions respectively, with $k$ degrees of freedom. Under $\psi$, the random variables $a_{\perp}$ and $a_{\parallel}$ are independent and satisfy

Indeed, we have $a_{\perp}^{2}\sim\chi^{2}\left(d-1\right)$ and $a_{\perp}\geq 0$. We give below some examples of kernels $Q_{x}(v,\text{d}v^{\prime})$ satisfying Equation (30).

Independent sampling : proposes using $Q_{x}\left(v,\text{d}v^{\prime}\right)\propto\psi\left(\text{d}v^{\prime}\right)\lambda\left(x,-v^{\prime}\right)\propto a^{\prime}_{\parallel+}\psi\left(\text{d}v^{\prime}\right)$ which satisfies (30) but a scheme to sample this distribution was not given. Using the parameterization (31)-(33), (34) shows this can be achieved by sampling $a_{\parallel}^{\prime}$ according to a density proportional to $a_{\parallel+}^{\prime}$ times the standard normal density, which is equivalent to sampling $a_{\parallel}^{\prime}\sim\chi\left(2\right)$. Finally, sample $v^{*}\sim\psi$ and set $a_{\perp}^{\prime}\thinspace n^{\prime}_{\perp}=v^{*}-\left\langle v^{*},n{}_{\parallel}\right\rangle n_{\parallel}$.

Autoregressive bounce: this is a new scheme where one samples $a_{\parallel}^{\prime}\sim\chi\left(2\right)$ with probability $p_{b}$ and $a_{\parallel}^{\prime}=-a_{\parallel}$ otherwise, sample $v^{*}\sim\psi$ and set $a_{\perp}^{*}\thinspace n{}_{\perp}^{*}=v^{*}-\left\langle v^{*},n{}_{\parallel}\right\rangle n_{\parallel}$. Finally, set $a_{\perp}^{\prime}\thinspace n^{\prime}_{\perp}=\rho\thinspace a_{\perp}\thinspace n{}_{\perp}+\sqrt{1-\rho^{2}}\thinspace a_{\perp}^{*}\thinspace n{}_{\perp}^{*}$ for $\rho\in$.

The properties of these randomized bounces are not yet well understood. In Section 6, we compare them experimentally on a variety of models.

4 Hamiltonian PD-MCMC

where $K(v)=v^{T}v/2$ and $\mu\left(x\right)\propto\exp(-V\left(x\right))$ is an auxiliary probability density ensuring $\Phi_{t}$ is analytically tractable, e.g., $V$ is quadratic or linear . For example if $\pi\left(x\right)$ is a posterior density arising from a Gaussian prior, then $\mu\left(x\right)$ could be this Gaussian prior. Alternatively, $\mu\left(x\right)$ can always be selected as a Gaussian approximation to $\pi\left(x\right)$. We can then rewrite the target as $\rho(z)=\exp\left(-H\left(z\right)\right)$ where

where $\widetilde{U}\left(x\right):=U\left(x\right)-V\left(x\right)$. This is the same rationale as in elliptical slice sampling-type algorithms : both schemes use an exact Hamiltonian dynamics associated with an approximation of $\pi$ to explore the space. The difference with these algorithms and the method proposed here is that we correct for the discrepancy between $\mu$ and $\pi$ by using a PDMP mechanism instead of slice sampling techniques.

The Hamiltonian flow $\Phi_{t}$ is induced by the ODE of drift $\phi=\left(\phi_{x},\phi_{v}\right)$ where $\phi_{x}=\nabla_{v}\widehat{H}\left(z\right)=v$ and $\phi_{v}=-\nabla_{x}\widehat{H}\left(z\right)=-\nabla V\left(x\right)$. Hence, we have $\nabla\cdot\phi\left(z\right)=0$ and

We can alternatively use the randomized bounces described in Section 2.3.3 substituting $\widetilde{U}$ for $U$. Figure 1 illustrates a sample path obtained from the resulting Hamiltonian BPS algorithm. Local and doubly stochastic versions of this algorithm as for BPS can also be directly developed.

In the big data examples considered in , one could for example use for $\mu$ a Gaussian approximation of $\pi$. A local algorithm can then be obtained using for $\nabla\widetilde{U}_{i}$ the difference of the gradient of the log-likelihood corresponding to data $i$ and the properly rescaled gradient of the log-approximate posterior, as in . If the terms $\nabla\widetilde{U}_{i}$ are locally bounded, we can simulate exactly the PDMP using thinning techniques which boil down to data subsampling . This provides an alternative to which also exploits Hamiltonian dynamics and subsampling but does not preserve $\pi$ as invariant distribution.

Finally, we also note that the methods introduced in this section can be combined with the HMC algorithm of proposed to perform exact simulation of constrained normal distributions. This extends significantly the applicability of the work in , which can be viewed as a special case where $\widetilde{U}=0$. An alternative approach to constrained problems is proposed in but it is limited to piecewise-linear dynamics.

5 Using generalized Metropolis–Hastings transitions at event times

All the algorithms we have considered so far are such that only a part of the state $z=\left(x,v\right)$ is updated at event times, i.e., the transition kernel is of the form $Q\left(z,\text{d}z^{\prime}\right)=\delta_{x}(\text{d}x^{\prime})Q_{x}\left(v,\text{d}v^{\prime}\right).$ We might be interested in designing more general transitions kernels satisfying Assumption A(A1).3 and similarly Assumption A(A2).3 or Assumption A(A3).3.

For sake of illustration, consider Assumption A(A1).3. This can be rewritten as

for the probability measure $\bar{\rho}\left(\text{d}z\right)\propto\rho\left(\text{d}z\right)\lambda\left(z\right)$ assuming that $\int\rho\left(\text{d}z\right)\lambda\left(z\right)<\infty$, a weak condition which we assume holds. If the mapping $\mathcal{S}$ is an involution, i.e., $\mathcal{S}^{-1}=\mathcal{S}$, and we can design a kernel $Q$ satisfying the so-called skewed detailed balance condition

then it follows directly by integrating both terms in this equality with respect to variable $z$ that it will satisfy (36).

We present here a generic mechanism which can be used to achieve this known as the Generalized Metropolis–Hastings (GMH) algorithm. The GMH algorithm is a simple extension of MH; see for example [31, pp. 74–77]. For a probability measure $\nu\left(\text{d}z\right)=\nu\left(z\right)\text{d}z$ on $\mathcal{Z}$, let us consider the following GMH kernel defined for a Markov proposal kernel $M$ by

Conditions on $\nu,$ $\mathcal{S},$ $M$ and $g$

The mapping $\mathcal{S}$ is an involution, i.e., $\mathcal{S}^{-1}=\mathcal{S}$.

is defined and positive for almost all $\left(z,z^{\prime}\right)\in\mathcal{Z\times Z}.$

Assumption A(A4).1 is satisfied for $g\left(r\right)=\min\left(1,r\right)$. For a deterministic proposal $M\left(z,\text{d}z^{\prime}\right)=\delta_{\Psi\left(z\right)}\left(\text{d}z^{\prime}\right),$ Assumption A(A4).3 is satisfied if $\Psi$ admits an inverse $\Psi^{-1}$ such that

and then the acceptance probability is given by

Assume (A(A4)). Then the GMH kernel $T$ defined by (38) satisfies the following skewed detailed balance condition

If additionally $\mathcal{S}$ is a $\nu$-preserving mapping then the GMH kernel is $\nu$-invariant.

The proof of this result follows from direct calculations given in the Appendix and can also be found in [31, pp. 74–77]. Using this result, it is possible to check easily Assumption A(A1).3 for the BPS and Zig-Zag processes. For example, for the BPS, $Q$ is of the form (38) with $\nu=\bar{\rho}$, $\mathcal{S}^{-1}=\mathcal{S}$, $g\left(r\right)=\min\left(1,r\right)$ as we use a deterministic proposal $\Psi\left(z\right)=\left(x,R_{\nabla U}\left(x\right)v\right)$ which verifies $\Psi^{-1}=\mathcal{S}\circ\Psi\circ\mathcal{S}$ so $\beta\left(z,z^{\prime}\right)=1$ for all $z,z^{\prime}$. Hence by Proposition 4, $Q$ satisfies the skewed detailed balance (37), hence it satisfies (36).

The benefit of the GMH approach is that it allows us to define much more general kernels at event times. For example one could use a deterministic proposal with $\Psi\left(z\right)=\left(x,R_{\nabla\widehat{U}}\left(x\right)v\right)$ where $\nabla\widehat{U}$ is a computationally cheap approximation of $\nabla U$. It is valid to use such a deterministic proposal at it satisfies $\Psi^{-1}\left(z\right)=\mathcal{S}\circ\Psi\circ\mathcal{S}(z)$. In this case, there is a probability of the bounce being rejected and setting $z^{\prime}\leftarrow\mathcal{S}\left(z\right)$. We can also use transition kernels which modify the component $x$ of $z$.

Discrete-time PDMP and PD-MCMC

We introduce here the class of discrete-time PDMP and present general conditions for such processes to ensure invariance w.r.t. a strictly positive density $\rho\left(z\right)=\exp(-H\left(z\right))$. These conditions parallel the conditions given Section 2.2 for continuous-time algorithms.

a diffeomorphism $\Phi:\mathcal{Z}\rightarrow\mathcal{Z}$ with the absolute value of the determinant of the Jacobian satisfying $|\nabla\Phi\left(z\right)|>0$ for all $z$,

an acceptance probability $\alpha:\mathcal{Z}\rightarrow$ with $1-\alpha\left(z\right)$ being the probability of having an event at the next time step when the current state is $z$, and

a Markov transition kernel $Q$ from $\mathcal{Z}$ to $\mathcal{Z}$ where the state at event time $t$ is given by $z_{t}\sim Q\left(z_{t-1},\cdot\right)$.

2 From discrete-time PDMP to PD-MCMC

Similarly to Section 2.2, assume we are interested in sampling a strictly positive density $\rho\left(z\right)$ given by (7) using a discrete-time PDMP process. Invariance of the kernel $K$ with respect to $\rho$ is satisfied if, by definition, one has

All the following developments could also be adapted to sample from distributions on discrete spaces but this will not be discussed here.

We provide here useful sufficient conditions on $\Phi,$ $\alpha$, and $Q$ to ensure $\rho$-invariance of the associated discrete-time PDMP, without making any structural assumption on these objects.

There exists a $\rho$-preserving mapping $\mathcal{S}:\mathcal{Z}\to\mathcal{Z}$.

The acceptance probability $\alpha$ satisfies

Conditions A(A5).1 to A(A5).3 parallel the conditions A(A1).1 to A(A1).3.

Assume (A(A5)). Then the discrete-time PDMP admits $\rho$ as invariant distribution.

When $\mathcal{S}$ is an involution so that $\rho(\mathcal{S}^{-1}\left(\text{d}z^{\prime}\right))=\rho(\mathcal{S}\left(\text{d}z^{\prime}\right))$, condition A(A5).3 can be interpreted as a “skewed” invariance condition on $\nu(\text{d}z)\propto\rho(\text{d}z)(1-\alpha(z))$. The quantity $\rho(\text{d}z)\left(1-\alpha\left(z\right)\right)$ is proportional to the invariant distribution of the “jump chain,” i.e. the distribution of those states where the proposal $\Phi\left(z\right)$ is rejected. It has a clear analogue in the continuous-time scenario where the jumps occur at states with distribution proportional to $\rho(\text{d}z)\lambda(z)$.

2.2 Sufficient conditions for local methods

In scenarios where $H\left(z\right)$ can be decomposed as in (11), it will prove convenient to consider an acceptance probability of the form

Based on these structural assumptions on $\alpha$ and $Q$, we can provide useful sufficient “local” conditions on $\Phi,$ $\{\alpha_{i}:i\in\left[n\right]\}$ and $\left\{Q_{B}:b\in\mathcal{B}\right\}$ to ensure invariance of the associated discrete-time PDMP w.r.t. $\rho$ is satisfied.

Conditions on $\phi,$ $\{\alpha_{i}:i\in\left[n\right]\}$, and $\left\{Q_{b}:b\in\mathcal{B}\right\}$

There exists a $\rho$-preserving mapping $\mathcal{S}:\mathcal{Z}\to\mathcal{Z}$.

The acceptance probabilities $\{\alpha{}_{i}:i\in\left[n\right]\}$ satisfy

For all $b\in\mathcal{B}$, the transition kernel $Q_{b}$ satisfies

For a mapping such that $\left|\nabla\Phi\right|=1$, then Assumption A(A6).2 is satisfied if for all $i\in\left[n\right]$

Assume (A(A6)). Then the discrete-time PDMP admits $\rho$ as invariant distribution.

2.3 Sufficient conditions for doubly stochastic methods

Consider finally the scenario where $H\left(z\right)$ is given by (17). In this context, we consider an acceptance probability of the form

Conditions on $\phi,$ $\{\alpha_{\omega}:\omega\in\Omega\}$ and $\left\{Q_{P}:P\in\mathcal{\mathcal{\mathscr{\mathcal{\mathcal{P}}}}}\right\}$

There exists a $\rho$-preserving mapping $\mathcal{S}:\mathcal{Z}\to\mathcal{Z}$.

The acceptance probabilities $\{\alpha_{\omega}:\omega\in\Omega\}$ satisfy

For all $P\in\mathcal{\mathcal{\mathscr{\mathcal{\mathcal{P}}}}}$, the transition kernel $Q_{P}$ satisfies

and the Radon-Nikodym derivative in the expression above is well-defined and strictly positive for $Q_{P}\left(z,\text{d}z^{\prime}\right)$ almost all $z^{\prime}$.

For a mapping such that $\left|\nabla\Phi\right|=1$, Assumption A(A7).2 is satisfied if for all $\omega\in\Omega$

Assume (A(A7)). Then the discrete-time PDMP admits $\rho$ as invariant distribution.

3 Existing PD-MCMC algorithms

This scheme satisfies Assumption A(A5).1 to Assumption A(A5).3 and is thus $\rho$-invariant. In particular Assumption A(A5).3 is satisfied as Steps 2 and 3 correspond to using for the event kernel $Q$ a GMH kernel satisfying the skewed-detailed balance condition (42) for $\nu\left(\text{d}z\right)\propto\rho\left(\text{d}z\right)\left(1-\alpha\left(z\right)\right)$.

Algorithm 3 can be alternatively viewed as a composition of reversible kernels. First, a delayed-rejection algorithm proposing $\Phi$ and, in case of rejection, then proposing $M(z,\mathcal{S}^{-1}(\cdot))$. Second, the involution $\mathcal{S}$ is applied unconditionally. In the delayed-rejection framework, we can view condition A(A5).3 as a condition on delayed-rejection kernels expressed in a sort of “remainder” form. While our algorithm uses two proposals, extending this remainder condition to multiple proposals would require that each $Q_{k}$ satisfies $\int\rho(\text{d}z)\prod_{i=1}^{k-1}(1-\alpha_{i}(z))Q_{k}(z,\text{d}z^{\prime})=\rho(\text{d}z^{\prime})\prod_{i=1}^{k-1}(1-\alpha_{i}(z^{\prime}))$.

This algorithm was proposed in . It is a special case of Algorithm 3 which uses $\Phi\left(z\right)=(x+v\epsilon,v)$ for some $\epsilon>0$ and a proposal $M\left(z,\text{d}z^{\prime}\right)=\delta_{\mathcal{S}\left(z\right)}\left(\text{d}z^{\prime}\right)$ which is accepted with probability 1.

3.2 Hamiltonian Monte Carlo

The celebrated HMC algorithm proposed in is also a special case of Algorithm 3 which uses a proposal $M\left(z,\text{d}z^{\prime}\right)=\delta_{\mathcal{S}\left(z\right)}\left(\text{d}z^{\prime}\right)$. However, contrary to guided random walk, it is using for $\Phi$ a symplectic integrator targeting the Hamiltonian $H$. This deterministic proposal satisfies indeed $\left|\nabla\Phi\right|=1$ and $\Phi^{-1}=\mathcal{S}\circ\Phi\circ\mathcal{S}$ (see, e.g., ). The resulting PD-MCMC kernel $K$ is usually combined with a momentum refreshment step $v\sim\psi$.

3.3 Reflective Slice Sampling: discrete-time BPS schemes

Several versions of slice sampling, known as reflective slice sampling, are based on bounces similar to the BPS and are also a special case of Algorithm 3; see [37, Section 7]. They rely $\Phi\left(z\right)=\left(x+v\epsilon,v\right)$ for some $\epsilon>0$ and a deterministic proposal $M\left(z,\text{d}z^{\prime}\right)=\delta_{\Psi\left(z\right)}\left(\text{d}z^{\prime}\right)$. Reflective slice sampling with inner reflections is using $\Psi\left(z\right)=\left(x^{*},v^{*}\right)=\left(x,R_{\nabla U}(x)v\right)$ while reflective slice sampling with outer reflections is using $\Psi\left(z\right)=\left(x^{*},v^{*}\right)=\left(x+v\epsilon+R_{\nabla U}(x+v\epsilon)v\epsilon,R_{\nabla U}(x+v\epsilon)v\right)$. Both proposals satisfy $\Psi^{-1}=\mathcal{S}\circ\Psi\circ\mathcal{S}$. The outer version of the algorithm has been recently proposed independently in ; see also for a related proposal in the context of nested sampling. In either case, the acceptance probability simplifies to

Under regularity conditions, reflective slice sampling with inner reflections converges weakly to the BPS for $\lambda_{\textrm{ref}}=0$ as $\epsilon\rightarrow 0$.

A precise mathematical statement, Theorem 12, and its proof are given in Appendix B. We can modify this algorithm to include a refreshment, i.e. by sampling $v^{\prime}\sim\psi$ with probability $\lambda_{\textrm{ref}}\epsilon$. This weak convergence result of Proposition 11 can be directly extended to this case to show that the resulting discrete-time process converges weakly to the BPS process with refreshment rate $\lambda_{\textrm{ref}}$. Note that the kernel $K$ would still be $\rho$-invariant if $\Phi$ were using a computationally cheap approximation $\nabla\widehat{U}$ of $\nabla U$ to bounce. However, this discrete-time algorithm does not converge to the BPS process as the probability of accepting $z^{\prime}=\mathcal{S}\left(z\right)$ does not vanish as $\epsilon\rightarrow 0$ in this scenario. Under regularity conditions, it will instead converge towards the algorithm described at the end of Section 2.5.

4 Extensions

For the kernel $M_{x}\left(v,\cdot\right)$, we can use the randomized bounces developed in Section 2.3.3 as well as $M_{x}\left(v,\cdot\right)=\psi\left(\cdot\right).$ The forward-event , generalized BPS , and autoregressive bouncing procedures discussed in Section 2.3.3 induce a transition kernel $M_{x}$ satisfying $\psi(v)\langle\nabla U(x),v\rangle_{+}M_{x}(v,v^{\prime})=\psi(-v^{\prime})\langle\nabla U(x),-v^{\prime}\rangle_{+}M_{x}(-v^{\prime},-v)$, for which we would expect that the acceptance ratio in Step 2.b of Algorithm 4 will be close to 1 for small $\epsilon$.

The invariance with respect to $\rho$ of the transition kernel is easy to check. Assumption A(A5).1 is clearly satisfied. Assumption A(A5).2 follows from direct calculations using $\left|\nabla\Phi\right|=1$ and $\Phi^{-1}=\mathcal{S}\circ\Phi\circ\mathcal{S}$. Finally Assumption A(A5).3 follows from the fact that the event kernel corresponding to steps 2.a to 2.c of Algorithm 4 is a GMH kernel with $\nu\left(z\right)\propto\rho\left(z\right)\left(1-\alpha\left(z\right)\right)$ with a proposal kernel $M_{x}\left(v,\text{d}v^{\prime}\right)$.

4.2 Discrete-time Hamiltonian BPS

The invariance with respect to $\rho$ of the transition kernel is easy to check. Assumption A(A5).1 is obviously satisfied. Assumption A(A5).2 follows from direct calculations using $\left|\nabla\Phi\right|=1$ and $\Phi^{-1}=\mathcal{S}\circ\Phi\circ\mathcal{S}$. Finally Assumption A(A5).3 follows from the fact that the event kernel corresponding to step (a) and (b) of Algorithm 5 is a GMH kernel with $\nu\left(z\right)\propto\rho\left(z\right)\left(1-\alpha\left(z\right)\right)$ with a deterministic transition kernel satisfying $\Psi^{-1}=\mathcal{S}\circ\Psi\circ\mathcal{S}$. If $\Phi$ is a leapfrog integrator of stepsize $\epsilon>0$ targeting the Hamiltonian $H\left(z\right)$, then the strategy described above is not directly applicable as $\widetilde{U}\left(x\right)=0$ for all $x$ so $R_{\nabla\widetilde{U}}\left(x\right)$ is not defined. However as $\Phi$ can be thought of as the exact time discretization of a shadow Hamiltonian of the form $\hat{H}_{\epsilon}\left(z\right)=H\left(z\right)-\epsilon^{2}\widetilde{H}\left(z\right)+\mathcal{O}\left(\epsilon^{4}\right)$ [30, p. 107], it may be possible to build bounces based on $\widetilde{H}\left(z\right)$ to correct for the discrepancy between the true Hamiltonian dynamics and its leapfrog approximation.

4.3 Discrete-time gradient-free BPS

The BPS-type algorithms given thus far all require computation of the gradient of the potential $\nabla U(x)$ in order to update the velocity $v$ when a bounce event occurs. However, we may wish to target potential functions where this gradient cannot be computed or is very expensive to compute. Additionally, the gradient may not be informative in some models, such as certain embeddings of discrete spaces where the gradient may be zero almost everywhere.

A scheme to approximate the gradient $\nabla U(x)$ by computing numerical differences was advanced in . Here, some number $n_{cpt}$ of orthogonal unit vectors $\zeta_{i},i\in[n_{cpt}]$ are selected, and the gradient approximated along each of these vectors by, e.g.,

for some small value $h$. The combination of these $n_{cpt}$ vectors yields an approximation to the gradient

which for $n_{cpt}=d$ is a typical numerical approximation to the gradient. The new velocity is found by a reversible map from the old velocity to the new velocity which preserves the magnitude of the velocity and maintains the projection of the velocity on the gradient vector.

We may derive an algorithm which operates in the same spirit as that of . By taking $n_{cpt}$ orthogonal unit vectors, here selected randomly and independently of $v$, we can achieve a reversible algorithm by simply taking the reflection off of the approximate gradient

and accepting this proposal in the same way we would accept a typical bounce in the discrete-time BPS algorithm; specifically, by accepting the bounce with probability

Alternatively, we propose an algorithm which is related to the continuous-time randomized bounces of Section 2.3.3. We had previously noted that the independent sampling algorithm proposed in consists of sampling from the distribution proportional to $\psi(v^{\prime})\lambda(x,-v^{\prime})$, independently of the current value of $v$. Based on the discrete-time invariance condition (50), we may analogously sample from the distribution proportional to $\psi(v^{\prime})\left[\pi(x)-\pi(x-v^{\prime}\epsilon)\right]_{+}$. This can be accomplished by using rejection sampling with instrumental distribution $\psi$, noting that the ratio between the densities is bounded above by $\pi(x)$; thus each rejection sampling proposal $v^{\dagger}$ is accepted with probability $\left[\pi(x)-\pi(x-v^{\dagger}\epsilon)\right]_{+}/\pi(x)$, and the first accepted proposal is also accepted as the new state $v^{\prime}$. See Algorithm 6 for details of this rejection-sampling scheme.

4.4 Efficient Implementation of Discrete-time PD-MCMC

Finally there are scenarios where it is possible to directly simulate an event time from (43). For example, assume that $\pi\left(x\right)=\exp(-U\left(x\right))$ where $U$ is strictly convex, $\Phi\left(z\right)=(x+v\epsilon,v)$ and $\alpha\left(z\right)=\min\left\{1,\rho\left(\Phi\left(z\right)\right)/\rho\left(z\right)\right\}=\min\left\{1,\exp\left(-\left(U\left(x+v\epsilon\right)-U\left(x\right)\right)\right)\right\}$ then it is easy to show that Algorithm 7 returns a sample from (43). This adapts the approaches developed in [10, Section 2.3.1] for the continuous-time BPS algorithm to the discrete-time case.

All these strategies can be easily combined. For example, we can use an upper bound $\bar{\alpha}\left(z,k\right)=\prod_{i=1}^{n}\bar{\alpha}_{i}\left(z,k\right)$ where $\rho_{i}\left(z\right)$ is strictly log-concave for some $i\in[n]$.

Discrete-time local PD-MCMC

Note that $\nabla\overline{U}\left(x\right)$ depends on both $u$, $v$ and $\epsilon$, we stress this dependence as it is omitted notationally.

Algorithms 8 and 9 might appear of limited interest as they require to sample $n$ Bernoulli random variables at each iteration. In the next sections, we show how we can propose implementations that parallel the priority queue implementation of the local BPS proposed in , see [10, Section 3.3.1] for a detailed description, as well as the subsampling algorithms proposed in [10, 6, 29, Section 3.3.2].

2 Prefetching implementation

We first describe a priority queue type implementation of Algorithm 9 based on parallel prefetching ideas in scenarios where

$x_{S_{i}}$ being a subset of the components of $x$ and $\pi_{i}(x)=\exp\left(-U_{i}\left(x_{S_{i}}\right)\right)$. There are many possible variations of this implementation.

The efficiency of Algorithm 10 relies on the capability of computing the $\tau_{i}$ efficiently. This may be possible when, for example, this is done in parallel or when we some property of $\pi_{i}$ allows it, such as in the case of log-concave targets detailed as in Algorithm 7 given above.

3 Subsampling implementations

For sufficiently small $\epsilon$, we might expect that in Step 1 of Algorithm 9 would yield very few indices for which $B_{i}=1$. This motivates an approach which can sample these variables more efficiently by finding an upper bound on the probability that $B_{i}=1$, essentially allowing us to bound the number of indices for which $B_{i}=1$. We present Algorithm 11; here, the acceptance of the bounce move (64) is computed in two stages: in Step 4.b we simulate events of probability $1-\min\left\{1,\frac{[\pi_{i}(x)-\pi_{i}(x-v^{*}\epsilon)]_{+}}{[\pi_{i}(x)-\pi_{i}(x+v\epsilon)]_{+}}\right\}$ for each $i$ where $B_{i}=1$, if these succeed then in Step 4.c we simulate events of probability $1-\min\left\{1,\frac{\min\left(\pi_{i}(x),\pi_{i}(x-v^{*}\epsilon)\right)}{\min\left(\pi_{i}(x),\pi_{i}(x+v\epsilon)\right)}\right\}$ for each $i$ where $B_{i}=0$. We suggest that one can make use of efficient procedures described in Algorithm 12 and Algorithm 13 to sample multiple Bernoulli random variables in both Steps 1 and 4.c; in both cases we expect few cases where the respective Bernoulli variables are 1. While Step 4.b also samples a set of Bernoulli variables, our assumption that $\epsilon$ is small suggests that the number of variables sampled here will be small; as such this step may be inexpensive and there is likely little to be gained by a more sophisticated simulation scheme.

In this case, we can determine the set $\left\{i:X_{i}=1\right\}$ using Algorithm 12. This incurs a computational complexity $O(1+|I|\bar{p})$ compared to $O(|I|)$ for the direct implementation . This implementation can be thought of as the discrete-time version of the thinning ideas leading to the “naive” subsampling techniques presented in .

Second, if we instead have access to local bounds $0\leq\bar{p}_{i}\leq 1$ such that

we could obviously use the previous strategy by setting $\bar{p}:=\max_{i\in I}\bar{p}_{i}$ but this strategy can be highly inefficient if, e.g., most bounds $\bar{p}_{i}$ are very close to zero and a few are close to 1. In this scenario, it is possible to use instead Algorithm 13 which relies on the simulation of Poisson random variables. This algorithm can be thought of as the discrete-time version of the thinning ideas leading to the “informed” subsampling techniques presented in .

For this algorithm to be of practical interest, the bounds $\overline{p}_{i}$ and the associated Poisson rates $\kappa_{i}$ should not have to be recomputed at each time step as for the examples considered in . In this scenario, it is then possible to use the alias method or ordered marginally uniform random variables on $\left[0,1\right]$ to sample efficiently from the multinomial distributions in complexity $O(S)$ .

The availability of an upper bound for Step 1, denoted here $\bar{p}_{i}(x,v)$, can be seen as equivalent to a lower bound on $\alpha_{i}(z)$ as discussed in Section 3.2.2, since

For Step 4.c, we would seek an upper bound

This bound may be achieved, for example, when $|\nabla U_{i}(x^{\prime})|<\delta$ for all $\left\{x^{\prime}:|x^{\prime}-x|<\epsilon|v^{*}|\right\}$. In this case, an upper bound can be derived using

Discrete-time doubly stochastic PD-MCMC

The kernel $Q_{P}$, which must satisfy (59), may be implemented using a scheme similar to the GMH. Using standard results on Poisson point processes and Assumption A(A7).3, the condition (61) can be simplified as

suggesting a GMH kernel with deterministic proposal $\Psi_{P}(z)$ satisfying $\Psi_{P}^{-1}=\mathcal{S}\circ\Psi_{P}\circ\mathcal{S}$ and acceptance probability

which arises by treating the integral terms and the product terms as two factors, each with its own acceptance probability. Based on this, we present Algorithm 14, wherein we sample an event of probability (68) using a two-stage acceptance procedure.

By selecting $\alpha_{\omega}(z)=\min\left(1,\frac{\rho_{\omega}(\Phi(z))}{\rho_{\omega}(z)}\right)=\exp\left(-\left[H_{\omega}(\Phi(z))-H_{\omega}(z)\right]_{+}\right)$ for $\rho_{\omega}\left(z\right)=\exp(-H_{\omega}\left(z\right))$, the acceptance probability (68) takes the form

Further allowing $\mathcal{S}(z)=(x,-v)$, $\Phi(z)=(x+v\epsilon,v)$ and $\Psi_{P}(z)=(x,R_{\nabla\overline{U}}(x)v)$ with $\nabla\overline{U}(x)=\sum_{\omega\in P}\nabla U_{\omega}(x)$ yields

The first term of this acceptance ratio, viewed as a void probability of a Poisson process, can be interpreted as the “excess” rate of $\alpha(x,-R_{\nabla\overline{U}}(x)v)$ over $\alpha(x,v)$; in other words, the probability that no extra points would be simulated for $P$ when in state $(x,-R_{\nabla\overline{U}}(x)v)$.

In either case, the simulation of Poisson processes $P$ and $P^{\prime}$ is possible when those rates can be bounded. If we have some lower bound $\underline{\alpha_{\omega}}(z)\leq\alpha_{\omega}(z)$ for which we can simulate a Poisson process of intensity $-\log\underline{\alpha_{\omega}}(z)\mu(\text{d}\omega)$, then we can recover $P$ by thinning this process. This condition is sufficient for simulation of $P^{\prime}$ as the corresponding intensity is bounded by $-\log\underline{\alpha_{\omega}}(\mathcal{S}\circ\Psi(z))$; however, it may be possible to bound the intensity of $P^{\prime}$ more tightly in some situations.

The idea of introducing a Poisson process so as to deal with the intractability of target distribution can also be exploited within a standard MCMC setting. For simplicity, assume a symmetric proposal density $q\left(\left.z^{\prime}\right|z\right)$ then it is easy to check that Algorithm 15 corresponds to a transition kernel which is reversible with respect to $\rho\left(z\right)=\exp(-\int H_{\omega}\left(z\right)\mu\left(\text{d}\omega\right))$.

2 For measures containing atoms

While it remains sufficient to use the acceptance probability (68), we note that a partition of $P^{*}$ into sets of equivalent $P_{\omega}$ (and therefore equivalent bounce proposals $\Psi_{P_{\omega}}(z)$) will yield a sufficient condition which is “integrated out” in the sense that the total density of the forward and reverse transitions are captured.

and that the Radon-Nikodym derivative above is well-defined and strictly positive for $Q_{P_{\omega}}^{*}(z,\text{d}z^{\prime})$-almost all $z^{\prime}$. The above implies an algorithm similar to Algorithm 14 but where $\Psi_{P_{\omega}}(z)$ would be accepted with a probability of

Numerical results

In , the local BPS algorithm was shown to outperform various state-of-the-art HMC algorithms in sparse precision Gaussian random field models with Poisson observations. In this section, we investigate the relative performance of local BPS and Hamiltonian BPS in the same setting. We find that Hamiltonian BPS has a modest advantage over local BPS when the number of observations is small but the dimensionality of the latent variables is high. On the other hand, when the number of observations is equal to the number of latent variables, the situation is reversed. However in both regimes Hamiltonian BPS outperforms global BPS, and it is worth keeping in mind that there are situations where Hamiltonian BPS is applicable while the local BPS is not computationally attractive, for example if a single variable is connected to all factors.

More generally, if $V$ is an arbitrary normal distribution, the situation considered here can be used after a change of variables. The computational trade-off results we present in this section are hence representative of situations where we have a high-dimensional Gaussian prior with a precision matrix admitting a Cholesky decomposition that can be computed in time $O(d)$, which arises for example in certain time series models and corresponds to a best case scenario for Hamiltonian BPS.

1.2 Exact simulation of bounce times

Let $j\in\left\{1,2,\dots,k\right\}$ index the observations. Assume that the negative log-likelihood $\widetilde{U}\left(x\right)$ can be decomposed as $\widetilde{U}(x)=\sum_{j=1}^{k}\widetilde{U}_{j}(x_{i(j)})$ for some function $i(\cdot)$ mapping observation indices to the latent variable indices. As a pre-processing step, we compute (numerically or analytically) a bound $B_{j}(b)\geq\sup\left\{|\nabla\widetilde{U}_{j}(x)|:|x|<b\right\}$.

Let $x=x_{i(j)}$ and $v=v_{i(j)}$ denote the initial position and velocity at the beginning of the current piecewise Hamiltonian segment for the latent variable $i(j)$ associated with observation $j$. From Section 2.3.3 of , it is enough to simulate the bounce time of a single factor $\widetilde{U}_{j}(x_{i(j)})$. Using the methodology developed in [10, Section 2.3.2], we simulate the bounce time of each factor using thinning and the following bound on the intensity $\chi(t)$:

where $\alpha=\arctan(-v/x)$, $\beta=\arctan(v/d)$.

1.3 Results

We consider a likelihood given by conditionally independent Poisson observations with observations $y_{i}$ having a natural exponential family parameter given by the latent random variable $x_{i}$:

We compare three algorithms: local and global BPS with piecewise linear trajectories, and Hamiltonian BPS. Computation of the bounce times for the piecewise linear trajectories is done as in . For the bounce times of Hamiltonian BPS, we use the result from Section 6.1.2 with $B_{j}(b)=\exp(b)+y_{j}$.

We show in Figure 2 the scaling of the CPU wall clock time required to obtain one effective sample size (ESS) as a function of the dimensionality $d$ (log-log scale). The wall clock time is measured in milliseconds on a 2.8 GHz Intel Core i7, and the ESS is computed using a batch mean estimator with a test function given by $f(x)=x_{1}^{2}$. Expectations from piecewise-deterministic trajectories are computed analytically as shown in and from piecewise Hamiltonian trajectories, using numerical integration. For each dimension and algorithm, we run $100$ independent chains and average the running times per ESS.

2 Empirical comparisons of local and global BPS to HMC and Standard and Elliptical Slice Sampling

We consider four models, built from two prior distributions: first, a Brownian bridge prior, and second, a diagonal precision prior. For each prior, we consider either a Poisson likelihood with synthetic observations (with the same structure as described in the previous section), or no likelihood function. We consider the following sampling methods: the Elliptical Slice Sampler , the “Standard” Slice Sampler (with exponential slice growing and slice shrinking) , HMC, or more precisely the NUTS algorithm implemented in Stan, the local and global BPS algorithm with linear trajectories, and the Hamiltonian BPS algorithm. For each combination, we run the algorithms on latent fields of dimensionality $\left\{2^{0},2^{1},2^{2},\dots,2^{7}\right\}$, and replicate the experiment 50 times with different random seeds. We measure ESS and wall clock time. ESS is computed using a batch mean estimator with a test function given by $f(x)=x_{1}^{2}$.

2.2 Results

We summarize the main results of this section in Figure 3, where the empirical computational complexity (wall clock time (ms) per ESS) is plotted in log-log scale against the dimensionality of the field for the four models. For sufficiently high-dimensional scenarios (> 10 dimensions), local BPS outperforms all other methods in 3 out of the 4 settings. In the fourth setting, (Diagonal Precision + Poisson Likelihood), NUTS (HMC) and Local BPS outperform the other methods, but neither strictly dominate the other. Elliptic Slice Sampling is competitive when there is no likelihood, but it is still not better than Local BPS, presumably because the latter can use the full trajectory when computing averages whereas Elliptical is discrete-time. However, once the Poisson Likelihood is added, Elliptical Sampling seems to have worse asymptotics, empirically roughly $O(n^{3/2})$ versus roughly $O(n^{1+\epsilon})$ for the best performing methods.

3 Randomized bounces

In this section, we compare the performance of several collision operators on two collections of problems of increasing dimensionality.

The first collection of target distributions we consider consists in funnel distributions from , namely multivariate normals of varying dimension $d$ with diagonal covariance matrix and standard deviations for each components given by $1,(d-1)/d,(d-2)/d,\dots,1/d$. Since the algorithms considered are rotationally invariant, this is representative of problems with averse conditioning. The second collection consists in isotropic multivariate normal of increasing dimensionality $d$. The isotropic examples are useful to identify cases where symmetries create a clear imperative for refreshment as discussed in . For each class of target distributions, we look at problems of dimensionality $2^{1},2^{2},\dots,2^{7}$.

We compare 8 algorithms, corresponding to $4$ different bounce operators and $2$ refreshment strategies (either independent refreshment at times determined by a unit rate homogeneous Poisson process, or no refreshment). The bounce operator labeled Flip corresponds to $Q_{x}(v,\text{d}v^{\prime})=\delta_{-v}(\text{d}v^{\prime})$, Det-Rand corresponds to the forward-event chain algorithm of , and Rand-Rand corresponds to the independent sampling algorithm of . We recorded the Monte Carlo averages $\hat{f_{i}}$ of the test function $f(x)=x_{1}^{2}$ for the trajectory up to event time index $i=2^{0},2^{1},\dots,2^{14}$ and computed the errors $e_{i}=|\hat{f}_{i}-1|$ . We then averaged the errors over $20$ independent executions of the algorithms using different random seeds. All experiments in this section are performed on a global (continuous-time) BPS algorithm. Both simulation of collision times and computation of Monte Carlo averaged are performed using closed form expressions that can be found in .

3.2 Results

We show in Figure 4 the average error as a function of the event index (log-log scale).

Our results show that in the low dimensional regime, at least two randomized bounce operators (Det-Rand and Rand-Rand) combined with no refreshment outperform the standard bounce with refreshment. However, this advantage asymptotically vanishes as the dimensionality of the problem increases. In fact, when refreshment is turned off, for all the operators but Rand-Rand, performance dramatically collapses with dimensionality. The performance drop-off is so pronounced that it may not be detected by conventional estimators of effective sample size. We can measure it here since the true value of the expectations are known.

We conjecture that this sharp drop in performance is due to a concentration of measure phenomenon making the variance of the randomized operators in the direction parallel to the gradient decrease with $d$, hence, informally speaking, making certain randomized operators such as Det-Rand more and more deterministic as $d$ increases. The lack of irreducibility of deterministic bounce operators without refreshment is shown formally in . This conjecture is also supported by the fact that reintroducing refreshment makes all methods behave similarly in high-dimensional settings (except for the cruder Flip operator).

This is noteworthy as one of the motivations for previous work on alternative bounce operators is that such operators may alleviate the need for refreshment in certain scenarios. Our results provide a cautionary example that in certain high-dimensional scenarios, it is still preferable to perform refreshment even when randomized bounces are used. Interestingly, this happens not only in the isotropic case but also in the non-isotropic, funnel distribution case, where one might expect refreshment to play a more minor role due to lack of symmetry.

Discussion

We have introduced a general framework which allows us to develop novel continuous-time and discrete-time PD-MCMC algorithms addressing some of the limitations of existing techniques. They allow to exploit dynamics dependent on the target distribution. Moreover, contrary to continuous-time algorithms, it is always possible to simulate exactly the event times.

There are many possible methodological extensions of these algorithms. To simplify presentation, we have presented our results for auxiliary distributions of the form $\psi\left(v\right)=g(|v|)$ but, as in the HMC context , it is possible to adapt these techniques to the scenario where $\rho\left(z\right)=\pi\left(x\right)\psi_{x}\left(v\right)$ with $\psi_{x}\left(v\right)=g(|v^{T}M\left(x\right)v|^{1/2})$ for $M\left(x\right)$ a positive definite matrix capturing the local curvature of $U$ around $x$. From preliminary experiments, we observe that using a position-dependent mass matrix $M(x)$ can provide significant gains in complex scenarios. Even selecting simply a suitable constant matrix $M$ can already improved substantially performance as already demonstrated for the BPS . Moreover, the proposed framework is very flexible but all the algorithms proposed so far in continuous-time are based on a divergence-free vector field and in discrete-time on a deterministic mapping with unit Jacobian determinant. There is conceptually no need to restrict ourselves to such scenarios and it would be interesting to come up with useful algorithms exploiting this degree of freedom.

From a theoretical point of view, PD-MCMC techniques appear to provide state-of-the-art performance on some interesting sampling problems but there are only few theoretical results available and there is much work to be done to better understand their properties.

References

Appendix A Proofs of invariance

Using Assumption A(A1).3 then Assumption A(A1).1, we obtain

under Assumption A(A1).2. This establishes the result. ∎

where we have used Assumptions A(A2).3 and A(A2).1. Hence, (8) is equal to

under Assumption A(A2).2. The result follows. ∎

The proof is similar to the proof of Proposition 2 and is therefore omitted. ∎

First notice that, if using Assumption A(A4).2, we define

then using the properties of the push-forward measure and Assumption A(A4).1, we have for any measurable function $h$

This establishes that the measure $\nu\left(\text{d}z\right)M\left(z,\text{d}z^{\prime}\right)$ is absolutely continuous w.r.t. $\nu\left(\mathcal{S}\left(\text{d}z^{\prime}\right)\right)M\left(\mathcal{S}\left(z^{\prime}\right),\mathcal{S}\left(\text{d}z\right)\right)$ with a Radon-Nikodym derivative given by

For the first term on the r.h.s. of (69), we have

where we have used Assumption A(A4).3 then Assumption A(A4).1.

The second term on the r.h.s. of (69) satisfies

using Assumption A(A4).1. The sum of the terms (70) and (71) is equal to $\nu\left(\mathcal{S}\left(\text{d}z^{\prime}\right)\right)T\left(\mathcal{S}\left(z^{\prime}\right),\mathcal{S}\left(\text{d}z\right)\right)$. Hence the GMH kernel satisfies the skewed detailed balance condition (37). ∎

The proof follows from simple manipulations. We have from Assumption A(A5).3 then Assumption A(A5).1 that the l.h.s. of (48) satisfies

Hence the condition (48) is satisfied if for all $z^{\prime}$

By rewriting this expression for $z^{\prime}=\Phi\left(z\right)$, and using the fact that $\left|\nabla\Phi^{-1}\left(z^{\prime}\right)\right|=\left|\nabla\Phi\left(\Phi^{-1}\left(z^{\prime}\right)\right)\right|^{-1}$ so $\left|\nabla\Phi^{-1}\left(z^{\prime}\right)\right|=\left|\nabla\Phi\left(z\right)\right|^{-1}$, we obtain Assumption A(A5).2. ∎

The proof follows from simple manipulations. We consider first the second term on the r.h.s. of (48). This satisfies

where we have used Assumption A(A6).3 then Assumption A(A6).1. The first term on the l.h.s. of (48) is given by

Hence the condition (48) is satisfied if for all $z^{\prime}$

By rewriting this expression for $z^{\prime}=\Phi\left(z\right)$, we obtain Assumption A(A6).2. which is also implied by (56) if $\left|\nabla\Phi\left(z\right)\right|=1$ for all $z$. ∎

The proof is very similar to the proof of Proposition 8. We similarly consider the second term on the r.h.s. of (48) which satisfies

where we have used Assumption A(A7).3 then Assumption A(A7).1. The first term on the l.h.s. of (48) is given by

Hence the condition (48) is satisfied if for all $z^{\prime}$

By rewriting this expression for $z^{\prime}=\Phi\left(z\right)$, we obtain Assumption A(A7).2.which is also implied by (62) if $\left|\nabla\Phi\left(z\right)\right|=1$ for all $z$. ∎

Appendix B Weak convergence of discrete-time BPS

for the infinitesimal generator of BPS where the domain will be discussed later on.

For any $\epsilon>0$, we write $K^{(\epsilon)}$ for the transition kernel of the discrete-time BPS, DPBS, with step size $\epsilon>0$. This kernel satisfies

To keep notation reasonably compact we will often write

with $p^{(i)}(x,v)$ for $i=1,2,3$ the probabilities appearing in (72).

We will write $\{Z^{(\epsilon)}(k);k\geq 0\}$ for the Markov chain generated by the transition kernel $K^{(\epsilon)}$, with $Z^{(\epsilon)}(0)\sim\rho$. We also define the càdlàg process $\{\zeta_{t}^{(\epsilon)}:t\in[0,\infty)\}$, through

For any $z=(x,v)\in\mathcal{Z}$ the function $t\mapsto\lambda(x+tv,v)$ is continuous.

where $\|\Delta f\|$ denotes the operator norm of the Hessian matrix of $\pi$.

and for some $\varepsilon,\delta>0$ we have for all $|v|=1$ and $r<\varepsilon$

Let $(\epsilon_{n};n\geq 1)$ be a positive sequence such that $\epsilon_{n}\to 0$ as $n\rightarrow\infty$. Under Assumptions 1, 2, 3 and 4 the law of $\{\zeta^{\left(\epsilon_{n}\right)}(\cdot)\}$ converges weakly to that of BPS as probability measures on $D_{\mathcal{Z}}[0,\infty)$ as $n\rightarrow\infty$.

Before we embark on the proof of Theorem 12 we prove some useful properties for the semigroup and the generator.

B.2 The Feller property

for all $t\geq 0$ and $f\in C_{0}(\mathcal{Z})$ we have $P^{t}f\in C_{0}(\mathcal{Z})$, and

$P^{t}f(z)\to f(z)$ as $t\to 0$ for $f\in C_{0}(\mathcal{Z})$ and $z\in\mathcal{Z}$.

Let Assumptions 1 and 2 hold. Then $\left\{P^{t};t\geq 0\right\}$ is a Feller semigroup and the martingale problem for $(\mathcal{L},\rho)$ admits a unique solution.

First we prove the uniqueness for the martingale problem assuming the Feller property. Then we will prove the Feller property.

Since the semigroup $\{P^{t}:t\geq 0\}$ is Feller it follows from [28, Theorem 19.6] that the semigroup is also strongly continuous, whence by the Hille-Yosida Theorem (see for example [19, Theorem 1.2.6]) if follows that $\mathcal{L}$ is dissipative, that is for any $f\in C_{0}(\mathcal{Z})$ we have

To complete the proof we now show that BPS is Feller. For $t\geq 0$ and $z=(x,v)\in\mathcal{Z}$, write $\Phi_{t}(z)=(x+tv,v)$. To prove (F2) notice that for any $f\in C_{0}(\mathcal{Z})$ and $z=(x,v)\in\mathcal{Z}$ we have

where it is clear that for any $z=(x,v)$ we have

as $t\downarrow 0$ and (F2) follows easily by continuity of $f$.

where $T_{1},T_{2},\dots$ are the event times of BPS. We can also write

Both integrals vanish by bounded convergence, since by continuity of $\lambda$ and $\phi$ the second integrand vanishes pointwise, while both integrands are bounded by boundedness of the flow $\Phi_{s}(z):s\in[0,t]\}$ and continuity of $\lambda$. On the other hand letting

We have thus shown that $Gg$ defined in (76) is continuous. In addition since $g$ is bounded, it follows that

Finally, recall from the proof of [15, Lemma 9.3] that

as $n\to\infty$, where $T_{n}$ is the time of $n$-th event, when BPS starts from $z$. Suppose now that $z=(x,v)\in B_{R^{\prime}}(0)\subset\mathcal{Z}$, the ball of radius $R^{\prime}$ around the origin. Then by construction of BPS there will be a compact set $K_{R^{\prime}}\subset\mathcal{Z}$, such that $\left\{Z_{s}=\left(X_{s},V_{s}\right);0\leq s\leq t\right\}\subset K_{R^{\prime}}$. Therefore, since $\lambda$ is locally bounded, we have that

as $n\to\infty$. It follows that $G^{n}g(t,z)\to P^{t}f(z)$ as $n\to\infty$ uniformly on compact sets. Since the functions $z\mapsto G^{n}g(t,z)$ are continuous, it follows that $P^{t}f(z)$ is continuous on every compact set and thus is continuous.

Then choose $K^{\prime}>K+Rt$ and $N$, such that for all $n\geq N$ we have $\left|x^{(n)}\right|\geq K^{\prime}$. Then since $X_{t}\in B\left(x^{(n)},Rt\right)$ it follows that $|X_{t}|>K$ and thus $|f(X_{t},V_{t})|\leq\epsilon$. Since $\epsilon>0$ is arbitrary the result follows. ∎

B.3 Preliminary calculations

We first need precise estimates for $p^{(i)}(x,v)$, $i=2,3$ and small $\epsilon$. We will often use the formula

Therefore if $\langle U(x),v\rangle\leq 0$, we will have $p_{\epsilon}^{(2)}(x,v)=0$, for all $\epsilon$ small enough. Thus we can assume that $\langle U(x),v\rangle>0$ in which case we also have $\langle U(x),-R(x)v\rangle>0$ and therefore for all $\epsilon>0$ small enough we have that $U(x+\epsilon v),U(x-R(x)v\epsilon)>U(x)$ and thus

Since we have assumed that for $\epsilon$ small enough we have $\pi(x-sR(x)v)<\pi(x)$, then we have for $\epsilon$ small enough

where we used Assumption 4. Overall we have that

B.4 Proof of Theorem 12

Let $\epsilon_{n}\to 0$. To ease notation we will write $\zeta^{(n)}$ rather than $\zeta^{(\epsilon_{n})}$. Define (see [19, Remark 8.3(b)])

where $\mathcal{G}_{t}^{n}:=\sigma\left(\zeta^{(n)}(s):s\leq t\right)$, the natural filtration of $\left\{\zeta^{(n)}(t):t\geq 0\right\}$. Recall that $\zeta^{(n)}(0)\sim\rho$ for all $n$.

To apply [19, Corollary 8.15 of Chapter 4] we need to check the following:

Compact Containment: For every $\eta>0$ and $T>0$ there is a compact set $\rho_{\eta,T}\subset\mathcal{Z}$ such that