On the efficiency of pseudo-marginal random walk Metropolis algorithms

Introduction

Markov chain Monte Carlo (MCMC) algorithms have proved particularly successful in statistics for investigating posterior distributions in Bayesian analysis of complex models; see, for example, RobSmi93 , Tierney1994 , MCMChandbook . Almost all MCMC methods are based on the Metropolis–Hastings (MH) algorithm which owes much of its success to its tremendous flexibility. However, in order to use the classical MH algorithm, it must be possible to evaluate the target density up to a fixed constant of proportionality. While this is often possible, it is increasingly common for exact pointwise likelihood evaluation to be prohibitively expensive, perhaps due to the sheer size of the data set being analysed. In these situations, classical MH is rendered inapplicable.

The pseudo-marginal Metropolis–Hastings algorithm (PsMMH) beau03 , AndrieuRoberts2009 provides a general recipe for circumventing the need for target density evaluation. Instead it is required only to be able to unbiasedly estimate this density. The target densities in the numerator and denominator of the MH accept/reject ratio are then replaced by their unbiased estimates. Remarkably, this yields an algorithm which still has the target as its invariant distribution. One possible choice of algorithm, the pseudo-marginal random walk Metropolis (PsMRWM), is popular in practice (e.g., GolightlyWilkinson2011 , KnapedeValpine2012 ) because it requires no further information about the target, such as the local gradient or Hessian, which are generally more computationally expensive to approximate than the target itself Poyiadjisetal2011 .

Broadly speaking, the mixing rate of any PsMMH algorithm decreases as the dispersion in the estimation of the target density increases AndrieuRoberts2009 . In particular, if the target density happens to be substantially over-estimated, then the chain will be overly reluctant to move from that state leading to a long run of successive rejections (a sticky patch). Now, in PsMMH algorithms, the target estimate is usually computed using an average of some number, $m$, of approximations; see Sections 1.1 and 3. This leads to a trade off, with increasing $m$ leading to better mixing of the chain, but also to larger computational expense. We shall consider the problem of optimising $m$.

It is well known (e.g., RobertsRosenthal2001 , SherlockFearnheadRoberts2010 ) that the efficiency of the random-walk Metropolis (RWM) algorithm varies enormously with the scale of the proposed jumps. Small proposed jumps lead to high acceptance rates but little movement across the state space, whereas large proposed jumps lead to low acceptance rates and again to inefficient exploration of the state space. The problem of choosing the optimal scale of the RWM proposal has been tackled for various shapes of target (e.g., RobertsGelmanGilks1997 , RobertsRosenthal2001 , BedardA2007 , BedardRosenthal2008 , BeskosRobertsStuart2009 , SherlockRoberts2009 , Sherlock2013 , BreRob00 ) and has led to the following rule of thumb: choose the scale so that the acceptance rate is approximately $0.234$. Although nearly all of the theoretical results are based upon limiting arguments in high dimension, the rule of thumb appears to be applicable even in relatively low dimensions (e.g., SherlockFearnheadRoberts2010 ).

This article focusses on the efficiency of the PsMRWM as the dimension of the target density diverges to infinity. For relatively general forms of the target distribution, under the assumption of additive independent noise in the log-target, we obtain (Theorem 1) expressions for the limiting expected squared jump distance (ESJD) and asymptotic acceptance rate. ESJD is now well established as a pragmatic and useful measure of mixing for MCMC algorithms in many contexts (see, e.g., PasGel10 ), and is particularly relevant when diffusion limits can be established; see, for example, the discussion in RobRos12 . We then prove a diffusion limit for a rescaling of the first component, in the case of a target with independent and identically distributed components (Theorem 2), the efficiency of the algorithm is then given by the speed of this limiting diffusion, which is equivalent to the limiting ESJD. We examine the relationship between efficiency, scaling, and the distributional form of the noise, and consider the joint optimisation of the efficiency of the PsMRWM algorithm (taking computational time into account) with respect to $m$, and the RWM scale parameter. Exact analytical results are obtained (Corollary 1) under an assumption of Gaussian noise in the estimate of the log-target, with a variance that is inversely proportional to $m$. In this case, we prove that the optimal noise variance is 3.283, and the corresponding optimal asymptotic acceptance rate is 7.001%, thus extending the previous 23.4% result of RobertsGelmanGilks1997 . Finally, we illustrate the use of these theoretical results in a simulation study (Section 4).

The PsMMH algorithm creates a Markov chain with a stationary density (since $c=1$) of

We are thus able to substitute the estimated density for the true density, and still obtain the desired stationary distribution for $\mathbf{x}$. Note that for symmetric proposals, this simplifies to $1\wedge[\hat{\pi}_{\mathbf{v}^{*}}(\mathbf{x}^{*})/\hat{\pi}_{\mathbf{v}}(\mathbf{x})]$.

Different strategies exist for producing unbiased estimators, for instance, using importance sampling or latent variable representations, as in MR2523903 , or using particle filters DelMoral2004 , GordonSalmondSmith1993 as in AndrieuDoucetHolenstein2010 . We shall illustrate our theory in the context of Bayesian analysis of a partially observed Markov jump process.

2 Previous related literature

Other theoretical properties of pseudo-marginal algorithms are considered in AndrieuVihola2014 , which gives qualitative (geometric and polynomial ergodicity) results for the method and some results concerning the loss in efficiency caused by having to estimate the target density.

3 Notation

In this paper, we follow the standard convention whereby capital letters denote random variables, and lower case letters denote their actual values. Bold characters are used to denote vectors or matrices.

Studying the pseudo marginal random walk Metropolis in high dimensions

We focus on the case where the proposal, $\mathbf{x}^{*}$, for an update to $\mathbf{x}$ is assumed to arise from a random walk Metropolis algorithm with an isotropic Gaussian proposal

and $\mathbf{I}$ is the $d\times d$ identity matrix, and $\lambda>0$ is the scaling parameter for the proposal. The results presented in this article extend easily to a more general correlation matrix by simply considering the linear co-ordinate transformation which maps this correlation matrix to the identity matrix and examining the target in this transformed space. In proving the limiting results we consider a sequence of $d$-dimensional target probabilities $\pi^{(d)}$. In dimension $d$ the proposal is $\mathbf{X}^{(d)*}\stackrel{{\scriptstyle\mathcal{D}}}{{\sim}}\mathbf{N}(\mathbf{x}^{({d})},\lambda^{(d)2}\mathbf{I}^{(d)})$.

2 Noise in the estimate of the log-target

We will work throughout with the log-density of the target, and it will be convenient to consider the difference between the estimated log-target [$\log\hat{\pi}_{V}(\mathbf{x})$] and the true log-target [$\log\pi(\mathbf{x})$] at both the proposed values ($\mathbf{x}^{*},V^{*}$) and the current values ($\mathbf{x},V$), as well as the difference between these two differences,

Throughout this article we assume the following.

The Markov chain $(\mathbf{X},W)=\{(\mathbf{X}_{k},W_{k})\}_{k\geq 0}$ is stationary, and the distribution of the additive noise in the estimated log-target at the proposal, $W^{*}$, is independent of the proposal itself, $\mathbf{X}^{*}$.

It is unrealistic to believe that the second part of Assumption 1 should hold in practice. Pragmatically, this assumption is necessary in order to make progress with the theory presented herein; however, in our simulation study in Section 4 we provide evidence that, in the scenarios considered, the variation in the noise distribution is relatively small.

This is Lemma $1$ of Pittetal2012 . Under Assumption 1, $W$ and $\mathbf{X}$ are therefore independent, and the stationary density of $W$ is $g^{*}(w)e^{w}$.

3 High-dimensional target distribution

We describe in this section conditions on the sequence of target densities $\pi^{(d)}$ that ensure that the quantity $\log[\pi^{(d)}(\mathbf{X}^{*})/\pi^{(d)}(\mathbf{X})]$ behaves asymptotically as a Gaussian distribution under an appropriate choice of jump scaling $\lambda^{(d)}$. The main assumption is that there exist sequences of scalings $s_{g}^{(d)}>0$ and $s_{L}^{(d)}>0$ for the gradient and the Laplacian of the log-likelihood $\log\pi^{(d)}$ such that the following two limits hold in probability:

for $\mathbf{X}^{({d})}\stackrel{{\scriptstyle\mathcal{D}}}{{\sim}}\pi^{(d)}$. In the rest of this article we assume that the sequence of densities $\pi^{(d)}$ is such that for each index $i\geq 1$, with all components of $\mathbf{x}$ fixed except the $i$th, the $i$th component satisfies

Under this regularity condition, an integration by parts shows that

Equation (5) thus yields $\lim_{d\to\infty}(s_{g}^{(d)})^{2}/s_{L}^{(d)}=1$. We will suppose from now on, without loss of generality, that $s_{g}^{(d)}=\sqrt{s_{L}^{(d)}}=:s^{(d)}$. We also require that no single component of the local Hessian $H^{(d)}(\mathbf{x}):=[\partial^{2}_{ij}\log\pi^{(d)}(\mathbf{x})]_{0\leq i,j\leq d}$ dominate the others in the sense that the limit

holds in probability. We also assume that the Hessian matrix is sufficiently regular so that for any $\sigma^{2},\varepsilon>0$ and $\mathbf{Z}^{(d)}\stackrel{{\scriptstyle\mathcal{D}}}{{\sim}}\mathbf{N}(0,\mathbf{I}^{(d)})$

These conditions are discussed in detail in Sherlock2013 where they are shown to hold, for example, when the target is the joint distribution of successive elements of a class of finite-order multivariate Markov processes. The targets considered in RobertsGelmanGilks1997 , RobertsRosenthal2001 and Section 2.5 all satisfy the conditions with $s^{(d)}\propto d^{1/2}$. We record the conditions formally as:

The sequence of densities $\pi^{(d)}$ satisfies equations (5), (7), (LABEL:eqn.hessian.regularity), and the regularity condition (6).

We shall show in next section that under these assumptions the choice of jump size

4 Expected squared jump distance

Clearly, (10) is satisfied when $\mathbf{T}^{({d})}=I_{d}$ (i.e., Euclidian ESJD).

with $B$ as in (3), where $\Phi$ is the cumulative distribution of a standard Gaussian distribution.

Expected squared jump distance. A rescaled expected squared jump distance converges as $d\to\infty$ to a related limit,

For these two examples, as might be expected, for any given scaling of the random walk proposal, the efficiency relative to the idealised algorithm decreases as the standard deviation of the noise increases, a phenomenon that is investigated more generally in AndrieuVihola2014 . Thus there is an implicit cost of having to estimate the target density. As a result of this, we should not expect the optimal acceptance probability for RWM of $0.234$ to hold here.

5 Diffusion limit

We next prove that PsMRWM in high dimensions can be well-approximated by an appropriate diffusion limit (obtained as $d\to\infty$). This provides further justification for measuring efficiency by the ESJD, as discussed in detail in RobRos12 . Briefly, the limiting ESJD (suitably scaled) is equal to the square of the limiting process’s diffusion coefficient, $h$ say. By a simple time change argument, the asymptotic variance of any Monte Carlo estimate of interest is inversely proportional to $h$. Minimising variance is thus equivalent to maximising $h$; that is, $h$ becomes (at least in the limit) unambiguously the right quantity to optimise. By constrast, MCMC algorithms which have nondiffusion limits can behave in very different ways, and ESJD may not be an appropriate way to compare algorithms in such cases.

We shall consider in this section the PsMRWM algorithm applied to a sequence of simple i.i.d. target densities

where $f$ is a one-dimensional probability density. We assume throughout this section that the following regularity assumptions hold.

The first four moments of the distribution with density $f$ are finite. The log-likelihood mapping $x\mapsto\log f(x)$ is smooth with second, third, and fourth derivatives globally bounded.

In the remainder of this article we consider the sequences of scaling functions $\sqrt{s_{L}^{(d)}}=s_{g}^{(d)}:=\sqrt{I\times d}$, with

Let $T>0$ be a finite time horizon. For all $d\geq 1$ let each Markov chain and the additive noise satisfy Assumption 1, let the sequence of product form densities $\pi^{(d)}$ satisfy the regularity Assumption 3 and set the scale of the jump proposals as in equation (15). Then, as $d\to\infty$,

in the Skorokhod topology on $D([0,T])$, where $V$ satisfies the Langevin SDE

with initial distribution $V_{0}\stackrel{{\scriptstyle\mathcal{D}}}{{\sim}}f$ and $B_{t}$ a standard Brownian motion. The speed function $h$ is proportional to the asymptotic rescaled ESJD function $J$,

with the constant of proportionality $I$ defined by equation (14).

Optimising the PsMRWM

We next consider the question of optimising the PsMRWM. Now, when examining the efficiency of a standard RWM, the expected computation (CPU) time is usually not taken into account since it is implicitly assumed to be independent of the choice of tuning parameter(s). This may indeed be approximately true for the RWM. However, for the PsMRWM the expected CPU time for a single iteration of the algorithm is usually approximately inversely proportional to the variance of the estimator $\hat{\pi}(x)$. For this reason, we measure the efficiency of the PsMRWM through a rescaled version of the ESJD,

Of course, for any increasing function $F$, the quantity $F\mbox{(ESJD)}/\mbox{(Expected}\break\mbox{one-step computing time)}$ is a possible measure of efficiency. However, the discussion at the start of Section 2.5 indicates that (17) is the appropriate measure of efficiency in the high-dimensional asymptotic regime considered in this article.

In the remainder of this section, we implicitly assume that the target distributions satisfy Assumption 2.

We shall restrict attention to the case in which the additive noise follows a Gaussian distribution. More precisely, we shall assume the following, which we shall refer to for brevity as “the standard asymptotic regime” (SAR):

For each $x\in\mathcal{X}$ and $\sigma^{2}>0$, we have an unbiased estimator $\hat{\pi}(x)$ of $\pi(x)$, such that $\log\hat{\pi}(x)$ follows a Gaussian distribution with variance $\sigma^{2}$. Furthermore, the expected one-step computing time is inversely proportional to $\sigma^{2}$.

Intuitively, Assumption 4 are designed to model the situation where $\pi(x)$ is estimated as a product of $n$ averages of $m$ i.i.d. samples in the limit as $n\rightarrow\infty$ and with $m\propto n$. For a fixed large $n$, approximate normality follows from the central limit theorem; moreover $\sigma^{2}\approx c/m$ for some $c>0$, and the computational time is proportional to $m$ and hence to $1/\sigma^{2}$. Assumption 4 have recently been shown to hold more generally, in the context of particle filtering for a hidden Markov model; see berard2013lognormal . There are other natural situations where multiplicative forms for the importance sampling estimator of the likelihood might make the estimator well-approximated as a log-Gaussian, for example, in correcting for a PAC likelihood approximation; see LiSte03 .

Under the SAR of Assumption 4, we will prove an optimality result in Section 3.2 which specifies a particular optimal variance for the estimate of the log-target.

2 Optimisation under the standard asymptotic regime

In this section we consider a sequence $\pi^{(d)}$ of target distributions satisfying Assumption 2 and assume that each unbiased estimator satisfies the independence in Assumption 1. Under these assumptions, the rescaled ESJD of the PsMRWM algorithm with jump size (9) is described by Theorem 1. Under the SAR, that is, Assumption 4, and with $\operatorname{Var}[\log\hat{\pi}(x)]=\sigma^{2}$, the noise difference is $B\stackrel{{\scriptstyle\mathcal{D}}}{{\sim}}\mathbf{N}(-\sigma^{2},2\sigma^{2})$. Since the mean one-step computing time is assumed to be inversely proportional to the variance, $\sigma^{2}$, the asymptotic efficiency, as $d\to\infty$, is proportional to

The point at which the maximal efficiency is achieved is detailed precisely in Corollary 1 below.

at which point the corresponding asymptotic acceptance rate is

For convenience, write $\tau^{2}:=2\sigma^{2}$, and introduce three independent standard Gaussian random variables $U,V,Z\stackrel{{\scriptstyle\mathcal{D}}}{{\sim}}\mathbf{N}(0,1)$. Notice that $B\stackrel{{\scriptstyle\mathcal{D}}}{{\sim}}-\tau^{2}/2+\tau U$ and

(3) In practice, $\sigma^{2}$ might be a function of a discrete number $m$ of samples or particles and hence only take a discrete set of values. In particular, if the variance in the noise using $m=1$ is already lower than $3.283$, then there can be little gain in increasing $m$.

Our theory applies in the limit when the dimension $d$ of the (marginal) target $\mathbf{X}$ goes to infinity. However, using a similar argument to that in SherlockRoberts2009 , when $\mathbf{X}\sim\mathbf{N}(0,\mathbf{I}_{d})$, it can be shown that under the SAR with the proposal as in (2) the ESJD and acceptance rate are

In the simulation study of Section 4 below, we find that Corollary 1 and its associated formulae provide a good description of the optimal settings for a particle filter with $T=50$ and $d=5$.

Simulation study

In practice the assumptions underlying this result may not hold: the dimension of the parameter space is finite, the distribution of the noise, $W^{*}$, may not be Gaussian, and it is likely to also vary with position, $\mathbf{x}^{*}$. We conduct a simulation study to provide an indication of both the extent of and the effect of such deviations.

We use the Particle Marginal RWM algorithm (PMRWM) of AndrieuDoucetHolenstein2010 to perform exact inference for the Lotka–Volterra predator-prey model; see GolightlyWilkinson2011 for a more detailed description of the PMRWM which focusses on this particular class of applications. Starting from an initial value, which is, for simplicity, assumed known, the two-dimensional latent variable $\mathbf{U}$ evolves according to a Markov jump process (MJP). Each component is observed at regular intervals with Gaussian error of an unknown variance. Appendix B provides details of the observation regime and of the transitions of the MJP and their associated rates. It also provides the parameter values, the priors and the lengths of the MCMC runs.

We perform three checks on our assumptions. The diagnostic runs provide a sample from the distribution of $W^{*}$, the estimate of the log-target at a proposed value; this allows us to investigate the second part of Assumption 1 and both parts of Assumption 4. We first examine the SAR Assumption 4. Figure 4 shows QQplots for $m=50$, $m=100$ and $m=400$ against a Gaussian distribution; it is clear that at $m=50$ the right-hand tail is slightly too light and the left-hand tail is much heavier than that of a Gaussian. Similar but much smaller discrepancies are present at $m=100$, whilst at $m=400$ the noise distribution is almost indistinguishable from that of a Gaussian.

The left-hand panel in Figure 5 plots $\log\operatorname{Var}[W^{*}]$ against $\log m$ and includes a line with the theoretical slope of $-1$ and passing through an additional point at $m=1600$. The heavy left-hand tail at $m=50$ leads to a considerably higher variance than that which would arise under the SAR; however, even by $m=80$ the fit is reasonably close.

We assess the degree of dependence of the distribution of $W^{*}$ on the position $\mathbf{x}$ by considering the joint distribution of $W^{*}$ and $L:=(\log\pi)(\mathbf{X})$, the true log-target evaluated at $\mathbf{X}$, where $\mathbf{X}$ is distributed according to the target. For a particular $m$, all of the runs with $\gamma>0$ provide a combined sample of size $n_{1}$ from the distribution of the estimate of the log-target at the current value, $\hat{L}=L+W$, whereas (after scaling so that $\frac{1}{n_{2}}\sum_{i=1}^{n_{2}}\exp{w^{*(i)}}=1$) each run with $\gamma=0$ provides a sample of size $n_{2}$ from the distribution of $W^{*}$ at $\mathbf{x}=\hat{\mathbf{x}}$. Equation (4) shows that subject to Assumption 1, $W$ and $L$ are independent and that the density of $W$ is an exponentially tilted version of the density of $W^{*}$. These two properties lead directly to the following.

The right-hand side of (20) is independent of the noise distribution, or equivalently of the number of particles, $m$. Moreover, if the noise is small enough then the ratio on the left-hand side should provide a good estimator of the true moment generating function (MGF) of $L$ even if there is dependence (since the impact of any dependence will be small).

In our scenario, realisations of $L$ are typically between $-385$ and $-375$ with a mode at approximately $-379$, so the MGFs of $L$ and $\hat{L}$ are dominated by the term $e^{-379t}$, whatever the noise distribution. To be able to discern any differences we therefore consider for each value of $m$, shifted estimators of the MGFs of $\hat{L}$ and of $L$

The central panel of Figure 5 shows $M_{1}(t)$ with a separate curve for each value of $m$; the lowest curve is our best estimate of the true MGF of $L$ ($M_{2}(t)$ from $m=1600$). The right-hand panel shows $M_{2}(t)$ for each value of $m$. Clearly the curves in the right-hand panel do not coincide, and so the assumption of independence does not hold precisely. However, it is clear from the very different vertical scales of the two figures that most of the difference between the distribution of $\hat{L}$ for any given $m$ and the distribution of $L$ can be explained by Assumption 1.

Proofs of results

Equation (4) yields that $B=W^{*}-W$ has density $\rho$ satisfying

This fact will be used in the proofs of Theorem 1 and Proposition 1.

For notational convenience, we drop the index $[\cdot]^{(d)}$ when the context is clear. As in Section 2.3, the Hessian matrix of the log-likelihood $L(\mathbf{x}):=\log\pi^{(d)}(\mathbf{x})$ at $\mathbf{x}\in\mathbf{R}^{d}$ is denoted by $H(\mathbf{x})=[\partial^{2}_{ij}L(\mathbf{x})]_{1\leq i,j\leq d}$.

Proof of equation (11). The mean acceptance probability equals

Equation (LABEL:eqn.hessian.regularity) shows that the remainder $R(\mathbf{X},\mathbf{Z},\lambda^{(d)})$ converges to zero in probability.

Proof of equation (12). The proof of equation (12) follows from equation (11). Note that we have

2 Proof of Proposition 1

This quantity is clearly strictly negative, completing the proof of Proposition 1.

3 Proof of Proposition 2

4 Proof of Theorem 2

The proof follows ideas from BedardA2007 , which itself is an adaptation of the original paper RobertsGelmanGilks1997 . It is based on ethier1986markov , Theorem $8.2$, Chapter $4$, which gives conditions under which the finite dimensional distributions of a sequence of processes converge weakly to those of some Markov process. ethier1986markov , Corollary $8.6$, Chapter $8$, provides further conditions for this sequence of processes to be relatively compact in the appropriate topology and thus establish weak convergence of the stochastic processes themselves.

The situation is slightly more involved than the one presented in RobertsGelmanGilks1997 , BedardA2007 ; the proof needs a homogenisation argument since the processes $\mathbf{X}^{(d)}$ and $W^{(d)}$ evolve on two different time scales. Indeed, it will become apparent from the proof that the process $\mathbf{X}^{(d)}$ takes $\mathcal{O}(d)$ steps to mix while the process $W^{(d)}$ takes $\mathcal{O}(1)$ steps to mix. In order to exploit this time-scales separation, we introduce an intermediary time scale $T_{d}=\lfloor{d^{\gamma}}\rfloor$ where $0<\gamma<1/4$ is an exponent whose exact value is not important to the proof. The intuition is that after $\mathcal{O}(T_{d})$ steps the process $W^{(d)}$ has mixed while each coordinate of $\mathbf{X}^{(d)}$ has only moved by an infinitesimal quantity. We introduce the subsampled processes $\widetilde{\mathbf{X}}^{(d)}$ and $\widetilde{W}^{(d)}$ defined by

One step of the process $\widetilde{\mathbf{X}}^{(d)}$ (resp., $\widetilde{W}^{(d)}$) corresponds to $T_{d}$ steps of the process $\mathbf{X}^{(d)}$ (resp., $W^{(d)}$). We then define an accelerated version $\widetilde{V}^{(d)}$ of the subsampled first coordinate process $k\mapsto\widetilde{X}^{(d)}_{k,1}$. In order to prove a diffusion limit for the first coordinate of the process $\mathbf{X}^{(d)}$, one needs to accelerate time by a factor of $d$; consequently, in order to prove a diffusion limit for the process $\widetilde{\mathbf{X}}^{(d)}$, one needs to accelerate time by a factor $d/T_{d}$, and thus define $\widetilde{V}^{(d)}$ by

The proof then consists of showing that the sequence $\widetilde{V}^{(d)}$ converges weakly in the Skorohod topology towards the limiting diffusion (16) and verifying that $\|\widetilde{V}^{(d)}-V^{(d)}\|_{\infty,[0,T]}$ converges to zero in probability; this is enough to prove that the sequence $V^{(d)}$ converges weakly in the Skorohod topology towards the limiting diffusion (16). The proof is divided into three main steps. First, we show that the finite dimensional marginals of the process $\widetilde{V}^{(d)}$ converge to those of the limiting diffusion (16). Second, we establish that the sequence $\widetilde{V}^{(d)}$ is weakly relatively compact. These two steps prove that the sequence $\widetilde{V}^{(d)}$ converges weakly in the Skorohod topology towards the diffusion (16). As a final step, we prove that the quantity $\|\widetilde{V}^{(d)}-V^{(d)}\|_{\infty,[0,T]}$ converges to zero in probability, establishing the weak convergence of the sequence $V^{(d)}$ towards the diffusion (16). Before embarking on the proof we define several quantities that will be needed in the sequel. We denote by $\mathcal{L}$ the generator of the limiting diffusion $(\ref{e.limiting.diffusion})$. Similarly, we define $\mathcal{L}^{(d)}$ and $\widetilde{\mathcal{L}}^{(d)}$ the approximate generators of the first coordinate process $\{X^{(d)}_{k,1}\}_{k\geq 0}$ and its accelerated version $\{\widetilde{X}^{(d)}_{k,1}\}_{k\geq 0}$; for any smooth and compactly supported test function $\varphi\dvtx\mathbf{R}\to\mathbf{R}$, vector $\mathbf{x}=(x_{1},\ldots,x_{d})\in\mathbf{R}^{d}$ and scalar $w\in\mathbf{R}$, we have

with $\delta=1/d$. Note that although $\varphi$ is a scalar function, the functions $\mathcal{L}^{(d)}\varphi$ and $\widetilde{\mathcal{L}}^{(d)}\varphi$ are defined on $\mathbf{R}^{d}\times\mathbf{R}$. In the sequel we sometimes write $\widetilde{\mathcal{L}}^{(d)}\varphi(x_{1},\ldots,x_{d},w)$ instead of $\widetilde{\mathcal{L}}^{(d)}\varphi(\mathbf{x},w)$.

In this section we prove that the finite dimensional distributions of the sequence of processes $\widetilde{V}^{(d)}$ converge weakly to those of the diffusion (16). Since the limiting process is a scalar diffusion, the set of smooth and compactly supported functions is a core for the generator of the limiting diffusion (ethier1986markov , Theorem $2.1$, Chapter $8$); in the sequel, one can thus work with test functions belonging to this core only. To prove the convergence of the finite dimensional marginals, one can apply ethier1986markov , Chapter $4$, Theorem $8.2$, Corollary $8.4$, to the pair $(\xi^{(d)},\varphi^{(d)})$ defined by

To establish that this result applies, we will concentrate on proving that for any smooth and compactly supported function $\varphi\dvtx\mathbf{R}\to\mathbf{R}$ the following limit holds:

Equation (24) follows from the telescoping expansion $\varphi(\mathbf{X}^{(d)}_{T_{d}})-\varphi(\mathbf{X}^{(d)}_{0})=\sum_{k=0}^{T_{d}-1}\varphi(\mathbf{X}^{(d)}_{k+1})-\varphi(\mathbf{X}^{(d)}_{k})$ and the law of iterated conditional expectations. The following lemma is crucial:

Let Assumptions 1 and 3 be satisfied. There exist two bounded and continuous functions $a,b\dvtx\mathbf{R}\to\mathbf{R}$ satisfying the following properties: {longlist}[(2)]

For any smooth and compactly supported function $\varphi\dvtx\mathbf{R}\to\mathbf{R}$ the averaged generator $\mathcal{G}\varphi$ defined for any $(x_{1},w)\in\mathbf{R}^{2}$ by

for an i.i.d. sequence $\{X_{k}\}_{k\geq 1}$ marginally distributed as $f(x)\,dx$ and constant $I$ defined by (14).

Let $h\dvtx\mathbf{R}\to\mathbf{R}$ be a bounded measurable function. Suppose that for any $d\geq 1$ the Markov chain $\{(\mathbf{X}^{(d)}_{k},W^{(d)}_{k})\}_{k\geq 0}$ is started at stationarity. The following limit holds:

with $W$ distributed according to the stationary distribution $W\stackrel{{\scriptstyle\mathcal{D}}}{{\sim}}e^{w}g^{*}(w)\,dw$.

The above lemma thus shows that $T_{d}=\lfloor{d^{\gamma}}\rfloor$ steps, with $0<\gamma<1/4$, are enough for the process $W^{(d)}$ to mix. The proof relies on a coupling argument and the ergodic theorem for Markov chains. Details can be found Appendix A.2. We now have all the tools in hands to prove equation (23). First, with the notation $\mathbf{X}^{(d)}=(X_{1},\ldots,X_{d})$, the telescoping expansion (24) and Jensen’s conditional inequality yields

One can then use the triangle inequality to obtain the bound

To complete the proof of the convergence of the finite dimensional distributions of $\widetilde{V}^{(d)}$ towards those of the limiting diffusion (16), it remains to prove that $E_{i}(d)\to 0$ as $d\to\infty$ for $i=1,2,3$:

The formula for the quantity $\mathcal{G}\varphi(x,w)$ shows that the expectation $E_{2}(d)$ also reads

Because the function $\varphi$ is smooth with compact support, it follows(Cauchy–Schwarz) that this quantity is less than a constant multiple of

The process $\widetilde{V}^{(d)}$ is started at stationarity and the space of smooth functions with compact support is an algebra that strongly separates points. Ethier and Kurtz (ethier1986markov , Chapter 4, Corollary 8.6) show that in order to prove that the sequence $\widetilde{V}^{(d)}$ is relatively weakly compact in the Skorohod topology it suffices to verify that equations (8.33) and (8.34) of ethier1986markov , Chapter 4, hold.

To prove (8.33) one needs to show that the expectation of $\sup\{|\xi_{d}(t)-\widetilde{V}^{(d)}(t)|\dvtx t\in[0,T]\}$ converges to zero as $d\to\infty$, where the process $\xi_{d}$ is defined in equation (LABEL:eq.xi.phi). Note that the supremum is less than

Let $\{\xi_{k}\}_{k\geq 1}$ an i.i.d. sequence of standard Gaussian random variables $\mathbf{N}(0,1)$. We have

Indeed, it suffices to prove that $\delta$ times the expectation of the supremum $\sup\{S(i,d)\dvtx i\leq d/T_{d}\}$, with $S(i,d)=\sum_{k=iT_{d}}^{(i+1)T_{d}}|\xi_{k}|$, converges to zero; this follows from Markov’s inequality and standard Gaussian computations.

This completes the proof of the relative weak compactness in the Skorohod topology. The sequence of processes $\widetilde{V}^{(d)}$ is weakly compact in the Skorohod topology, and the finite dimensional distributions of $\widetilde{V}^{(d)}$ converge to the finite dimensional distribution of the diffusion (16). Consequently, the sequence of processes $\widetilde{V}^{(d)}$ converges weakly in the Skorohod space $D([0,T])$ to the diffusion (16). The next section shows that the discrepancy between $V^{(d)}$ and $\widetilde{V}^{(d)}$ is small and thus proves that the sequence of processes $V^{(d)}$ also converges to the diffusion (16).

Since $\sup_{t\leq T}|V^{(d)}_{t}-\widetilde{V}^{(d)}_{t}|$ is less than the supremum of equation (27), Lemma 3 yields that $\|\widetilde{V}-\break V^{(d)}\|_{\infty,[0,T]}$ converges to zero in probability. This ends the proof of Theorem 2.

Discussion

We have examined the behaviour of the pseudo-marginal random walk Metropolis algorithm in the limit as the dimension of the target approaches infinity, under the assumption that the noise in the estimate of the log-target at a proposed new value, $\mathbf{x}$, is additive and independent of $\mathbf{x}$.

Subject to relatively general conditions on the target, limiting forms for the acceptance rate and for the efficiency, in terms of expected squared jump distance (ESJD), have been obtained. We examined two different noise distributions (Gaussian and Laplace), and found that the optimal scaling of the proposal is insensitive to the variance of the noise and to whether the noise has a Gaussian or a Laplace distribution.

We then examined the behaviour of the Markov chain on the target, $\mathbf{x}$, and the noise, obtaining a limiting diffusion for the first component of a target with independent and identically distributed components. The efficiency function in this case is proportional to the speed of the diffusion, thus further justifying the use of ESJD in this context.

We identified a “standard asymptotic regime” under which the additive noise is Gaussian with variance inversely proportional to the number of unbiased estimates that are used. In this regime the efficiency function is especially tractable, and we showed that it is maximised when the acceptance rate is approximately 7.0% and the variance of the Gaussian noise is approximately 3.3. We noted that in this regime the optimal noise variance is also insensitive to the choice of scaling.

A detailed simulation study on a Lotka–Volterra Markov jump process using a particle filter suggested that in the scenario considered the assumptions of the standard asymptotic regime are reasonable provided the number of particles is not too low. Furthermore, whilst the assumption that the distribution of the noise does not depend on the current position is not true, variations in the distribution have a small effect on the distribution of the estimates of the log-target compared with the effect of the noise itself. The optimal scaling was found to be insensitive to the noise variance (or equivalently the number of particles), and the optimal noise variance was relatively insensitive to the choice of scaling. The overall optimal scaling was consistent with the theoretical value obtained; however the optimal variance was a little lower than the theoretically optimal value. Investigations showed that this discrepancy can be explained by the differences between our theoretical measure of efficiency (ESJD) and empirical measures used in the simulation study (ESS).

The results from the simulation study suggest that in low dimension a safer option than tuning to a particular variance and acceptance rate might be to take advantage of the insensitivity of the optimal scaling to the variance and vice versa and optimise scaling and variance independently.

The diffusion limit provides strong support for the optimisation strategies suggested by the ESJD criterion. However, in an ideal world it would be good to show that the sequence of algorithms which achieves the minimal optimal integrated autocorrelation time for a given functional might converge to the optimal diffusion. This is a generic question which is relevant to all diffusion limits for MCMC algorithms, and there are still important open questions regarding the relationships between ESJD, diffusion limits, and limiting optimal integrated autocorrelation. In this direction, a recent paper RRcomplexity has shown that diffusion limits can be translated into complexity results, thus demonstrating that at least the order of magnitude of the number of iterations to “converge” can be read off from the diffusion limit.

Appendix A Proof of technical lemmas

and $(X^{(d)}_{k+1},W^{(d)}_{k+1})=(X^{(d)}_{k},W^{(d)}_{k})$ otherwise. We define $\mathbf{X}^{(d)}=(X^{(d)}_{k,1},\ldots,X^{(d)}_{k,d})$. For any dimension $d\geq 1$ the process $\{\mathbf{X}^{(d)}_{k},W^{(d)}_{k})\}_{k\geq 1}$ is a Metropolis–Hastings Markov chain started at stationarity, that is, $(\mathbf{X}^{(d)}_{0},W^{(d)}_{0})=(X_{1},\ldots,\break X_{d},W)\stackrel{{\scriptstyle\mathcal{D}}}{{\sim}}\pi^{(d)}$, targeting the distribution $\pi^{(d)}$.

In this section, for notational convenience, we write $Z_{j}$ instead of $Z_{0,j}$ and $W^{*}$ instead of $W^{*}_{0}$. We set

converges in $L^{2}$ to $a(W)A^{\prime}(X_{1})$. By the Cauchy–Schwarz inequality, this reduces to proving that

By the Portmanteau’s theorem, the dominated convergence theorem, and the definition of $\Omega^{(d)}$, this reduces to proving that for almost every realisation $\{x_{i}\}_{i\geq 1}$ of the i.i.d. sequence $\{X_{i}\}_{i\geq 1}$ the following limit holds in distribution:

A.2 Proof of Lemma 2

For convenience, we first give a high-level description of the reasoning. We construct processes $\{\widehat{W}^{(d)}_{k}\}_{k\geq 0}$, $\{\widehat{Y}^{(d)}_{k}\}_{k\geq 0}$, and $\{Y_{k}\}_{k\geq}$ satisfying the following:

With high probability $\widehat{W}^{(d)}_{k}=W^{(d)}_{k}$ for all $k\leq T_{d}$.

The process $\{\widehat{Y}^{(d)}_{k}\}_{k\geq 0}$ has the same law as the process $\{\widehat{W}^{(d)}_{k}\}_{k\geq 0}$.

With high probability $\widehat{Y}^{(d)}_{k}=Y_{k}$ for all $k\leq T_{d}$.

The process $\{Y_{k}\}_{k\geq 0}$ is a Markov chain that is ergodic with invariant distribution $e^{w}g^{*}(w)\,dw$.

One can thus use an approximation of the type

and the usual ergodic theorem gives the conclusion. We use at several places the following elementary lemma.

Let $T_{d}=\lfloor{d^{\gamma}}\rfloor$ with $0<\gamma<\frac{1}{4}$. Let $\{P^{(d)}_{k}\}_{k,d\geq 0}$ and $\{Q^{(d)}_{k}\}_{k,d\geq 0}$ be two arrays of $(0,1)$-valued random variables. Let $\{U_{k}\}_{k\geq 0}$ be a sequence of random variables uniformly distributed on the interval $(0,1)$. We suppose that for all dimension $d\geq 1$ the random variable $U_{k}$ is independent from $\{P^{(d)}_{j}\}_{j=0}^{k-1}$ and $\{Q^{(d)}{j}\}_{j=0}^{k-1}$. Consider the event

and $\widehat{W}^{(d)}_{k+1}=\widehat{W}^{(d)}_{k}$ otherwise. We define $\widehat{Y}^{(d)}_{k+1}=W^{*}_{k}$ if

and $\widehat{Y}^{(d)}_{k+1}=\widehat{Y}^{(d)}_{k}$ otherwise. We define $Y_{k+1}=W^{*}_{k}$ if

Under Assumption 3 the second and third derivatives of $A$ are bounded so that bound (30) follows from a second-order Taylor expansion,

$\widehat{W}^{(d)}$ and $\widehat{Y}^{(d)}$ have same law. It is straightforward to verify that the processes $\{\widehat{W}^{(d)}_{k}\}_{k\geq 0}$ and $\{\widehat{Y}^{(d)}_{k}\}_{k\geq 0}$ have the same law.

We now show that the Markov chain $\{Y_{k}\}_{k\geq 0}$ is a Markov chain that is reversible with respect to the distribution $e^{w}g^{*}(w)\,dw$,

which is indeed symmetric. Note that this Markov chain corresponds to the penalty method of ceperley1999penalty ; see also nicholls2012coupled for a discussion of this algorithm. The ergodic theorem for Markov chains applies; for any bounded and measurable function $h\dvtx\mathbf{R}\to\mathbf{R}$ we have

One can thus use the triangle inequality several times to see that for any bounded and measurable function $h\dvtx\mathbf{R}\to\mathbf{R}$, we have

Appendix B Details of the Lotka Volterra model

the rate for any other transition is zero. Observations of the Markov chain, when they occur, are subject to Gaussian error,

Using $\mathbf{x}=(0.006,0.6,0.3,25,49)$, a realisation of the stochastic process was simulated from initial value $\mathbf{u}(0)=(70,70)$ for $T=50$ time units. The state, perturbed with Gaussian noise, $\mathbf{y}(t)$, was recorded at $t=1,2,\dots,T$. For inference, $X_{1},\ldots,X_{5}$ were assumed to be independent, a priori with $\log X_{i}\sim\operatorname{Unif}$, ($i=1,\ldots,5$).

The initial value for each chain was a vector of estimates of the posterior median for each parameter, obtained from the initial run; hence no “burn-in” was required. Each algorithm was run for $2.5\times 10^{5}$ iterations, except with $m=50$ and $m=80$, where $10^{6}$ iterations were used. Output was thinned by a factor of $10$ for storage.

Acknowledgements

We are grateful to the Associate Editor and three referees for their comments, which helped improve both the presentation and the content of this article. Gareth Roberts and Jeffrey Rosenthal are grateful for financial support in carrying out this research from, respectively, EPSRC of the UK, through the CRiSM (EP/D002060/1) and iLike (EP/K014463/1) projects, and NSERC of Canada.

References