Concave-Convex PDMP-based Sampling

Abstract

Recently nonreversible samplers based on simulating piecewise deterministic Markov processes (PDMPs) have shown potential for efficient sampling in Bayesian inference problems. However, there remains a lack of guidance on how to best implement these algorithms. If implemented poorly, the computational costs of simulating event times can outweigh the statistical efficiency of the nonreversible dynamics. Drawing on the adaptive rejection literature, we propose the concave-convex adaptive thinning approach for simulating a piecewise deterministic Markov process, which we call CC-PDMP. This approach provides a general guide for constructing bounds that may be used to facilitate PDMP-based sampling. A key advantage of this method is its additive structure—adding concave-convex decompositions yields a concave-convex decomposition. This makes the construction of bounds modular, as given a concave-convex decomposition for a class of likelihoods and a family of priors, they can be combined to construct bounds for the posterior. We show that constructing our bounds is simple and leads to computationally efficient thinning. Our approach is well suited to local PDMP simulation where conditional independence of the target can be exploited for potentially huge computational gains. We provide an R package and compare with existing approaches to simulating events in the PDMP literature. Supplementary materials for this article are available online.

Keywords:

• Algorithms
• Bayesian methods
• Markov chain Monte Carlo (MCMC)
• Piecewise deterministic Markov process (PDMP)

1 Introduction

Monte Carlo methods based on continuous-time Markov processes have shown potential for efficient sampling in Bayesian inference challenges (Fearnhead et al. 2018; Goldman and Singh 2021). These new sampling algorithms are based on simulating piecewise-deterministic Markov processes (PDMPs). Such processes evolve deterministically between random event times with possibly random dynamics at event times (see Davis 1993, for an introduction to PDMPs). Monte Carlo methods based on simulating a PDMP with a desired invariant distribution were first proposed in the computational physics literature where they were motivated as examples of nonreversible Markov processes (Peters and de With 2012; Michel, Kapfer, and Krauth 2014). These ideas transitioned to the statistics community as an alternative to traditional Markov chain Monte Carlo (MCMC) samplers. Popular algorithms in this development include the bouncy particle sampler (Bouchard-Côté, Vollmer, and Doucet 2018) and the Zig-Zag sampler (Bierkens, Fearnhead and Roberts 2019) amongst others (Vanetti et al. 2017; Wu and Robert 2020; Michel, Durmus, and Sénécal 2020; Bierkens et al. 2020).

There is substantial theoretical and empirical evidence to support the claim that nonreversible MCMC samplers offer more efficient sampling (Diaconis, Holmes, and Neal 2000; Bierkens and Roberts 2017; Bouchard-Côté, Vollmer, and Doucet 2018) as nonreversibility reduces the chance of moving back to areas that have recently been visited. Samplers based on PDMPs simulate from a density $\pi (\mathit{\theta })$ by introducing a velocity component v which is used to evolve the position of the state, $\mathit{\theta }$, with deterministic dynamics. At stochastic event times the velocity component is updated using a Markov kernel after which the process continues with the new velocity. The main barrier to implementing a PDMP-based sampler is the simulation of the random event times. This involves simulating the next event time in a time-inhomogeneous Poisson process with rate function that depends on the current position of the sampler.

Simulation of event times in a Poisson process is often facilitated through a method known as thinning (Lewis and Shedler 1979). The aim of this article is to provide a general framework to aid practitioners in implementing PDMP-based samplers. Specifically we introduce the concave-convex thinning approach to facilitate simulating a PDMP-based sampler, which we call CC-PDMP. This method can be applied whenever the rate function of the PDMP can be decomposed into the sum of concave and convex functions, and can be used to facilitate thinning from polynomial rate functions, allowing for broad applicability of the method. We discuss the efficiency of the method compared to alternative approaches in the literature for exact sampling of event times.

Related to our work is the concave-convex adaptive rejection sampler (CC-ARS) of Görür and Teh (2011) which uses the same upper-bounding approach to construct bounds on the log density within a rejection sampler that draws independent samples from univariate density functions. In contrast, by applying this technique within PDMP sampling, the bounds are constructed on the gradient of the log density (i.e., the rate) and used within the nonreversible Markov process framework to generate samples from multivariate densities. Moreover, the CC-PDMP framework allows for both subsampling and local updating schemes to facilitate these methods for high dimensional sampling. A short comparison between concave-convex adaptive rejection sampling and CC-PDMP for sampling from a univariate density is given in Appendix C, supplementary materials.

Some authors have proposed using automatic but approximate methods for simulating a PDMP-based sampler. Among these approaches, Cotter, House, and Pagani (2020) propose numerical integration and root-finding algorithms to facilitate simulation. Others have considered using approximate local bounds which can be simulated exactly (Pakman et al. 2017; Goldman and Singh 2021; Goan et al. 2023). Both of these approaches sacrifice exact sampling of the posterior and involve a tradeoff between the computational cost and the approximation of the sampling distribution.

The article is organized as follows. In Section 2, we introduce the technical details of simulating from a PDMP-based sampler and the details of some popular samplers. Section 3 reviews the literature on thinning for PDMP-based samplers. We introduce the concave-convex PDMP approach for adaptive thinning in Sections 4 and 5. Empirical evaluation of CC-PDMP, and comparison to existing methods, is given in Section 6. We conclude with discussion of limitations and extensions that are possible for the method. Code for implementing our method and replicating our results is available at https://github.com/matt-sutton/ccpdmp.

2 Sampling using a PDMP

2.1 Piecewise Deterministic Markov Processes

Before we explore PDMP-based samplers, we first review PDMPs. A PDMP is a stochastic process, and we will denote its state at time t by ${\mathit{z}}_t$. There are three key components that define the dynamics of a PDMP:

1. A set of deterministic dynamics defined as ${\mathit{z}}_{t+s}=\psi ({\mathit{z}}_t,s)$.

2. Event times that are driven by an inhomogeneous Poisson process with rate $\lambda ({\mathit{z}}_t)$.

3. A Markov transition kernel $q(\mathit{z}|{\mathit{z}}_{t-})$ where ${\mathit{z}}_{t-}:=\underset{s\uparrow t}{\lim }{\mathit{z}}_s$.

A PDMP with these specifications will evolve according to its deterministic dynamics between event times and at event times will update according to the Markov transition kernel $q(\mathit{z}|{\mathit{z}}_{t-})$. The event times must be simulated from a Poisson process with rate $\lambda ({\mathit{z}}_t)$, which depends on the state of the process at time t. The result of simulating this process is a PDMP skeleton ${\{(t_k,{\mathit{z}}_{t_k})\}}_{k=1}^n$ where for $k=1,\dots ,n$, t_k denotes the kth event time and we have stored the state of the process immediately after this event time, ${\mathit{z}}_{t_k}$. Given this skeleton it is possible to fill in the values of the state between the events using the deterministic dynamics.

2.2 PDMP-based Samplers

Assume we wish to construct a PDMP process to sample from a target distribution of interest, $\pi (\mathit{\theta })$. Current PDMP-based samplers are defined on an augmented space $\mathit{z}\in \mathcal{E}$, where $\mathcal{E}=\mathcal{X}\times \mathcal{V}$, that can be viewed as having a position, ${\mathit{\theta }}_t$, and a velocity, ${\mathit{v}}_t$. As described below, the dynamics of the PDMP can be chosen so that the PDMP’s invariant distribution has the form $\nu (\mathit{z})=\pi (\mathit{\theta })p(\mathit{v}|\mathit{\theta })$, for some conditional distribution of the velocities, $p(\mathit{v}|\mathit{\theta })$. As a result, the $\mathit{\theta }$-marginal of the invariant distribution is the target distribution we wish to sample from. We will consider target distributions of the form$$\pi (\mathit{\theta })\propto \hspace{0.17em}\exp \hspace{0.17em}(-U(\mathit{\theta }))$$where $\mathit{\theta }\in {\mathbb{R}}^p$ and $U(\mathit{\theta })$ is the associated potential.

There are two ways to evaluate Monte Carlo estimates of $\int g(\mathit{z})\nu (\mathit{z})\mathrm{d}\mathit{z}$ using the output of a PDMP. The first estimator averages over the continuous time path of the PDMP,$$\frac{1}{t_n-t_1}\sum _{k=1}^{n-1}{\int }_{t_k}^{t_{k+1}}g(\psi ({\mathit{z}}_{t_k},s))ds.$$

This requires that the integral of $g(\psi ({\mathit{z}}_{t_k},s))$ with respect to s may be computed. The second, simpler, approach is to “discretise” the PDMP path. This involves taking an integer M > 0, defining $s=(t_n-t_1)/M$, calculating the value of the state at discrete-time points, ${\mathit{z}}_{t_1+s},\dots ,{\mathit{z}}_{t_1+Ms}$, and then using the familiar Monte Carlo estimator$$\frac{1}{M}\sum _{j=1}^{M}g({\mathit{z}}_{t_1+js}).$$

The points ${\mathit{z}}_{t_1+js}$ may be found by evolving the process along its deterministic dynamics for the appropriate event time interval.

2.3 Bouncy Particle Sampler

The bouncy particle sampler (BPS) takes the velocity space $\mathcal{V}$ to be either ${\mathbb{R}}^p$ or the unit sphere ${\mathcal{S}}^{p-1}$ and targets an invariant distribution $\nu (\mathit{z})=\pi (\mathit{\theta })p(\mathit{v})$ where $p(\mathit{v})$ is either the standard p-dimensional Normal distribution, or a uniform distribution on ${\mathcal{S}}^{p-1}$. The sampler moves with linear deterministic dynamics $\psi (({\mathit{\theta }}_t,{\mathit{v}}_t),s)=({\mathit{\theta }}_t+s{\mathit{v}}_t,{\mathit{v}}_t)$ until a random time t occurs according to a Poisson process with canonical event rate given by$${\lambda }^c({\mathit{z}}_t)=\max \left\{0,〈{\mathit{v}}_t,{\nabla }_{\mathit{\theta }}U({\mathit{\theta }}_t)〉\right\}.$$

At an event time t the velocity is updated as ${\mathit{v}}_t=R_{{\mathit{\theta }}_{t-}}({\mathit{v}}_{t-})$ where$$R_{\mathit{\theta }}(\mathit{v})=\mathit{v}-2\frac{〈\mathit{v},{\nabla }_{\mathit{\theta }}U(\mathit{\theta })〉}{\Vert {\nabla }_{\mathit{\theta }}U(\mathit{\theta }){\Vert }^2}{\nabla }_{\mathit{\theta }}U(\mathit{\theta }).$$

This update can be interpreted as reflecting the velocity off the hyperplane tangential to the gradient of $U(\mathit{\theta })$. Additionally, at random times t an additional event occurs according to a homogeneous Poisson process with rate ${\lambda }^{\text{ref}}>0$ in which the velocity is refreshed, drawing ${\mathit{v}}_t\sim p(\mathit{v})$. Refreshment ensures that the Bouncy Particle Sampler samples the invariant distribution and does not get “stuck” on contours of the potential (Bouchard-Côté, Vollmer, and Doucet 2018). Formally the Bouncy Particle Sampler has event rate $\lambda ({\mathit{z}}_t)={\lambda }^c({\mathit{z}}_t)+{\lambda }^{\text{ref}}$ with Markov transition kernel$$q(\mathit{z}|{\mathit{z}}_{t-})=\frac{{\lambda }^c({\mathit{z}}_{t-})}{\lambda ({\mathit{z}}_{t-})}{\delta }_{{\mathit{\theta }}_{t-}}(\mathit{\theta }){\delta }_{R_{{\mathit{\theta }}_{t-}}({\mathit{v}}_{t-})}(\mathit{v})+\frac{{\lambda }^{\text{ref}}}{\lambda ({\mathit{z}}_{t-})}{\delta }_{{\mathit{\theta }}_{t-}}(\mathit{\theta })p(\mathit{v}),$$where ${\delta }_a(A)$ denotes the Dirac delta function with point mass at a. The BPS is an example of a global PDMP as all components of the velocity v are changed according to the Markov transition kernel.

2.4 Zig-Zag Sampler

The Zig-Zag sampler takes the space of velocities, $\mathcal{V}$, to be ${\{-1,1\}}^p$ with $p(\mathit{v})$ uniform on this space. Similarly to BPS, the Zig-Zag sampler uses linear dynamics $\psi (({\mathit{\theta }}_t,{\mathit{v}}_t),s)=({\mathit{\theta }}_t+s{\mathit{v}}_t,{\mathit{v}}_t)$ between events, but unlike the BPS only a single element of the velocity vector is updated at event times. The overall Zig-Zag process can be viewed as a composition of p local event rates where each event rate has a corresponding Markov transition kernel which operates on a single component of the velocity. Specifically, the ith component of the velocity is updated with local event rate$${\lambda }_i({\mathit{z}}_t)=\max \left\{0,v_i{\partial }_{{\theta }_i}U({\mathit{\theta }}_t)\right\}$$for $i=1,\dots ,p$, where ${\partial }_{{\theta }_i}$ denotes the partial derivative of $U({\mathit{\theta }}_t)$ with respect to θ_i. Simulating an event time in this local framework consists of simulating an event time for each local event rate and then taking the minimum time as the event time t for the overall process. At an event time t triggered by event rate λ_i the velocity is updated as ${\mathit{v}}_t=F_i({\mathit{v}}_{t-})$ where$$F_i(\mathit{v})={(v_1,\dots ,v_{i-1},-v_i,v_{i+1},\dots ,v_p)}^T.$$

Once the velocity is updated for component i local event times must be re-simulated only for rates ${\lambda }_j({\mathit{z}}_t)$ which are functions of θ_i. This fact can induce massive computational savings when there is conditional independence between the components of $\mathit{\theta }$. To see this, consider the extreme case where the target is fully independent across dimensions. In this case each event rate, ${\lambda }_i({\mathit{z}}_t)$, will only depend on the ith components of ${\mathit{\theta }}_t$ and ${\mathit{v}}_t$. Thus, at each event, it is only the time of the next event for the component whose velocity changes that needs to be recalculated. The computational cost of simulation in this case will scale as O(1) per event as opposed to O(d) for the global BPS. This computational advantage was explored in the context of simulating from diffusion bridges in Bierkens et al. (2021).

2.5 Global and Local PDMP Methods

As described above, we can divide PDMP samplers into two main classes: global methods and local methods. These PDMPs differ in how they operate on the velocity vector when an event is triggered. Global methods are PDMP-based methods where the transition kernel acts on the entire velocity vector. By contrast, local methods are defined so that the transition kernel only affects a subset of the velocity. In this section we will give a general algorithm for constructing a local PDMP-based sampler (though see also Bouchard-Côté, Vollmer, and Doucet 2018). Let ${\mathit{v}}^{[S]}$ denote the sub-vector of v with elements indexed by $S\subseteq \{1,\dots ,p\}$ and ${\mathit{v}}^{[-S]}$ denote the sub-vector of v without the elements indexed in S. For some set S define the local kernel as one that satisfies$$q_{\mathit{\theta }}^{\left[S\right]}(\mathit{v}|{\mathit{v}}_{t-})={\delta }_{{\mathit{v}}_{t-}^{[-S]}}\left({\mathit{v}}^{[-S]}\right)q_{\mathit{\theta }}^{\left[S\right]}(\mathit{v}|{\mathit{v}}_{t-}),$$that is, $q_{\mathit{\theta }}^{[S]}(\mathit{v}|{\mathit{v}}_{t-})$ only updates the sub-vector ${\mathit{v}}^{[S]}$. Suppose we have a partition of the components of ${\theta }_1,\dots ,{\theta }_p$ given by $\mathcal{S}=\{S_1,S_2,\dots ,S_F\}$. Subset S_f of v is updated with local kernel $q_{\mathit{\theta }}^{[S_f]}(\mathit{v}|{\mathit{v}}_{t-})$ and associated event rate$${\lambda }_f({\mathit{z}}_t)=\max \left\{0,〈{\mathit{v}}^{\left[S_f\right]},{\nabla }_{{\mathit{\theta }}^{\left[S_f\right]}}U({\mathit{\theta }}_t)〉\right\},$$ for $f=1,\dots ,F$. Simulating a local PDMP with this structure involves simulating an event time for each of the F rates and applying the local transition kernel corresponding to the smallest simulated event time. For the partition $\mathcal{S}$ define $\mathcal{N}=\{N_1,\dots ,N_F\}$ with$$N_f=\{m\in \{1,\dots ,F\left\}\right|{\mathit{\theta }}^{\left[S_m\right]}\mathrm{}/\perp \text{}\perp {\mathit{\theta }}^{\left[S_f\right]}\}$$ for $f=1,\dots ,F$ where ${\mathit{\theta }}^{[S_m]}\mathrm{}/\perp \text{}\perp {\mathit{\theta }}^{[S_f]}$ denotes that the rate of events associated with S_m depends on some element of ${\mathit{\theta }}^{[S_f]}$. When an event is triggered for a local rate f the velocity ${\mathit{v}}^{[S_f]}$ is updated and new event times are simulated for rates which depend on the value of one or more θ_i for $i\in S_f$. This simulation process is summarized in alg1.

In this framework the Zig-Zag sampler is a local PDMP with partition $\mathcal{S}=\left\{\right\{1\},\{2\},\dots ,\{p\left\}\right\}$ and local kernels $q_{\mathit{\theta }}^{[S_f]}(\mathit{v}|{\mathit{v}}_{t-})={\delta }_{F_f({\mathit{v}}_{t-})}(\mathit{v})$. The local BPS uses kernel $q_{\mathit{\theta }}^{[S]}(\mathit{v}|{\mathit{v}}_{t-})={\delta }_{{\mathit{v}}_{*}}(\mathit{v})$ where ${\mathit{v}}_{*}^{[-S]}={\mathit{v}}_{t-}^{[-S]}$ and ${\mathit{v}}_{*}^{[S]}=R_{\mathit{\theta }}^{[S]}({\mathit{v}}_{t-}^{[S]})$ with$$R_{\mathit{\theta }}^{\left[S\right]}({\mathit{v}}^{\left[S\right]})={\mathit{v}}^{\left[S\right]}-2\frac{〈{\mathit{v}}^{\left[S\right]},{\nabla }_{{\mathit{\theta }}^{\left[S\right]}}U(\mathit{\theta })〉}{\Vert {\nabla }_{{\mathit{\theta }}^{\left[S\right]}}U(\mathit{\theta }){\Vert }^2}{\nabla }_{{\mathit{\theta }}^{\left[S\right]}}U(\mathit{\theta }).$$

If there is a single factor $\mathcal{S}=\left\{S\right\}$ containing all indices $S=\{1,\dots ,p\}$ we recover the global BPS.

Algorithm 1

: Simulating a PDMP with local structure

Inputs: $t_1=0,{\mathit{\theta }}_{t_1},{\mathit{v}}_{t_1}$, factorization $(\mathcal{S},\mathcal{N})$ with $\mathcal{S}=\{S_1,S_2,\dots ,S_F\}$ and $\mathcal{N}=\{N_1,N_2,\dots ,N_F\}$.

Sample τ_f for $f=1,\dots ,F$ where$$\mathbb{P}({\tau }_f\ge t)=\hspace{0.17em}\exp \hspace{0.17em}\left(-{\int }_0^t{\lambda }_f({\mathit{\theta }}_0+{\mathit{v}}_{0}s,{\mathit{v}}_0)\mathrm{d}s\right)$$

For $k=1,\dots ,n$

1. Let $f^{*}={\text{argmin}}_f({\tau }_f)$ and update $t_{k+1}=t_k+{\tau }_{f^{*}}$

2. Update states: ${\mathit{\theta }}_{t_{k+1}}={\mathit{\theta }}_{t_k}+{\tau }_{f^{*}}{\mathit{v}}_{t_k}$

3. Update velocity: ${\mathit{v}}_{t_{k+1}}=q_{{\mathit{\theta }}_{t_{k+1}}}^{[S_{f^{*}}]}({\mathit{v}}_{t_k}|\mathit{v})$

4. For $f\in N_{f^{*}}$ update times, that is, resample τ_f where$$\mathbb{P}({\tau }_f\ge t)=\hspace{0.17em}\exp \hspace{0.17em}\left(-{\int }_0^t{\lambda }_f({\mathit{\theta }}_{t_{k+1}}+{\mathit{v}}_{t_{k+1}}s,{\mathit{v}}_{t_{k+1}})\mathrm{d}s\right)$$

Outputs: PDMP skeleton ${\{(t_k,{\mathit{\theta }}_{t_k},{\mathit{v}}_{t_k})\}}_{k=1}^n$

3 Algorithms for Event-Time Simulation

Simulating the first event time, τ, from a Poisson process with event rate $\lambda (t)$ is equivalent to simulating $u\sim \text{Uniform}[0,1]$ and solving$$\tau ={\text{argmin}}_t\left\{-\hspace{0.17em}\log \hspace{0.17em}(u)={\int }_0^t\lambda ({\mathit{z}}_s)ds\right\}.$$

If the event rate is sufficiently simple this can be done exactly. For example if the rate is constant or linear, then there are analytical solutions for τ in terms of u. If the rate function is convex a solution may be found using numerical methods (Bouchard-Côté, Vollmer, and Doucet 2018). For more complex functions, the rate can be simulated by a process known as thinning. This process makes use of Theorem 1.

Theorem 1

(Lewis and Shedler (1979)). Let $\lambda (t)$ and $\overline{\lambda }(t)$ be continuous functions where $\lambda (t)\le \overline{\lambda }(t)$ for all $0\le t_{\max }$. Let ${\tau }_1,{\tau }_2,\dots ,{\tau }_n$ be event times of the Poisson process with rate $\overline{\lambda }(t)$ over the interval $[0,t_{\max })$. For $i=1,\dots ,n$ retain the event times τ_i with probability $\lambda ({\tau }_i)/\overline{\lambda }({\tau }_i)$ to obtain a set of event times from a nonhomogeneous Poisson process with rate λ over the interval $(0,t_{\max }]$.

The efficiency of simulating via thinning depends on how tightly the rate $\overline{\lambda }$ upper-bounds λ and how costly it is to simulate from $\overline{\lambda }$. Another key tool for enabling the simulation of event times is known as super-positioning and the general process is outlined in Theorem 2.

Theorem 2

(Kingman (1992)). Suppose that ${\Lambda }_1,\dots ,{\Lambda }_n$ are a set of independent Poisson processes with rates ${\lambda }_1(t),\dots ,{\lambda }_n(t)$ with first arrival times ${\tau }^{[1]},\dots ,{\tau }^{[n]}$. The Poisson process with rate $\lambda (t)={\sum }_{i=1}^n{\lambda }_i(t)$ may be simulated by returning the first arrival time as $\tau ={\min }_{i=1,\dots ,n}{\tau }^{[i]}$.

Both super-positioning and thinning provide useful tools in the simulation of event times from a nonhomogeneous Poisson process. However, constructing the upper-bounds required for thinning is largely an ad hoc and time consuming process for the practitioner. A common desire for a practitioner is to reuse simulation knowledge. For example if an event rate can be written as the sum of several sub-event rates, that can be exactly simulated, ideally this knowledge should be used to simulate from the sum. We refer to this class of event rates as additive rates. Specifically, we will refer to a Poisson process with a rate$$\lambda ({\mathit{z}}_t)=\max \{0,f_1(t)+f_2(t)\}$$as process with an additive rate. The super-positioning method from Theorem 2 facilitates additive event-time simulation by thinning using an upper-bounding event rate$$\overline{\lambda }(t)=\max \{0,f_1(t)\}+\max \{0,f_2(t)\},$$ which satisfies $\lambda ({\mathit{z}}_t)\le \overline{\lambda }(t)$. We can simulate the first event from a process with this rate as $\tau =\min ({\tau }_1,{\tau }_2)$, where τ₁ is simulated with rate ${\overline{\lambda }}_1(t)=\max \{0,f_1(t)\}$ and τ₂ with rate ${\overline{\lambda }}_2(t)=\max \{0,f_2(t)\}$. This simulated time will be accepted with probability$$\frac{\max \{0,f_1(\tau )+f_2(\tau )\}}{\max \{0,f_1(\tau )\}+\max \{0,f_2(\tau )\}},$$ which will be one when both $f_1(\tau )>0$ and $f_2(\tau )>0$. This procedure is useful since it allows the reuse of thinning procedures for f₁ and f₂. Thinning via super-positioning is often recommended in the literature; notable examples include Bouchard-Côté, Vollmer, and Doucet (2018) who illustrate this approach when simulating from an exponential family and Wu and Robert (2020) who discuss the approach generally. The approach was also recommended broadly by Sen et al. (2020) where it was suggested for practical implementation of the Zig-Zag where f₁ and f₂ correspond to the terms from the likelihood and prior, respectively.

4 Concave-Convex Adaptive Thinning

Our general proposal for simulating events in a PDMP is based on concave-convex adaptive thinning. As the dynamics between events are deterministic, conditional on the current state of the PDMP, we can rewrite the rate of the next event as a function of time. With slight abuse of notation we will denote this rate as $\lambda (t)$ with the argument of the function making it clear whether we are viewing it as a function of the state, or as here, as a function of time. Suppose $\lambda (t)=\max \{0,f(t)\}$ where f(t) can be decomposed as(1) $$f(t)=f_{\cup }(t)+f_{\cap }(t)$$(1) on a finite interval $t\in [0,{\tau }_{\max })$ where $f_{\cup }(t)$ is a convex function and $f_{\cap }(t)$ is a concave function in time. We will later discuss general conditions where such a decomposition is possible. The problem of upper-bounding f(t) is recast to finding upper-bounding piecewise linear functions ${\mathcal{l}}_u(t)$ and ${\mathcal{l}}_n(t)$ such that $f_{\cup }(t)\le {\mathcal{l}}_u(t)$ and $f_{\cap }(t)\le {\mathcal{l}}_n(t)$. We may then apply the bound $f(t)\le \mathcal{l}(t)$ where $\mathcal{l}(t)={\mathcal{l}}_n(t)+{\mathcal{l}}_u(t)$. Since $\mathcal{l}(t)$ is the sum of piecewise linear functions it will also be a piecewise linear function and direct simulation from a Poisson process with rate $\mathcal{l}(t)$ is possible (see Appendix A, supplementary materials). Concave-convex adaptive thinning proceeds by constructing a piecewise linear function $\mathcal{l}(t)$ over a set of m abscissae $t_1=0<t_2<\dots <t_m={\tau }_{\max }$ at which the evaluation of $f_{\cup }(t_i),f_{\cap }(t_i)$ and derivative $f_{\cap }^{\prime }(t_i)$ are known. An event time τ is simulated from the Poisson process with rate $\overline{\lambda }(t)=\max \{0,\mathcal{l}(t)\}$ and is accepted with probability$$\frac{\max \{0,f_{\cup }(\tau )+f_{\cap }(\tau )\}}{\overline{\lambda }(\tau )}.$$

If the event is rejected, we can reuse the information from the evaluation of $f_{\cup }(\tau )$ and $f_{\cap }(\tau )$ with an additional evaluation of $f{\prime }_n(\tau )$ to refine the simulation on the abscissae $t_1=\tau$ to ${\tau }_{\max }$ where all existing evaluations for $\tau <t_k\le {\tau }_{\max }$ may be reused. If an event does not occur on the range $t\in [0,{\tau }_{\max })$ the PDMP process is evolved by ${\tau }_{\max }$ and the thinning process is repeated. We give the general construction of the piecewise linear function $\mathcal{l}(t)$ for $t\in [t_1,t_2)$ and note this construction may be applied iteratively over the abscissae. This construction is also depicted visually in Figure 1. A convex function $f_{\cup }(t)$, is by definition a function that can be upper-bounded by the line segment connecting any two function evaluations. So, for any $t\in [t_1,t_2)$ we have the upper-bound $f_{\cup }(t)\le {\mathcal{l}}_u(t)$ where$${\mathcal{l}}_u(t)=f_{\cup }(t_1)+\frac{f_{\cup }(t_2)-f_{\cup }(t_1)}{t_2-t_1}(t-t_1).$$

Fig. 1 Upper bounds for a rate function based on concave (left) and convex (right) information on three abscissae ($t_i=0,1,2$). Concave bounds are formed using linear segments from the gradient of $f_{\cap }(t_i)$. Convex bounds are constructed using linear segments connecting evaluations of $f_{\cup }(t_i)$.

For a concave function $f_{\cap }(t)$ with derivative $f_{\cap }^{\prime }(t)$ the line segment corresponding to the tangent of $f_{\cap }(t)$ will upper-bound the function. Thus, for $t\in [t_1,t_2),f_{\cap }(t)$ will be upper-bounded by$${\mathcal{l}}_n(t)=\min \left\{f_{\cap }(t_1)+f_{\cap }^{\prime }(t_1)(t-t_1),f_{\cap }(t_2)+f_{\cap }^{\prime }(t_2)(t_2-t)\right\}.$$

This minimum will switch at the point of intersection between the lines $f_{\cap }(t_1)+f_{\cap }^{\prime }(t_1)(t-t_1)$ and $f_{\cap }(t_2)+f_{\cap }^{\prime }(t_2)(t_2-t)$. It is simple to find this intersection point (Görür and Teh 2011),$$t^{*}=\frac{f_{\cap }(t_2)-f_{\cap }^{\prime }(t_2)t_2-f_{\cap }(t_1)+f_{\cap }^{\prime }(t_1)t_1}{f_{\cap }^{\prime }(t_1)-f_{\cap }^{\prime }(t_2)}$$which will be a point on the interval $[t_1,t_2]$ provided the derivatives $f_{\cap }^{\prime }(t_1)$ and $f_{\cap }^{\prime }(t_2)$ are not equal. If $f_{\cap }^{\prime }(t_1)=f_{\cap }^{\prime }(t_2)$, the linear segment will not change over the interval and we take $t^{*}=t_2$. So we take$${\mathcal{l}}_n(t)=\{\begin{array}{cc}{f}_{\cap }(t_1)+f_{\cap }^{\prime }(t_1)t\hfill & t\in [t_1,t^{*})\hfill \\ f_{\cap }(t_2)+f_{\cap }^{\prime }(t_2)(t_2-t)\hfill & t\in [t^{*},t_2)\hfill \end{array},$$ and combining these bounds we have $\mathcal{l}(t)={\mathcal{l}}_u(t)+{\mathcal{l}}_n(t)$ which is a piecewise linear function upper-bounding $f(t)\le \mathcal{l}(t)$ for $t\in [t_1,t_2)$.

4.1 Concave-Convex Decompositions

The proposition below gives simple conditions for the class of nonhomogeneous Poisson processes that admit thinning using the concave-convex approach.

Proposition 1.

If a Poisson process with rate function $\lambda (t)$ can be written as $\lambda (t)=\max \{0,f(t)\}$ where f(t) is continuous then the process admits thinning using concave-convex adaptive thinning.

Proof.

Consider f(t) on the interval $[0,{\tau }_{\max }]$ where ${\tau }_{\max }$ is some arbitrary value ${\tau }_{\max }>0$. On this closed compact interval f(t) is continuous and by the Stone-Weierstrass theorem (Stone 1948), for any $ϵ>0$ there exists a polynomial g(t) on the interval with $\Vert f-g\Vert <ϵ$. The function g is a polynomial so it admits the concave-convex decomposition$$g(t)=\sum _{\{i:a_i>0\}}{a}_i{t}^i+\sum _{\{i:a_i<0\}}{a}_i{t}^i$$with convex function $g_u(t)=\sum _{\{i:a_i>0\}}{a}_i{t}^i$ and concave function $g_n(t)=\sum _{\{i:a_i<0\}}{a}_i{t}^i$. The Process with rate $\lambda (t)$ admits thinning on the interval $[0,{\tau }_{\max })$ using the Poisson process with rate $\widehat{\lambda }(t)=\max \{0,g(t)+ϵ\}$ which has a convex-concave decomposition. If an event does not occur on the interval $[0,{\tau }_{\max })$ then the process is evolved to ${\tau }_{\max }$ and the thinning process repeated. □

Remark 1.

While continuity of f(t) is required for the Stone-Weirestrass theorem, the concave-convex process may also be applied more generally. If f(t) is continuous for all $0<t<t^{*}$ and $\underset{t\to (t^{*})-}{\lim }f(t)=\infty$ then concave-convex thinning can be applied. Consider f(t) using the abscissae $t_0=0,t_1=\frac{t^{*}}{2},\text{tmax}=t^{*}$. Proposition 1 ensures a concave-convex decomposition on the interval $[0,t^{*}/2)$. On $[t^{*}/2,t^{*})$ we upper bound the rate by infinity. So if a time is not simulated on $[0,t^{*}/2)$, we will propose an event at t₁, and this event will be rejected. The concave-convex bound will then be updated based on $f(t^{*}/2)$, and used to simulate an event on $[t^{*}/2,t^{*})$. The abscissae on this increment will be $t^{*}/2, 3t^{*}/4$, and $t^{*}$, and the process will be repeated. Importantly, as $\underset{t\to (t^{*})-}{\lim }f(t)=\infty$ there will be an actual event before $t^{*}$ with probability one, thus, repeating this process will yield an (accepted) event time.

The proof of Proposition 1 relies on finding a polynomial that can upper-bound the process over a finite interval. If a bound can be given for higher order derivatives of f(t) then there are two main approaches for finding an upper-bounding polynomial. Using Taylor’s expansion of f(t) about zero, we get$$\begin{array}{c}f(t)=f(0)+f\prime (0)t+\frac{f″(0)}{2!}{t}^2+\cdots \\ \hspace{1em}+\frac{f^{(k)}(0)}{k!}{t}^k+{\int }_0^t{f}^{(k+1)}(s)\frac{t^k}{k!}ds.\end{array}$$

If there is a known constant M such that $|f^{(k+1)}(t)|\le M$ we can employ thinning based on the polynomial bound $g(t)=f(0)+f\prime (0)t+\frac{f″(0)}{2!}{t}^2+\cdots +\frac{f^{(k)}(0)}{k!}{t}^k+\frac{M}{(k+1)!}{t}^{(k+1)}$. Alternatively polynomial upper-bounds can be constructed based on interpolation where the error is controlled. For example using Lagrange polynomial interpolation on $[0,{\tau }_{\max }]$ with k + 1 evaluations of the function has error bounded by ${\tau }_{{\max }^{k+1}}\frac{M}{(k+1)!}$. Adding this bound as a constant yields an upper-bounding polynomial without requiring evaluation of the derivatives of f(t) required for the Taylor series expansion.

4.2 Concave-Convex Adaptive Thinning and PDMP-based Samplers

A key advantage to CC-PDMP is that it can easily and efficiently simulate from rates with an additive structure. That is, if $f_1(t)$ and $f_2(t)$ both have concave-convex decompositions then $f(t)=f_1(t)+f_2(t)$ will also have a concave-convex decomposition. The rates associated to a PDMP-based sampler with $\pi (\mathit{\theta })\propto \hspace{0.17em}\exp \hspace{0.17em}(-U(\mathit{\theta }))$ as its target will have an additive structure. In particular, a local PDMP sampler as described in Section 2.4 will have rates$${\lambda }_f({\mathit{z}}_t)=\max \left\{0,\sum _{j\in S_f}{f}_j(t)\right\}$$where $f_j(t)=v_j{\partial }_{{\theta }_j}U(\mathit{\theta }+t\mathit{v})$. This means that the concave-convex decompositions for the partial derivatives $f_j(t)$ can be trivially combined to simulate from new local or global implementations of the PDMP. There is a subtle tradeoff between computational and statistical efficiency offered by global and local PDMP-based samplers which we explore in Section 6.

Additive rates are also present when using a PDMP to simulate from a Bayesian posterior distribution. For a model with parameter $\mathit{\theta }\in {\mathbb{R}}^p$, the target has the form$$\pi (\mathit{\theta })\propto p(\mathit{\theta })p(y_{1:n}|\mathit{\theta })$$where $p(y_{1:n}|\mathit{\theta })$ is the likelihood of the observed data $y_{1:n}=(y_1,\dots ,y_n)$, and $p(\mathit{\theta })$ is the prior for the parameters of the model. The function associated to the jth partial derivative of $U(\mathit{\theta })$ can be written as $f_j(t)=f_j^p(t)+f_j^{\mathcal{l}}(t)$ where $f_j^p(t)=v_j{\partial }_{{\theta }_j}\hspace{0.17em}\log \hspace{0.17em}p(\mathit{\theta }+t\mathit{v})$ depends on the log prior and $f_j^{\mathcal{l}}(t)=v_j{\partial }_{{\theta }_j}\hspace{0.17em}\log \hspace{0.17em}p(y_{1:n}|\mathit{\theta }+t\mathit{v})$ depends on the log-likelihood. This allows for the rate function to be split into more manageable pieces. Once a decomposition has been found for a choice of prior or likelihood, it can be reused to facilitate simulating from that prior or likelihood in future models. Specific examples are given in Section 6.

5 Tuning Parameters

5.1 Choice of Abscissae

A key user choice in CC-PDMP sampling is that of the position and number of abscissae on the interval $[0,{\tau }_{\max })$. Suppose that CC-PDMP simulation is implemented with abscissae $t_0=0<t_1<\cdots <t_n={\tau }_{\max }$ if an event does not occur on the interval $[0,{\tau }_{\max })$ the PDMP will be evolved by ${\tau }_{\max }$ and the thinning process will repeat where the function evaluations at ${\tau }_{\max }$ can be reused at $t_0=0$. Thus, one simulation with n abscissae would be computationally equivalent to n evolutions of the CC-PDMP approach using abscissae with only two points $t_0=0$ and $t_1={\tau }_{\max }$. This encourages choosing a minimal number of abscissae and more carefully choosing the length of the interval. There are two main issues that can occur when tuning the parameter ${\tau }_{\max }$. If ${\tau }_{\max }$ is too small, then the proposal frequently lie outside the interval $[0,{\tau }_{\max })$ and simulating an event will requiring many iterations of bounding the event rate. While if ${\tau }_{\max }$ is too large, then the bound on the rate may be poor, leading to many events that are rejected. In this article, we refer to an iteration of the CC-PDMP that does not generate an event time (i.e., a rejected proposal event or not generating on the interval) as a shadow event. The total number of iterations will be the sum of the number of events and shadow events. Efficiency is measured as the proportion of iterations that are events. Adapting the value of ${\tau }_{\max }$ will not change the sampling dynamics. Consequently, unlike in adaptive MCMC, this parameter may be adapted throughout the course of simulating the PDMP. Ideally this parameter should be a little larger than the average event time so that events are proposed on the interval. A simple automatic approach for choosing this parameter is to set ${\tau }_{\max }$ equal to the qth percentile of previous simulated event times, and update ${\tau }_{\max }$ every 100 iterations of the sampler. We found that setting q = 80th percentile worked well for automatically selecting ${\tau }_{\max }$. We investigate the dependency on this tuning parameter by implementing the Zig-Zag sampler on the two-dimensional Banana distribution, with$$U(\mathit{\theta })={({\theta }_1-1)}^2+\kappa {({\theta }_2-{\theta }_1^2)}^2,$$where κ controls how much mass concentrates around the region ${\theta }_2\approx {\theta }_1^2$. Since the potential is polynomial in both parameters it is trivial to see that the Zig-Zag rate functions will also be polynomial. Further implementation details may be found in the supplementary material and associated GitHub.

Figure 2 shows the efficiency of the Zig-Zag sampler on the Banana distribution with κ = 1 when changing the interval length ${\tau }_{\max }$ from 0.1 to 4 with two abscissae. The Zig-Zag algorithm was simulated for 10,000 events at each value of ${\tau }_{\max }$. The effect of selecting ${\tau }_{\max }$ too small or too large is clearly seen and the red line shows the performance using the adaptive approach. The main computational cost is due to an artifact of the otherwise fast C++ code having to call an R function each time the concave-convex decomposition is evaluated. For simple event rates such as polynomials this computation can be implemented in C++ as well. Appendix D, supplementary materials illustrates the effect of ${\tau }_{\max }$ on computation time when the convex-concave decomposition is evaluated in C++.

Fig. 2 Efficiency, defined as the proportion of iterations (events + shadow events) that are events, of Zig-Zag using the CC-PDMP thinning with varying choices for ${\tau }_{\max }$ on the Banana target. Plot on the right shows the actual computation time for running the sampler with different choices of ${\tau }_{\max }$. The red line shows the performance of an automatic choice for this parameter.

5.2 Choice of Decomposition

The choice of concave-convex decomposition is not unique (Görür and Teh 2011). Arbitrary convex functions can be added to $f_{\cup }(t)$ and subtracted from $f_{\cap }(t)$ to give a new valid decomposition. A concave-convex decomposition is minimal if there is no non-affine convex function that can be added to $f_{\cap }(t)$ and subtracted from $f_{\cup }(t)$ while preserving the convexity of the decomposition (Hartman 1959). For example, consider the function$$f(t)=-t^3+3t^2-3t+3,$$where $f″(t)=-6t+6$. It is simple to see that the function is convex for t < 1 and concave for t > 1. A minimal decomposition is given by the piecewise functions$$\begin{array}{c}{f}_{\cup }^1(t)=\{\begin{array}{cc}-t^3+3t^2-3t+3\hfill & t\le 1\hfill \\ 0\hfill & t>1\hfill \end{array}\text{ and }\\ f_{\cap }^1(t)=\{\begin{array}{cc}0\hfill & t\le 1\hfill \\ -t^3+3t^2-3t+3\hfill & t>1\hfill \end{array}.\end{array}$$

Görür and Teh (2011) give a general construction for a minimal concave-convex decomposition. However, this construction relies on finding all points of inflection, which are points where the function changes convexity. Alternatively, Proposition 1 gives a simple decomposition $f_{\cup }^2(t)=3t^2+3$ and $f_{\cap }^2(t)=-t^3-3t$. While the decomposition from Proposition 1 is not minimal, it does not require finding all points of inflection. These two decompositions are shown in Figure 3 using abscissae at 0 and 1.

Fig. 3 The upper-bounds resulting from two different concave-convex decompositions of the function $f(t)=-t^3+3t^2-3t+3$, based on abscissae at 0 and 1. The top row corresponds to the decomposition from Proposition 1 and the bottom row corresponds to an optimal decomposition which recognizes f(t) as a convex function on $[0,1)$. The columns show the functions $f_{\cup },f_{\cap }$ and f as well as their piece-wise linear upper-bounds.

The optimal decomposition gives a tighter bound, here reducing the bounding rate by approximately 0.33 on average across the interval. However, using this bound comes with additional computation cost of finding inflection points. This has the potential to reduce the overall efficiency of the method. Our experience is that the efficiency gains for using the optimal polynomial are generally not sufficient for it to be beneficial. See Appendix D, supplementary materials for further discussion and simulations.

6 Experiments

We now present empirical evaluation of CC-PDMP and comparison with other approaches to simulate PDMPs. Our experiments were implemented using the R package ccpdmp available at https://github.com/matt-sutton/ccpdmp.

The package enables you to simulate a PDMP provided one specifies the concave and convex decomposition of the rate function. If the rate function is polynomial (or bounded by a polynomial) the practitioner may parse this instead and the concave-convex decomposition will be handled internally. The package contains some basic example use cases and code to reproduce the experiments.

6.1 Application in Generalized Linear Models

Generalized Linear Models (GLMs) provide a rich class of models that are frequently used in Bayesian inference. Let ${\{(y_i,{\mathit{x}}_i)\}}_{i=1}^n$ be the observed data, where y_i is an observed response and ${\mathit{x}}_i\in {\mathbb{R}}^p$ is a vector of associated covariates for $i=1,\dots ,n$. The expected value of y_i is modeled by $g^{-1}({\mathit{x}}_i^T\mathit{\theta })$ where $g^{-1}:\mathbb{R}\to \mathbb{R}$ is the inverse link function. The potential of the likelihood has the form$$U^{(\mathcal{l})}(\mathit{\theta })=\sum _{i=1}^n\varphi ({\mathit{x}}_i^T\mathit{\theta },y_i)$$where $\varphi :{\mathbb{R}}^2\to \mathbb{R}$ is the GLM mapping function $\varphi (a,y)=-\hspace{0.17em}\log \hspace{0.17em}p(y|g^{-1}(a))$ returning the negative log-likelihood for observation y given a. The partial derivative with respect to θ_k has the form$${\partial }_{{\theta }_k}{U}^{(\mathcal{l})}(\mathit{\theta })=\sum _{i=1}^n\varphi \prime ({\mathit{x}}_i^T\mathit{\theta },y_i)x_{ik}$$ where $\varphi \prime (a,y)={\partial }_a\varphi (a,y)$ and higher order derivatives are defined similarly. Over the time interval $t\in [0,\tau ]$ let $f_k(t)=v_k{\partial }_{{\theta }_k}{U}^{(\mathcal{l})}(\mathit{\theta }+t\mathit{v})$ which is$$f_k(t)=v_k\sum _{i=1}^n\varphi \prime (a_i(t),y_i)x_{ik}$$ where $a_i(t)={\mathit{x}}_i^T(\mathit{\theta }+t\mathit{v})$. We can use this to define local rates as ${\lambda }_k(t)=\max \{0,f_k(t)\}$ for $k=1,\dots ,p$. For GLMs, repeated application of the chain rule yields:$$f_k^{(j)}(t)=v_k\sum _{i=1}^n{\varphi }^{(j+1)}(a_i(t),y_i){\left(\frac{\partial a_i}{\partial t}\right)}^j{x}_{ik}$$ where $\frac{\partial a_i}{\partial t}={\mathit{x}}_i^T\mathit{v}$ and $f_k^{(j)}$ denotes the jth derivative of $f_k(t)$. When there are bounds on ${\varphi }^{(j)}$ we can use the upper-bounding Taylor polynomial. In Appendix B, supplementary materials we provide bounds for several modeling choices. Here we consider logistic regression, which has the mapping function $\varphi (a,y)=\hspace{0.17em}\log \hspace{0.17em}(1+\hspace{0.17em}\exp \hspace{0.17em}(a))-ya$ where $y\in \{0,1\}$. We look at the efficiency of CC-PMDP thinning for a five dimensional logistic regression problem with n = 200 observations. The covariates were generated from a multivariate normal, ${\mathit{x}}_i\sim \mathcal{N}(\mathbf{0},V^{-1})$ for $i=1,\dots ,200$, with mean zero, precision matrix$$V=\left(\begin{array}{ccccc}1& \rho & 0& 0& 0\\ \rho & 1& 0& 0& 0\\ 0& 0& 1& 0& 0\\ 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 1\end{array}\right)$$ and data generated using $\theta ={(-1.25,0.5,-0.4,-0.4,-0.4)}^T$. We take a Gaussian prior ${\theta }_j\sim \mathcal{N}(0,1)$ for $j=1,\dots ,5.$ As the correlation is increased the thinning becomes more challenging. We investigate the performance of CC-PDMP for increasing polynomial bounds with increasing correlation ρ.

To date only linear bounds (polynomial order 1) for logistic regression have been used for thinning (Bierkens et al. 2018; Bouchard-Côté, Vollmer, and Doucet 2018). These are based on the bound $|\varphi ″(a,y)|\le 1/4$. Table 1 shows the efficiency for thinning using polynomials of order 1–3 (using bounds $|\varphi ″(a,y)|\le 1/4,|\varphi ″\prime (a,y)|\le 1/\left(6\sqrt{3}\right)$ and $|{\varphi }^{(4)}(a,y)|\le 1/8$). The linear 1st order bound matches the bound used in Bierkens et al. (2018) for logistic regression using the Zig-Zag sampler. Table 1 shows that higher order polynomial thinning facilitated through the CC-PDMP allows more efficient thinning. However, these bounds can incur higher computational cost as they require evaluation of the higher order derivatives. For this example an order 2 polynomial appears to have the tradeoff in computation and thinning efficiency giving the best performance. For the linear bounds we used a large fixed ${\tau }_{\max }$ to ensure the linear rate proposal is simulated on the interval $[0,{\tau }_{\max })$. For other polynomial orders we used the adaptive procedure described in the tuning section to select ${\tau }_{\max }$.

Table 1 Efficiency of thinning in Zig-Zag sampling of logistic regression for increasing order of Taylor series polynomial thinning bound.

Polynomial	Correlation
Order	0.00	0.25	0.50	0.65	0.75	0.85	0.95
1	0.64 (4.13)	0.58 (4.41)	0.49 (4.87)	0.41 (5.57)	0.34 (6.47)	0.26 (8.06)	0.13 (14.80)
2	0.77 (4.13)	0.76 (4.12)	0.75 (4.16)	0.72 (4.26)	0.68 (4.42)	0.62 (4.75)	0.39 (6.56)
3	0.80 (5.56)	0.80 (5.56)	0.79 (5.58)	0.78 (5.67)	0.75 (5.76)	0.70 (6.03)	0.47 (7.91)

NOTE: Efficiency is measured as the average proportion of iterations that are events over 20 repetitions of the sampler. Average runtime, in seconds, of each sampler is shown in parentheses.

6.2 Comparison to Thinning via Super-Positioning

In this example we demonstrate the advantages of CC-PDMP when the event rate is additive. In particular, we compare thinning using the CC-PDMP approach with thinning using super-positioning as outlined in Section 3. Let ${\left\{y_i\right\}}_{i=1}^n$ be the observed data with the following model$$\begin{array}{c}{y}_i\mid {\theta }_i\sim \text{Poisson}(\hspace{0.17em}\exp \hspace{0.17em}({\theta }_i)),\\ {\theta }_i\sim \mathcal{N}(0,1)\end{array}$$independently for $i=1,\dots ,n$. The derivative of the potential is ${\partial }_{{\theta }_k}U(\mathit{\theta })={\theta }_k-y_k+\hspace{0.17em}\exp \hspace{0.17em}({\theta }_k)$ and$$f_k(t)=v_k{\partial }_{{\theta }_k}U(\mathit{\theta }+t\mathit{v})=v_k({\theta }_k+v_{k}t)-y_k{v}_k+v_k\hspace{0.17em}\exp \hspace{0.17em}({\theta }_k+v_{k}t).$$

This has the concave-convex decomposition $f_{\cup }(t)=f_k(t)$ when $v_k>0$ and $f_{\cup }(t)=v_k({\theta }_k+v_{k}t)-y_k{v}_k,f_{\cap }(t)=v_k\hspace{0.17em}\exp \hspace{0.17em}({\theta }_k+v_{k}t)$ when $v_k<0$. The global BPS rate can be defined as$$\lambda (t)=\max \left\{0,\sum _{k=1}^n{v}_k{\partial }_{{\theta }_k}U(\mathit{\theta }+t\mathit{v})\right\}=\max \left\{0,\sum _{k=1}^n{f}_k(t)\right\},$$which has concave-convex decomposition defined by the decompositions of the individual f_k functions. In contrast, thinning via super-positioning involves simulating an event time for the linear and exponential terms of the rate individually. The proposed event time is the minimum of all simulated times. This approach to simulating event times for exponential likelihood and Poisson-Gaussian Markov random fields has been previously used in implementing PDMPs (Vanetti et al. 2017; Bouchard-Côté, Vollmer, and Doucet 2018).

Figure 4 compares the thinning and overall computational efficiency for the BPS applied to the Poisson likelihood problem. In the CC-PDMP approach we used the adaptive update for ${\tau }_{\max }$ as described in the tuning section. As the dimension increases we see the proportion of iterations that are events using super-positioning drops quickly, consequently the overall computation time for this approach scales poorly. The poor performance of the super-positioning approach with increasing dimension is the result of the large number of exponential terms, $v_k\hspace{0.17em}\exp \hspace{0.17em}({\theta }_k+v_{k}t)$, in the event rate. The thinning acceptance rate for this proposal is$$\frac{\max \left\{0,\sum _{k=1}^n{v}_k({\theta }_k+v_{k}t-y_k+\hspace{0.17em}\exp \hspace{0.17em}({\theta }_k+v_{k}t))\right\}}{\begin{array}{c}\max \left\{0,\sum _{k=1}^n(v_k({\theta }_k+v_{k}t)-y_k{v}_k)\right\}\hfill \\ \hspace{1em}+\sum _{k=1}^n\max \left\{0,v_k\hspace{0.17em}\exp \hspace{0.17em}({\theta }_k+v_{k}t)\right\}\hfill \end{array}}.$$

Fig. 4 Efficiency of thinning for the global BPS using thinning via super positioning and CC-PDMP thinning with increasing dimension. Efficiency of the thinning (left) and total computation time (right) is averaged over 20 repeated runs. The BPS samplers were simulated for 1000 event times on each run.

When $v_k<0$ these exponential terms contribute zero to the denominator of the super-positioning approach. In comparison the CC-PDMP approach uses a linear upper-bound on these exponential terms which can be negative and reduce the denominator term, leading to more efficient thinning proposals. This is seen in Figure 4 where the thinning efficiency remains roughly constant with increasing dimension.

6.3 Local Methods

Local PDMP methods take advantage of the conditional independence between parameters. Consider the following extension to the previous Poisson example$$\begin{array}{c}{y}_i\mid {\theta }_i\sim \text{Poisson}(\hspace{0.17em}\exp \hspace{0.17em}({\theta }_i))\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{independently for }i=1,\dots ,n\\ {\theta }_1\sim \mathcal{N}(0,1/(1-{\rho }^2))\\ {\theta }_i\mid {\theta }_{i-1}\sim \mathcal{N}(\rho {\theta }_{i-1},1),\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\text{for }i=2,\dots ,n\end{array}$$where we fix $\rho =0.5$. The prior on $\mathit{\theta }$ corresponds to an AR(1) process. The partial derivative of the potential is$${\partial }_{{\theta }_k}U(\mathit{\theta })=\{\begin{array}{cc}(1+{\rho }^2){\theta }_k-\rho {\theta }_{k+1}-y_k+\hspace{0.17em}\exp \hspace{0.17em}({\theta }_k)\hfill & k=1\hfill \\ (1+{\rho }^2){\theta }_k-\rho {\theta }_{k-1}-\rho {\theta }_{k-1}\hfill & \hfill \\ \hspace{1em}-y_k+\hspace{0.17em}\exp \hspace{0.17em}({\theta }_k)\hfill & 1<k<n\hfill \\ {\theta }_k-\rho {\theta }_{k-1}-y_k+\hspace{0.17em}\exp \hspace{0.17em}({\theta }_k)\hfill & k=n\hfill \end{array}$$ which has the same linear and exponential form as the rate in Section 6.2 so we can use an analogous concave-convex decomposition for this rate. In this section we consider local PDMP implementations. The CC-PDMP implementation facilitates simple construction of thinning bounds for local PDMP factorizations. For the partition $\mathcal{S}=\{S_1,\dots ,S_F\}$ the local rate for factor f is$${\lambda }_f({\mathit{z}}_t)=\max \left\{0,\sum _{k\in S_f}{v}_k{\partial }_{{\theta }_k}U({\mathit{\theta }}_t)\right\}=\max \left\{0,\sum _{k\in S_f}{f}_k(t)\right\},$$ for $f=1,\dots ,F$. In this section we consider a number of local decompositions of the form ${\mathcal{S}}^{(j)}=\{S_1,\dots ,S_F\}$ where $S_1=\{1,2,\dots ,j\},S_2=\{j+1,j+2,\dots ,2j\},\dots ,S_F=\{n-j,\dots ,n-1,n\}$. We assume n is divisible by F so each S_f contains j elements. For these decompositions the conditional independence between parameters gives the following neighbors $\mathcal{N}=\{N_1,\dots ,N_F\}$ where$$N_f=\{\begin{array}{cc}\{f,f+1\}\hfill & f=1\hfill \\ \{f-1,f,f+1\}\hfill & 1<f<F\hfill \\ \{f-1,f\}\hfill & f=F\hfill \end{array}.$$

This local decomposition offers computational advantages since updating a single rate requires recalculation of at most three rates regardless of the dimension of the problem. However, this computational efficiency comes at a loss of statistical efficiency. In particular, a global PDMP will be able to move further in stochastic time than local methods before requiring a new event simulation as the global rate function will lower bound the local rate. The overall efficiency of the PDMP method should therefore be considered in terms of both its computational and statistical efficiency.

In Figure 5 we measure overall efficiency as the effective sample size (ESS) per second, the computational efficiency as the number of events or MCMC samples per second and the statistical efficiency as the ESS per MCMC sample or PDMP event. This allows a clear description of the scaling performance of local methods in comparison with alternative state-of-the-art MCMC methods. Linear trends on the log-log scale are stated in the figure legend and used to give approximate asymptotic scaling rates with dimension. We compare local PDMP methods with well-tuned implementations of Metropolis adjusted Langevin algorithm (MALA) and Hamiltonian Monte Carlo (HMC). MALA is tuned to obtain acceptance rate approximately equal to 0.5 and HMC tuned to scale with acceptance probability approximately equal to 0.6. These are close to the optimal values from the MCMC optimal scaling literature (Roberts and Rosenthal 1998; Beskos et al. 2013). For MALA this involves scaling the variance in the proposal and for HMC this involves adjusting the number of leapfrog steps per iteration. These implementations are in line with theoretical scaling results when the parameters of the distribution are independent.

Fig. 5 Log-log breakdown of computational, statistical and overall empirical efficiency scaling with dimension. The ESS values are calculated with respect to the first coordinate θ₁ using the coda R package on a discretized trajectory of the PDMPs. Plotted are the average rates and error bars of all methods calculated from 50 repeated runs of the methods. The legend shows the slope and intercept fitted for each method giving empirical evidence for scaling rates.

The computational efficiency of the local PDMP methods remains (approximately) constant with increasing dimension, while for MALA and BPS the computational efficiency scales approximately as $O(d^{-1})$. The computational efficiency for HMC scales at a rate worse than $O(d^{-1})$ due to the increase in the number of leapfrog iterations. The statistical efficiency for the local methods improves with larger local factors and is highest for the Global BPS. The statistical efficiency for MALA drops at a rate roughly equal to $O(d^{-1/3})$. The Zig-Zag and local BPS methods attain the best overall efficiency scaling rates around O(d) or better. However, it is clear that there is a tradeoff between the statistical efficiency of larger local factors and the increased computational cost. This decomposition can be thought of as an additional choice for PDMP implementation that can easily be tuned using the CC-PDMP approach for thinning. Despite a poorer scaling with dimension, HMC remains competitive with local methods to a very high dimension.

7 Conclusion

PDMP-based samplers have shown advantages over traditional MCMC samplers, but their use has been limited by the perceived difficulty in simulating the event times. We have introduced CC-PDMP as a general approach that can simulate the PDMP provided we can specify a concave-convex decomposition of the event rate. This method has broad applicability, enables simple implementation of local PDMP methods, and empirically outperforms alternative simulation approaches. Additional generalizations of the CC-PDMP approach are also possible by making use of other ideas for adaptive rejection sampling. For example, to bound the concave term, $f_{\cap }(t)$, of the rate function we require that the derivative $f_{\cap }^{\prime }(t)$ is known on the interval $t\in [0,{\tau }_{\max })$. However, we could instead use a looser bound that does not require this derivative information based on the adaptive rejection bounds proposed in Gilks (1992). If the rate function is particularly computationally intensive, a lower bound based on the abscissae where $f_{\cup },f_{\cap }$ and $f{\prime }_n$ and $f_{\cup }^{\prime }$ are evaluated may be used to perform early rejection in the thinning. Additional generalizations may also be found in the difference of convex functions programming literature (Le Thi and Pham Dinh 2018).

A general approach to facilitate thinning using a concave-convex decomposition was described in Section 4. This approach advocates finding an upper-bounding polynomial approximation for the rate function. This upper-bounding polynomial may be thinned using the concave-convex adaptive thinning as described in Proposition 1. It is natural to consider an approximate version of our algorithm where this polynomial is estimated via interpolation. Without controlling for the potential error in the polynomial the overall algorithm would be biased. Recently Corbella, Spencer, and Roberts (2022) have proposed an automatic Zig-Zag implementation which uses optimization methods to obtain a valid linear thinning bound. Future research may consider how to automate the concave-convex thinning procedure for higher order polynomial thinning based on this new work. This may enable more efficient automatic PDMP-based samplers. Finally while our CC-PDMP approach has been illustrated only on PDMPs with linear dynamics the approach could be used more generally on PDMPs with nonlinear dynamics, such as the Boomerang sampler of Bierkens et al. (2020). The only requirement is that the rate function can be bounded by a function with a concave-convex decomposition.

Supplementary Materials

Appendix: Appendix containing additional coding and mathematical details for the work (appendix.pdf)

R-package: R-package “ccpdmp” containing code to perform methods described in the article. The package also contains all data generating code to generate the examples in the article. (ccpdmp.zip, zipped file)

Readme: A markdown file giving instructions for installing the package and replicating the results (README.md)

Additional information

Funding

This research was supported by EPSRC grant EP/R018561 and EP/R034710/1