Accelerating Metropolis-Hastings algorithms by Delayed Acceptance

Introduction

When running an MCMC sampler such as Metropolis–Hastings algorithms (Robert and Casella 2004), the complexity of the target density required by the acceptance ratio may lead to severe slow-downs in the execution of the algorithm. A direct illustration of this difficulty is the simulation from a posterior distribution involving a large dataset of $n$ points for which the computing time is at least of order $\text{O}(n)$. Several solutions to this issue have been proposed in the recent literature (Korattikara et al. 2013, Neiswanger et al. 2013, Scott et al. 2013, Wang and Dunson 2013), taking advantage of the likelihood decomposition

to handle subsets of the data on different processors (CPU), graphical units (GPU), or even computers. However, there is no consensus on the method of choice, some leading to instabilities by removing most prior inputs and others to approximations delicate to evaluate or even to implement.

Our approach here is to delay acceptance (rather than rejection as in Tierney and Mira 1998) by sequentially comparing parts of the MCMC acceptance ratio to independent uniforms, in order to stop earlier the computation of the aforesaid ratio, namely as soon as one term is below the corresponding uniform.

More formally, consider a generic Metropolis–Hastings algorithm where the acceptance ratio $\pi(y)q(y,x)/\pi(x)q(x,y)$ is compared with a $\,\mathcal{U}(0,1)$ variate to decide whether or not the Markov chain switches from the current value $x$ to the proposed value $y$ (Robert and Casella 2004). If we now decompose the ratio as an arbitrary product

where the only constraint is that the functions $\rho_{k}$ are all positive and satisfy the balance condition $\rho_{k}(x,y)=\rho_{k}(y,x)^{-1}$ and then accept the move with probability

i.e. by successively comparing uniform variates $u_{k}$ to the terms $\rho_{k}(x,y)$, the motivation for our delayed approach is that the same target density $\pi(\cdot)$ is stationary for the resulting Markov chain.

In practice, sequentially comparing those probabilities with $d$ uniform variates means that the comparisons stop at the first rejection, implying a gain in computing time if the most costly items in the product (2) are saved for the final comparisons.

Examples of the specific two-stage Delayed Acceptance as defined by Christen and Fox 2005 can be found in Golightly et al. 2014, in the pMCMC context, and in Shestopaloff and Neal 2013.

However, even in these simple cases, it is possible to find examples where Delayed Acceptance may be profitable. Consider for instance resorting to a costly non-informative prior distribution (as illustrated in Section 5.3 in the case of mixtures); here the first acceptance step can be solely based on the ratio of the likelihoods and the second step, which involves the ratio of the priors, does not require to be computed when the first test leads to rejection. Even more often, the converse decomposition applies to complex or just costly likelihood functions, in that the prior ratio may first be used to eliminate values of the parameter that are too unlikely for the prior density. As shown in Figure 1, a standard normal-normal example confirms that the true posterior and the histogram resulting from such a simulated sample are in agreement.

In more complex settings, as for example in “Big Data” settings where the likelihood is made of a very large number of terms, the above principle also applies to any factorisation of the like of (1) so that each individual likelihood factor can be evaluated separately. This approach increases both the variability of the evaluation and the potential for rejection, but, if each term of the factored likelihood is sufficiently costly to compute, the decomposition brings some improvement in execution time. The graphs in Figure 2 illustrate an implementation of this perspective in the Beta-binomial case, namely when the binomial $\mathcal{B}(N,p)$ observation $x=32$ is replaced with a sequence of $N$ Bernoulli observations. The fit is adequate on $10^{5}$ iterations, but the autocorrelation in the sequence is very high (note that the ACF is for the 100 times thinned sequence) while the acceptance rate falls down to 9%. (When the original $y=32$ observation is (artificially) divided into 10, 20, 50, and 100 parts, the acceptance rates are 0.29, 0.25, 0.12, and 0.09, respectively.) The gain in using this decomposition is only appearing when each Bernoulli likelihood computation becomes expensive enough.

On one hand, the order in which the product (3) is explored determines the computational efficiency of the scheme, while, on the other hand, it has no impact on the overall convergence of the resulting Markov chain, since the acceptance of a proposal does require computing all likelihood values. It therefore makes sense to try to optimise this order by ranking the entries in a way that improves the execution speed of the algorithm (see Section 3.2).

We also stress that the Delayed Acceptance principle remains valid even when the likelihood function or the prior are not integrable over the parameter space. Therefore, the prior may well be improper. For instance, when the prior distribution is constant, a two-stage acceptance scheme reverts to the original Metropolis–Hastings one.

Finally, while the Delayed Acceptance methodology is intended to cater to complex likelihoods or priors, it does not bring a solution per se to the “Big Data” problem in that (a) all terms in the product must eventually be computed; and (b) all terms previously computed (i.e., those computed for the last accepted value of the parameter) must be either stored for future comparison or recomputed. See, e.g., Scott et al. 2013, Wang and Dunson 2013, for recent entries on different parallel ways of handling massive datasets.

The plan of the paper is as follows: in Section 2, we validate the decomposition of the acceptance step into a sequence of decisions, arguing about the computational gains brought by this generic modification of Metropolis–Hastings algorithms and further analysing the relation between the proposed method and the Metropolis–Hastings algorithm in terms of convergence properties and asymptotic variances of statistical estimates. In Section 4 we briefly state the relations between Delayed Acceptance and other methods present in the literature. In Section 3 we aim at giving some intuitions on how to improve the behaviour of Delayed Acceptance by ranking the factors in a given decomposition to achieve optimal computational efficiency and finally give some preliminary results in terms of optimal scaling for the proposed method. Then Section 5 studies Delayed Acceptance within three realistic environments, the first one made of logistic regression targets, the second one alleviating the computational burden from a Geometric Metropolis adjusted Langevin algorithm and a third one handling an original analysis of a parametric mixture model via genuine Jeffreys priors. Section 6 concludes the paper.

Validation and convergence of Delayed Acceptance

In this section, we establish that Delayed Acceptance is a valid Markov chain Monte Carlo scheme and analyse on a theoretical basis the differences with the original version.

We assume for simplicity that the target distribution $\pi$ and the proposal distributions $Q(x,\cdot)$ all admit densities w.r.t. the Lebesgue or counting measures. We also denote by $\pi$ the target density and let $q(x,y)$ denote the proposal density.

Let $(X_{n})_{n\geq 1}$ be a Markov chain evolving on $\mathsf{X}$ with Metropolis–Hastings Markov transition kernel $P$ associated with $q$ and $\pi$, i.e. for $A\in\mathcal{B}(\mathsf{X})$

Above, $\alpha(x,y)$ is known as the Metropolis–Hastings acceptance probability and $r(x,y)$ as the Metropolis–Hastings acceptance ratio.

We consider the class of “Delayed acceptance” Markov kernels associated with $P$, which are defined by factorisations of the function $r$ as

with all components in the product satisfying $\rho_{k}(x,y)=\rho_{k}(y,x)^{-1}$. The associated Delayed Acceptance Markov kernel is then defined as

The order in which the sequence of functions $\rho_{k}$ appears in the factorisation (4) is important for algorithmic specification, as can be seen in Algorithm 1. It means that $\rho_{1}(x,Y)$ is evaluated first, $\rho_{2}(x,Y)$ second, and so on until $\rho_{d}(x,Y)=r(x,Y)/\prod_{k=1}^{d-1}\rho_{k}(x,Y)$ which is last, with the motivation that “early rejection” can allow computational savings by avoiding the computation of the subsequent $\rho_{k}(x,Y)$.

2 Validation

The first lemma is a standard representation leading to the validation of the Delayed Acceptance Markov chain:

For any Markov chain with transition kernel $\Pi$ of the form

and satisfying detailed balance, the function $a(\cdot)$ satisfies (for $\pi$-a.a. $x,y$)

This follows immediately from the detailed balance condition

Given a probability measure $\mu$, let us denote

For a generic Markov kernel $\Pi:\mathsf{E}\times\mathcal{B}(\mathsf{E})$ with unique invariant probability measure $\mu$, we define

where $(X_{n})_{n\geq 1}$ is a Markov chain with Markov kernel $\Pi$ initialised with $X_{1}\sim\mu$.

For any $f\in L^{2}(\mathsf{E},\mu)$ we define the Dirichlet form associated with a $\mu$-reversible Markov kernel $\Pi:\mathsf{E}\times\mathcal{B}(\mathsf{E})$ as

The (right) spectral gap of a generic $\mu$-reversible Markov kernel has the following variational representation

which leads to the following domination lemma:

Let $\Pi_{1}$ and $\Pi_{2}$ be $\mu$-reversible Markov transition kernels of $\mu$-irreducible and aperiodic Markov chains, and assume that there exists $\varrho>0$ such that for any $f\in L_{0}^{2}\bigl(\mathsf{E},\mu\bigr)$

We will need the following condition in the sequel:

Defining $A:=\{(x,y)\in\mathsf{X}^{2}:r(x,y)\geq 1\}$, there exists a $c$ such that

$\inf_{(x,y)\in A}\min_{k\in\{1,\ldots,d\}}\rho_{k}(x,y)\geq c$.

Since $r(x,y)\geq 1$, we have $\alpha(x,y)=1$. On the other hand, for $(x,y)\in A$

since at least one $\rho_{k}(x,y)\geq 1$ whenever $r(x,y)\geq 1$.

4 Modification of a given factorisation

The easiest way to use the above result is to modify any candidate factorisation. Given a factorisation of the function $r$

satisfying the balance condition, we can define a sequence of functions $\rho_{k}$ such that both $r(x,y)=\prod_{k=1}^{d}\rho_{k}(x,y)$ and Condition 1 holds. To that effect, take an arbitrary $c\in(0,1]$ and define $b:=c^{\frac{1}{d-1}}$. Then, if we set

From this modification, we deduce the following result:

We note that $\inf_{(x,y)\in\mathsf{X}^{2}}\prod_{k=1}^{d-1}1\wedge\rho_{k}(x,y)\geq b^{d-1}=c$ and that

5 Example: unmodified surrogate targets

One common usage (Christen and Fox 2005) of Delayed Acceptance is to substitute a surrogate target $\bar{\pi}$ for $\pi$ in $\rho_{1}(x,y)$. We consider the case $d=2$ and a random walk proposal to examine Condition 1 in this context. Here we have $q(x,y)=q(y,x)$ and so

Considering $(x,y)\in A=\{(x,y)\in\mathsf{X}^{2}:r(x,y)\geq 1\}$ we require $c>0$ satisfying simultaneously

The first of these says that $\bar{\pi}(y)/\bar{\pi}(x)$ cannot be too small when $\pi(y)\geq\pi(x)$. The second says that $\bar{\pi}(y)/\bar{\pi}(x)$ should not be a large multiple of $\pi(y)/\pi(x)$. There are a large variety of choices of $\bar{\pi}$ that allow one to take $c=1$. For example, $\bar{\pi}(x)=\pi(x)+C$ for some constant $C\geq 0$ and $\bar{\pi}(x)\propto\pi(x)^{\beta}$ for some $\beta\in$. Note that $\beta=0$ corresponds to $\bar{\pi}$ being a constant function and $\beta=1$ corresponds to $\bar{\pi}\propto\pi$. In between, one can think of $\bar{\pi}$ as being a flattened version of $\pi$.

6 Counter-example: failure to reproduce geometric ergodicity

Mengersen and Tweedie 1996 showed that a random-walk Metropolis–Hastings chain for targets with super-exponential tails is geometrically ergodic. We now exploit this result to derive that, if $\sigma^{2}<1$, then the unmodified delayed acceptance Markov chain is not geometrically ergodic.

The unmodified Delayed Acceptance Markov chain using the factorisation into $\rho_{1}$ and $\rho_{2}$ as above is not geometrically ergodic when $\sigma<1$.

since $\log(\tau/3)<0$ and $\sigma^{2}<1$. Therefore,

Optimisation

When considering Markov Chain Monte Carlo methods in practice, their efficiency as measured by mixing properties and computational cost is a fundamental issue. This section addresses both perspectives in connection with Delayed Acceptance. Section 3.1 examines the proposal distribution and derives its optimal scaling from standard random–walk Metropolis–Hastings theory. Then Section 3.2 covers the ranking of the factors $\rho_{i}$, which drives the total computational cost of the procedure.

The explorative performances of a random–walk MCMC are strongly dependent on its proposal distribution. As exemplified in Roberts et al. 1997, finding the optimal scale parameter does lead to efficient ‘jumps’ in the state space and the overall acceptance rate of the chain is directly connected to the average jump distance and to the asymptotic variance of ergodic averages. This provides practitioners with an approach to ‘auto-tune’ the resulting random–walk MCMC algorithm. Extending this calibration to the Delayed Acceptance scheme is equally important, on its own ground towards finding a reasonable scaling for the proposal distribution and to avoid comparisons with the standard Metropolis–Hastings version.

The original framework of Roberts et al. 1997 is cantered on estimating a collection of expected functionals, say $g$, where a plausible criterion for the performances of the MCMC is the minimisation of the stationary integrated auto-correlation time (ACT) of the Markov chain, defined as

where the index $g$ stresses the dependence on the considered functional, which is connected to the asymptotic variance through $\text{var}(P,g)=\tau_{g}\times\text{var}_{\pi}(g)$ whenever the chain is $\phi$-irreducible, aperiodic, and reversible, $\text{var}_{\pi}(g)$ is finite and $g\in L^{2}(\pi)$.

Research on this optimisation focus on two main cases:

Sherlock and Roberts 2009 focus on unimodal elliptically symmetric targets and show that a proxy for the ACT in finite dimensions is the Expected Square Jumping Distance (ESJD), defined as

where $X$ and $X^{\prime}$ are two successive points in the chain and $\|\cdot\|_{\beta}$ represent the norm on the principal axes of the ellipse rescaled by the coefficients $\beta_{i}$ so that every direction contributes equally.

An interesting result in Sherlock and Roberts 2009 is that, as $d\rightarrow\infty$, the ESJD on one marginal component of the chain converges with the same speed as the diffusion process described in Roberts et al. 1997. These authors furthermore show the asymptotic result holds for rather small dimension, roughly starting from $d=5$.

Moreover, when considering efficiency for Delayed Acceptance, which is a technique tailored on costly computations, we need to focus on the execution time of the algorithm as well. We then proceed to define our measure of efficiency as

similarly to Sherlock et al. 2013 for Pseudo-Marginal MCMC.

Due to the complex acceptance ratio in Delayed Acceptance, an extension of the previous results requires rather stringent assumptions, albeit providing a proper guideline in practice. Section 5 will further demonstrate optimality extends beyond those conditions. Note that our assumptions are quite standard in the literature on the subject.

(H1) We assume for simplicity’s sake that the Delayed Acceptance procedure operates on two factors only, i.e., that $r(x,y)=\rho_{1}(x,y)\times\rho_{2}(x,y)$. The acceptance probability of the scheme is thus

(H2) We further assume that the target distribution satisfies (A1) and (A2) in Roberts et al. 1997, which are regularity conditions on $\pi$ and its first and second derivatives, and that $\pi(x)=\prod\limits_{i=1}^{n}f(x_{i})$.

The second part of the lemma follows directly from Theorem 1.1 in Roberts et al. 1997. ∎

Let us stress that almost all assumptions in the above Lemma can be relaxed and that performances are robust against small deviances from those assumptions, as shown by the literature on standard Metropolis–Hastings. Obtaining analytical results without such conditions, while possible, requires however a considerable mathematical effort.

We now state the main practical implication of Lemma 4.

If the conditions of Lemma 4 holds, the optimal average acceptance rate $\alpha^{*}(\delta)$ is independent of $I$.

2 Ranking the Blocks

Consider an i.i.d. framework; the target (in $x$) is given by

where $Z=(z_{1},\dots,z_{n})$ is an i.i.d. sample from $f(z|x)$ and $p(x)$ is the prior distribution for $x$, we can always consider the decomposition

where each $\xi_{i}(x,y)$ is made of a small number of density ratio terms, with one including the prior and proposal ratios. In the limit, it is feasible if not necessarily efficient to consider the case $K=n+1$ with

Assuming the computing cost is comparable for all terms, a solution for optimising the order of these factors ranks the entries according to the success rates observed so far, starting with the least successful values. Alternatively, the factorisation can start with the ratio that has the highest variance, since it is the most likely to be rejected. (Note however that poor factorisations (6) lead to very low acceptance rates, as for instance when picking only outliers in a given group of observations.) Lastly, we can rank factors by their correlation with the full Metropolis–Hastings ratio; taking the argument to the limit, if the first factor has a perfect correlation with $r(\cdot,\cdot)$ then all the successive terms must be accepted and their order is hence of no interest.

As all these features can be evaluated for each subsample while running a chain with acceptance ratio factored as in (6), an implementation based on this intuition is then to take

A relevant warning is that if we rearrange terms during the run, not only reordering but also merging them, in accordance to their correlation with the unmodified ratio, the resulting method has no theoretical guarantee since the kernel is potentially changing at each iteration depending on properties of previous samples (Gelfand and Sahu 1994).

As with standard adaptive MCMC (Roberts and Rosenthal 2005) we resort thus to a finite adaptation scheme; we start with a fixed number of iterations to rank and rearrange the factors, followed by a fixed ordering to achieve ergodicity of the chain. We test this procedure in Section 5.1 on a simulated example.

Finally note that while we focused on the i.i.d. setting, in more complex cases where the ratio is factored and Delayed Acceptance can be applied, it is often the case that the optimal ordering of such factors is already known.

Relation with other methods

Prefetching, as defined by Brockwell 2006, is a programming method that accelerates the convergence of a single MCMC chain by guessing future states in the path of a random walk Metropolis–Hastings Markov chain in order to use any additional computing power available, in the form of extra parallel processors, to calculate in advance necessary quantities (like the Metropolis–Hastings ratio) so that when the chain reaches a given state the computationally-heavy part of that iteration are ready.

Since an in-depth exploration of prefetching is outside the scope of this work the reader is referred to Strid 2010 and citations therein for a complete discussion of the argument.

2 Alternative procedure for Delayed Acceptance

In the case that every factor $\rho_{i}(x,y)$ has roughly the same computational cost, Philip Nutzman suggested (personal communication) that Delayed Acceptance can be slightly modified by taking the overall acceptance probability

Such a decomposition follows from the same idea that one would like to compute as few factors as possible once one realizes that the proposal is likely to be rejected. Under this modification the associated Markov chain still achieves the correct target in the stationary regime and the procedure satisfies detailed balance, provided the ordering of the terms is uniformly random.

An interesting consequence of this modification is that, as the number of factor increases, the acceptance rate eventually stabilises, while for the method described in Section 1 the acceptance rate decreases to zero. This property is indeed appealing, even thought this procedure logically takes longer to complete when compared with the standard Delayed Acceptance (albeit less than the reference Metropolis–Hastings procedure).

The evident disadvantage of the modification in a general setting is that detailed balance implies that the factors are computed in a random order at each iteration, making vain any attempt to adapt in terms of the ordering (Section 3.2) or to set the order based on respective computational costs.

This drawback can be somewhat reduced by combining the above two approaches; consider the decomposition

where $d_{1}+d_{2}=d$ and the factors $\rho_{i}$ and $\phi_{j}$ represent respectively cheap factors and costly factors. By taking now

the algorithm computes cheap factors first and expensive factors last, applying the symmetry requirement to satisfy detail balance inside each of both subsets. Clearly the above can be generalised to a larger number of subsets, each with $d_{i}$ factors in it. Intuitively, this last modification can be explained as an early rejection of each of the intermediate acceptance/rejection steps inside a Delayed Acceptance scheme.

Interestingly if $d_{l}=1\ \forall l$ ($l$ being the number of subsets considered) this procedure reduces to Delayed Acceptance, and for $l$ that increases and $d_{l}>1\ \forall l$ this combined technique will have a even lower overall acceptance rate than standard Delayed Acceptance.

3 Delayed Acceptance and Slice Sampling

As a final remark, we stress another analogy between our Delayed Acceptance algorithm and slice sampling (Neal 1997, Robert and Casella 2004). Based on the same decomposition (1), slice sampling proceeds as follows

to compare with Delayed Acceptance which conversely

simulate $\theta^{\prime}\sim q(\theta^{\prime}|\theta)$;

The differences between both schemes are thus that (a) slice sampling always accepts a move, (b) slice sampling requires the simulation of $\theta^{\prime}$ under the constraints, which may prove infeasible, and (c) Delayed Acceptance re-simulates the uniform variates in the event of a rejection. In this respect, Delayed Acceptance appears as a “poor man’s” slice sampler in that values of $\theta^{\prime}s$ are proposed until one is accepted.

Examples

To illustrate the improvement brought by Delayed Acceptance, we study three different realistic settings that reflect on the generality of the method. First, in Section 5.1 we consider a Bayesian analysis of a logistic regression model, to assess the computational gain brought by our approach in a “Big-Data” environment where obtaining the likelihood is the main computational burden. Secondly (Section 5.2) we examine a high dimensional toy Normal-Normal model, sample with a geometric Metropolis adjusted Langevin algorithm where the main computational cost comes from the proposal distribution which is position specific and involves derivatives of the density up till the third level, which are computed numerically at each iteration. Finally in Section 5.3 we investigate a mixture model where a formal Jeffreys prior is used, as it is not available in closed-form and does require an expensive approximation by numerical or Monte Carlo means. The latter example comes as a realistic setting where the prior itself is a burdensome object, even for small datasets.

While a simple model, or due to its simplicity, logistic regression is widely used in applied statistics, especially in classification problems. The challenge in the Bayesian analysis of this model is not generic, since simple Markov Chain Monte Carlo techniques providing satisfactory approximations, but stems from the data-size itself. This explains why this model is used as a benchmark in some of the recent accelerating papers (Korattikara et al. 2013, Neiswanger et al. 2013, Scott et al. 2013, Wang and Dunson 2013). Indeed, in “big Data” setups, MCMC is deemed to be progressively inefficient and researchers are striving to keep simulation effective, focusing mainly on parallel computing and on sub-sampling but also on replacing the classic Metropolis scheme itself.

We tested the proposed method against the standard Metropolis–Hastings algorithm on $10^{6}$ simulated data with a $100$-dimensional parameter space. The proposal distribution is Gaussian: $y|x\sim\mathcal{N}(x,\Sigma)$ with $\Sigma$ initialised to be $0.2\times I_{d}$ ($d$ being the dimension of the parameter space) and then adapted. The Metropolis–Hastings benchmark was made adaptive by targeting the asymptotic optimal acceptance rate of $\alpha^{*}=0,234$ (Roberts et al. 1997).

We collected $100$ repetitions of the experiment and the results are presented in Table 1. Before commenting the results we highlight the fact that this situation may seem not particularly appealing for Delayed Acceptance and in fact straight application of the method by randomly choosing the composition of $\rho_{1}$ and $\rho_{2}$ may lead to variable results. Further coding effort is required here in order to choose adaptively how to split the MH ratio. Borrowing from both Section 3.2 and the end of Section 3.1, i.e. by choosing during the brief burn-in of the chain which subset best represents the whole likelihood and then, based on how populated that subset is, targeting a specific acceptance ratio, produces both a completely automated MCMC version for this kind of data (iid) and better results under a time constraint.

As shown in Table 1, while the assumption made in Section 3 not completely satisfied, the relative efficiency of Delayed Acceptance is higher that for MH by a factor of almost $6$. We measured efficiency trough effective sample size (ESS, from the coda R package (Plummer et al. 2006)) or expected square jumping distance (ESJD). By choosing the first subsample to be informative on the whole ratio there is practically no loss on ESS (while the estimated ESJD actually increased) and, given the significantly reduced acceptance rate, the computing time is usually less then a fourth of the computing time of the corresponding optimal MH, taking into account the first part of chain used to determine the blocks ranking.

2 G-MALA with Delayed Acceptance

Random walk Metropolis–Hastings, while generic and popular, can struggle with posterior distributions in high dimensions or in the presence of high correlation between some components. In such cases it is inefficient, with low acceptance rate, poor mixing and highly correlated samples. Metropolis adjusted Langevin algorithm (MALA, see for instance Roberts and Stramer 2002) has been devised to overcome these difficulties by taking advantage of the gradient of the target distribution in the proposal mechanism, making the Markov chain more robust with respect to the dimension of the problem and proposing broader moves with higher probability. MALA is based on a Langevin diffusion, with the target (the posterior distribution $\pi(\theta|y)$ in our case) as a stationary distribution, defined by the SDE

where $B$ is a Brownian motion. Using a first-order discretisation the diffusion gives the following proposal mechanism:

where $\varepsilon$ is the step-size for the Euler’s integration and $Z\sim\mathcal{N}(\mathbf{0},\mathbf{I})$. This discretisation is then compensated by introducing an accept/reject probability similar to a Metropolis–Hastings algorithm.

This diffusion is isotropic and will hence still be inefficient for highly correlated components or with very different scales, as the step size $\varepsilon$ is fixed across dimensions. Roberts and Stramer 2002 propose to alleviate the issue using a pre-conditioning matrix $A$ so that the proposal becomes

Christensen et al. 2003 demonstrate however that defining this matrix in general can be difficult and that tuning on the go may result in an inappropriate asymptotic behaviour.

In a recent work Girolami and Calderhead 2011 propose the Geometric-MALA in order to overcome this difficulty, advising the use of a position specific metric for the matrix $A$, which takes advantage of the geometry of the target space that the chain is exploring. They suggest in particular the Fisher-Rao metric tensor. In terms of Bayesian inference, where the target distribution is the posterior density, this choice translates into $A^{T}A$ being the expected Fisher information matrix plus the negative Hessian of the log-prior.

This theoretically efficient solution also performs well in practice but comes with a serious computational burden in the fact that at every evaluation of the Metropolis–Hastings ratio derivatives up till the third order of our log-target distribution are needed and, in the event of them being analytically not available, expensive numerical approximations are to be computed (see equation (10) of Girolami and Calderhead 2011).

2.2 Sampling with Delayed Acceptance and GMALA:

Geometric-MALA represent a perfect application for Delayed Acceptance since we can naturally divide its acceptance ratio into the product of the posterior ratio and the ratio of the proposals, the latter to be only computed when the proposed point is associated with a relatively large posterior probability.

As described above, the computational bottleneck of the G-MALA lays in the computation of the third derivative of our log-target at the proposed point, while the computation of the posterior itself has usually a low relative cost. Moreover even with a non-symmetric efficient proposal mechanism (the discretised Langevin diffusion) G-MALA is still close to a random walk and we expect the ratio of the proposal to be near 1, especially at equilibrium especially when $\varepsilon$ is small. Therefore, the first ratio is inexpensive, relative to the second one, while the decision reached at the first stage should be consistent with the overall acceptance rate.

2.3 Simulation study:

To test the above assumptions we ran a toy MALA example where we drew $100$ samples from a $\mathcal{N}_{d}(\theta,I)$, with $d=10$; $\pi(\theta)$ was set to be $\mathcal{N}_{d}(0,100)$. Figure 5 presents an example run. We then repeated the experiment $100$ times and computed an average efficiency gain, defined either as ESS or as the ESJD, over the computing time. We computed $\delta$ at each run by averaging a few computed derivatives, required by the proposal ratio. We then adapt $\varepsilon$ to get the optimal acceptance rate, being conservative in order to avoid overflow issues with the first-order numerical integrator. Results are presented in Table 2. Delayed Acceptance exhibits improvement by a factor of $10$ in this example, obtained almost for free in terms of to coding time.

2.4 HMC with Delayed Acceptance:

As a side note, while the reasoning applied to MALA does theory apply to Hamiltonian Monte Carlo (HMC), the computational gain obtained through Delayed Acceptance is only connected with avoiding some proposal computations. In a general HMC though (with both point-dependent and independent pre-conditioning matrices), proposing a new value still involves the computation of $L-1$ (with $L$ the number of steps in the discretised–Hamiltonian integration) derivatives, as only the starting point is recovered from the previous iteration, while computing the final Metropolis–Hastings ratio involves just the extra computation at the end point. Therefore, in this setting, the computational gain is much reduced.

3 Mixture Model

Consider a standard mixture model (MacLachlan and Peel 2000) with a fixed number of components

This standard setting nonetheless offers a computational challenge in that the reference objective Bayesian approach based on the Fisher information and the associated Jeffreys prior (Jeffreys 1939, Robert 2001) is not readily available for computational reasons and has thus not been implemented so far. Proxys using Jeffreys priors on the components of (7) have been proposed instead, with the drawback that since they always lead to improper posteriors, ad hoc corrections have to be implemented (Diebolt and Robert 1994, Roeder and Wasserman 1997, Stephens 1997).

When relying instead on dependent improper priors, it is not always the case that the posterior distribution is improper. For instance, Robert and Titterington 1998 provide a location-scale representation that allows for some improper prior. In the current paper, we consider instead the genuine Jeffreys prior for the complete set of parameters in (7), derived from the Fisher information matrix for the whole model. While establishing the analytical properness of the associated posterior is beyond the goal of the current paper, we handle large enough samples to posit that a sufficient number of observations is allocated to each component and hence the likelihood function dominates the prior distribution. (In the event the posterior remains improper, the associated MCMC algorithm should exhibit a transient behaviour.)

Therefore, this is an appropriate and realistic example for evaluating Delayed Acceptance since the computation of the prior density is clearly costly, relying on many integrals of the form:

3.2 Simulation study:

An experiment comparing a standard Metropolis–Hastings implementation with a Metropolis–Hastings version relying on Delayed Acceptance is summarised in Table 3. Data were simulated from the following Gaussian mixture model:

Both the standard Metropolis–Hastings and the Delayed Acceptance version are adapted against their respective optimal acceptance rate, which is computed to be $2\%$, given that $\delta$ is empirically established to be $0.01$ using $500$ samples for the Monte Carlo estimation of the prior. As a consequence the MH+DA algorithm will produce less unique samples in the total $10^{5}$ iterations of the chain, as reflected in the lesser ESS in Table 3, but this is counterbalanced by the impressive decrease in computing time, leading again to an overall gain in terms of $\text{ESS}/t$ of about $9$.

Conclusion

We introduced in this paper Delayed Acceptance, a generic and easily implemented modification of the standard Metropolis–Hastings algorithm that splits the acceptance rate into more than one step in order to increase the computational efficiency of the resulting MCMC, under the sole condition that the Metropolis–Hastings ratio can be factorised this way.

The choice of splitting the target distribution into parts ultimately depends on the respective costs of computing the said parts and of reducing theoretically the overall acceptance rate and expected square jump distance (ESJD). Still, this generic alternative to the standard Metropolis–Hastings approach can be considered on a customary basis, since it both requires very little modification in programming and can be easily tested against the basic version, both empirically and theoretically by the results of (2). The Delayed Acceptance algorithm presented in (1) can significantly decrease the computational time per se as well as the overall acceptance rate. Nevertheless, the examples presented in Section 5 suggest that the gain in terms of computational time is not linear in the reduction of the acceptance rate, especially in the presence of optimisation techniques like (3).

Furthermore, our Delayed Acceptance algorithm does naturally merge with the widening range of prefetching techniques, in order to make use of parallelisation and reduce the overall computational time even more significantly. Most settings of interest are open to take advantage of the proposed method, if mostly in the situation of Bayesian statistics where the target density and/or the Metropolis–Hastings ratio always allow for a natural factorisation. The case when the likelihood function can be factorised in an useful way represents the best gain brought by our solution, in terms of computational time, and it may easily improve even more by exploiting parallelisation techniques.

Thanks to Christophe Andrieu for a very helpful discussion on an earlier version of the manuscript. The massive help provided by Jean-Michel Marin and Pierre Pudlo towards an implementation on a large cluster has been fundamental in the completion of this work. Thanks to Samuel Livingstone for suggesting the Geometric MALA example and finally thanks to Philip Nutzman for the interesting conversation and for the suggestion of the method proposed in (4.2). Christian P. Robert research is partly financed by Agence Nationale de la Recherche (ANR, 212, rue de Bercy 75012 Paris) on the 2012–2015 ANR-11-BS01-0010 grant “Calibration” and by a 2010–2015 senior chair grant of Institut Universitaire de France. Marco Banterle PhD is funded by Université Paris Dauphine.

References