A Discrete Bouncy Particle Sampler

Introduction

Markov Chain Monte Carlo (MCMC) algorithms provide Monte Carlo approximations to expectations with respect to a given probability distribution, $\pi$, via an ergodic Markov chain whose invariant distribution is $\pi$. Non-reversible Markov Chaine Monte Carlo samplers, of which the Hamiltonian Monte Carlo algorithm (Duane et al. 1987) is perhaps one of the most successful and widely-used examples, are known to enjoy desirable mixing properties in several contexts. Indeed, several theoretical results quantify the advantages of non-reversible samplers. For example, Diaconis et al. 2000 obtains rates of convergence for a non-reversible version of the random walk algorithm; subsequently and inspired by Diaconis et al. 2000, Chen et al. 1999 describes the best acceleration achievable through the idea of lifting. On a different note, Hwang et al. 2015, Lelièvre et al. 2013, Rey-Bellet & Spiliopoulos 2015 and Duncan et al. 2016 investigate and quantify the advantages offered by leveraging (a discretization of) a non-reversible diffusion process for computing Monte-Carlo averages: in many settings, it can be proved that the standard reversible Langevin dynamics is the worst in terms of asymptotic variances among a large class of diffusion processes that are ergodic with respect to a given target distribution. More recently, different designs of non-reversible MCMC sampler have been proposed. The Zig-Zag sampler, an instance of the large class of Piecewise-Deterministic-Markov-Processes, first obtained as a scaling limit of a lifted Metropolis–Hastings Markov chain (Bierkens & Roberts 2017), is a continuous-time non-reversible Markov process that can be used for computing ergodic averages, and can be used for efficiently exploring Bayesian posterior distributions in the Big-Data regime (Bierkens et al. 2019). Inspired from the physics literature (Peters & de With 2012), the Bouncy Particle Sampler (Bouchard-Côté et al. 2017) is another continuous-time non-reversible Monte-Carlo sampler that demonstrates state-of-the-art performance when used to explore certain Bayesian posterior distributions. Fearnhead et al. 2018 reviews the Zig-Zag sampler and the Bouncy Particle Sampler and describes some extensions. The Bouncy Particle Sampler or the Zig-Zag Sampler requires more than simple point-wise evaluations of the log-target density and its gradient: one typically needs local upper bounds on derivatives of the log-target density. Unfortunately, those bounds are unavailable or difficult to compute in many applied situations. Consequently, such continuous-time samplers cannot be directly used in these settings. This article presents a discrete-time MCMC sampler, inspired by the Bouncy Particle Sampler, that can be implemented when only point-wise evaluations of the target-density and its gradient are available.

Our algorithm, the Discrete Bouncy Particle Sampler, is described in detail in Section 2.1. It extends the statespace from a position to a position and a direction in the same way that the Zig-Zag and Bouncy Particle samplers do; however, our algorithm operates in discrete time and is based upon the guided random walk of Gustafson 1998. The guided random walk combines two reversible kernels to create a non-reversible kernel which heads in a specific direction until a rejection occurs, at which point it reverses direction. Our key addition is a particular delayed-rejection proposal (Tierney & Mira 1999), used after any initial rejection, potentially avoiding many inefficient direction reversals. The delayed-rejection move is analogous to the bounce in the Bouncy Particle Sampler and we show that the Discrete Bouncy Particle Sampler can be viewed as a time discretization of the Bouncy Particle Sampler. An alternative discretization of the Bouncy Particle Sampler, based on the reflective slice sampler (Neal 2003) is described and extended in the independent work of Vanetti et al. 2017. Importantly, several interesting extensions have recently been proposed (Vanetti et al. 2017; Wu & Robert 2017; Wu & Robert 2020; Monmarché 2019; Michel et al. 2020) to scale and enhance this class of Piecewise-Deterministic-Markov-Processes MCMC samplers.

As with the Bouncy Particle Sampler, our algorithm can be reducible. To solve this issue, we perturb the direction vector at the end of every iteration. The size of this per-iteration perturbation has a substantial impact on the performance of the algorithm, in a similar way to the occasional, complete direction refresh of the Bouncy Particle Sampler (Bierkens et al. 2018). Analogously to the independent investigations for the Bouncy Particle Sampler in Bierkens et al. 2018, for the Discrete Bouncy Particle Sampler exploring a Gaussian target we use diffusion-approximation arguments to describe the limit of the radial process as dimension increases to infinity. We then leverage this to obtain the theoretical efficiency of the Discrete Bouncy Particle Sampler as a function of the partial-refreshment parameter, $\kappa$. This leads to a simple and robust tuning criterion that is described in Section 3.2. In more practical developments, we show that a surrogate may be substituted for the gradient of the target density when it is computationally expensive or impossible to obtain. Our final contribution is a construction that allows the user to choose to only calculate a fixed number of orthogonal components of the gradient vector.

The Discrete Bouncy Particle Sampler

The Discrete Bouncy Particle Sampler operates on the extended state space $\mathcal{X}\times\mathcal{S}$, where $\mathcal{S}\subseteq\mathcal{X}$, and explores the extended target distribution

where $\rho(\mathsf{d}u)$ is an auxiliary spherically symmetric distribution with support $\mathcal{S}\subset\mathcal{X}$. Section 2.2 describes several standard choices of auxiliary distributions. Henceforth, we will refer to the variable $x\in\mathcal{X}$ as the position of a particle and the variable $u\in\mathcal{S}$ as its direction. The bounce after which the Discrete Bouncy Particle Sampler is named enters through the operator $u\mapsto\mathscr{R}_{v}(u)$ that reflects the vector $u\in\mathcal{S}$ with respect to the hyperplane orthogonal to the vector $v\in\mathcal{X}\setminus\{0\}$,

For any vector $v\in\mathcal{X}\setminus\{0\}$, the reflection operator $\mathscr{R}_{v}$ is an involution $\mathscr{R}_{v}\circ\mathscr{R}_{v}(u)=u$ that preserves norms. This remark underlies the proof of correctness of the Discrete Bouncy particle Sampler whose details are presented in the Supplementary Material. As discussed in Michel et al. 2020, it is possible to rely on more general reflection operators. Most of the methods developed in this text extends to these variants, although we concentrate on Equation (1) for ease of exposition. The Discrete Bouncy Particle Sampler relies on a non-vanishing vector field

In practice, this vector field is either chosen as $\mathscr{F}(x)=\nabla\log\pi(x)$, replaced by an arbitrary modification when the gradient vanishes, or as an approximation of it, as described in Section 2.4. The quantity $\mathscr{R}_{\mathscr{F}(x)}(u)$ represents the resulting direction when a particle with incoming direction $u$ performs an elastic bounce off the hyperplane orthogonal to the vector $\mathscr{F}(x)$.

The Discrete Bouncy Particle Sampler deterministically cycles through two Markovian transitions that leave the extended target distribution $\widetilde{\pi}$ invariant, (1) a Position Update, with a possible Direction Reflection (2) a Direction Refreshment. The resulting scheme is, in general, non-reversible. For the direction refreshment, one can choose any Markov transition kernel $\mathsf{P}_{\rho}(u,\mathsf{d}u^{\prime})$ that leaves the auxiliary distribution $\rho$ invariant. For a discretization parameter $\delta>0$, and a current state $(x_{k},u_{k})\in\mathcal{X}\times\mathcal{S}$, the algorithm proceeds as follows.

Position Update: Generate a proposal $(x^{\prime},u^{\prime})=(x_{k}+\delta\,u_{k},u_{k})$. With position update probability

set $(\hat{x}_{k},\hat{u}_{k})=(x^{\prime},u^{\prime})$ and go to Step 3. Otherwise, proceed to Step 2.

Direction Reflection: consider $u^{\prime\prime}=\mathscr{R}_{\mathscr{F}(x^{\prime})}(u^{\prime})$ and $x^{\prime\prime}=x^{\prime}+\delta\,u^{\prime\prime}$. With direction reflection probability

set $(\hat{x}_{k},\hat{u}_{k})=(x^{\prime\prime},u^{\prime\prime})$. Otherwise, negate the direction by setting $(\hat{x}_{k},\hat{u}_{k})=(x_{k},-u_{k})$.

Direction Refreshment: Set $(x_{k+1},u_{k+1})=(\hat{x}_{k},U)$ where $\operatorname{P}(U\in A)=\mathsf{P}_{\rho}(\hat{u}_{k},A)$.

By construction, the Direction Refreshment step preserves the extended target distribution. The fact that the combination of the Position Update and Direction Reflection steps also preserves the extended distribution is discussed in the Supplementary Material. This can be understood as a slight generalization of the standard delayed rejection mechanism (Tierney & Mira 1999) when applied to a deterministic and volume preserving proposal. A similar scheme was proposed independently in Vanetti et al. 2017. Furthermore, and importantly for Section 2.4, the algorithm remains valid if the deterministic vector field is replaced by a randomized version of it. The proof of correctness is identical to the deterministic case and is briefly discussed in the Supplementary Material.

A simple thought experiment, such as considering a target density $\pi$ with spherically-symmetric contours and $\mathsf{P}_{\rho}(u,\mathsf{d}u^{\prime})=\delta_{u}(\mathsf{d}u^{\prime})$, shows that, as with the Bouncy Particle Sampler, the Discrete Bouncy Particle Sampler can be reducible. The choice and tuning of the direction refreshment operator $\mathsf{P}_{\rho}$ is consequently important in practice and is discussed at length in the sequel.

2 Direction dynamics

Full Refresh: for $\rho=\rho_{S}$ or $\rho=\rho_{G}$ and an update rate $\kappa>0$ the Markov process with generator $\mathcal{L}^{(V)}\varphi(u)=\kappa\,{\left(\rho(\varphi)-\varphi(u)\right)}$ completely refreshes the direction at rate $\kappa$ and has a mixing time of $\mathcal{O}(1/\kappa)$. For a time discretization parameter $0<\delta<1$, set

where $\xi\sim\rho$ and $B_{\kappa}^{\delta}$ is a Bernoulli random variable with $\operatorname{P}(B_{\kappa}^{\delta}=1)=\exp(-\kappa\,\delta)$.

Ornstein-Uhlenbeck refresh: for $\rho=\rho_{G}$ and an update rate $\kappa>0$, the Ornstein-Uhlenbeck process $\mathsf{d}V_{t}=-(\kappa/2)\,V_{t}\,\mathsf{d}t+(\kappa/d)^{1/2}\,\mathsf{d}W$ leaves $\rho_{G}$ invariant and has a mixing time of $\mathcal{O}(1/\kappa)$. Set

for $\xi\sim\rho_{G}$ and $\alpha=\exp(-\kappa\,\delta/2)$.

for $\xi\sim\rho_{G}$ and $\alpha=\exp(-\kappa\,\delta/2)$.

3 Continuous-time limit

for any integer $0\leq k\leq T/\delta$ and linearly interpolation in between. The main result of this section, Proposition 1, whose proof is presented in the Supplementary Material, relies on the following regularity assumptions.

The function $x\mapsto\log\pi(x)$ is twice differentiable with a bounded second derivative.

The vector field $\mathscr{F}:\mathcal{X}\to\mathcal{X}\setminus\{0\}$ is continuous.

There exists a continuous time Markov process $\{V_{t}\}_{t\geq 0}$ with generator $\mathcal{L}^{(V)}$ such that, for any time discretization parameter $\delta>0$, the transition kernel $\mathsf{P}^{\delta}_{\rho}$ describes the transition of the Markov process $V$ in the sense that $\mathsf{P}^{\delta}_{\rho}(u,\mathsf{d}u^{\prime})=\operatorname{P}(V_{t+\delta}\in\mathsf{d}u^{\prime}|V_{t}=u)$. We assume that the trajectories of the Markov process $\{V_{t}\}_{t\geq 0}$ are almost surely continuous.

Let Assumptions A(1-2-3) hold and consider a fixed time horizon $T>0$. As $\delta\to 0$, the sequence of continuous time processes $\overline{z}^{\delta}_{t}=(\overline{x}^{\delta}_{t},\overline{u}^{\delta}_{t})$ converges weakly in the Skorokhod topology to the bivariate Markov process $\overline{Z}_{t}=(\overline{X}_{t},\overline{U}_{t})$ with generator

with rate $\lambda(x,u)\equiv\left\langle-\nabla\log\pi(x),u\right\rangle_{+}$ and acceptance probability

The limiting Markov process $\overline{Z}_{t}=(\overline{X}_{t},\overline{U}_{t})$ with generator (6) evolves according to the dynamics

in between events that arrive at rate $\lambda(\overline{X}_{t},\overline{U}_{t})$. When such an event is triggered, the direction is reflected, i.e. $\overline{U}_{t}=\mathscr{R}_{\mathscr{F}(X_{t^{-}})}(\overline{U}_{t^{-}})$, with probability $\mathcal{A}(\overline{X}_{t^{-}},\overline{U}_{t^{-}})$, and completely reversed, i.e. $\overline{U}_{t}=-\overline{U}_{t^{-}}$, with probability $1-\mathcal{A}(\overline{X}_{t^{-}},\overline{U}_{t^{-}})$. Possible choices of Markovian dynamics with generator $\mathcal{L}^{(V)}$ in Assumption A2 are detailed in Section 2.2. In the case when $\mathscr{F}(x)=\nabla\log\pi(x)$ and $\mathcal{L}^{(V)}\varphi(u)=\kappa\,{\left\{\rho(\varphi)-\varphi(u)\right\}}$, for a fixed refreshment rate $\kappa>0$, the limiting Markov process is the standard Bouncy Particle Sampler (Bouchard-Côté et al. 2017). The interested reader is referred to Vanetti et al. 2017; Wu & Robert 2017 for other interesting generalizations

In order to understand the influence of the vector field $\mathscr{F}$, it is instructive to study the limiting acceptance probability (7). The limiting process $\overline{Z}$ is rejection free, i.e. never backtracks, if for any $(x,u)\in\mathcal{X}\times\mathcal{S}$ we have that $\lambda\{x,-\mathscr{R}_{\mathscr{F}(x)}(u)\}=\lambda(x,u)$. It is readily seen that this condition is equivalent to choosing $\mathscr{F}(x)$ proportional to $\nabla\log\pi(x)$ for any $x\in\mathcal{X}$ where this quantity does not vanish. In other words, any other choice of vector field $\mathscr{F}$ leads to a limiting process that is not rejection-free. Section 2.4 describes ways to efficiently approximate this optimal choice when evaluating $\nabla\log\pi(x)$ is not computationally efficient.

4 Approximate reflections

As described in the previous section, vector fields that lead to a rejection-free algorithm in the limit $\delta\to 0$ are such that $\mathscr{F}(x)$ is proportional to $\nabla\log\pi(x)$ for all $x\in\mathcal{X}$ where this quantity is non-zero. When computing the gradient of the log-density is not computationally feasible, one can instead use a vector field $\mathscr{F}$ that only approximates $\nabla\log\pi$, necessarily paying the price of a non-zero probability for the limiting algorithm to backtrack. Another strategy, similar to the one presented in Fielding et al. 2011, consists in choosing $\mathscr{F}(x)$ as the gradient of an approximate surrogate target distribution.

One can also completely reflect the component of $u\in\mathcal{S}$ that is orthogonal to the plane $V(\underline{\zeta})$. In other words, the updated direction $u^{\prime}$ is defined as

While both reflection operators are valid, we have empirically found that the operator (9) leads to better mixing properties.

5 Preconditioning

Algorithm tuning through diffusion approximation

Note that, in the definition of the process $R_{t}^{d}$, time has been accelerated by a factor of $d$. Proposition 2 stated below shows that, in order to observe a non-degenerate scaling limit as $d\to\infty$, this acceleration factor is the correct one. As will be demonstrated, the mixing properties of $R^{d}_{t}$ are closely related to the mixing properties of the scalar jump-diffusion $\{\theta^{\kappa}_{t}\}_{t\geq 0}$ with generator

The operator $\mathcal{L}^{(K)}$ is the generator of an Ornstein-Uhlenbeck process that is reversible with the standard Gaussian density. Similarly, $\mathcal{L}^{(J)}$ is the generator of the Markov process with unit drift and reflections $\theta\mapsto-\theta$ that occur at rate $\theta_{+}=\max(0,\theta)$. It can readily be checked that this process also leaves the standard Gaussian distribution invariant. Combining these two facts show that the process $\theta^{\kappa}_{t}$ also leaves the standard Gaussian distribution invariant. We denote by $V_{\sigma}(\kappa)$ the asymptotic variance of ergodic averages along $\theta^{\kappa}_{t}$ defined as

There is no closed form expression for the quantity $V_{\sigma}(\kappa)$ but it can easily be approximated numerically, as displayed in Figure 1. Note that $V_{\sigma}(\kappa)\to 0$ as $\kappa\to 0$ and $\kappa\to\infty$. Proposition 2, whose proof can be found in the Supplementary Material, shows that the asymptotic variance $V_{\sigma}(\kappa)$ dictates the mixing rate of the log-target process. The higher the asymptotic variance $V_{\sigma}(\kappa)$, the faster the mixing of the radial process.

The velocity function $V_{\sigma}(\kappa)$ is defined in Equation (12).

The process (13) is an Ornstein-Uhlenbeck that is reversible with respect to the centred Gaussian distribution with variance $\sigma^{2}/2$. Since the Ornstein-Uhlenbeck (13) has a mixing time of order $\mathcal{O}(1)$, this indicates that in the high-dimensional regime $d\to\infty$ and $\delta\to 0$, one can expect (Roberts & Rosenthal 2016) the log-target process to mix on a time scale of order $\mathcal{O}(d/\delta)$. When implementing the Discrete Bouncy Particle Sampler in practice, the parameter $\delta$ should be chosen small enough to guarantee that the acceptance rate remains high-enough, but not smaller. Similar guidelines for the Hamiltonian Monte-Carlo method are described in Betancourt et al. 2014. The optimal tuning of the parameter $\delta$, and study of its dependence with respect to the dimensionality of the target distribution, is beyond the scope of this article. Instead, we concentrate on the tuning of the refreshment parameter $\kappa$. When optimising the mixing of the log-target process, we observe empirically that the tuning of the parameter $\kappa$ is insensitive to the value of $\delta$. This is in part because whatever the value of $\delta$, as long as it is sufficiently small, the log-target process is an approximation of the limiting diffusion (13) (see also Section 3.2); empirical evidence that this insensitivity continues to hold for large $\delta$ is provided in Section 4.1.

2 Tuning of the refreshment parameter κ\kappa

The limiting diffusion (13) obtained in Section 3.1 indicates that, for an isotropic Gaussian target distribution with marginal variance $\sigma^{2}$ and in the regime $d\to\infty$, optimising the efficiency of the Discrete Bouncy Particle Sampler is achieved by choosing a refreshment parameter $\kappa>0$ that maximizes the velocity $V_{\sigma}(\kappa)$. In other words, for a given marginal variance $\sigma^{2}>0$, the optimal refreshment rate $\kappa_{\star}(\sigma)$ is given by

Furthermore, a change of time argument immediately shows that $\kappa_{\star}(\sigma)=\kappa_{\star}/\sigma$ with $\kappa_{\star}\equiv\kappa_{\star}(\sigma=1)$. In practice, the variance parameter $\sigma^{2}>0$ is not known so that the optimal refreshment parameter is not directly accessible. To make progress, denote by $\{\tau_{j}\}_{j\geq 1}$ the (strictly increasing) sequence of time indices at which Direction Reflection events are attempted (and always accepted in the Gaussian setting). We denote by $u_{\tau^{-}_{j}}\in\mathcal{S}$ the direction right before a reflection event, and by $u_{\tau^{+}_{j}}\in\mathcal{S}$ the direction right after the reflection. For tuning purposes, we propose to monitor the dot product $\beta\in$ between the direction vectors right after and before the Direction Reflection attempts,

For a Discrete Bouncy Particle Sampler evolving at stationarity, as $\delta\to 0$ and for any fixed dimension $d\geq 2$, consider the distribution $\mu^{d}(\mathsf{d}\beta;\kappa,\sigma)$ of these dot products. One can readily check that if $\{x^{d,\delta}_{k},u^{d,\delta}_{k}\}_{k\geq 0}$ is Discrete Bouncy Particle Sampler chain with parameters $\kappa,\delta>0$ exploring the centred $d$-dimensional Gaussian with marginal standard deviation $\sigma$ then, for any scaling factor $s>0$, the Markov chain defined as $\{s\times x^{d,\delta}_{k},u^{d,\delta}_{k}\}_{k\geq 0}$ is also Discrete Bouncy Particle Sampler chain, with refreshment parameter $\kappa/s$ and time discretization parameter $s\,\delta>0$, exploring the centred $d$-dimensional Gaussian with marginal standard deviation $s\,\sigma$. It follows that $\mu^{d}(\mathsf{d}\beta,\kappa,\sigma)=\mu^{d}(\mathsf{d}\beta;\kappa/s,s\,\sigma)$ for any scaling factor $s>0$. Consequently, since $\kappa_{\star}(\sigma)=\kappa_{\star}/\sigma$, the distribution $\mu^{d}(\mathsf{d}\beta;\kappa_{\star}(\sigma),\sigma)\equiv\mu_{\star}^{d}(\mathsf{d}\beta)$ does not depend on the standard deviation $\sigma$. It is straightforward to numerically estimate the average dot product at optimality,

Figure 1 illustrates this optimality result. Very low values $\beta\ll\beta_{\star}$ indicate that the directions are updated too frequently, leading to an inefficient random-walk behaviour. High values $\beta\approx 1$ indicate that the directions are not updated frequently enough, leading to an inefficient exploration of the state space. The case $\beta=1$ corresponds to the case when the direction are not updated at all, which is known in the Gaussian case to lead to a reducible Markov Chain with an incorrect invariant distribution. For tuning the refreshment parameter $\kappa>0$ of a general Discrete Bouncy Particle Sampler, we consequently propose to estimate empirically the expectation at stationarity of the quantity $\beta$ in (14). Let $\{\tau_{j}\}_{j\geq 1}$ now denote the realization of the sequence of indices at which direction reflection attempts occur. A direction reflection attempt is accepted with probability (3), otherwise the direction is negated. We define

and choose $\kappa>0$ so that this quantity approximately equals its optimal value $\beta_{\star}\approx 0.2$. Just as with the estimation of acceptance rates when tuning the scaling of various algorithms (Roberts & Rosenthal 2001), and unlike the Effective Sample Size itself that is notoriously difficult to reliably estimate, the mean dot product can be estimated accurately from short MCMC runs. Importantly, and as described in Section 4, we have found this tuning procedure to be robust with respect to departure from Gaussianity and to approximately hold in non-isotropic and relatively low-dimensional settings with $\delta\gg 0$.

3 Non-isotropic target and non-zero δ\delta

The diffusion limit in Section 3.1 was obtained as $\delta\downarrow 0$ and for an isotropic Gaussian target where the direction reflection proposals are always accepted. In this Section, we consider the non-isotropic case of a $d$-dimensional target distribution defined as

where $\Phi(t)=(2\pi)^{-1/2}\,\int_{0}^{t}e^{-t^{2}/2}\,dt$ is the standard Gaussian cumulative function.

Simulation studies

2 Robustness of advice to departures from Gaussianity and isotropy

We next investigate the robustness of our tuning advice to departures of the target from Gaussianity and isotropy.

Consider three scenarios: an isotropic multivariate logistic distribution $\pi^{\textrm{logis}}$ with density $\prod_{i=1}^{d}\,\exp(x_{i})/[1+\exp(x_{i})]^{2}$, an isotropic Gaussian distribution $\pi^{\textrm{iso}}$ with density proportional to $\prod_{i=1}^{d}\,\exp(-x_{i}^{2}/2)$ and a non-isotropic multivariate Gaussian distribution $\pi^{\textrm{aniso}}$ with density proportional to $\prod_{i=1}^{d}\,\exp\{-x_{i}^{2}/(2\,\sigma_{i}^{2})\}$. In this section, we choose $d=100$ for both isotropic targets, and $d\in\{20,50,200\}$ for the anisotropic target. In order to test the robustness of our tuning guidelines to non-isotropic distributions, we chose the scales $\sigma_{1}<\ldots<\sigma_{d}$ linearly separated between $\sigma_{1}=1$ and $\sigma_{d}=10$. The scaled effective sample size curves as a function of the mean dot-product are displayed in Figure 3. There is extremely good agreement with the theory developed in Section 4 for the isotropic distribution $\pi^{\textrm{iso}}$ and the approximately isotropic distribution $\pi^{\textrm{logistic}}$. Not surprisingly, mild departure from the theory is observed for strongly non-isotropic distributions such as $\pi^{\textrm{aniso}}$. However, for dimension $d=200$ the departure, especially in terms of the optimal dot product, is barely noticeable. When $d=20$ and $d=50$, however, our proposed guideline, i.e. tune the refreshment rate $\kappa$ such that the mean dot-product $\beta\approx 0.2$, leads only to a loss of efficiency of approximately $10\%$ and $5\%$ respectively.

3 Convergence and tail behaviour

One of the most commonly used algorithms for inference in high-dimensional scenarios is Hamiltonian Monte Carlo Duane et al. 1987. However, it is well known (Livingstone et al. 2019) that due to the dependence of the Leapfrog step on $\|\nabla\log\pi\|$, Hamiltonian Monte Carlo is not geometrically ergodic on targets with tails lighter than those of a Gaussian. By contrast, the Discrete Bouncy Particle Sampler depends on $\nabla\log\pi$ only through the equivalent normalized vector. The Supplementary Material details a simulation study on a non-isotropic, light-tailed target where a tuned Hamiltonian Monte Carlo algorithm is nearly $50\%$ more efficient than a tuned discrete bouncy particle sampler when started from stationarity. However, with the same tunings, but when started from a random point in the tail of the distribution, Hamiltonian Monte Carlo does not even move, whereas the discrete bouncy particle sampler quickly converges to the centre of the posterior mass.

4 The Markov modulated Poisson process

A simulation study on the eight-dimensional posterior of a non-trivial statistical model, the Markov modulated Poisson process, is detailed in the Supplementary Material. We find that mixing efficiency of $\log\pi$ is optimized at a dot product statistic of $\beta\approx 0.4$, tuning to $\beta\approx 0.2$ would lead to only a $10\%$ reduction in efficiency. When either preconditioning or using only $n_{\text{cpt}}=3$ random components of the gradient vector, the optimal efficiency is achieved for a dot product statistics of $\beta\approx 0.2$.

Discussion

The key advantage of the Discrete Bouncy Particle Sampler over its continuous-time counterpart is that posterior and gradient evaluations can be treated as a “black box” with no requirement to bound the gradient so as to apply Poisson thinning. Unlike the continuous-time algorithm, there is a chance that an attempted bounce will be rejected and the particle will approximately backtrack, however for sensible tunings, the back-tracking probability converges towards zero as the dimension of the problem increases.

The average computational cost per iteration is insensitive to the choice of $\kappa$ since the direction is updated every iteration, and $\kappa$ has little effect on the number of steps between potential bounces. Thus it is sufficient for the purposes of this article to describe and empirically record efficiency in terms of effective sample size rather than effective samples per unit of time. The only exception to this is when we compare against Hamiltonian Monte Carlo in Appendix C.1.

The dot-product tuning diagnostic maps $\kappa$ to an absolute scale via the properties of the posterior, just as the acceptance rate diagnostic does for the scaling in the random-walk Metropolis algorithm. When $\kappa=0$ the velocity direction just before the next bounce is identical to that just after the previous bounce, and the dot product is unity. The mixing time of the refreshment process is $1/\kappa$; when this is small compared with the time between bounces or, equivalently, the length scale of the target, the two velocity directions bear little relation to each other, and the dot product is small.

Whilst this article offers theory-based practical advice on the tuning of the refreshment parameter, $\kappa$, it does not tackle the choice of the discretization parameter, $\delta$. In contrast to the insensitivity of computational cost to the choice of $\kappa$, increasing $\delta$ increases the frequency of potential bounces and hence of expensive gradient calculations, and this would need to be accounted for in any analysis.

Acknowledgements

Work by CS was supported by EPSRC grant EP/P033075/1. AHT acknowledges support from the Singapore Ministry of Education Tier 2 Grant (MOE2016-T2-2-135) and a Young Investigator Award Grant (NUSYIA FY16 P16; R-155-000-180-133).

References

Appendix A Correctness and proofs of propositions

The mapping $u\mapsto\mathsf{R}(u,x)$ is volume preserving.

The mapping $u\mapsto\mathsf{R}(u,x)$ preserves norms, $\|\mathsf{R}(u,x)\|=\|u\|$.

Generate a proposal $(x^{\prime},u^{\prime})=(x_{k}+\delta\,u_{k},-u_{k})$. With probability

set $(\hat{x}_{k},\hat{u}_{k})=(x^{\prime},-u^{\prime})$ and go to Step 3. Otherwise, proceed to Step 2.

consider $u^{\prime\prime}=-\mathsf{R}(u,x^{\prime})$ and $x^{\prime\prime}=x^{\prime}-\delta\,u^{\prime\prime}$. With probability

set $(\hat{x}_{k},\hat{u}_{k})=(x^{\prime\prime},u^{\prime\prime})$. Otherwise, set $(\hat{x}_{k},\hat{u}_{k})=(x_{k},u_{k})$.

Reverse the direction: $(x_{k+1},u_{k+1})=(\hat{x}_{k},-\hat{u}_{k})$

Consider any spherically symmetric probability density $\rho(u)$. Under Assumptions B1-2-3, the Markov kernel described by Step 1-2-3 leaves the density $\widetilde{\pi}(x,u)=\pi(x)\,\rho(u)$ invariant.

with $\alpha_{1}(z)=1\wedge\mu\{T_{1}(z)\}/\mu(z)$ and $\alpha_{3}(z)=1-\alpha_{1}(z)-\{1-\alpha_{1}(z)\}\,\alpha_{2}(z)$ and

Algebra shows that this is equivalent to proving that

Since $T_{1}$ and $T_{2}$ are involutions that preserve volume, a change of variable $z\mapsto T_{1}(z)$ shows that the first integral in Equation (17) also equals its negation, and hence vanishes. And similarly, the change of variable $z\mapsto T_{2}(z)$ shows that the second integral in Equation (17) also vanishes. This concludes the proof of the lemma. ∎

In Section 2.1, the combination of the Position Update and Direction Update is equivalent to Step 1-2-3 with the operator $\mathsf{R}(u,x)=\mathscr{R}_{\mathscr{F}(x)}(u)$. Since algebra shows that the Conditions B1-2-3 are satisfied, Lemma 1 thus shows the correctness of the Discrete Bouncy Particle Sampler as described in Section 2.1.

A.2 Proof of Proposition 1

Recall that the quantity $\lambda(x,u)$ is defined as $\lambda(x,u)\equiv\left\langle-\nabla\log\pi(x),u\right\rangle_{+}$. Under Assumption (A1) and a discretization parameter $\delta>0$, the acceptance probability $\alpha^{\delta}(x,u)$that the proposal $(x,u)\mapsto(x+\delta\,u,u)$ is accepted reads

where $\mathcal{O}(\delta^{2})$ is a quantity whose absolute value is less than a constant times $\delta^{2}$. For $t>0$, the probability that the Discrete Bouncy Particle Sampler algorithm accepts $\lfloor t/\delta\rfloor+1$ consecutive proposals $(x,u)\mapsto(x+\delta\,u,u)$ without reflection attenpts equals $\prod_{k=0}^{\lfloor t/\delta\rfloor}\alpha^{\delta}(x^{\delta}_{k},u^{\delta}_{k})$. Under Assumption (A3), one can condition upon a fixed trajectory of the Markov process $V$, i.e. $V_{t}=v_{t}$ for all $0\leq t\leq T$ and $u^{\delta}_{k}=\overline{u}^{\delta}_{k\delta}=v_{k\delta}$, not depending on the parameter $\delta$, so that $\overline{x}^{\delta}_{k\delta}=\overline{x}^{\delta}_{0}+\delta\sum_{j=0}^{k-1}v_{j\delta}$. Equation (18), the continuity of the rate function $\lambda$ as well as the continuity of the trajectories of the Markov process $V$, show that

where $\overline{x}_{s}=x_{0}+\int_{0}^{s}v_{t}\,dt$. This means that, in the limit $\delta\to 0$, bounce attempts arrive at rate $\lambda(x,u)$ and, in between the bounces, the limiting process simply evolves according to the dynamics (8).

Finally, once a proposal $(x,u)\mapsto(x+\delta\,u,u)$ is rejected, the second proposal $(x,u)\mapsto(x^{\prime\prime},u^{\prime\prime})$, i.e. the bounce, is accepted with probability $\alpha_{\textrm{DR}}(x,u)$ described in Equation (3). By continuity of the density $x\mapsto\pi(x)$, we have that $\pi(x^{\prime\prime})/\pi(x)\to 1$ as $\delta\to 0$. Furthermore, the Taylor expansion (18) gives that $1-\alpha^{\delta}(x,u)=\delta\,\lambda(x,u)+\mathcal{O}(\delta^{2})$. Under Assumption, the vector field $\mathscr{F}$ is continuous, which implies that

where we have dropped the dependence on $\delta$ from the notation $x^{\prime\prime}=x+\delta\,u+\delta\,\mathscr{R}_{\mathscr{F}(x+\delta u)}(u)$ and $u^{\prime\prime}=\mathscr{R}_{\mathscr{F}(x+\delta u)}(u)$. It follows from (19) that, in the limit as $\delta\to 0$, a proposed bounce $(x,u)\mapsto(x^{\prime\prime},u^{\prime\prime})$ is accepted with probability $\mathcal{A}(x,u)$ described in Equation (7), with $(x^{\prime\prime},u^{\prime\prime})\to(x,\mathscr{R}_{\mathscr{F}(x)}(u))$ as $\delta\to 0$. This completes the proof of Proposition 1.

A.3 Proof of Proposition 2

where $\mathcal{L}^{(\kappa,B)}$ is the generator of the Brownian motion on the united sphere (5) and $\mathcal{F}$ is the flip operator defined as $\mathcal{F}\varphi(x,u)=\varphi(x,-u)-\varphi(x,u)$. In order to obtain the limit of the process defined in Equation (10), set

Note that time has been accelerated by a factor $d$. The process $\theta^{d}_{t}$ describes the dot product between the position $X^{d}$ and the direction $U^{d}$, scaled by a factor $d^{1/2}$ in order to observe a non-degenerate limiting process. Itô’s lemma, neglecting terms of order $1/d$, directly shows (after straightforward algebra) that the Markov process $(R^{d}_{t},\theta^{d}_{t})$ has a generator $\mathcal{L}^{\varepsilon}$ that reads

with the standard multiscale expansion notation $\varepsilon=1/\sqrt{d}$, generators $\mathcal{L}^{(J)}$ and $\mathcal{L}^{(K)}$ defined in Equation (11) and

where $\theta^{\kappa}_{t}$ is the Markov process with generator $\mathcal{L}^{(\text{Fast})}\equiv\frac{1}{\sigma}\,\mathcal{L}^{(J)}+\frac{\kappa}{2}\,\mathcal{L}^{(K)}$. Indeed, Equation (22) is the Kolmogorov backward equation associated to the Ornstein-Uhlenbeck

which concludes the proof of Proposition 2.

Appendix B High-dimensional behaviour for finite δ\delta and non-isotropic target

where $|T_{1}^{(d)}|\leq\delta^{3}/d\times L\sum_{i=1}^{d}\gamma_{i}^{3}|Z_{i}|^{3}\rightarrow\delta^{3}L\times\operatorname{E}[\gamma^{3}]\operatorname{E}[|Z_{i}|^{3}]<\infty$. Hence as $d\rightarrow\infty$, the Central Limit Theorem gives

Also, by (24) and the central limit theorem,

where $|T_{2}^{(d)}|\leq\delta^{3}L\sum_{i=1}^{d}\gamma_{i}^{3}V_{i}^{2}|Z_{i}|=\mathcal{O}(1)$ by the Lipschitz condition on $f^{\prime\prime}$, the boundedness of $\operatorname{E}[\gamma^{3}]$ and because $\|U-V\|^{2}=\mathcal{O}(1/d)$ and $U=Z/\sqrt{d}$. Further, the quantity $-B(X,U)$ also reads

Since $\|U\|=1$ and $\|U-V\|=\mathcal{O}(1/\sqrt{d})$, the term to control in (26) is

in probability. If there is a delayed-rejection event then the standard move must have been rejected and, for example, $B$ must be negative. Let $\mathsf{DR}$ be the event that the standard move has been rejected and so a delayed-rejection step is being attempted. Let $f_{B}(b)$ be the a priori density for $B$ at stationarity, and let $f_{B|\mathsf{DR}}(b)$ be the density conditional on there being a delayed-rejection event. Then

which is well-behaved and has no mass where $\alpha_{\text{dr}}(X^{(d)},U^{(d)})$ is undefined. In the limit as $d\rightarrow\infty$, $f_{B}(b)$ is the density of the Gaussian distribution in (23), and (27) gives the limiting conditional density. It follows from the Bounded Convergence Theorem that

B.2 Simulation study varying dd for fixed δ\delta

Appendix C Further simulation studies

We reran Hamiltonian Monte Carlo for $10^{6}$ iterations $40$ additional times with $X_{0}=\gamma(\sigma_{1}z_{1},\dots,\sigma_{d}z_{d})\times r_{\star}$ for each $\gamma\in\{1.5,2.0.2.5,3.0\}$, with a new, independent $z$ vector on each of the $160$ occasions. On each occasion we counted the fraction of times where, by iteration $10^{6}$ the algorithm had ever had a value with $\|x\|_{\text{M}}\leq r_{\star}$; i.e., the algorithm had reached the main posterior mass. The number of runs which converged by this measure were: $40/40~(\gamma=1.5)$, $36/40~(\gamma=2.0)$, $4/40~(\gamma=2.5)$ and $0/40~(\gamma=3.0)$; indeed, for every run with $\gamma=3.0$ the empirical acceptance rate was exactly zero. This fits with the known lack of geometric ergodicity of Hamiltonian Monte Carlo on light-tailed targets. By contrast, for the discrete bouncy particle sampler with $\gamma=3.0$, all $40$ runs converged within $1000$ iterations, and, indeed, $26$ of the runs converged within $300$ iterations.

In summary, on this occasion, when both algorithms were started from the main posterior mass, the discrete bouncy particle sampler was competetive with Hamiltonian Monte Carlo, though less efficient. However, because our algorithm depends on $\nabla\log\pi$ only through the unit vector, it is robust to large $\|\nabla\log\pi\|$, unlike Hamiltonian Monte Carlo.

C.2 The Markov modulated Poisson process

Finally, we consider a $k$-state, continuous-time Markov chain $Z_{t}$ started from state $1$, and a Poisson process $N_{t}$ whose rate $\lambda_{t}$ is a fixed function of $Z_{t}$. The doubly-stochastic process is parameterized by the rate matrix for the Markov chain, $Q$, and a vector of rates for the Poisson process, $\lambda$, where $\lambda_{i},~(i=1,\dots,k)$ is the rate of $N_{t}$ when $Z_{t}=i$.

Code was written in C++ where auto-differentiation was not available for general matrix exponentials, and so numerical differentiation via centred differences was used (the cheaper, first-order Euler approximation led to precision problems). We applied the Discrete Bouncy Particle Sampler for $2\times 10^{5}$ iterations for a number of $\kappa$ values and repeated this but evaluating only $n_{\text{cpt}}=3$ randomly-orientated components of the eight-dimensional gradient vector on each delayed-rejection step. Figure 5 plots scaled effective sample size against $\kappa$ and suggests that the optimal mean dot product is around $0.5$ when all gradient components are used and around $0.2$ when three random components are used. The optimal effective sample size in the latter case is around $3/8$ of the former; since the number of gradient calculations during a bounce is also reduced by $3/8$ this suggests no loss in overall efficiency. When all components are used, tuning to a dot product of $0.2$ brings only a $10\%$ reduction in effective sample size. The posterior variance matrix has a condition number of $49.2$, so, following typical practice, for each $\kappa$ the Discrete Bouncy Particle Sampler was rerun using a crude preconditioning matrix (Section 2.5) of $M=\mbox{diag}(1/2,2,1,1,2,2,2,2)$. Although the effective variance matrix still has a condition number of $\approx 4.3$ the right panel of Figure 5 shows that the optimal choice of $\kappa$ is now $\approx 0.2$.