Accelerating Metropolis-Hastings algorithms by Delayed Acceptance

MCMC algorithms such as Metropolis-Hastings algorithms are slowed down by the computation of complex target distributions as exemplified by huge datasets. We offer in this paper a useful generalisation of the Delayed Acceptance approach, devised to reduce the computational costs of such algorithms by a simple and universal divide-and-conquer strategy. The idea behind the generic acceleration is to divide the acceptance step into several parts, aiming at a major reduction in computing time that out-ranks the corresponding reduction in acceptance probability. Each of the components can be sequentially compared with a uniform variate, the first rejection signalling that the proposed value is considered no further. We develop moreover theoretical bounds for the variance of associated estimators with respect to the variance of the standard Metropolis-Hastings and detail some results on optimal scaling and general optimisation of the procedure. We illustrate those accelerating features on a series of examples


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 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 (1) n i=1 (θ|x i ) 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 π(y)q(y, x)/π(x)q(x, y) is compared with a 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 (2) π(y) q(y, x) π(x)q(x, y) where the only constraint is that the functions ρ k are all positive and satisfy the balance condition ρ k (x, y) = ρ k (y, x) −1 and then accept the move with probability (3) min {ρ k (x, y), 1} , i.e. by successively comparing uniform variates u k to the terms ρ k (x, y), the motivation for our delayed approach is that the same target density π(•) is stationary for the resulting Markov chain.
The mathematical validation of this simple if surprising result can be seen as a consequence of Christen and Fox (2005).This paper reexamines Fox and Nicholls (1997), where the idea of testing for acceptance using an approximation and before computing the exact likelihood was first suggested.In Christen and Fox (2005), the original proposal density q is used to generate a value y that is tested against an approximate target π.If accepted, y can be seen as coming from a pseudo-proposal q that simply is formalising the earlier preliminary step and is then tested against the true target π.The validation in Christen and Fox (2005) follows from standard detailed balance arguments; we will focus formally on this point in Section 2.
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).The major drawback of the scheme is that Delayed Acceptance efficiently reduces the computing cost only when the approximation π is "good enough" or "flat enough", since the probability of acceptance of a proposed value will always be smaller than in the original Metropolis-Hastings scheme.In other words, the original Metropolis-Hastings kernel dominates the new one in Peskun's (Peskun, 1973a) sense.The most relevant question raised by Christen and Fox (2005) is thus how to achieve a proper approximation; note in fact that while in Bayesian statistics a decomposition of the target is always available, by breaking original data in subsamples and considering the corresponding likelihood parts or even by just separating the prior, proposal and likelihood ratio into different factors, these decompositions may just lead to a deterioration of the algorithm properties without impacting the computational efficiency.
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 normalnormal 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 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.

The general scheme
We assume for simplicity that the target distribution π and the proposal distributions Q(x, •) all admit densities w.r.t. the Lebesgue or counting measures.We also denote by π the target density and let q(x, y) denote the proposal density.
Above, α(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 ρ k (x, y) = ρ k (y, x) −1 .The associated Delayed Acceptance Markov kernel is then defined as P (x, A) := ˆA q(x, y)α(x, y)dy + 1 − ˆX q(x, y)α(x, y)dy 1 A (x), where α(x, y) := We will denote by ( Xn ) n≥1 the Markov chain associated with P .
The order in which the sequence of functions ρ k appears in the factorisation ( 4) is important for algorithmic specification, as can be seen in Algorithm 1.It means that ρ 1 (x, Y ) is evaluated first, ρ 2 (x, Y ) second, and so on until ) which is last, with the motivation that "early rejection" can allow computational savings by avoiding the computation of the subsequent ρ k (x, Y ).

Validation
The first lemma is a standard representation leading to the validation of the Delayed Acceptance Markov chain: Lemma 1.For any Markov chain with transition kernel Π of the form Π(x, A) = ˆA q(x, y)a(x, y)dy + 1 − ˆX q(x, y)a(x, y)dy 1 A (x), and satisfying detailed balance, the function a(•) satisfies (for π-a.a.x, y) a(x, y) a(y, x) = r(x, y).
Proof.This follows immediately from the detailed balance condition π(x)q(x, y)a(x, y) = π(y)q(y, x)a(y, x).
The Delayed Acceptance Markov chain ( Xn ) n≥1 is then associated with the intended target: Lemma 2. ( Xn ) n≥1 is a π-reversible Markov chain.
Proof.From Lemma 1 it suffices to verify that α(x, y)/α(y, x) = r(x, y).Indeed, we have 2.3 Comparisons of the kernels P and P Given a probability measure µ, let us denote For a generic Markov kernel Π : E × B(E) with unique invariant probability measure µ, we define var(f, Π) := lim where (X n ) n≥1 is a Markov chain with Markov kernel Π initialised with X 1 ∼ µ.
For any f ∈ L 2 (E, µ) we define the Dirichlet form associated with a µ-reversible Markov kernel Π : E × B(E) as The (right) spectral gap of a generic µ-reversible Markov kernel has the following variational representation which leads to the following domination lemma: Andrieu et al., 2013, Lemma 34)).Let Π 1 and Π 2 be µ-reversible Markov transition kernels of µ-irreducible and aperiodic Markov chains, and assume that there exists > 0 such that for any f ∈ L 2 0 E, µ We will need the following condition in the sequel: Condition 1. Defining A := {(x, y) ∈ X 2 : r(x, y) ≥ 1}, there exists a c such that inf (x,y)∈A min k∈{1,...,d} ρ k (x, y) ≥ c.
The implication of this result is that, if P admits a right spectral gap, then so does P , whenever Condition 1 holds.Furthermore, and irrespective of whether or not P admits a right spectral gap, quantitative bounds on the asymptotic variance of MCMC estimates using ( Xn ) n≥1 in relation to those using (X n ) n≥1 are available.

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 r(x, y) = d k=1 ρk (x, y) , satisfying the balance condition, we can define a sequence of functions ρ k such that both r(x, y) = d k=1 ρ k (x, y) and Condition 1 holds.To that effect, take an arbitrary c ∈ (0, 1] and define b := c it then follows that one must define .
Proof.We note that inf (x,y)∈X
While this modification ensures that one can take = c 2 in Proposition 1, it is too general to suggest that using P can be more computationally efficient than using P when the cost of evaluating each ρ k is taken into account.Indeed, Proposition 2 holds when the functions ρk are chosen completely arbitrarily.Of course in practice, one should choose ρk and hence ρ k so that they are in some sense in agreement with r.
We will show in the next example that a certain class of ρk 's are beneficial, namely those which correspond to Metropolis-Hastings acceptance ratios with "flattened" surrogate target densities.On the other hand, it is far from difficult to come up with surrogate target densities for which unmodified use of ρk can lead to disastrous performance.

Example: unmodified surrogate targets
One common usage (Christen and Fox, 2005) of Delayed Acceptance is to substitute a surrogate target π for π in ρ 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 The first of these says that π(y)/π(x) cannot be too small when π(y) ≥ π(x).The second says that π(y)/π(x) should not be a large multiple of π(y)/π(x).There are a large variety of choices of π that allow one to take c = 1.For example, π(x) = π(x) + C for some constant C ≥ 0 and π(x) ∝ π(x) β for some β ∈ [0, 1].Note that β = 0 corresponds to π being a constant function and β = 1 corresponds to π ∝ π.In between, one can think of π as being a flattened version of π.
2.6 Counter-example: failure to reproduce geometric ergodicity Consider the case π(x) = N (x; 0, 1) and π(x) = N (x; 0, σ 2 ) with Q(x, •) a normal distribution with mean x and fixed variance for each x ∈ R.Here we have 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 σ 2 < 1, then the unmodified delayed acceptance Markov chain is not geometrically ergodic.
Proposition 3. The unmodified Delayed Acceptance Markov chain using the factorisation into ρ 1 and ρ 2 as above is not geometrically ergodic when σ < 1.
Intuitively, when x is large P (x, (−∞, x)) ≈ 1 2 but lim x→∞ P (x, {x}) = 1 because ρ 1 (x, y) takes on smaller and smaller values for y > x and ρ 2 (x, y) takes on smaller and smaller values for y < x.
Proof.From Roberts and Tweedie (1996, Theorem 5.1), it suffices to show that π−ess inf x∈X P (x, {x} ) = 0, i.e. that for any τ ∈ (0, 1) we can find A ⊆ X such that π(A) > 0 and sup x∈A P (x, {x} ) ≤ τ .Let B s (z) denote the ball of radius s around z. Given τ ∈ (0, 1), we define Now let x ∈ A, y ∈ B r (x) ∩ R + and assume y < x.It follows that ρ 2 (x, y) attains its maximum when y = x − r and therefore Similarly, let x ∈ A, y ∈ B r (x) ∩ R + and assume y > x.It follows that ρ 1 (x, y) attains its maximum when y = x + r and therefore since log(τ /3) < 0 and σ 2 < 1.Therefore, so P (x, {x} ) ≤ τ and we conclude.
The same argument can be made for much more general targets and proposals, albeit at the expense of brevity and clarity.We refrain from such a generalisation as our purpose here is to demonstrate that the DA chain may fail to inherit geometric ergodicity and that the simple proposed modification of the Delayed Acceptance kernel provided in Section 2.4 allows one to avoid this.

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 ρ i , which drives the total computational cost of the procedure.

Optimising the proposal mechanism
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 autocorrelation 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 var(P, g) = τ g × var π (g) whenever the chain is φ-irreducible, aperiodic, and reversible, var π (g) is finite and g ∈ L 2 (π).
Research on this optimisation focus on two main cases: • Consider the limit in the dimension of the target distribution toward ∞, where Roberts et al. (1997) gave conditions under which each marginal chain converges toward a Langevin diffusion.Maximising the speed of that diffusion, say h( ) where is a parameter of the scale of the proposal, implies a minimisation of the ACT and also that τ is free from the dependence on the functional, defining thus an independent measure of efficiency for the algorithm; • 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 are two successive points in the chain and • β represent the norm on the principal axes of the ellipse rescaled by the coefficients β i so that every direction contributes equally.
An interesting result in Sherlock and Roberts (2009) is that, as d → ∞, 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 (5) Eff := ESJD cost per iteration 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.
(H2) We further assume that the target distribution satisfies (A1) and (A2) in Roberts et al. (1997), which are regularity conditions on π and its first and second derivatives, and that π(x) = n i=1 f (x i ).
(H3) We consider a random walk proposal y = x + 2 /d Z, where Z ∼ N (0, I d ).Note that Gaussianity can be easily relaxed to distributions with finite fourth moment and similar results are available for more heavy-tailed distributions (Neal and Roberts, 2011).
(H4) Finally we assume that the cost of computing f (•), say c, is proportional to the cost of computing π(•), named C, with c = δC.
Normalising costs by setting C = 1, the average total cost per iteration of the Delayed Acceptance chain is δ + E [α] and the efficiency of the proposed method under the above conditions is Lemma 4.Under the above conditions (H1-H4) on the target π(x), on the proposal q(x, y) and on the factorised acceptance probability α(x, y) and that as d → ∞ where as defined in Roberts et al. (1997).
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. Proposition 4. If the conditions of Lemma 4 holds, the optimal average acceptance rate α * (δ) is independent of I.
Proof.Consider Eff (δ, ) in terms of (δ, a( )): where we dropped the dependence on in a for notation's sake.It is now evident that to maximise Eff (δ, a) in a we only need maximise The optimal scale of the proposal * (δ) and the optimal acceptance rate a * (δ) are thus given as functions of δ.In particular, as the relative cost of computing ρ 1 (x, y) with respect to ρ 2 (x, y) decreases, the proposed moves become bolder, in that * increases and a * decreases, since rejecting costs the algorithm little in terms of time, while every accepted move results in an almost independent sample.On the contrary when δ grows larger the chain rapidly approaches a Metropolis-Hastings behaviour, as it is no longer convenient to reject early.Figure 3 helps visualise the result.

Ranking the Blocks
As mentioned at the end of Section 1, the order in which the factors ρ i (x, y) are tested has a strong influence on the performance of the algorithm.Delayed Acceptance was first developed in Christen and Fox (2005), Fox and Nicholls (1997) to speed up computations using a cheap approximation π(•) of the target distribution π(•) as a first step before computing the actual, and costly, Metropolis-Hastings ratio π(y) π(x) only in the cases where the acceptance test based on the approximation π was satisfied.The main idea, namely to avoid the computation of the most costly parts as often as possible, remains relevant even for factorisations composed of more than two terms.
Consider an i.i.d.framework; the target (in x) is given by where Z = (z 1 , . . ., 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 ξ 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 . ., n and ξ n+1 (x, y) = p(y)q(y, x) p(x)q(x, y) .
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(•, •) then all the successive terms must be accepted and their order is hence of no interest.
This later setting is akin to considering the hypothetical optimal solution introduced in Section 3.1 with only two terms in the decomposition.Let a small number of the best scoring terms be merged to form ρ 1 and let the remaining factors become ρ 2 .ρ 1 (x, y) = π(y)/π(x) is then highly correlated to r(x, y), ρ 2 (x, y) ∼ 1 for every (x, y) and hence π(x) is a close approximation of the target, albeit probably flattened, which is exactly what we want (see Section 2).
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 with m < n, where Z * = (z * 1 , . . ., z * m ) is a subsample of Z.At each iteration t of the Markov chain we compute all the ξ i (x (t−1) , y) and Z * (t) is chosen as the subset that maximise the observed correlation between the values of πZ * (x (j) ) πZ * (x (j) ) (j = 1, . . ., t − 1) and the full Metropolis-Hastings ratio (or whatever other selected criterion).As computing the real arg max Z * ⊂Z is expensive, in our practical implementation we resort to a forward selection scheme; starting with the factor ξ i with the maximal correlation we build Z * (t) merging one term at a time until a desired correlation level is achieved, the observed correlation after including another term does not grow more than a small > 0 or the size of Z * (t) has reached a critical point for computational purposes (e.g.10% of the whole sample Z).
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.

Delayed Acceptance and Prefetching
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.
Clearly the usefulness of this technique depends on our ability to guess the path of the chain correctly and hence many advanced prefetching strategies make use of the observed acceptance rate of the chain or even of a fast approximation π of the target distribution to select the most likely future outcomes.
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.
As mentioned above and demonstrated in Angelino et al. (2014), Strid (2010) if a cheap approximation π of the target density is available, it can be used to select more likely future paths of the chain and this results in an efficient prefetching algorithm.
In our case the master process sequentially samples from the (Delayed Acceptance) chain by checking only the (assumed) inexpensive first approximation ρ 1 (x, y) = π(y)/π(x) while the other additional processors provide him the more expensive ρ 2 (x, y) = π(y)π(x)/π(x)π(y) computed beforehand thanks to prefetching.The theoretical properties of the chain are unchanged while the achievable speed-up may be substantial, especially for the first few additional processors.

Alternative procedure for Delayed Acceptance
In the case that every factor ρ 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 π(y) q(y, x) π(x)q(x, y) = 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.
Remark 3. Interestingly if d l = 1 ∀l (l being the number of subsets considered) this procedure reduces to Delayed Acceptance, and for l that increases and d l > 1 ∀l this combined technique will have a even lower overall acceptance rate than standard Delayed Acceptance.
to compare with Delayed Acceptance which conversely 1. simulate θ ∼ q(θ |θ); 2. simulate u 1 , . . ., u n ∼ U(0, 1) and set The differences between both schemes are thus that (a) slice sampling always accepts a move, (b) slice sampling requires the simulation of θ 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 θ 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.

Logistic Regression
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 subsampling 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 ∼ N (x, Σ) with Σ initialised to be 0.2 × 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 α * = 0, 234 (Roberts et al., 1997).
Delayed Acceptance was optimised first against the ordering of the factors as explained in Section 3; we split the data into subsamples of 10 elements and computed their empirical correlation with the full Metropolis-Hastings ratio as a criterion.Once these estimates were stable we merged into the surrogate target f the smallest number of subsamples needed to achieve a ≥ 0, 85 correlation with r(x, y).As soon as the ordering was fixed we computed δ, the relative cost of the obtained ρ 1 with respect to the full ratio, and run the chain for the remaining iterations optimising Σ against the optimal acceptance rate found through (5).
We also added the modification explained in Section 2.4 with c set such that b was slightly lower than the optimal acceptance rate above.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 ρ 1 and ρ 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.
5.2 G-MALA with Delayed Acceptance 5.2.1 MALA and Geometric MALA: 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 π(θ|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 ε is the step-size for the Euler's integration and Z ∼ N (0, 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 ε is fixed across dimensions.Roberts and Stramer (2002) propose to alleviate the issue using a pre-conditioning matrix A so that the proposal becomes θ = θ (i−1) + ε 2 A T A∇ θ log(π(θ (i−1) |y))/2 + εAZ.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)).

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 ε 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.
Given that optimal scaling for MALA in terms of the dimension d of the target differs from the random-walk setting (see Roberts and Rosenthal, 2001), we set the variance of the random-walk normal component as σ 2 d = 2 d 1/3 .Borrowing from Section 3.1, we can obtain the optimal acceptance rate for the DA-MALA, through Equation ( 5), by maximising In the above the computational cost per iteration is taken to be c = δC for the posterior ratio, C = 1 for the proposal ratio (and hence c + E [α] (C − c) for the whole kernel), h( ) is again the speed of the limiting diffusion process and K is a measure of roughness of the target distribution, depending on its derivatives.Since the optimal a * is independent from K, we do not define it more rigorously, referring to Roberts and Rosenthal (2001).Figure 4 shows that a * decreases with δ, as is the case with random-walk Metropolis-Hastings.It reaches the known optimum for the standard MALA when δ = 1.

Simulation study:
To test the above assumptions we ran a toy MALA example where we drew 100 samples from a N d (θ, I), with d = 10; π(θ) was set to be 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 δ at each run by averaging a few computed derivatives, required by the proposal ratio.We then adapt ε 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.

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.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: Indeed, these integrals cannot be computed analytically and thus their derivation involve numerical or Monte Carlo integration.This setting is such that the prior ratio-as opposed to the more common case of the likelihood ratio-is the costly part of the target evaluated in the Metropolis-Hastings acceptance ratio.Moreover, since the Jeffreys prior involves a determinant, there is no easy way to split the computation in more parts than "prior × likelihood".Hence, the Delayed Acceptance algorithm can be applied by simply splitting between the prior p J (ψ) and the likelihood (ψ|x) ratios, the later being computed first.Moreover, since the proposed prior is "non informative", its influence on the definition of the posterior distribution should be small with respect to the likelihood function and, then, computing the likelihood ratio first should not have a substantial impact on the acceptance rate.However, the improper nature of the prior means using a second acceptance ratio solely based on the prior can create trapping states in practice, even though the method remains theoretically valid.We therefore opted for stabilising this second step by saving a small fraction of the likelihood, corresponding to 5% of the sample, to regularise this second acceptance ratio.This choice translates into Algorithm 2.
Algorithm 2 Metropolis-Hastings with Delayed Acceptance for Mixture Models  3 Comparison using different performance indicators in the example of mixture estimation, based on 100 replicas of the experiments according to model (9) with a sample size n = 500, 10 5 MH simulations and 500 samples for the prior estimation.("ESS" is the effective sample size,"time" is the computational time).The actual averaged gain ( ESS DA /ESS M H time DA /time M H ) is 9.58, higher than the "double average" that the table above suggests as being around 5.
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 δ 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 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.

Figure 3 .
Figure 3. Two top panels: behaviour of * (δ) and α * (δ) as the relative cost varies.Note that for δ >> 1 the optimal values converges towards the values computed for the standard Metropolis-Hastings (dashed in red).Two bottom panels: close-up of the interesting region for 0 < δ < 1.

1 +
d 2 = d and the factors ρ i and φ j represent respectively cheap factors and costly factors.By taking now min (m=1,... ,d 1 ) cheap factors first and expensive factors last, applying the symmetry requirement to satisfy detail balance inside each of both subsets.

Figure 4 .
Figure 4. Optimal acceptance rate for the DA-MALA algorithm as a function of δ.In red, the optimal acceptance rate for MALA obtained by Roberts and Rosenthal (2001) is met for δ = 1.

Figure 5 .
Figure 5.Comparison between geometric MALA (top panels) and geometric MALA with Delayed Acceptance (bottom panels): marginal chains for two arbitrary components (left), estimated marginal posterior density for an arbitrary component (middle), 1D chain trace evaluating mixing (right).

5. 3
Mixture Model 5.3.1 Non-Informative inference on a Mixture Model: Consider a standard mixture model (MacLachlan and Peel, 2000) with a fixed number of components (7) k i=1 w i f (x|θ i ) , with k i=1 w i = 1 .

Table 1
Comparison between MH and MH with Delayed Acceptance on a logistic model.ESS is the effective sample size, ESJD the expected square jumping distance, time is the computation time.

Table 2
Comparison between standard geometric MALA and geometric MALA with Delayed Acceptance, with ESS the effective sample size, ESJD the expected square jumping distance, time the computation time and a the observed acceptance rate.