The Bouncy Particle Sampler: A Non-Reversible Rejection-Free Markov Chain Monte Carlo Method

Introduction

However, the implementation of the BPS proposed in is not applicable to most target distributions arising in statistics. In this article we make the following contributions:

by reformulating explicitly the bounces times of the BPS as the first arrival times of inhomogeneous Poisson Processes (PP), we leverage standard sampling techniques [8, Chapter 6] and methods from chemical kinetics to obtain new computationally efficient ways to simulate the BPS process for a large class of target distributions.

when the target distribution can be expressed as a factor graph , a representation generalizing graphical models where the target is given by a product of factors and each factor can be a function of only a subset of variables, we adapt a physical multi-particle system method discussed in [27, Section III] to achieve additional computational efficiency. This local version of the BPS only manipulates a restricted subset of the state components at each bounce but results in a change of all state components, not just the one being updated contrary the Gibbs sampler.

we present a proof of the ergodicity of BPS when the velocity of the particle is additionally refreshed at the arrival times of an homogeneous PP. When this refreshment step is not carried out, we exhibit a counter-example where ergodicity does not hold.

we propose alternative refreshment schemes and compare their computational efficiency experimentally.

Empirically, these new MCMC schemes compare favorably to state-of-the-art MCMC methods on various Bayesian inference problems, including for high-dimensional scenarios and large data sets. Several additional original extensions of the BPS including versions of the algorithm which are applicable to mixed continuous-discrete distributions, distributions restricted to a compact support and a method relying on the use of curved dynamics instead of straight lines can be found in . For brevity, these are not discussed here.

The rest of this article is organized as follows. In Section 2, we introduce the basic version of the BPS, propose original ways to implement it and prove its ergodicity under weak assumptions. Section 3 presents a modification of the basic BPS which exploits a factor graph representation of the target distribution and develops computationally efficient implementations of this scheme. In Section 4, we demonstrate this methodology on various Bayesian models. The proofs are given in the Appendix and the Supplementary Material.

The bouncy particle sampler

2 Algorithm description

When the particle bounces, its velocity is updated in the same way as a Newtonian elastic collision on the hyperplane tangential to the gradient of the energy. Formally, the velocity after bouncing is given by

3 Algorithms for bounce time simulation

Implementing BPS requires sampling the first arrival time $\tau$ of a one-dimensional inhomogeneous PP $\Pi$ of intensity $\chi(t)=\lambda(x+vt,v)$ given by (1). Simulating such a process is a well-studied problem; see [8, Chapter 6, Section 1.3]. We review here three methods and illustrate how they can be used to implement BPS for examples from Bayesian statistics. The first method described in Section 2.3.1 will be particularly useful when the target is log-concave, while the two others described in Section 2.3.2 and Section 2.3.3 can be applied to more general scenarios.

If we let $\varXi\left(t\right)=\int_{0}^{t}\chi\left(s\right){\rm d}s$, then the PP $\Pi$ satisfies

and therefore $\tau$ can be simulated from a uniform variate $V\sim\mathcal{U}\left(0,1\right)$ via the identity

where $\varXi^{-1}$ denotes the quantile function of $\varXi,$ \varXi^{-1}(p)=\inf\mbox{\left\{t:p\leq\varXi(t)\right\}}. Refer to Figure 1, top right for a graphical illustration. This identity corresponds to the method proposed in to determine the bounce times and is also used in to simulate related processes.

In general, it is not possible to obtain an analytical expression for $\tau$. However, when the target distribution is strictly log-concave and differentiable, it is possible to solve Equation (5) numerically (see Example 1 below).

Log-concave densities. If the energy is strictly convex (see Figure 1, bottom right), we can minimize it along the line specified by $(x,v)$

where $\tau_{*}$ is well defined and unique by strict convexity. On the interval $\left[0,\tau_{*}\right)$, which might be empty, we have ${\rm d}U\left(x+vt\right)/{\rm d}t<0$ and ${\rm d}U\left(x+vt\right)/{\rm d}t\geq 0$ on \left[\tau_{*},\text{\infty}\right). The solution $\tau$ of (5) is thus necessarily such that $\tau\geq\tau_{*}$ and (5) can be rewritten using the gradient theorem as

Even if we only compute $U$ pointwise through a black box, we can solve (6) through line search within machine precision.

We note that (6) also provides an informal connection between the BPS and MH algorithms. Exponentiating this equation, we get indeed

Hence, in the log-concave case, and when the particle is climbing the energy ladder (i.e., $\tau_{*}=0$), BPS can be viewed as “swapping” the order of the steps taken by the MH algorithm. In the latter, we first sample a proposal and second sample a uniform $V$ to perform an accept-reject decision. With BPS, $V$ is first drawn then the maximum distance allowed by the same MH ratio is travelled. As for the case of a particle going down the energy ladder, the behavior of BPS is simpler to understand: bouncing never occurs. We illustrate this method for Gaussian distributions.

Multivariate Gaussian distributions. Let $U\left(x\right)=\left\|x\right\|^{2}$, then simple calculations yield

3.2 Simulation using adaptive thinning

When it is difficult to solve (5), the use of an adaptive thinning procedure provides an alternative. Assume we have access to local-in-time upper bounds $\bar{\chi}{}_{s}\left(t\right)$ on $\chi(t)$, that is

where $\triangle$ is a positive function (standard thinning corresponds to $\Delta=+\infty$). Assume additionally that we can simulate the first arrival time of the PP $\bar{\Pi}_{s}$ with intensity $\bar{\chi}_{s}(t)$. Such bounds can be constructed based on upper bounds on directional derivatives of $U$ provided the remainder of the Taylor expansion can be controlled. Algorithm 2 shows the pseudocode for the adaptive thinning procedure.

The case $V>\{\chi\left(\tau\right)/\bar{\chi}_{s}\left(\tau\right)\}$ corresponds to a rejection step in the thinning algorithm but, in contrast to rejection steps that occur in standard MCMC samplers, in the BPS algorithm this means that the particle does not bounce and just coasts. Practically, we would like ideally $\triangle$ and the ratio $\chi\left(\tau\right)/\bar{\chi}_{s}\left(\tau\right)$ to be large. Indeed this would avoid having to simulate too many candidate events from $\bar{\Pi}_{s}$ which would be rejected as these rejection steps incur a computational cost.

3.3 Simulation using superposition and thinning

Assume that the energy can be decomposed as

where $\chi^{[j]}(t)=\max\left(0,\left\langle\nabla U^{[j]}(x+tv),v\right\rangle\right)$ for $j=1,...,m$. It is therefore possible to use the thinning algorithm of Section 2.3.2 with $\bar{\chi}_{0}(t)=\sum_{j=1}^{m}\chi^{[j]}\left(t\right)$ for $t\geq 0$ (and $\Delta=+\infty$), as we can simulate from $\bar{\Pi}_{0}$ via superposition by simulating the first arrival time $\tau^{[j]}$ of each PP with intensity $\chi^{[j]}\left(t\right)\geq 0$ then returning

Exponential families. Consider a univariate exponential family with parameter $x$, observation $y,$ sufficient statistic $\phi(y)$ and log-normalizing constant $A(x)$. If we assume a Gaussian prior on $x,$ we obtain

The time $\tau^{}$ is computed analytically in Example 1 whereas the times $\tau^{}$ and $\tau^{}$ are given by

Using the superposition and thinning method (Section 2.3.3), simulation of the bounce times can be broken into subproblems corresponding to $R+1$ factors: one factor coming from the prior, with corresponding energy

and $R$ factors coming from the likelihood of each datapoint, with corresponding energy

Simulation of $\tau^{[R+1]}$ is covered in Example 1. Simulation of $\tau^{[r]}$ for $r\in\left\{1,2,\dots,R\right\}$ can be approached using thinning. In Appendix C.1, we show that

Since the bound is constant for a given $v$, we sample $\tau^{[r]}$ by simulating an exponential random variable.

4 Estimating expectations

see, e.g., . When $\varphi\left(x\right)=x_{k}$, $k\in\left\{1,2,\dots,d\right\}$, we have

When the above integral is intractable, we may just discretize $x\left(t\right)$ at regular time intervals to obtain an estimator

where $\delta>0$ is the mesh size and $L=1+\left\lfloor T/\delta\right\rfloor$. Alternatively, we could approximate these univariate integrals through quadrature.

5 Theoretical results

If we add the condition ${{\lambda^{\text{ref}}}}>0$, we get the following stronger result.

The local bouncy particle sampler

In numerous applications, the target distribution admits some structural properties that can be exploited by sampling algorithms. For example, the Gibbs sampler takes advantage of conditional independence properties. We present here a “local” version of the BPS introduced in [27, Section III] which can similarly exploit these properties and, more generally, any representation of the target density as a product of positive factors

where $x_{f}$ is a restriction of $x$ to a subset N_{f}\subseteq\mbox{\lx@text@lbrace 1,2,\dots$,$d\}} of the components of $x$, and $F$ is an index set called the set of factors. Hence the energy associated to $\pi$ is of the form

with $\partial U_{f}\left(x\right)/\partial x_{k}=0$ for any variable absent from factor $f$, i.e. for any $k\in\left\{1,2,\dots,d\right\}\backslash N_{f}$.

Such a factorization of the target density can be formalized using factor graphs (Figure 2, top). A factor graph is a bipartite graph, with one set of vertices $N$ called the variables, each corresponding to a component of $x$ ($|N|=d$), and a set of vertices $F$ corresponding to the local factors $\left(\gamma_{f}\right)_{f\in F}$. There is an edge between $k\in N$ and $f\in F$ if and only if $k\in N_{f}.$ This representation generalizes undirected graphical models [32, Chap. 2, Section 2.1.3] as, for example, factor graphs can have distinct factors connected to the same set of components (i.e. $f\neq f^{\prime}$ with $N_{f}=N_{f^{\prime}}$) as in the example of Section 4.6.

2 Local BPS: algorithm description

Similarly to the Gibbs sampler, each step of the local BPS manipulates only a subset of the $d$ components of $x$. Contrary to the Gibbs sampler, the local BPS does not require sampling from any full conditional distribution and each local calculation results in a change of all state components, not just the one being updated—how this can be done implicitly without manipulating the full state at each iteration is described below. Related processes exhibiting similar characteristics have been proposed in .

We can check that $R_{f}\left(x\right)$ satisfies

We define a collection of PP intensities based on the previous event position $x^{(i-1)}$ and velocity $v^{(i-1)}$: $\chi_{f}(t)=\lambda_{f}(x^{\left(i-1\right)}+v^{\left(i-1\right)}t,v^{\left(i-1\right)})$. In the local BPS, the next bounce time $\tau$ is the first arrival of a PP with intensity $\chi(t)=\sum_{f\in F}\chi_{f}(t)$. However, instead of modifying all velocity variables at a bounce as in the basic BPS, we sample a factor $f$ with probability $\chi_{f}(\tau)/\chi(\tau)$ and modify only the variables connected to the sampled factor. More precisely, the velocity $v_{f}$ is updated using $R_{f}\left(x_{f}\right)$ defined in (18). A generalization of the proof of Proposition 1 given in the Supplementary Material shows that the local BPS algorithm results in a $\pi-$invariant kernel. In the next subsection, we describe various computationally efficient procedures to simulate this process.

For all these implementations, it is useful to encode trajectories in a sparse fashion: each variable $k\in N$ only records information at the times $t_{k}^{(1)},t_{k}^{(2)},\dots$ where an event (a bounce or refreshment) affected it. By (19), this represents a sublist of the list of all event times. At each of those times $t_{k}^{(i)},$ the component’s position $x_{k}^{(i)}$ and velocity $v_{k}^{(i)}$ right after the event is stored. Let $L_{k}$ denote a list of triplets $(x_{k}^{(i)},v_{k}^{(i)},t_{k}^{(i)})_{i\geq 0}$, where $x_{k}^{(0)}$ and $v_{k}^{(0)}$ denote the initial position and velocity and $t_{k}^{(0)}=0$ (see Figure 2, where the black dots denote the set of recorded triplets). This list is sufficient to compute $x_{k}(t)$ for $t\leq t_{k}^{(|L_{k}|+1)}$. This procedure is detailed in Algorithm 3 and an example is shown in Figure 2, where the black square on the first variable’s trajectory shows how Algorithm 3 reconstructs $x_{1}(t)$ at a fixed time $t$: it identifies $i(t,1)=3$ as the index associated to the largest event time $t_{1}^{(3)}$ before time $t$ affecting $x_{1}$ and return $x_{1}\left(t\right)=x_{1}^{(3)}+v_{1}^{(3)}(t-t_{1}^{(3)})$.

3 Local BPS: efficient implementations

We can sample arrivals from a PP with intensity $\chi(t)=\sum_{f\in F}\chi_{f}(t)$ using the superposition method of Section 2.3.3, the thinning step therein being omitted. To implement this technique efficiently, we store potential future bounce times (called “candidates”) $t_{f}$, one for each factor, in a priority queue $Q$: only a subset of these candidates will join the lists $L_{k}$ which store past, “confirmed” events. We pick the the smallest time in $Q$ to determine the next bounce time and the next factor $f$ to modify. The priority queue structure ensures that finding the minimum element of $Q$ or inserting/updating an element of $Q$ can be performed with computational complexity $O(\log|F|)$. When a bounce occurs, a key observation behind efficient implementation of the local BPS is that not all the other candidate bounce times need to be resimulated. Suppose that the bounce was associated with factor $f$. In this case, only the candidate bounce times $t_{f^{\prime}}$ corresponding to factors $f^{\prime}$ with $N_{f^{\prime}}\cap N_{f}\neq\emptyset$ need to be resimulated. For example, consider the first bounce in Figure 2 (shown in purple), which is triggered by factor $f_{\text{a}}$ (rectangles represent candidate bounce times $t_{f}$; dashed lines connect bouncing factors to the variables that undergo an associated velocity change). Then only the velocities for the variables $x_{1}$ and $x_{2}$ need to be updated. Therefore, only the candidate bounce times for factors $f_{\text{a}}$ and $f_{\text{b}}$ need to be re-simulated while the candidate bounce time for $f_{\text{c}}$ stays constant (this is shown by an exclamation mark in Figure 2).

The method is detailed in Algorithm 4. Several operations of the BPS such as step 4, 6.iii, 6.iv and 7.ii can be easily parallelized.

3.2 Implementation via thinning

When the number of factors involved in Step 6(d)iv is large, the previous queue-based implementation can be computationally expensive. Implementing the local BPS in this setup is closely related to the problem of simulating stochastic chemical kinetics and innovative solutions have been proposed in this area. We adapt here the algorithm proposed in to the local BPS context. For ease of presentation, we present the algorithm without refreshment and only detail the simulation of the bounce times. This algorithm relies on the ability to compute local-in-time upper bounds on $\lambda_{f}$ for all $f\in F$. More precisely, we assume that given a current position $x$ and velocity $v$, and $\Delta\in(0,\infty]$, we can find a positive number $\bar{\chi}_{f}$, such that for any $t\in[0,\Delta)$, we have $\bar{\chi}_{f}\geq\lambda_{f}(x+vt,v)$. We can also use this method on a subset $G$ of $F$ and combine it with the previously discussed techniques to sample candidate bounce times for factors in $F$\$G$ but we restrict ourselves to $G=F$ to simplify the presentation.

Algorithm 6 will be particularly useful in scenarios where summing over the bounds (Step 5a) and sampling a factor (Step 5(c)ii) can be performed efficiently. A scenario where it is possible to implement these two operations in constant time is detailed in Section 4.6. Another scenario where sampling quickly from $\mathcal{F}$ is feasible is if the number of distinct upper bounds is much smaller than the number of factors. For example, we only need to sample a factor $\mathcal{F}$ uniformly at random if $\Lambda=\bar{\chi}_{f}=\bar{\chi}_{f^{\prime}}$ for all $f,f^{\prime}$ in $F$ and $\bar{\chi}=\left|F\right|\cdot\Lambda$ (that is no factor needs to be inspected in order to execute Algorithm 5(c)ii) and the thinning procedure in Step 5(c)iii boils down to

A related approach has been adopted in for the analysis of big data. In this particular scenario, an alternative local BPS can also be implemented where $s>1$ factors $\mathcal{F}=\left(\mathcal{F}_{1},\dots,\mathcal{F}_{s}\right)$ are sampled uniformly at random without replacement from $F$, the thinning occurs with probability

and the components of $x$ belonging to $N_{\mathcal{F}}$ bounce based on $\sum_{j=1}^{s}\nabla U_{\mathcal{F}_{j}}(x)$. One can check that the resulting dynamics preserves $\pi$ as an invariant distribution. In contrast to $s=1$, this is not an implementation of local BPS described in Algorithm 6, but instead this corresponds to a local BPS update for a random partition of the factors.

Numerical results

We consider an isotropic multivariate Gaussian target distribution, $U\left(x\right)=\left\|x\right\|^{2}$, to illustrate the need for refreshment. Without refreshment, we obtain from Equation (7)

In this scenario, we show that BPS without refreshment admits a countably infinite collection of invariant distributions. Let us define $r\left(t\right)=\left\|x\left(t\right)\right\|$ and $m\left(t\right)=\left\langle x\left(t\right),v\left(t\right)\right\rangle/\left\|x\left(t\right)\right\|$ and denote by $\chi_{k}$ the probability density of the chi distribution with $k$ degrees of freedom.

For any dimension $d\geq 2$, the process $\left(r\left(t\right),m\left(t\right)\right)_{t\geq 0}$ is Markov and its transition kernel is invariant with respect to the probability densities $\left\{f_{k}(r,m)\propto\chi_{k}(\sqrt{2}r)\cdot(1-{{m}}^{2})^{(k-3)/2};k\in\left\{2,3,\ldots\right\}\right\}$.

Next, we look at the scaling of the Effective Sample Size (ESS) per CPU second of the basic BPS algorithm for $\varphi\left(x\right)=x_{1}$ when $\lambda_{\text{ref}}=1$ as the dimension $d$ of the isotropic normal target increases. The ESS is estimated using the R package mcmcse by evaluating the trajectory on a fine discretization of the sampled trajectory. The results in log-log scale are displayed in Figure 3. The curve suggests a decay of roughly $d^{-1.47}$, slightly inferior to the $d^{-1.25}$ scaling for an optimally tuned Hamiltonian Monte Carlo (HMC) algorithm [6, Section III], [24, Section 5.4.4]. It should be noted that BPS achieves this scaling without varying any tuning parameter, whereas HMC’s performance critically depends on tuning two parameters (leap-frog stepsize and number of leap-frog steps). Both BPS and HMC compare favorably to the $d^{-2}$ scaling of the optimally tuned random walk MH .

2 Comparison of the global and local schemes

We compare the basic “global” BPS of Section 2 to the local BPS of Section 3 on a sparse Gaussian field. We use a chain-shaped undirected graphical model of length $d=1000$ and perform separate experiments for various pairwise precision parameters for the interaction between neighbors in the chain. Both methods are run for 60 seconds. We compare the Monte Carlo estimate of the variance of $x_{500}$ to its true value. The results are shown in Figure 4. The smaller computational complexity per local bounce of the local BPS offsets significantly the associated decrease in expected trajectory segment length. Moreover, both versions appear insensitive to the pairwise precision used in this sparse Gaussian field.

3 Comparisons of alternative refreshment schemes

In Section 2, the velocity was refreshed using a Gaussian distribution. We compare here this global refreshment scheme to three alternatives:

if the local BPS is used, the factor graph structure can be exploited to design computationally cheaper refreshment operators. We pick one factor $f\in F$ uniformly at random and resample only the components of $v$ with indices in $N_{f}$. By the same argument used in Section 3, each refreshment requires bounce time recomputation only for the factors $f^{\prime}$ with $N_{f}\cap N_{f^{\prime}}\neq\emptyset$.

the velocities are refreshed according to $\phi\left(v\right)$, the uniform distribution on $\mathcal{S}^{d-1}$, and the BPS admits now $\rho\left(z\right)=\pi\left(x\right)\phi\left(v\right)$ as invariant distribution.

a variant of restricted refreshment where we sample an angle $\theta$ by multiplying a Beta($\alpha$, $\beta$)-distributed random variable by $2\pi.$ We then select a vector uniformly at random from the unit length vectors that have an angle $\theta$ from $v$. We used $\alpha=1,\beta=4$ to favor small angles.

One limitation of the results in this section is that the optimal refreshment scheme and refreshment rate will in general be problem dependent. Adaptation methods used in the HMC literature could potentially be adapted to this scenario , but we leave these extensions to future work.

4 Comparisons with HMC methods on high-dimensional Gaussian distributions

Next, we compare the local BPS to NUTS (“adapt=true,fit_metric=true,nuts=true”) as the dimension $d$ increases. Experiments are performed on the chain-shaped Gaussian Random Field of Section 4.2 with the pairwise precision parameter set to 0.5. We vary the length of the chain (10, 100, 1000), and run Stan’s implementation of NUTS for 1000 iterations + 1000 iterations of adaptation. We measure the wall-clock time (excluding the time taken to compile the Stan program) and then run our method for the same wall-clock time 40 times for each chain size. The absolute value of the relative error averaged on 10 equally spaced marginal variances is measured as a function of the percentage of the total computational budget used; see Figure 7. The gap between the two methods widens as $d$ increases. To visualize the different behavior of the two algorithms, three marginals of the Stan and BPS paths for $d=100$ are shown in Figure 8. Contrary to Section 4.1, BPS outperforms here HMC as its local version is able to exploit the sparsity of the random field.

5 Poisson-Gaussian Markov random field

6 Bayesian logistic regression for large data sets

Consider the logistic regression model introduced in Example 3 when the number of data $R$ is large. In this context, standard MCMC schemes such as the MH algorithm are computationally expensive as they require evaluating the likelihood associated to the $R$ observations at each iteration. This has motivated the development of techniques which only evaluate the likelihood of a subset of the data at each iteration. However, most of the methods currently available introduce either some non-vanishing asymptotic bias, e.g. the subsampling MH scheme proposed in , or provide consistent estimates converging at a slower rate than regular MCMC algorithms, e.g. the Stochastic Gradient Langevin Dynamics introduced in . The only available algorithm which only requires evaluating the likelihood of a subset of data at each iteration yet provides consistent estimates converging at the standard Monte Carlo rate is the Firefly algorithm .

In this context, we associate $R+1$ factors to the target posterior distribution: one for the prior and one for each data point with $x_{f}=x$ for all $f\in F$. As a uniform upper bound on the intensities of these local factors is available for restricted refreshment, see Appendix C.1, we could use (21) in conjunction with Algorithm 6 to provide an alternative to the Firefly algorithm which selects at each bounce a subset of $s$ data points uniformly at random without replacement. For $s=1$, a related algorithm has been recently explored in . In presence of outliers, this strategy can be inefficient as the uniform upper bound becomes very large, resulting in a computationally expensive implementation. After a pre-computation step of complexity $O(R\log R)$ only executed once, we show here that Algorithm 6 can be implemented using data-dependent bounds mitigating the sensitivity to outliers while maintaining the computational cost of each bounce independent of $R$. We first pre-compute the sum of covariates over the data points, $\iota_{k}^{c}=\sum_{r=1}^{R}\iota_{r,k}{\mathbf{1}}[y_{r}=c]$, for $k\in\left\{1,\ldots,d\right\}$ and class label $c\in\left\{0,1\right\}$. Using these quantities, it is possible to compute

with $\bar{\chi}^{[r]}$ given in (12). If $d$ is large, we can keep the sum $\bar{\chi}$ in memory and add-and-subtract any local updates to it. The implementation of Step 5(c)ii relies on the alias method [8, Section 3.4]. A detailed description of these derivations and of the algorithm is presented in Appendix C.

We compare the local BPS with thinning to the MAP-tuned Firefly algorithm implementation provided by the authors. This version of Firefly outperforms experimentally significantly the standard MH in terms of ESS per-datum likelihood . The two algorithms are here compared in terms of this criterion, where the ESS is averaged over the $d=5$ components of $x$. We generate covariates $\iota_{rk}\overset{\text{i.i.d.}}{\sim}\mathcal{U}(0.1,1.1)$ and data $y_{r}\in\{0,1\}$ for $r=1,\dots,R$ according to (9) and set a zero-mean normal prior of covariance $\sigma^{2}I_{d}$ for $x$. For the algorithm, we set $\lambda_{\text{ref}}=0.5,\,\sigma^{2}=1$ and $\Delta=0.5$, which is the length of the time interval for which a constant upper bound for the rate associated with the prior is used, see Algorithm 7. Experimentally, local BPS always outperforms Firefly, by about an order of magnitude for large data sets. However, we also observe that both Firefly and local BPS have an ESS per datum likelihood evaluation decreasing in approximately $1/R$ so that the gains brought by these algorithms over a correctly scaled random walk MH algorithm do not appear to increase with $R$. The rate for local BPS is slightly superior in the regime of up to $10^{4}$ data points, but then returns to the approximate $1/R$ rate. To improve this rate, one can adapt the control variate ideas introduced in for the MH algorithm to these schemes. This has been proposed in for a related algorithm and in for the local BPS.

7 Bayesian inference of evolutionary parameters

We compare against a state-of-the-art HMC sampler that uses Bayesian optimization to adapt the key parameters of HMC, the leap-frog stepsize and the number of leap-frog steps, while preserving convergence to the target distribution. Both our method and this HMC method are implemented in Java and share the same gradient computation code. Refer to the Supplementary Material for additional background and motivation behind this adaptation method.

We first perform various checks to ensure that both BPS and HMC chains are in close agreement given a sufficiently large number of iterations. After 20 millions HMC iterations, the credible intervals estimates from the HMC method are in close agreement with those obtained from BPS (result not shown) and both methods pass the Geweke diagnostic .

To compare the effectiveness of the two samplers, we look at the ESS per second of the model parameters. We show the maximum, median, and maximum over the 10 parameter components for 10 runs, for both BPS and HMC in Figure 11. We observe a speed-up by a factor two for all statistics considered (maximum, median, minimum). In the supplement, we show that the HMC chain displays much larger autocorrelations than the BPS chain.

Discussion

Most MCMC methods currently available, such as the MH and HMC algorithms, are discrete-time reversible processes. There is a wealth of theoretical results showing that non-reversible Markov processes mix faster and provide lower variance estimates of ergodic averages . However, most of the non-reversible processes studied in the literature are diffusions and cannot be simulated exactly. The BPS is an alternative continuous-time Markov process which, thanks to its piecewise deterministic paths, can be simulated exactly for many problems of interest in statistics.

As any MCMC method, the BPS can struggle in multimodal scenarios and when the target exhibits very strong correlations. However, for a range of applications including sparse factor graphs, large datasets and high-dimensional settings, we have observed empirically that BPS is on par or outperforms state-of-the art methods such as HMC and Firefly. The main practical limitation of the BPS compared to HMC is that its implementation is model-specific and requires more than knowing $\nabla U$ pointwise. An important open problem is therefore whether its implementation, and in particular the simulation of bouncing times, can be fully automated. However, the techniques described in Section 2.3 are already sufficient to handle many interesting models. There are also numerous potential methodological extensions of the method to study. In particular, it has been shown in how one can exploit the local geometric structure of the target to improve HMC and it would be interesting to investigate how this could be achieved for BPS. More generally, the BPS is a specific continuous-time piecewise deterministic Markov process . This class of processes deserves further exploration as it might provide a whole new class of efficient MCMC methods.

Acknowledgements

Alexandre Bouchard-Côté’s research is partially supported by a Discovery Grant from the National Science and Engineering Research Council. Arnaud Doucet’s research is partially supported by the Engineering and Physical Sciences Research Council (EPSRC), grants EP/K000276/1, EP/K009850/1 and by the Air Force Office of Scientific Research/Asian Office of Aerospace Research and Development, grant AFOSRA/AOARD-144042. Sebastian Vollmer’s research is partially supported by the EPSRC grants EP/K009850/1 and EP/N000188/1. We thank Markus Upmeier for helpful discussions of differential geometry as well as Nicholas Galbraith, Fan Wu, and Tingting Zhao for their comments.

References

Appendix A Proofs of Section 2

The BPS process is a specific piecewise-deterministic Markov process so the expression of its generator follows from [7, Theorem 26.14]. Its adjoint is given in and a derivation of this expression from first principles can be found in the Supplementary Material. To establish the invariance with respect to $\rho$, we first show that $\mathcal{\int L}h(z)\rho\left(z\right){\rm d}z=0$. We have

As $\rho\left(z\right)=\pi\left(x\right)\psi\left(v\right)$, the term (24) is trivially equal to zero, while a change-of-variables shows that

as $R^{-1}\left(x\right)v=R\left(x\right)v$ and$\left\|R\left(x\right)v\right\|=\left\|v\right\|$ implies $\psi\left(R\left(x\right)v\right)=\psi\left(v\right)$. Additionally, by integration by parts, we obtain as $h$ is bounded

Substituting (25) and (26) into (22)-(23)-(24), we obtain

where we have used $\left\langle\nabla U(x),R\left(x\right)v\right\rangle=-\left\langle\nabla U(x),v\right\rangle$ and $\max\{0,-f\}-\max\{0,f\}=-f$ for any $f$. The result now follows by [7, Proposition 34.7].

A.2 Proof of Theorem 1

We can now define our tractable set and establish its key properties.

Let $t>0$, and assume the initial point of the BPS, satisfies $\|x_{0}\|\leq t$, $\|v_{0}\|\leq 1$. If $\|\nabla U\|^{*}=\sup\left\{\|\nabla U(x)\|:\|x\|\leq 3t\right\}$ then the event

On the event $\mathscr{E}$, we have $\|v(t^{\prime})\|\leq 1$ and $\|x(t^{\prime})\|\leq 2t$ for all $t^{\prime}\in\left[0,t\right]$,

On the event $\mathscr{E}$, there are exactly two refreshments and no bouncing in the interval $(0,t)$, i.e. $\tau_{1}=\tau_{1}^{\text{(ref)}}$, $\tau_{2}=\tau_{2}^{\text{(ref)}},$ and $\tau_{1}+\tau_{2}\leq t\leq\tau_{1}+\tau_{2}+\tau_{3}$,

$v_{i}|\mathscr{E}\overset{\text{i.i.d.}}{\sim}\psi_{\leq 1}(0,I)$ for $i\in\left\{1,2\right\}$, where $\psi_{\leq 1}$ denotes the truncated Gaussian distribution, with $\|v_{i}\|\leq 1$,

$\left(\tau_{1}^{\text{(ref)}},\tau_{2}^{\text{(ref)}}\right)|\mathscr{E}\sim\mathcal{U}\left(\left\{(\tau_{1},\tau_{2})\in(0,t)^{2}:\tau_{1}+\tau_{2}\leq t\right\}\right)$.

To prove Part 1 and 2, we will make use of this preliminary result: on $\mathscr{E}_{3},$ $\|v_{i-1}\|\leq 1,\|x_{i-1}\|\leq 2t$ implies $\tau_{i}^{\text{(bounce)}}\geq t$ for $i\in\left\{1,2,3\right\}$. Indeed, $\|v_{i-1}\|\leq 1$ and $\|x_{i-1}\|\leq 2t$ imply that $\chi_{z_{i-1}}(t^{\prime})\leq\|\nabla U\|^{*}$ for all $t^{\prime}\in\left[0,t\right]$. It follows that $\Xi_{z_{i-1}}(t)\leq\|\nabla U\|^{*}t$. Hence, by the continuity of $\Xi_{z_{i-1}}$ and standard properties of the quantile function, $\tau_{i}^{\text{(bounce)}}=\Xi_{z_{i-1}}^{-1}(e_{i}^{\text{(bounce)}})\geq t$.

Part 1 and 2: by the assumption on $x_{0}$ and $v_{0}$ and our preliminary result, $\tau_{1}^{\text{(bounce)}}\geq t$, and hence, combining with $\mathscr{E}_{1}$ and $\mathscr{E}_{2}$ we have $\tau_{1}=\tau_{1}^{\text{(ref)}}\leq t$ and $\|v_{1}\|=\|n_{1}\|\leq 1$. Also, by the triangle inequality, $\|x_{1}\|\leq\|x_{0}\|+\|x_{1}-x_{0}\|\leq t+\tau_{1}^{\text{(ref)}}\leq 2t$. We can therefore apply our preliminary result again and obtain $\tau_{2}^{\text{(bounce)}}\geq t$, and hence, combining again with $\mathscr{E}_{1}$ and $\mathscr{E}_{2}$, we have $\tau_{2}=\tau_{2}^{\text{(ref)}}$, $\tau_{1}+\tau_{2}\leq t$, $\|v_{2}\|=\|n_{2}\|\leq 1$. Applying the triangle inequality a second time yields $\|x_{2}\|\leq t+\tau_{1}^{\text{(ref)}}+\tau_{2}^{\text{(ref)}}\leq 2t$. We apply our preliminary result one last time to obtain $\tau_{3}^{\text{\text{(bounce)}}}\geq t$. Hence, if $\tau_{3}^{\text{(ref)}}>\tau_{3}^{\text{(bounce)}}$, $\tau_{1}+\tau_{2}+\tau_{3}=\tau_{1}^{\text{(ref)}}+\tau_{2}^{\text{(ref)}}+\tau_{3}^{\text{(bounce)}}\geq t$, while if $\tau_{3}^{\text{(ref)}}\leq\tau_{3}^{\text{(bounce)}}$, we can use $\mathscr{E}_{1}$ to conclude that $\tau_{1}+\tau_{2}+\tau_{3}=\tau_{1}^{\text{(ref)}}+\tau_{2}^{\text{(ref)}}+\tau_{3}^{\text{(ref)}}\geq t$. It follows from the triangle inequality that $\|x(t^{\prime})\|\leq 2t$ for all $t^{\prime}\in\left[0,t\right]$.

Part 3, 4 and 5: these follow straightforwardly from the construction of $\mathscr{E}$. ∎

Note that the statement and proof of Part 4 is simple because $\mathscr{E}_{3}\in\sigma(e_{i}^{\text{(bounce)}}:i\in\left\{1,2,3\right\})$. In contrast, conditioning on conceptually simpler events of the form $(\tau_{i}^{\text{(bounce)}}>t)\in\sigma\left(e_{i}^{\text{(bounce)}},z_{i-1}\right)$ leads to conditional distributions on $v_{k}$ which are harder to characterize.

In the following, $B_{R}\left(x\right)$ denotes the $d$-dimensional Euclidean ball of radius $R$ centered at $x$.

For all $\epsilon,t>0$ such that $\epsilon\leq t/6$, and $v,v_{2}\in B_{1}(0)$, $x,x^{\prime}\in B_{\epsilon}(0)$, $0\leq\tau_{1}\leq\frac{t}{6}$, and $\frac{2t}{3}\leq\tau_{2}\leq\frac{5t}{6}$, we have $\left\|v_{1}\right\|\leq 1$, where $v_{1}$ is defined by:

Using Lemma 1, Parts 4 and 5, we can rewrite the above conditional expectation as:

We will use the coarea formula to reorganize the order of integration, see e.g. Section 3.2 of . We start by introducing some notation.

For $C^{1}$ Riemannian manifolds $M$ and $N$ of dimension $m$ and $n$, a differentiable map $F:M\rightarrow N$ and $h$ a measurable test function, the coarea formula can be written as:

Here $\mathcal{H}_{m}$, $\mathcal{H}_{n}$ and $\mathcal{H}_{m-n}$ denote the volume measures associated with the Riemannian metric on $M$, $N$ and $F^{-1}(u)$ (with the induced metric of $M$). In the above equations, $JF$ is a generalization of the determinant of the Jacobian $JF:=\sqrt{\det g(\nabla f_{i},\nabla f_{j})}$ where $g$ is the corresponding Riemannian metric and $f$ is the representation of $F$ in local coordinates, see . Here JF\text{=\sqrt{\det DF\,DF^{\top}}} where $DF$ is defined in equation (A.2) below.

We apply the coarea formula to $M=\left\{(\tau_{1},\tau_{2})\in(0,t)^{2}:\tau_{1}+\tau_{2}\leq t\right\}\times B_{1}(0)\times B_{1}(0)$, $N=B_{t}(x)\times B_{1}(0)$, F_{z}(\tau_{1},{{v}}_{1},\tau_{2},{\hyperlink{v}{\color[rgb]{0,0,0}{v}}}_{2})=\left(x+\tau_{1}v+\tau_{2}v_{1}+(t-\tau_{1}-\tau_{2})v_{2},v_{2}\right), $m=2d+2$, $n=2d$ and obtain:

where $p(\tau_{1},\tau_{2})$ denotes the joint conditional density of $\tau_{1},\tau_{2}|\mathscr{E}$ described in Part 5 of Lemma 1, and:

We define $\delta^{\prime}=\inf\left\{I_{z}(z^{\prime}):z,z^{\prime}\in B\right\}$, and obtain the following inequality

It is therefore enough to show that $\delta^{\prime}>0$. To do so, we will derive the following bounds related to the integral in $I_{z}(z^{\prime})$:

Its domain of integration $F_{z}^{-1}(z^{\prime})$ is guaranteed to contain a set of positive $\mathcal{H}_{2}$ measure.

Its integrand is bounded below by a strictly positive constant.

To establish 1, we let $z^{\prime}=(x^{\prime},v^{\prime})=(x^{\prime},v_{2})$, and notice that rearranging

yields an expression for $v_{1}$ given in (27). From Lemma 2, it follows that

and $\mathcal{H}_{2}(C_{z}(z^{\prime}))\geq(t/6)^{2}$ since the surface of the graph of a function is larger than the surface of the domain.

To establish 2, we start by analyzing $JF_{z}$. Exploiting its block structure, we obtain:

Moreover, it follows from basic properties of the truncated Gaussian distribution $\psi_{\leq 1}$ and of $p(\tau_{1},\tau_{2})$ that

To prove Part 2 of the lemma, we divide the trajectory of length $nt_{0}$ into three “phases” namely a deceleration, travel, and acceleration phases, or respective lengths $t_{\text{d}}+t_{\text{t}}+t_{\text{a}}=nt_{0}$ defined below (see also Supplement for a figure illustrating the notation used in this part of the lemma). This allows us to use Part 1 of the present lemma which requires velocities bounded in norm by one. We require an acceleration phase since $W$ may not necessarily include velocities of norms bounded by one.

First, we show that we decelerate with positive probability by time $t_{\text{\text{d}}}=t_{0}$. Let $\mathscr{R}$ denote the event that there is exactly one refreshment in the interval $(0,t_{0})$, and that the refreshed velocity has norm bounded by one. Define also $r_{0}=t_{0}\max\left\{1,\|v\|\right\}$, which bounds the distance travelled in $(0,t_{0})$ for outcomes in $\mathscr{R}$, since bouncing does not change the norm of the velocity. We have:

Next, to prepare applying the first part of the lemma, set

Informally, $\epsilon$ is selected so that the ball of radius $\epsilon$ around the origin contains both any position attained after deceleration, as well as ball around a point $x^{\star}$ in $W$. Indeed, since $W$ is open and that $\epsilon>r^{\prime}$, there exists some $r>0,\left(x^{\star},v^{\star}\right)\in W$ such that $B_{r}(x^{\star})\times B_{r}(v^{\star})\subseteq W$ and $B_{\epsilon}(0)\supseteq B_{r_{0}}(x)\cup B_{r}(x^{\star})$. Let also $n=2+\left\lceil\frac{6\epsilon}{t_{0}}\right\rceil$ .

Of the total time $nt_{0}$, we reserve time $t_{\text{a}}=\min\left\{t_{0},(2(\|v^{\star}\|/r+1))^{-1},r/2\right\}$ to accelerate. This time is selected so that (a) $t_{0}-t_{\text{a}}\geq 0$, and (b), if we start with a position in $B_{r/2}(x^{\star})$, move with a velocity bounded in norm by $\bar{v}=\max\left\{1,\|v^{\star}\|+r\right\}$ for a time $\Delta t\leq\min\left\{(2(\|v^{\star}\|/r+1))^{-1},r/2\right\}$, we have that the final position is in $B_{r}(x^{\star})$. This holds since the distance travelled is bounded by $\bar{v}\Delta t\leq r/2$. Hence, by a similar argument as used for deceleration, we have, for all $z^{\prime\prime}\in B_{r/2}(x^{\star})\times B_{1}(0)$,

We can now exploit this Lemma to prove Theorem 1.

This contradicts that $\mu_{i}(A_{i})=0$ for $i\in\left\{1,2\right\}$.

The law of large numbers then follows by Birkhoff’s pointwise ergodic theorem; see e.g. [9, Theorem 2.30, Section 2.6.4]. ∎

Appendix B Proof of Proposition 2

The dynamics of the BPS can be lumped into a two-dimensional Markov process involving only the radius $r\left(t\right)=\left\|x\left(t\right)\right\|$ and $m\left(t\right)=\left\langle x\left(t\right),v\left(t\right)\right\rangle/\left\|x\left(t\right)\right\|$ for any dimensionality $d\geq 2$. The variable $m\left(t\right)$ can be interpreted (via $\arccos(m\left(t\right))$) as the angle between the particle position $x\left(t\right)$ and velocity $v\left(t\right)$. Because of the strong Markov property we can take $\tau_{1}=0$ without loss of generality and let $t$ be some time between the current event and the next, yielding:

If there is a bounce at time $t$, then $r\left(t\right)$ is not modified but $m\left(t\right)=-m\left(t\right)$.The bounce happens with intensity ${{\lambda}}(x+tv,v)=\max\left(0,\left\langle x+vt,v\right\rangle\right)$. These processes can also be written as an Stochastic Differential Equation (SDE) driven by a jump process whose intensity is coupled to its position. This is achieved by writing the deterministic dynamics given in (35) between events as the following Ordinary Differential Equation (ODE):

Taking the bounces into account turns this ODE into an SDE with

where $N_{t}$ is the counting process associated with a PP with intensity $\max\left(0,r\left(t\right)m\left(t\right)\right)$.

Now consider the push forward measure of $\mathcal{N}\left(0,\frac{1}{2}I_{k}\right)\otimes\mathcal{U}(\mathcal{S}^{k-1})$ under the map $\left(x,v\right)\mapsto\left(\left\|x\right\|,\left\langle x,v\right\rangle/\left\|x\right\|\right)$ where $\mathcal{U}(\mathcal{S}^{k-1})$ is the uniform distribution on $\mathcal{S}^{k-1}$. This yields the collection of measures with densities $f_{k}(r,m)$. One can check that $f_{k}(r,m)$ is invariant for (36) for all $k\geq 2$.

Appendix C Bayesian logistic regression for large datasets

We derive here a datapoint-specific upper bound $\bar{\chi}^{[r]}$ to $\chi^{[r]}(t)$. First, we need to compute the gradient for one datapoint:

We then consider two sub-cases depending on $y_{r}=0$ or $y_{r}=1$. Suppose first $y_{r}=0$, and let $x(t)=x+tv$

When implementing Algorithm 6, we need to bound $\sum_{r=1}^{R}\chi^{[r]}(t)$. We have

The bound $\bar{\chi}$ is constant between bounce events and only depends on the magnitude of $v$. If we further assume that we use restricted refreshment then this bound is valid for any $t>0$ allowing us to implement Algorithm 6 using (20) or (21).

C.2 Sampling the thinned factor

We show here how to implement Step 5(c)ii of Algorithm 6 without enumerating over the $R$ datapoints. We begin by introducing some required pre-computed data structures. The pre-computation is executed only once at the beginning of the algorithm, so its running time, $O(R\log R)$ is considered negligible (the number of bouncing events is assumed to be greater than $R$). For each dimensionality $k$ and class label $c$, consider the categorical distribution with the following probability mass function over the datapoints:

This is just the distribution over the datapoints that have the given label, weighted by the covariate $k$. An alias sampling data-structure [8, Section 3.4] is computed for each $k$ and $c$. This pre-computation takes total time $O(R\log R)$. This allows subsequently to sample in time $O(1)$ from the distributions $\mu_{k}^{(c)}$.

We now show how this pre-computation is used to to implement Step 5(c)ii of Algorithm 6. We denote the probability mass function we want to sample from by

To sample this distribution efficiently, we construct an artificial joint distribution over both datapoints and covariate dimension indices

We denote by $q_{\textrm{k}}(k)$, respectively $q_{\textrm{r}|\textrm{k}}(r|k)$, the associated marginal, respectively conditional distribution. By construction, we have

It is therefore enough to sample $(r,k)$ and to return $r$. To do so, we first sample (a) $k\sim q_{\textrm{k}}(\cdot)$ and then (b) sample $r|k\sim q_{\textrm{r}|\textrm{k}}(\cdot|k)$.

so this sampling step again does not require looping over the datapoints thanks the pre-computations described earlier.

and therefore this sampling step can be computed in $O(1)$ thanks to the pre-computed alias sampling data structure.

C.3 Algorithm description

Algorithm 7 contains a detailed implementation of the local BPS with thinning for the logistic regression example (Example 4.6).

Supplemental Material: The Bouncy Particle Sampler A Non-Reversible Rejection-Free Markov Chain Monte Carlo Method

Appendix D Illustration for Lemma 3

The following figure illustrates the different phases considered in the proof of Lemma 3.

Appendix E Direct proof of invariance

Let $\mu_{t}$ be the law of $z\left(t\right)$. In the following, we prove invariance by explicitly verifying that the time evolution of the density $\frac{{\rm d}\mu_{t}}{{\rm d}t}=0$ is zero if the initial distribution $\mu_{0}$ is given by $\rho(z)=\pi\left(x\right)\psi\left(v\right)$ in Proposition 1. This is achieved by deriving the forward Kolmogorov equation describing the evolution of the marginal density of the stochastic process. For simplicity, we start by presenting the invariance argument when $\lambda^{\text{{ref}}}=0$.

It follows that the probability of having no bounce in the interval $[0,t]$ is given by:

and the density of the random variable $T_{1}$ is given by:

If a bounce occurs, then the algorithm follows a translation path for time $T_{1}$, at which point the velocity is updated using a bounce operation $C(z)$, defined as:

The algorithm then continues recursively for time $t-T_{1}$, in the following sense: a second bounce time $T_{2}$ is simulated by adding to $T_{1}$ a random increment with density $q(\cdot;C\circ\Phi_{t_{1}}(z))$. If $T_{2}>t$, then the output of the algorithm is $\Phi_{t-t_{1}}\circ C\circ\Phi_{t_{1}}(z)$, otherwise an additional bounce is simulated, etc. More generally, given an initial point $z$ and a sequence $\mathbf{t}=(t_{1},t_{2},\dots)$ of bounce times, the output of the algorithm at time $t$ is given by:

where $(\thinspace)$ denotes the empty list and $\mathbf{t}^{\prime}$ the suffix of $\mathbf{t}$: $\mathbf{t}^{\prime}=(t_{2},t{}_{3},\dots)$. As for the bounce times, they are distributed as follows:

where $T_{1:i-1}=\left(T_{1},T_{2},\ldots,T_{i-1}\right).$

From this, we get the following decomposition:

Indeed, on the event that $n\geq 1$ bounces occur, the random variable $h(\Phi_{t}(z))$ only depends on a finite dimensional random vector, $(T_{1},T_{2},\dots,T_{n})$, so we can write the expectation as an integral with respect to the density $\widetilde{q}(t_{1:n};t,z)$ of these variables:

Marginal density. Let us fix some arbitrary time $t>0$. We seek a convenient expression for the marginal density at time $t$, $\mu_{t}(z)$, given an initial vector $Z\sim\rho$, where $\rho$ is the hypothesized stationary density $\rho(z)=\pi\left(x\right)\psi\left(v\right)$ on $Z$. To do so, we look at the expectation of an arbitrary non-negative measurable test function $h$:

We used the following in the above derivation successively the law of total expectation, equation (S12), equation (S18), Tonelli’s theorem and the change of variables, $z^{\prime}=\Psi_{t_{1:n},t}(z)$, justified since for any fixed $0<t_{1}<t_{2}<\dots<t_{n}<t<t_{n+1}$, $\Psi_{t_{1:n},t}(\cdot)$ is a bijection (being a composition of bijections). Now the absolute value of the determinant is one since $\Psi_{\mathbf{t},t}\left(z\right)$ is a composition of unit-Jacobian mappings and, by using Tonelli’s theorem again, we obtain that the expression above the brace is necessarily equal to $\mu_{t}(z^{\prime})$ since $h$ is arbitrary.

Derivative. Our goal is to show that for all $z^{\prime}\in\mathcal{Z}$

Since the process is time homogeneous, once we have computed the derivative, it is enough to show that it is equal to zero at $t=0$. To do so, we decompose the computation according to the terms $I_{n}$ in Equation (S22):

The categories of terms in Equation (S23) to consider are:

No bounce: $n=0$, $\Psi_{t_{1:n},t}(z)=\Phi_{t}(z)$, or,

Exactly one bounce: $n=1$, $\Psi_{t_{1:n},t}(z)=F_{t,t_{1}}:=\Phi_{t-t_{1}}\circ C\circ\Phi_{t_{1}}(z)$ for some $t_{1}\in(0,t)$, or,

Two or more bounces: $n\geq 2$, $\Psi_{t_{1:n},t}(z)=\Psi_{t-t_{2}}\circ C\circ F_{t_{2},t_{1}}(z)$ for some $0<t_{1}<t_{2}<t$

In the following, we show that the derivative of the terms in the third category, $n\geq 2$, are all equal to zero, while the derivative of the first two categories cancel each other.

No bounce in the interval. From Equation (S14):

We now compute the derivative at zero of the above expression:

The first term in the above equation can be simplified as follows:

using Equation (S4). In summary, we have:

Exactly one bounce in the interval. From Equation (S16), the trajectory consists in a bounce at a time $T_{1}$, occurring with density (expressed as before as a function of the final point $z^{\prime}$) $q(t_{1};F_{t,t_{1}}^{-1}(z^{\prime}))$, followed by no bounce in the interval $(T_{1},t]$, an event of probability:

where we used that $C^{-1}=C$. This yields:

To compute the derivative of the above equation at zero, we use again Leibniz’s rule:

Two or more bounces in the interval. For a number of bounce, we get:

and hence, using Leibniz’s rule on the integral over $t_{1}$:

Putting all terms together. Putting everything together, we obtain:

where we used that $\rho(z^{\prime})=\rho(C(z^{\prime}))$, $\left\langle\nabla U(x^{\prime}),R\left(x^{\prime}\right)v^{\prime}\right\rangle=-\left\langle\nabla U(x^{\prime}),v^{\prime}\right\rangle$ and $-\max\{0,f\}+\max\{0,-f\}=-f$ for any function $f$. Hence we have $\left.\frac{{\rm d}\mu_{t}(z^{\prime})}{{\rm d}t}\right|_{t=0}=0$, establishing that that the bouncy particle sampler $\lambda^{\text{ref}}=0$ admits $\rho$ as invariant distribution. The invariance for $\lambda^{\text{ref}}>0$ then follows from Lemma 4 given below.

Suppose $P_{t}$ is a continuous time Markov kernel and $Q$ is a discrete time Markov kernel which are both invariant with respect to $\mu.$ Suppose we construct for $\lambda^{\text{{ref}}}>0$ a Markov process $\hat{P}_{t}$ as follows: at the jump times of an independent PP with intensity $\lambda^{\text{{ref}}}$ we make a transition with $Q$ and then continue according to $P_{t}$, then $\hat{P}_{t}$ is also $\mu$-invariant.

Hence $\hat{P}_{t}$ is $\mu$-invariant. ∎

Appendix F Invariance of the local sampler

The generator of the local BPS is given by

The proof of invariance of the local BPS is very similar to the proof of Propostion 1. We have

where the term (S42) is straightforwardly equal to 0 while, by integration by parts, the term (S40) satisfies

as $h$ is bounded. Now a change-of-variables shows that for any $f\in F$

as $R_{f}^{-1}\left(x\right)v)=R\left(x\right)v$ and $\left\|R_{f}\left(x\right)v\right\|=\left\|v\right\|$ implies $\psi\left(R_{f}\left(x\right)v\right)=\psi\left(v\right)$. So the term (S41) satisfies

where we have used $\left\langle\nabla U_{f}(x),R_{f}\left(x\right)v\right\rangle=-\left\langle\nabla U_{f}(x),v\right\rangle$ and $\max\{0,-f\}-\max\{0,f\}=-f$ for any $f$. Hence, summing (S43)-(S45)-(S42), we obtain$\mathcal{\int L}h(z)\rho\left(z\right){\rm d}z=0$ and the result follows by [7, Proposition 34.7].

Appendix G Calculations in the isotropic normal case

As we do not use refreshment, it follows from the definition of the collision operator that

It follows that $\left\langle x^{(j)},v^{(j)}\right\rangle\leq 0$ for $j>0$ if $\left\langle x^{(0)},v^{(0)}\right\rangle\leq 0$ so, in this case, we have

Appendix H Supplementary information on the evolutionary parameters inference experiments

We consider an over-parameterized generalized time reversible rate matrix with $d=10$ corresponding to 4 unnormalized stationary parameters $x_{1},\dots,x_{4}$, and 6 unconstrained substitution parameters $x_{\{i,j\}}$, which are indexed by sets of size 2, i.e. where $i,j\in\left\{1,2,3,4\right\},\thinspace i\neq j$. Off-diagonal entries of $Q$ are obtained via $q_{i,j}=\pi_{j}\exp\left(x_{\{i,j\}}\right)$, where

We assign independent standard Gaussian priors on the parameters $x_{i}.$ We assume that a matrix of aligned nucleotides is provided, where rows are species and columns contains nucleotides believed to come from a shared ancestral nucleotide. Given $x=\left(x_{1},\dots,x_{4},x_{\{1,2\}},\dots,x_{\{3,4\}}\right),$ and hence $Q$, the likelihood is a product of conditionally independent continuous time Markov chains over $\{$A, C, G, T$\}$, with “time” replaced by a branching process specified by the phylogenetic tree’s topology and branch lengths. The parameter $x$ is unidentifiable, and while this can be addressed by bounded or curved parameterizations, the over-parameterization provides an interesting challenge for sampling methods, which need to cope with the strong induced correlations.

H.2 Baseline

We compare the BPS against a state-of-the-art HMC sampler that uses Bayesian optimization to adapt the the leap-frog stepsize $\epsilon$ and trajectory length $L$ of HMC. This sampler was shown in to be comparable or better to other state-of-the-art HMC methods such as NUTS. It also has the advantage of having efficient implementations in several languages. We use the author’s Java implementation to compare to our Java implementation of the BPS. Both methods view the objective function as a black box (concretely, a Java interface supporting pointwise evaluation and gradient calculation). In all experiments, we initialize at the mode and use a burn-in of 100 iterations and no thinning. The HMC auto-tuner yielded $\epsilon=0.39$ and $L=100$. For our method, we use the global sampler and the global refreshment scheme.

H.3 Additional experimental results

To ensure that BPS outperforming HMC does not come from a faulty auto-tuning of HMC parameters, we look at the ESS/s for the log-likelihood statistic when varying the stepsize $\epsilon$. The results in Figure S3(right) show that the value selected by the auto-tuner is indeed reasonable, close to the value 0.02 found by brute force maximization. We repeat the experiments with $\epsilon=0.02$ and obtain the same conclusions. This shows that the problem is genuinely challenging for HMC.