MSC 2020 subject classifications: Primary 60J22; secondary 65C05.

Keywords: Markov chain Monte Carlo; Metropolis–Hastings; optimal scaling; Hamiltonian dynamics; Langevin algorithms.

1 Introduction↩︎

This paper presents a unified approach to the investigation of the scaling properties of different MCMC sampling algorithms. Hundreds or perhaps thousands of MCMC algorithms based on Metropolisation have been suggested to sample from probability densities \[\label{eq:target} \pi(q) \propto \exp(-V(q)), \qquad q\in\mathbb{R}^m.\tag{1}\] Two of the best known proposals are given by the Random Walk Metropolis (RWM) formula \[\label{eq:RWMproposal} q^\star = q+\sqrt{\delta} p,\qquad p\sim N(0,I_m),\tag{2}\] and the MALA formula [1] \[\label{eq:malaproposal} q^\star = q - \frac{\delta}{2} \nabla V(q)+ \sqrt{\delta} p, \qquad p \sim N(0,I_m).\tag{3}\] The latter provides a consistent discretisation of the Langevin equation \[\label{eq:langevin} dq(t) = -\frac{1}{2}\nabla V(q(t))\,dt+\, dW(t),\tag{4}\] which preserves the target 1 .

A criterion to choose between all the different algorithms is to study how their parameters have to be varied as the dimensionality of the target increases. For targets consisting of \(d\) independent copies of a univariate distribution, it was proved in the pioneering contribution [2] that the variance of the RWM proposal should be scaled as \(\delta_d= \ell^2/d\), with a constant \(\ell\), to ensure that, as \(d\rightarrow \infty\), the acceptance probability converges to a nontrivial limit not equal to zero or one. For MALA [3], the scaling is less demanding: \(\delta_d= \ell^2/d^{1/3}\). The best result available for a MALA-like sampler corresponds to fMALA [4] where \(\delta_{d}=\ell^2/d^{1/5}\), but this comes at the cost of calculating higher order derivatives of the potential \(V\). These studies also identify optimal values of the acceptance probability, so that the algorithms are most efficient if \(\ell\) is tuned to achieve those optimal values. It is remarkable that, at least in the product of identical copies scenario, the optimal values of the acceptance probability may be proved to be independent of the target. The literature on optimal scaling is by now substantial; a selection of useful references include [4]–[15] among others.

In this paper we provide a simple, yet general approach to deriving scaling results that, on the one hand, makes it possible to recover the results in [2]–[4], [6] (with fewer hypotheses) and, on the other, provides, in an easy way, optimal scaling results for many other proposals. These include implicit proposals that need the solution of systems of nonlinear equations and proposals based on a variety of integrators taken from the Hamiltonian literature. Our approach allows to easily construct MALA-like gradient-based proposals with scaling \(\delta_d= \ell^2/d^\mu\) with arbitrarily small \(\mu>0\) in the i.i.d.product scenario. The analysis applies to products of a given, not necessarily univariate, distribution, and also to cases where the different components are scaled differently. In the latter case, we investigate the effect of preconditioning the algorithms.

The initial motivation for our research came when trying to study the scaling for the following implicit counterpart of 3 \[\label{eq:maimlaproposal} q^\star = q - \frac{\delta}{2} \nabla V\left(\frac{1}{2}(q^\star+q)\right)+ \sqrt{\delta} p, \qquad p \sim N(0,I_m),\tag{5}\] where the Langevin equation is discretised via the implicit midpoint rule. For the corresponding Metropolis MCMC algorithm we use the acronym MAIMLA, Metropolis Adjusted Implicit Midpoint Algorithm. When compared with MALA, MAIMLA is of interest in some circumstances see [16]: it is geometrically ergodic in a rather general setting and, in addition, regardless of the choice of \(\delta\), the acceptance rate is 100% for Gaussian targets (and accordingly may be expected to be very high in the neighbourhood of a mode).

The key to the derivation of our results is the exploitation of symmetry inherent in the Metropolis-Hastings formula. For clarity in our presentation, we first focus on proposals that can be interpreted as approximations to Hamiltonian deterministic dynamics. For integrators that are volume preserving and reversible, the symmetries of the Hamiltonian formalism [17] can be used to show the average of the negative logarithm of the acceptance probability is equal to half of its variance. This relationship between the mean and variance is key in identifying the specific form of the limiting acceptance rate. After having discussed the Hamiltonian set-up we consider completely general Metropolis-Hastings algorithms, where it turns out that remarkably the same relation between the mean and the variance holds. This explains why all the scaling results in the literature lead to similar expressions for the limiting acceptance probability.

The paper is structured in six sections. In Section 2 we present a general setup to view proposals as discretisations of deterministic Hamiltonian dynamics; we then study the energy error, which determines the acceptance probability. The main scaling result in the Hamiltonian setting is given in Section 3. That section also contains the analysis of how best to tune the parameters to maximize the progress of the Markov chain. Also studied is the effect of preconditioning both the dynamics and the algorithms in situations where the target includes components whose scales are widely different. Section 4 leaves the Hamiltonian scenario and discusses the general Metropolis-Hastings case, while Section 5 contains the more technical parts of the proofs of the results. The final section provides an outlook on how to use discretisations available in the Hamiltonian integration literature to obtain samplers with arbitrarily favourable scalings. A simple numerical illustration is also provided.

2 Preliminaries↩︎

2.1 MCMC algorithms↩︎

We introduce an auxiliary variable \(p\) (the momentum) and the extended target in \((q,p)\) space \[\label{eq:extendedtarget} \Pi(q,p) \propto \exp(-H(q,p)),\qquad H(q,p) = \frac{1}{2}\| p\|^2+V(q),\qquad (q,p)\in\mathbb{R}^m\times\mathbb{R}^m.\tag{6}\] (\(\|\cdot\|\) always denotes the standard Euclidean norm.) The marginal on \(q\) is the target 1 and \(p\sim N(0,I_m)\).

The Hamiltonian function \(H\) (energy) gives rise to the Hamiltonian differential equations \[\label{eq:hamiltonianequations} \frac{d}{dt} q = p,\qquad\qquad \frac{d}{dt}p = -\nabla V(q).\tag{7}\] For each \(h>0\), the flow (solution map) \((q(0),p(0))\mapsto (q(h),p(h))\) of this system preserves the extended target 6 .

We consider MCMC algorithms with proposals generated with the help of one time step of a numerical integrator for 7 . The integrator is described by a transformation: \(\psi_h(q,p):\mathbb{R}^m\times\mathbb{R}^m \rightarrow \mathbb{R}^m\times\mathbb{R}^m\), parameterized by a steplength \(h>0\), such that for any solution of 7 , the value \((q(h),p(h))\in\mathbb{R}^m\times\mathbb{R}^m\) may be approximated by \((q^\star,p^\star) = \psi_h(q(0),p(0))\) in the sense that \[\label{eq:consistency} (q^\star,p^\star) -(q(h),p(h)) = o(h^\eta), \qquad h \rightarrow 0,\tag{8}\] for some \(\eta\geq 0\). Of course, the transformation \(\psi_h\) changes when the potential \(V\) is changed, but the dependence of \(\psi_h\) on \(V\) is not incorporated to this notation; we will take this point up below.

We are interested in integrators with the following properties, that are discrete counterparts of properties of the Hamiltonian dynamics:

Preservation of volume: \[\forall (q,p)\in\mathbb{R}^m\times\mathbb{R}^m\qquad |{\rm det}(\psi_h^\prime(q,p))| = 1,\]

where \(\psi_h^\prime(q,p)\) is the Jacobian matrix of \(\psi_h(q,p)\) with respect to \((q,p)\).

Reversibility: \[\forall (q,p)\in\mathbb{R}^m\times\mathbb{R}^m \qquad \psi_h(q,p)= (q^\star,p^\star) \Leftrightarrow \psi_h(q^\star,-p^\star) = (q,-p) .\]

Each transformation \(\psi_h\) gives rise to a Metropolis–Hastings MCMC algorithm as follows. If \(q\) is the current state of the Markov chain, draw \(p\sim N(0,I_m)\) (independent from \(q\) and from previous history), compute \((q^\star,p^\star) = \psi_h( q, p)\), and take \(q^\star\) as the proposal. If the integrator is volume preserving and reversible and the acceptance probability for the accept/reject mechanism is taken as \[\label{eq:acceptanceprobability} a(q,p) = \min(1,\exp(-\Delta(q,p,h))), \qquad \Delta(q,p,h) = H(\textcolor{blue}{q^\star},p^\star)-H(q,p),\tag{9}\] then this procedure defines a Markov chain in \(\mathbb{R}^m\) reversible with respect to the target 1 .

Some examples:

• Random Walk Metropolis (RWM): The transformation \((q^\star,p^\star) = \psi_h(q,p)\) reads \[q^\star = q+hp,\qquad p^\star = p.\] This is volume preserving and reversible. It satisfies 8 with \(\eta=0\), so that from a numerical analysis point of view the discretisation is not consistent with 7 (which is not surprising as it does not use the potential \(V\)). This lack of consistency poses no problem for the analysis below.

• MALA: The transformation \((q^\star,p^\star) = \psi_h(q,p)\) is given by the Störmer–Verlet integrator: \[p_{1/2} = p -\frac{h}{2}\nabla V(q), \qquad q^\star = q+hp_{1/2},\qquad p^\star = p_{1/2} -\frac{h}{2}\nabla V(q^\star).\] It is trivial to show that \(\psi_h\) is volume preserving and reversible; for sufficiently smooth potentials, property 8 holds with \(\eta = 2\) (one says that Störmer–Verlet is a consistent second order integrator).

  After eliminating \(p_{1/2}\) and setting \[\label{eq:hdelta} \delta = h^2\tag{10}\] the formula to compute the proposal may be written in the standard form 3 . Note that \(\delta\) is the variance of the random term \(\sqrt{\delta} p\) in the definition of \(q^\star\).

• MAIMLA: Now, \(\psi_h\) is defined by the implicit midpoint discretisation: \[q^\star -q = \frac{h}{2} (p^\star+p), \qquad p^\star-p = -h \nabla V\Big(\frac{1}{2}(q^\star+q)\Big).\] Finding \(q^\star\) requires the solution of a system of \(m\) nonlinear equations. We assume that, for each \((q,p)\) the system has a unique solution for sufficiently small \(h\). Again, it is trivial to check that the transformation is volume preserving and reversible and 8 holds with \(\eta=2\) (second order).

Many other integrators are possible, as discussed in Section 6.

2.2 Hamiltonian equations: properties of the energy error↩︎

For solutions of 7 , \(H(q(t),p(t))\) does not vary with \(t\) (conservation of energy). Therefore, the value \(H(q,p)\) in 9 coincides with the value \(H(q(h),p(h))\), where \((q(t),p(t))\) denotes the solution with initial value \((q(0),p(0))=(q,p)\) and thus \(\Delta(q,p,h)\) is the energy error \(H(q^\star,p^\star)-H(q(h),p(h))\) resulting from approximating the true value \((q(h),p(h))\) by means of the numerical solution \((q^\star,p^\star) =\psi_h(q,p)\).

We will consider the following (mild) conditions on \(\Delta\).

Condition 1. There exist real functions \(\alpha(q,p)\), \(\rho(q,p,h)\) and a constant \(\nu>0\) such that \[\Delta(q,p,h) = h^{\nu}\alpha(q,p)+h^\nu\rho(q,p,h),\qquad \lim_{h\rightarrow 0}\rho(q,p,h) = 0.\] A Taylor expansion of \(\Delta\) as a function of \(h\) shows that this condition will typically hold with \(\nu = \eta+1\) (see 8 ) provided that \(V\) is sufficiently smooth. To be precise, for the examples considered above one may easily prove:

• RWM: For continuously differentiable \(V\), \(\nu = 1\) and (superscripts denote Cartesian components of the vectors \(q\) and \(p\)) \[\label{eq:alpharwm} \alpha(q,p) = \sum_{i=1}^m \frac{\partial V}{\partial q^i}p^i.\tag{11}\]

• MALA: This has \(\nu=3\) and \[\label{eq:alphamala} \alpha(q,p) = \frac{1}{4} \sum_{i,j=1}^m \frac{\partial^2 V}{\partial q^i \partial q^j}\frac{\partial V}{\partial q^i} p^j - \frac{1}{12} \sum_{i,j,k=1}^m \frac{\partial^3 V}{\partial q^i \partial q^j\partial q^k} p^ip^jp^k,\tag{12}\] provided that \(V\) is three times continuously differentiable.

• MAIMLA: For \(V\) three times continuously differentiable, \(\nu = 3\) and \[\label{eq:alphamaimla} \alpha(q,p) = \frac{1}{24}\sum_{i,j,k=1}^m \frac{\partial^3 V}{\partial q^i \partial q^j\partial q^k} p^ip^jp^k.\tag{13}\]

The three integrators we have just described have odd values of \(\nu\). This happens for all reversible algorithms [17]. Algorithms with \(\nu>3\) are discussed in Section 6.

The summations that feature in the formulas 11 –13 involving derivatives of the potential and components of the momentum are studied in the theory of symplectic integration, where they are called elementary Hamiltonians [17]. Their structure may be easily described by using suitable graphs.

While Condition 1 describes the behaviour of \(\Delta\) as a function of \(h\), the following condition controls the variation of \(\Delta\) as a function of \((q,p)\).

Condition 2. There exists a real function \(D(q,p)\) such that \[\textcolor{black}{\sup_{0<h\leq 1} \frac{|\Delta(q,p,h)|^3}{h^{3\nu}} \leq D(q,p)}\] and \(D\) is integrable with respect to the probability density 6 .

For typical integrators, this condition will hold in particular whenever \(\exp(-V(q))\) decays exponentially as \(\|q\|\rightarrow \infty\) and the derivatives of \(V\) grow at most polynomially. For RWM, MALA and MAIMLA, Taylor expansions in powers of \(h\) with integral remainder show that it is sufficient that \(V\) has polynomially growing continuous derivatives of order less than or equal to \(2, 4, 4\) respectively.

The properties of reversibility and volume preservation impose remarkable symmetry on the geometric behaviour of \(\psi_h\) (see [18]). That symmetry has important implications for the mean, second moment and variance of \(\Delta(q,p,h)\) \[\label{eq:mean} \mu(h) = \mathbb{E}(\Delta(q,p,h)),\qquad s(h) = \mathbb{E}(\Delta(q,p,h)^2),\qquad \sigma^2(h) = s(h)-\mu(h)^2.\tag{14}\] (Expectations \(\mathbb{E}\) are taken with respect to the extended target 6 .) We have the following key result, which is essentially [18], a generalization of [7]. These references operate with a version of Condition 2 slightly weaker than that just considered. Here we have strengthened this hypothesis in order to avoid the somewhat obscure use of the dominated convergence theorem made in the proofs provided in those references.

2.3 Scaling the potential↩︎

We now investigate the effect of scaling the potential. Let \(\lambda>0\) be a scaling factor, and consider the scaled potential \(V^\lambda(q) = V(q/\lambda)\), where \(V\) is as above. We shall associate the superindex \(\lambda\) with objects pertaining to the scaled potential. The scaled probability density is \(\pi^\lambda(q) \propto \exp( -V^\lambda(q))\), which leads to the scaled Hamiltonian \(H^\lambda(q,p)= (1/2)\|p\|^2+V(q/\lambda)\) with Hamiltonian equations \[\frac{d}{dt} q = p,\qquad\qquad \frac{d}{dt}p = -\frac{1}{\lambda}\nabla V(q/\lambda).\]

After the change of variables \[\bar{q} = q/\lambda,\qquad \bar{p}= p,\qquad \bar{t} = t/\lambda\] the scaled equations of motion become \[\frac{d}{d\bar{t}} \bar{q} = \bar{p},\qquad\qquad \frac{d}{d\bar{t}}\bar{p} = -\nabla V(\bar{q}),\] a system that coincides with the original, unscaled 7 . Thus, for the differential equations, the effect of scaling the potential may be reproduced by using the unscaled dynamics as applied with a scaled time and scaled dependent variables. We now demand that the same happens at the level of the integrator, i.e. if \(\psi_h(q,p)\) is the transformation associated with the original potential \(V\) and \(\psi_h^\lambda(q,p)\) the transformation when the integrator is used for the rescaled differential equations: \[\label{eq:equivariance} (q^\star,p^\star) = \psi_h^\lambda(q,p)\Leftrightarrow (q^\star/\lambda,p^\star) = \psi_{h/\lambda}(q/\lambda,p).\tag{16}\] We will say that an integrator is equivariant if the two following processes lead to the same result: (i) rescaling the variables in the differential system and then applying the integrator, (ii) applying the integrator to the original differential equations and then rescaling variables in the discrete equations. The notion of equivariance is akin to the notion of dimensional correctness in physics. All integrators of any interest are equivariant; the performance of a nonequivariant integrator would depend on the particular units used to measure the quantities involved.

For an equivariant integrator the energy error \(\Delta^\lambda(q,p,h)\) for the scaled potential satisfies \[\Delta^\lambda(q,p,h)= \Delta(q/\lambda,p,h/\lambda),\] where \(\Delta\) is the energy error for the original potential. Furthermore we have the following elementary result:

In addition, with the notation \[\label{eq:momentsoldandnew} \mu^\lambda(h) = \mathbb{E}^\lambda(\Delta^\lambda(q,p,h)),\qquad s^\lambda(h) = \mathbb{E}^\lambda(\Delta^\lambda(q,p,h)^2),\qquad \sigma^{2,\lambda}(h) = s^\lambda(h)-\mu^\lambda(h)^2,\tag{17}\] where \(\mathbb{E}^\lambda\) refers to expectations with respect to the density \(\Pi^\lambda(q,p)\propto\exp(-H^\lambda(q,p))\). Equation 16 and Lemma 1 show that, for an equivariant integrator, \[\label{eq:scaledmoments} \mu^\lambda(h) = \mu(h/\lambda),\qquad s^\lambda(h)=s(h/\lambda), \qquad \sigma^{2,\lambda}(h) = \sigma^2(h/\lambda),\tag{18}\] where the quantities in the right hand-sides refer to the unscaled potential and were defined in 14 .

3 Optimal scaling using the Hamiltonian framework↩︎

For \(d\geq 1\), we set \(N=dm\) and write vectors in \(\mathbb{R}^N\) as \(Q=(q_1,\dots,q_d)\) with \(q_i\in\mathbb{R}^m\). We fix the potential \(V\) in \(\mathbb{R}^m\) as in 1 and consider the target probability density in \(\mathbb{R}^N\) defined by \[\label{eq:product32target} \pi^\lambda(Q)\propto \prod_{i=1}^d \exp(-V(q_i/\lambda_i)),\tag{19}\] where \(\lambda_i>0\) are scaling constants. After defining the potential \[{\cal V}(Q) = {\cal V}(q_1,\dots,q_d) = \sum_{i=1}^d V(q_i/\lambda_i),\] the target is of the form \(\pi^\lambda(Q)\propto \exp(-{\cal V}(Q))\) and the algorithms in Section 2 are applicable. The Hamiltonian that features in the extended target \(\propto \exp(-{\cal H}(Q,P))\) is \[{\cal H}(Q,P) = {\cal H}(q_1,\dots,q_d,p_1\dots,p_d) = \frac{1}{2} \sum_{i=1}^d \|p_i\|^2+\sum_{i=1}^d V(q_i/\lambda_i),\] and the Hamiltonian equations in \(\mathbb{R}^N\times\mathbb{R}^N\) are obtained by simply juxtaposing the \(d\) Hamiltonian systems in \(\mathbb{R}^m\times\mathbb{R}^m\) corresponding to the different \((q_i,p_i)\). There is no coupling between the differential equations for the different components and the same is true for the numerical integrator. However in the accept/reject mechanism the components come together because the acceptance probability is \(\min(1,\exp(-\Delta_d))\) with (stars denote proposals) \[\Delta_d (Q,P,h) = {\cal H}(q_1^\star,\dots,q_d^\star,p_1^\star,\dots,p_d^\star) - {\cal H}(q_1,\dots,q_d,p_1\dots,p_d).\] With the notation as in the previous section, \[\Delta_d (Q,P,h) = \sum_{i=1}^d \Delta^{\lambda_i}(q_i,p_i,h).\]

3.1 A general scaling result↩︎

We are interested in finding conditions that ensure that \(\Delta_d\) has a distributional limit as \(d\rightarrow \infty\). It is reasonable to assume that the time-step \(h\) in the integration has to shrink as \(d\) increases, so that the larger number of terms being summed in the preceding display is compensated by a decrease in the size of the individual terms. We denote by \(h_d\) the integration time-step to be used when there are \(d\) components.

Here is the main result of this paper:

A simple outline of the proof of this result will be given presently; the missing technical details are provided in Section 5 below. Since the integrator is assumed to be equivariant, according to 16 , \[\label{eq:Deltad} \Delta_d = \sum_{i=1}^d \Delta(q_i/\lambda_i,p_i,h/\lambda_i).\tag{23}\] We will use a Central Limit Theorem to identify the distributional limit of \(\Delta_d\) as \(d\rightarrow \infty\). We consider the triangular array of random variables whose \(d\)-row, \(d=1,2,\dots\), is \[\label{eq:row} \Delta(q_1/\lambda_1,p_1,h_d/\lambda_1), \Delta(q_2/\lambda_2,p_2,h_d/\lambda_d),\dots, \Delta(q_d/\lambda_d,p_1,h_d/\lambda_d).\tag{24}\] Each \(q_i\) is distributed as \(\propto \exp(-V(q_i/\lambda_i))\) and each \(p_i\) is normal with zero mean and covariance matrix \(I_m\); all these variables are mutually independent. According to 17 , the variances of the variables in 24 are \[\sigma^2(h_d/\lambda_1), \sigma^2(h_d/\lambda_2),\dots, \sigma^2(h_d/\lambda_d),\] and by Theorem 1, the variance of \(\Delta_d\) will approximately be \(\sum_{i=1}^d {h_d^{2\nu}}/{\lambda_i^{2\nu}} \Sigma\) ; the corresponding expectation will be approximately a half of this value. Using the definition of \(h_d\): \[\sum_{i=1}^d \frac{h_d^{2\nu}}{\lambda_i^{2\nu}} \Sigma = \ell^{2\nu} \: d^{-\gamma} \sum_{i=1}^d \frac{1}{\lambda_i^{2\nu}}\Sigma,\] a quantity that, according to the hypotheses on the \(\lambda_i\), tends to \(\ell^{2\nu} K \Sigma\) as \(d\rightarrow \infty\). By invoking the Central Limit Theorem for triangular arrays (21 is required to check the Lindeberg condition), we conclude that \(\Delta_d\) converges in distribution to a normal variable \(\Delta_\infty\) of mean \((1/2) \ell^{2\nu}K\Sigma\) and variance \(\ell^{2\nu}K\Sigma\). From the boundedness of \(x\mapsto \min(1,\exp(-u))\), \(\mathbb{E}(a_d)\) converges to \(\mathbb{E}(\min(1,\exp(-\Delta_\infty))\); the last expectation may be computed analytically and has the value 22 . This concludes the outline of the proof.

Here are some useful choices of scaling constants:

1. When \(\lambda_i=1\) for each \(i\) (the different \(q_i\) share the common density \(\propto\exp(-V)\)), \(\gamma=1\), \(K=1\) (21 trivially holds) and \(h_d = \ell/d^{1/(2\nu)}\).

2. For \(\lambda_i = 1/i^\kappa\), with \(\kappa>0\), \[\sum_{i=1}^d \frac{1}{\lambda_i^{2\nu}} = \sum_{i=1}^d i^{2\nu\kappa} \sim \frac{d^{2\nu\kappa+1}}{2\nu\kappa+1},\qquad d\rightarrow \infty,\] and therefore \(\gamma =2\nu\kappa+1\), \(K = 1/(2\nu\kappa+1)\) leading to \(h_d = \ell/d^{\kappa+1/(2\nu)}\). Again 21 holds.

Note that for choice 2 the exponent of \(d\) in the denominator of \(h_d\) exceeds the exponent \(1/(2\nu)\) of the \(\lambda_i = 1\) case by an amount \(\kappa\). In fact, for choice 2, the time step \(h_d\) has to be selected to cater for the small \(\lambda_i\) components for which \(q_i\) undergoes fast changes. As a result the components \(q_i\) with large \(\lambda_i\) are integrated with a time step that is too small relative to their own rate of change. The need to use very small values of \(h_d\) may be circumvented by using preconditioning (see Section 3.3). On the other hand, if \(\kappa\) increases, then \(K = 1/(2\nu\kappa+1)\) decreases and the acceptance rate in 22 increases; this was of course to be expected because \(h_d = \ell/d^{\kappa+1/(2\nu)}\) is a decreasing function of \(\kappa\).

Let us finish this section by applying Theorem 2 to the algorithms considered in Section 2:

• RWM: This has \(\nu=1\). For the i.i.d. case with \(\lambda_i=1\), where \(\gamma=1\), \(K=1\), the scaling is \(h_d = \ell/\sqrt{d}\) (the variance of the proposal will then be \(\delta_d=\ell^2/d\)). Under this scaling, the limit expected acceptance rate will be \[A = 2 \Phi\left(-\frac{\ell\sqrt{\Sigma}}{2}\right)\] with (see 11 ) \[\Sigma = \mathbb{E}\big( \left(\sum_{i=1}^m \frac{\partial V}{\partial q^i}p^i\right)^2 \big) = \mathbb{E}\big( \sum_{i=1}^m \left(\frac{\partial V}{\partial q^i}\right)^2 \big)=\mathbb{E}\big(\|\nabla V(q)\|^2\big),\] where we have taken into account that the \(p^i\) are uncorrelated and have unit variance. In the particular case \(m=1\), the scaling \(\delta_d = \ell^2/d\) and the formula for \(A\) were first presented in the pioneering paper [2].

• MALA: This has \(\nu=3\). For the i.i.d. case with \(\lambda_i=1\), where \(\gamma=1\), \(K=1\), the scaling is therefore \(h_d = \ell/d^{1/6}\), i.e.\(\delta_d = \ell^2/d^{1/3}\). Under this scaling, the limit expected acceptance rate is \[A = 2 \Phi\left(-\frac{\ell^3\sqrt{\Sigma}}{2}\right),\] where \(\Sigma = \mathbb{E}(\alpha^2)\) with \(\alpha\) given in 12 .

  For the particular case \(m=1\), the scaling of \(\delta_d\) and the formula for \(A\) were first derived in [3] (see also [5]). In that particular case, noting that \(p\sim N(0,1)\) and integrating by parts in the integral for the expectation, one may show that \[\Sigma = \frac{1}{48} \mathbb{E}\big(5V^{\prime\prime\prime}(q)^2+3 V^{\prime\prime}(q)^3\big).\]

  For the scaling \(\lambda_i = 1/i^\kappa\), the theorem gives \(h_d = \ell/d^{\kappa+1/6}\) or \(\delta_d = \ell^2/d^{2\kappa+1/3}\), a result that, for the particular case \(m=1\), was derived in [6]; this reference does not identify the value of \(A\), which, as we shall discuss later, is required to derive optimal scaling results and was given in [19].

• MAIMLA: The only difference with MALA is in the value of \(\Sigma\). A comparison between 12 and 13 suggests that, in general, \(|\alpha_{MALA}|> |\alpha_{MAIMLA}|\) so that MAIMLA enjoys larger acceptance rates for a given steplength. Numerical experiments confirm this [16].

3.2 Squared Jumping Distance↩︎

According to 22 , the expected acceptance rate approaches \(100\%\) as \(\ell\rightarrow 0\). However very small values of \(\ell\) entail proposals that are too close to the current state and result in large correlations in the Markov chain. To analyse the optimal choice of \(\ell\), it is useful to consider the squared jumping distance defined as \[\label{eq:sjd} {\cal SJD}_d = \frac{1}{d} \sum_{i=1}^d \mathbb{E}(\|q_i^\star-q_i\|^2 \mathbb{I} _{\{a_d(Q,P,h_d) \leq U\}}),\qquad U\sim{\rm Uniform}(0,1);\tag{25}\] here, a star denotes proposal and the indicator function determines whether the proposal is accepted or otherwise. Note the averaging over the different \(q_i\).

For the RWM algorithm applied to the target 1 , \(q^\star-q = h p\). For algorithms (including MALA and MAIMLA) based on integrators consistent with the Hamiltonian system 7 , \(q^\star-q = h p+ o(h)\). More generally, we introduce the following condition:

Condition 3. When the integrator is applied to the target 1 , as \(h\rightarrow 0\), \[r(h): = \mathbb{E}\left(\|\frac{1}{h}(q^\star-q)\|^2-\|p\|^2 \right) \rightarrow 0.\] (Here \(q \sim \exp(-V(q))\) and \(p \sim N(0,I_m)\).)

This condition is easily checked by Taylor expansion of \(q^\star\) as a function of \(h\). In particular, it will hold if \(\nabla V(q)\) grows at most as a polynomial and \(\exp(-V(q))\) has tails that decrease exponentially.

The next result identifies the limit of the squared jumping distance:

The proof may be seen in Section 5.

The theorem suggests that \(\ell\) should be chosen so as to maximize \[E = \ell^2 A(\ell) = 2 \ell^2 \Phi\left(-\frac{\ell^{\nu}\sqrt{K\Sigma}}{2}\right).\] In practice \(\Sigma\) will be unavailable; however it is well known and remarkable that, in situations like this, the optimization problem for \(E\) may be solved by using \(A\) as an independent variable in lieu of \(\ell\). In terms of \(A\), \[\begin{align} \label{eq:EA} E = \frac{2^{2/\nu}}{(K\Sigma)^{1/\nu}}A \left(\Phi^{-1}(1-\frac{A}{2})\right)^{2/\nu} , \end{align}\tag{26}\] and clearly the value \(A_{\rm opt}\) that maximizes \(E\), while changing with \(\nu\), is independent of \(\Sigma\) and \(K\) and therefore of the target and the scaling constants. Values of \(A_{\rm opt}\) are provided in Table 1. Recall that necessarily reversible integrators cannot have an even value of \(\nu\).

With implicit algorithms, like MAIMLA, enlarging the value of the step length \(h\) may result in the need for more iterations to solve the nonlinear equations at each step. In that scenario, maximizing the squared jumping distance per step as above will be different from maximizing the squared jumping distance per unit of computational work. It is likely that in practice for implicit algorithms the acceptance rate should be chosen to be slightly higher than the values we just quoted.

3.3 Preconditioned dynamics↩︎

To sample from the target 1 , one may consider, rather than 6 , the extended target \[\label{eq:extendedtargetmass} \Pi^M(q,p)\propto \exp(-H^M(q,p)),\qquad H^M(q,p) = \frac{1}{2}p^TM^{-1} p+V(q),\qquad (q,p)\in\mathbb{R}^m\times\mathbb{R}^m,\tag{27}\] where \(M\) is a positive definite symmetric \(m\times m\) matrix, referred to as mass matrix. The auxiliary momentum variable \(p\) now has distribution \(N(0,M)\) and Hamilton’s equations read \[\label{eq:hamiltonianequationsM} \frac{d}{dt} q = M^{-1}p,\qquad\qquad \frac{d}{dt}p = -\nabla V(q).\tag{28}\] The integrators described in Section 2 are easily extended to this more general setting and, after defining the acceptance probability \[\label{eq:acceptanceprobabilityM} a(q,p,h) = \min(1,\exp(-\Delta^M(q,p,h))), \qquad \Delta^M(q,p,h) = H^M(\textcolor{blue}{q^\star},p^\star)-H^M(q,p),\tag{29}\] and considering the marginal on \(q\), they generate Markov chains reversible with respect to the target 1 .

For reasons of brevity, we only write down the formulas for the Störmer–Verlet integrator. This reads: \[p_{1/2} = p -\frac{h}{2}\nabla V(q), \qquad q^\star = q+hM^{-1}p_{1/2},\qquad p^\star = p_{1/2} -\frac{h}{2}\nabla V(q^\star).\] Elimination of \(p_{1/2}\) and 10 yield \[q^\star = q - \frac{\delta}{2} M^{-1}\nabla V(q)+ \sqrt{\delta} M^{-1}p, \qquad p \sim N(0,M),\] or \[q^\star = q - \frac{\delta}{2} M^{-1}\nabla V(q)+ \sqrt{\delta} M^{-1/2}z, \qquad \textcolor{black}{z} \sim N(0,I_m),\] a consistent discretisation of the preconditioned Langevin equation \[dq(t) = -\frac{1}{2}M^{-1}\nabla V(q(t))\,dt+ M^{-1/2}\, dW(t),\] which preserves the target 1 . Thus the introduction of mass matrices different from \(I_m\) in the Hamiltonian formulation provides preconditioning for the dynamics.

Let us study the effect of scaling the potential as in Section 2.3. Assume that to sample from \(V(q/\lambda)\), we choose the mass matrix \(M\) to be \(\lambda^{-2}I_m\). From a physical point of view, this means increasing the mass as \(\lambda\) decreases and the potential becomes steeper; this avoids large values of the velocity \((d/dt) q\), see [20], [21] for further discussion of this point. For this choice of \(M\) the differential equations 28 are \[\frac{d}{dt} q = \lambda^2 p,\qquad\qquad \frac{d}{dt}p = -\frac{1}{\lambda}\nabla V(q/\lambda),\] and the scaling of variables \[\bar{q} = q/\lambda,\qquad \bar{p} = \lambda p\] transforms them into \[\frac{d}{dt} \bar{q} = \bar{p},\qquad\qquad \frac{d}{dt}\bar{p} = - \nabla V(\bar{q}),\] i.e.the Hamiltonian differential equations for the unscaled potential when the mass matrix is chosen to be the identity. Therefore, for an equivariant integrator and \(M = \lambda^{-2}I_m\): \[\Delta^M(q,p,h) = \Delta(q/\lambda,\lambda p, h),\] where \(\Delta\) corresponds to the integrator as applied to 6 . In addition we have the following analogue of Lemma 1, whose elementary proof will not be given:

Therefore, if to sample from the product target 19 in \(\mathbb{R}^N\), we use the Hamiltonian \[\label{eq:HM} {\cal H}^M(Q,P) = {\cal H}^M(q_1,\dots,q_d,p_1\dots,p_d) = \frac{1}{2} \sum_{i=1}^d \|\lambda_ip_i\|^2+\sum_{i=1}^d V(q_i/\lambda_i),\tag{30}\] then the energy error to accept/reject will satisfy \[\Delta_d^M(Q,P,h) = \sum_{i=1}^d \Delta^M(q_i,p_i,h) = \sum_{i=1}^d \Delta(q_i/\lambda,\lambda p_i,h).\] This is very similar to 23 and reduces the study of \(\Delta_d^M\) to the study of the energy errors for the potential \(V\) and unit mass matrix. However, we note that here and as distinct from 23 , a common value of \(h\) features in the different \(\Delta(q_i/\lambda,\lambda p_i,h)\) being summed: preconditioning ensures that the different \(q_i\) evolve in the same time scale. Therefore the present situation is analogous to the setting where in Theorem 2 all the \(\lambda_i\) take the value 1. By mimicking the proof of Theorem 2, one derives the following result:

Note that there are no hypotheses on the \(\lambda_i\) and that, when preconditioning is used, the scaling for \(h_d\) coincides with the one found in the nonpreconditioned case when all the \(q_i\) share a common distribution. In the particular case of the MALA algorithm with \(m=1\) and \(\lambda_i = 1/i^\kappa\) the formula for \(h_d\) was presented in [6], and the value of \(A\) was given in [19].

Let us now discuss the squared jumping distance for preconditioned algorithms. With preconditioning, due to the presence of the mass matrix, \(q_i^\star-qi \approx h_d \lambda_i^2 p_i\) and \(p_i\sim N(0,(1/\lambda_i^2)I_m)\). Therefore, we would expect that \[\mathbb{E}(\|q_i^\star-q_i\|^2) \approx m h_d^2 \lambda_i^2;\] in this way the mean quadratic displacement from current state to proposal varies with \(i\) and is proportional to \(\lambda_i^2\), and therefore to the second moment of \(q_i\), which has distribution \(\exp(-V(q_i/\lambda_i))\). (Achieving this proportionality is precisely the motivation for preconditioning.) The last display should be compared with the non-preconditioned case where (cf.Condition 3) regardless of the value of \(i\) \[\mathbb{E}(\|q_i^\star-q_i\|^2) \approx m h_d^2.\] Due to this difference, the metric defined in 25 , based on an average with respect to \(i\), makes little sense for preconditioned dynamics. One may instead use, for \(i\leq d\), the alternative metric \[{\cal SJD}_d^i = \mathbb{E}(\|q_i^\star-q_i\|^2 \mathbb{I} _{\{a_d(Q,P,h_d) \leq U\}}),\qquad U\sim{\rm Uniform}(0,1).\]

With the scaling \(h_d= \ell/d^{1/(2\nu)}\) from Theorem 5, one has \[\begin{align} {\cal SJD}_d^i \sim m \lambda_i^2 d^{-1/\nu}\ell^2 A, \quad d \rightarrow \infty. \end{align}\] The main difficulty in proving this estimate stems from \(\|q_i^\star-q_i\|^2\) and \(\mathbb{I} _{\{\Delta_d(Q,P,h_d) \leq U\}}\) not being independent. This obstacle may be circumvented by using the technique in the proof of [7], that takes into account that, for \(d\) large, \(\mathbb{I} _{\{\Delta_d(Q,P,h_d) \leq U\}}\) is ’almost’independent of the single component \((q_i,p_i)\). Details will not be given. The proof used in this paper for Theorem 3 does not have to address the dependence issue because it deals with an average \((1/d)\sum_i\|q_i^\star-q_i\|^2\approx (1/d)h_d^2\sum_i\|p_i\|^2\) whose limit, being constant by the weak law of large numbers, is stochastically independent from the indicator function.

In any case, as in the non-preconditioned case, the best value of \(\ell\) is determined by maximising \(\ell^2A\). As discussed above, the maximum is achieved for the values reported in Table 1.

4 General Metropolis–Hastings algorithms↩︎

We now leave the Hamiltonian framework and study general Metropolis–Hastings algorithms. For the target 1 , consider a proposal of \(q^\star\) from a position \(q\) according to the probability density function \(\wp_h(q,q^\star)\). Here \(h\) is a small parameter for the algorithm, perhaps related to the standard deviation of the proposal. Acceptance is based on the ratio \[\frac{\exp(-V(q^\star))\: \wp_h(q^\star,q)}{\exp(-V(q))\:\wp_h(q,q^\star)};\] swapping \(q\) and \(q^\star\) takes this ratio into its inverse. To parallel the notation of the preceding sections, we introduce \[\Delta(q,q^\star,h) = V(q^\star)-V(q)-\log\left(\frac{\wp_h(q^\star,q)}{\wp_h(q,q^\star)}\right)\] and then the acceptance probability is \(\min\{1,\exp(-\Delta(q,q^\star,h))\}\). Note that \[\label{eq:oppositeqstar} \Delta(q,q^\star,h) = -\Delta(q^\star,q,h);\tag{31}\] this just rephrases the symmetry of the Metropolis-Hastings formula, but is the foundation of all subsequent developments. We encountered a similar relation in 15 when working in the Hamiltonian framework.

The relation 31 implies that when computing \(\mathbb{E}(\Delta)\) much cancellation takes place. More precisely, after defining \[\mu(h) = \mathbb{E}[\Delta(q,q^\star,h)], \qquad s(h) = \mathbb{E}[\Delta(q,q^\star,h)^2], \qquad \sigma^2(h)=s(h)-\mu(h)^2,\] (expectations \(\mathbb{E}\) are with respect to the probability measure \(\Pi(q,q^\star)\propto\exp(-V(q))\wp(q,q^\star)dq dq^\star\)) we have the following analogue to Theorem 1. Note, once more, the property ‘’mean equals half the variance’’.

We now study the scaling for products as in 19 and, to simplify the exposition, only consider the case where the target is given as an i.i.d. product, with \(\lambda_i=1\) for all \(i\). (The scaled case may be dealt with via dimensionality arguments as we did in the Hamiltonian framework.) We only give the main idea. The acceptance probability may be written as \[a_d(Q,Q^\star,h) = \min\{1,e^{-\Delta_d(Q,Q^\star,h)}\}, \qquad \Delta_d = \sum_{i=1}^d \Delta(q_i,q_i^\star,h).\] As a consequence of Theorem 6, \(\{\Delta(q_i,q_i^\star,h)\}_{i=1}^d\) is a sequence of i.i.d. \(d\) random variables with mean approximately \((1/2) h^{2\nu}\Sigma\) and variance approximately \(h^{2\nu}\Sigma\). Setting \(h_d = \ell d^{-1/(2\nu)}\) and applying the Central Limit Theorem, the distributions of \(\Delta_d\) will approach the distribution of a normal variable \(\Delta_\infty\) of mean \((1/2)\ell^{2\nu}\Sigma\) and variance \(\ell^{2\nu}\Sigma\) and, by arguing as in Theorem 2, the limit acceptance probability will have the value \[A = 2\Phi\left(-\frac{\ell^\nu\sqrt{\Sigma}}{2}\right).\] Of course this is exactly 22 (for the case at hand with \(K=1\)). Note that the expression for \(A\) is a consequence of the ‘’mean equal half the variance’’property of the limit normal variable \(\Delta_\infty\). It is clear that the values of the optimal acceptance probability reported in Table 1 are valid in the present scenario.

As a first, simple illustration, we apply these considerations to RWM, a sampler we already discussed within the Hamiltonian framework. This has \[\Delta(q,q^\star,h) = V(q^\star)-V(q).\] Conditional on \(q\), the proposal has a distribution \(N(q, h^2I_m)\) and then, with \(p\sim N(0,I_m)\), \[s(h) = \mathbb{E}[\Delta(q,q^\star)^2] = \mathbb{E}[(V(q+hp)-V(q))^2].\] Therefore, for suitably smooth \(V\), \(s(h)/h^{2}\) approaches \(\Sigma = \mathbb{E}\big(\sum_i (\partial V/\partial q^ip^i)^2\big)\) as \(h\rightarrow 0\); so \(\nu=1\) in Theorem 6 and we then recover the RWM results in Section 3.

As a second illustration, we consider fMALA, a modification of MALA constructed in [4] to achieve \(\delta_d = \ell^2/d^{1/5}\) scaling in the i.i.d. product case. This scaling for fMALA (and for its variants) is attained at the price of using, in addition to the gradient, the Hessian and the tensor of third derivatives of the potential. We have not been able to cast fMALA in the Hamiltonian framework, but the material in this section may be applied. By expanding \(s(h)\), one shows that the algorithm satisfies the requirements in Theorem 6 with \(\nu=5\); in fact fMALA was constructed by annihilating the lower order terms in the expansion of \(s(h)\) in powers of \(h\). The material in this section then implies that the variance of the proposal has to be scaled as \(\delta_d = \ell^2/d^{1/5}\) and the optimal acceptance rate is 0.704. It is perhaps of interest to observe that [4] implicitly proves the ‘’mean equal half the variance’’property (see formulas (47)–(48) in that paper), but does so by finding the expansion of \(\mathbb{E}(\Delta)\) in powers of \(h\) to order 10, a task that required the use of a symbolic algebra package. Also [4] was limited to the univariate case \(m=1\).

5 Missing proofs↩︎

Proof of Theorem 2:

We begin by proving that for the triangular array in 24 , the variance \({\rm Var}_d\) of the sum \(\Delta_d\) of the independent variables in the \(d\)-th row converges to \(\ell^{2\nu}K\Sigma\) as \(d\rightarrow\infty\). We decompose as \[{\rm Var}_d = \sum_{i=1}^d \sigma^2(h_d/\lambda_i) = \sum_{i=1}^d \left(\sigma^2(h_d/\lambda_i)- \frac{h_d^{2\nu}}{\lambda_i^{2\nu}} \Sigma\right)+ \sum_{i=1}^d \frac{h_d^{2\nu}}{\lambda_i^{2\nu}} \Sigma.\] and recall from the discussion following Theorem 2, that the last of these sums converges to \(\ell^{2\nu}K\Sigma\). On the other hand \[\begin{align} \sum_{i=1}^d \left(\sigma^2(h_d/\lambda_i)- \frac{h_d^{2\nu}}{\lambda_i^{2\nu}} \Sigma\right)& = & h_d^{2\nu} \sum_{i=1}^d\frac{1}{\lambda_i^{2\nu}} \left(\frac{\sigma^2(h_d/\lambda_i)}{ \frac{h_d^{2\nu}}{\lambda_i^{2\nu}} }-\Sigma\right)\nonumber\\ &=& \ell^{2\nu}d^{-\gamma}\sum_{i=1}^d\frac{1}{\lambda_i^{2\nu}} \left(\frac{\sigma^2(h_d/\lambda_i)}{ \frac{h_d^{2\nu}}{\lambda_i^{2\nu}} }-\Sigma\right)\nonumber\\ &\leq&\ell^{2\nu}\left( d^{-\gamma} \sum_{i=1}^d\frac{1}{\lambda_i^{2\nu}}\right) \max_{i=1,\dots,d}\left|\frac{\sigma^2(h_d/\lambda_i)}{\frac{h_d^{2\nu}}{\lambda_i^{2\nu}}} -\Sigma\right|.\label{eq:boundvariance} \end{align}\tag{33}\] Now, according to 21 , the quotients \(h_d/\lambda_i\), \(i=1,\dots,d\) converge uniformly to zero: \[\label{eq:uniformity} \lim_{d\rightarrow \infty}\max_{1,\dots,d}\frac{h_d}{\lambda_i} = \lim_{d\rightarrow \infty} \max_{1,\dots,d} \ell \left(d^{-\gamma}\frac{1}{\lambda_i^{2\nu}}\right)^{1/(2\nu)} = 0.\tag{34}\] Therefore, from Theorem 1, \[\lim_{d\rightarrow\infty}\max_{i=1,\dots,d}\left|\frac{\sigma^2(h_d/\lambda_i)}{\frac{h_d^{2\nu}}{\lambda_i^{2\nu}}} -\Sigma\right| = 0,\] which, in tandem with 20 , implies that 33 converges to zero. We conclude that \({\rm Var}_d\rightarrow \ell^{2\nu}K\Sigma\), as we had announced. In a similar way, one proves that the expectation of \(\Delta_d\) converges to \(\ell^{2\nu}K\Sigma/2\).

It remains to show that the Lindeberg condition is fulfilled. Fix \(\epsilon>0\) and consider the centered variables \(\Delta(q_i/\lambda_i,p_i,h_d/\lambda_i)-\mu(h_d/\lambda_i\)), \(i=1,\dots,d\), and the quotient \[\frac{\sum_{i=1}^d\mathbb{E}\Big( \big(\Delta(q_i/\lambda_i,p_i,h_d/\lambda_i)-\mu(h_d/\lambda_i)\big)^2 \mathbb{I} _{\{|\Delta(q_i/\lambda_i,p_i,h_d/\lambda_i)-\mu(h_d/\lambda_i)| > \epsilon{\rm Var}_d^{\frac{1}{2}}\}} \Big)}{{\rm Var}_d}.\] As we have just proved, the denominator has a finite limit and therefore our task is to prove that numerator converges to zero. By using Lemma 1, this numerator may be written as \[\sum_{i=1}^d\mathbb{E}\Big( \big(\Delta(q,p,h_d/\lambda_i)-\mu(h_d/\lambda_i)\big)^2 \mathbb{I} _{\{|\Delta(q,p,h_d/\lambda_i)-\mu(h_d/\lambda_i)| > \epsilon{\rm Var}_d^{\frac{1}{2}}\}} \Big),\] where now \((q,p\)) has density \(\propto\exp(-H(q,p))\), or, after rearrangement

\[\sum_{i=1}^d \sigma^2(h_d/\lambda_i) \mathbb{E}\Big( \Big(\frac{\Delta(q,p,h_d/\lambda_i)-\mu(h_d/\lambda_i)}{\sigma(h_d/\lambda_i)}\Big)^2 \mathbb{I} _{\{|\Delta(q,p,h_d/\lambda_i)-\mu(h_d/\lambda_i)| > \epsilon{\rm Var}_d^{\frac{1}{2}}\}} \Big).\] Since we know that \(\sum_i \sigma^2(h_d/\lambda_i)\) has a finite limit, the proof will be complete if we show that \[\lim_{d\rightarrow \infty} \max_{i=1,\dots,d} \mathbb{E}\Big( \Big(\frac{\Delta(q,p,h_d/\lambda_i)-\mu(h_d/\lambda_i)}{\sigma(h_d/\lambda_i)}\Big)^2 \mathbb{I} _{\{|\Delta(q,p,h_d/\lambda_i)-\mu(h_d/\lambda_i)| > \epsilon{\rm Var}_d^{\frac{1}{2}}\}} \Big)= 0\] or, recalling the uniform convergence in 34 , if we prove that \[\lim_{t\rightarrow 0} \mathbb{E}\Big( \Big(\frac{\Delta(q,p,t)-\mu(t)}{\sigma(t)}\Big)^2 \mathbb{I} _{\{|\Delta(q,p,t)-\mu(t)| > \epsilon{\rm Var}_d^{\frac{1}{2}}\}} \Big)= 0.\] In order to do so, apply the dominated convergence theorem, after noting that from Condition 1 and Theorem 1, \[\lim_{t\rightarrow 0} \frac{\Delta(q,p,t)-\mu(t)}{\sigma(t)} = \frac{\alpha(q,p)}{\sqrt{\Sigma}},\] and \[\lim_{t\rightarrow 0} \mathbb{I} _{\{|\Delta(q,p,t)-\mu(t)| > \epsilon{\rm Var}_d^{\frac{1}{2}}\}} = \mathbb{I} _{\{\frac{|\Delta(q,p,t)-\mu(t)|}{\sigma(t)} > \epsilon\frac{{\rm Var}_d^{\frac{1}{2}}\}}{\sigma(t)}} = 0,\] because \(\sigma(t)\rightarrow 0\). To dominate the integrand, observe that \(\Delta^2/\sigma^2\) behaves like \((\Delta/t^\nu)^2\) and is therefore upper bounded by the integrable function \(D^{2/3}\), where \(D\) is from Condition 2.

Proof of Theorem 3:

By integrating with respect to \(U\): \[{\cal SJD}_d = \frac{1}{d} \sum_{i=1}^d \mathbb{E}(\|q_i^\star-q_i\|^2 a_d(Q,P,h_d)).\] Then, we have \[d^{\gamma/\nu}{\cal SJD}_d = \ell^2\frac{1}{d} \sum_{i=1}^d \mathbb{E}\left(\|\frac{1}{h}(q_i^\star-q_i)\|^2 a_d(Q,P,h_d)\right),\] an expression that may be decomposed as the sum of \[\label{eq:aux1} \ell^2\frac{1}{d} \sum_{i=1}^d \mathbb{E}(\|p_i\|^2a_d(Q,P,h_d) ) = \ell^2 \mathbb{E}\Big(\frac{1}{d}\sum_{i=1}^d \|p_i\|^2a_d(Q,P,h_d)\Big)\tag{35}\] and a remainder \[\ell^2\frac{1}{d} \sum_{i=1}^d \mathbb{E}\left(\Big(\|\frac{1}{h}(q_i^\star-q_i)\|^2-\|p_i\|^2\Big) a_d(Q,P,h_d)\right).\] In 35 , we know from the proof of Theorem 2 and the continuous mapping theorem that \(a_d(Q,P,h)\) \(= \min(1,\exp(-\Delta_d))\) converges in distribution to \(\min(1,\textcolor{black}{\exp(-\Delta_\infty)})\). The random variable \((1/d)\sum_i \|p_i\|^2\) has a Gamma distribution with parameters \(\alpha= dm/2\), \(\theta = 2/d\). It has mean \(m\) and variance \(2m/d\) and therefore converges to the constant \(m\) (as one may alternatively have proved by the weak law of large numbers). By Slutsky’s theorem the variables \((1/d)\sum_i \|p_i\|^2a_d\) have a distributional limit and, since their second moments are bounded, we conclude that 35 converges to the limit \[m \ell^2 \lim_{d\rightarrow \infty} \mathbb{E}(a_d(Q,P,h_d)) = m\ell^2 A = 2 m \ell^2 \Phi\left(-\frac{\ell^{\nu}\sqrt{K\Sigma}}{2}\right).\]

The proof concludes by checking that the remainder converges to zero. In fact, from equivariance and Condition 3, \[\mathbb{E}\left(\|\frac{1}{h}(q_i^\star-q_i)\|^2-\|p_i\|^2 \right) = r(h_d/\lambda_i)\] and it is therefore enough to recall the uniform convergence in 34 .

6 Outlook: high order algorithms for better scaling↩︎

In Section 2 we presented three volume preserving, reversible integrators that generate MCMC algorithms. The associated values of \(\nu\) were 1, 3, 3. We now discuss alternative integrators with higher values of \(\nu\) and therefore better scaling properties.

Many symplectic, reversible integrators for Hamiltonian problems have been suggested in the literature [17]. Since they are automatically volume preserving, all of them may be used as described in Section 2. The class of reversible splitting integrators [18] is of particular interest. Those algorithms are simple, easy to implement generalizations of the Verlet scheme that use several evaluations of \(\nabla V\) to generate a single proposal. As an illustration consider as in [22] the family of splitting integrators \((q^\star,p^\star) = \psi_h(q,p)\) described by the equations \[\begin{align} p^{(1)} &=& p- h (1/2-b)\nabla V(q),\\ q^{(1)} &=& q+ h a p^{(1)},\\ p^{(2)} &=& p^{(1)}- h b\nabla V(q^{(1)}),\\ q^{(2)} &=& q^{(1)}+ h (1-2a)p^{(2)},\\ p^{(3)} &=& p^{(2)}- h b\nabla V(q^{(2)}),\\ q^\star &=& q^{(2)}+ h ap^{(3)},\\ p^\star &=& p^{(3)}- h (1/2-b)\nabla V(q^\star), \end{align}\] where \(a\) and \(b\) are real parameters. On acceptance, the first evaluation of \(\nabla V\) to be used for the next proposal coincides with the last evaluation \(\nabla V(q^\star)\) in the current step. On rejection, the first evaluation of \(\nabla V\) to generate the next proposal coincides with the first evaluation \(\nabla V(q)\) in the current step. Thus each step of the Markov chain, except the very first, uses three gradient evaluations.

For the choice \[a = \frac{1}{2-\sqrt[3]{2}},\qquad b = \frac{1-a}{2}\] sometimes attributed to Yoshida, the integrator is fourth order accurate with \(\nu= 5\). For the i.i.d.case with \(\lambda_i=1\), the scaling is \(h_d = \ell/ d^{1/10}\) or \(\delta_d = \ell^2/d^{1/5}\), a marked improvement on \(\delta_d = \ell^2/d\) for RWM and \(\delta_d = \ell^2/d^{1/3}\) for MALA. By writing \(q^\star\) as a function of \(q\) and \(p\) and assuming \(p\sim N(0,I_m)\), it is easy to check that the Yoshida algorithm provides an approximation to the Langevin dynamics of weak order exactly one: approximating to high order the Hamiltonian dynamics does not imply approximating to high order the Langevin dynamics. In addition, this also implies that higher order algorithms for the Langevin equations are not required for better scaling, thus answering negatively one of the conjectures in the discussion of [23] All other choices of \(a\) and \(b\) in the family above lead to integrators with \(\nu=3\) and some of them show very good performance in the context of HMC algorithms [24], [25]. It may be proved [18] that by increasing the number of gradient evaluations per proposal it is possible to have splitting integrators which attain arbitrarily high values of \(\nu\) and therefore arbitrarily better scaling properties of the sampler. (The relation between a target order \(\nu\) and the required number of gradient evaluations is known, but not simple.)

It is also possible to reach arbitrarily high values of \(\nu\) by means of implicit integrators that generalize the implicit midpoint rule used in MAIMLA. The Gauss method of \(s\) stages [17], \(s = 1,2,\dots\), achieves \(\nu=2s+1\). The implicit midpoint corresponds to \(s=1\). These integrators have the property of exactly conserving energy for quadratic potentials and therefore possess an acceptance rate of \(100\%\) for Gaussian distributions. When applied to the potential 1 in \(\mathbb{R}^m\) each proposal requires the solution of an \(s\times m\) dimensional system of nonlinear equations.

A comparison between the merits of the different proposals is beyond the scope of this paper. However we would like to illustrate the theory with a numerical example. We have considered as a target a product of \(d\) independent \(N(0,1)\) univariate distributions, \(d= 10^{3}, 10^{4}, 10^{5}\), and sampled, when the chain is started at stationarity, with MALA \(\delta_d = \ell^2/d^{1/3}\), fMALA \(\delta_d = \ell^2/d^{1/5}\) and the Yoshida algorithm presented above, \(\delta_d = \ell^2/d^{1/5}\). Recall that fMALA involves much more complexity than MALA and the splitting algorithm. For each sampler, simulations were carried out with a range of values of \(\ell\) and we recorded the empirical Squared Jumping Distance and acceptance probability. The data reported are averages over \(10^{4}\) realizations. Figure 1 displays the Squared Jumping Distance \(\times\: d^{1/5}\) as a function of acceptance rate. According to our analysis, the \(d^{1/5}\) scaling is appropriate for fMALA and Yoshida. We see that indeed for those two algorithms the (scaled) jumping distance does not degrade as \(d\) increases. However, the MALA results move downwards as \(d\) increases, because for this algorithm a scaling \(d^{1/3}\) would be required. The figure also displays, for each sampler, the relation in equation 26 between \(\mathcal{SJD}_d\) and the acceptance rate; since the prefactor \(2^{2/\nu}/(K\Sigma)^{1/\nu}\) was not known (but could have been computed as the target is Gaussian), we scaled the graphs by matching the maximum of the empirical results to equal the maximum of the corresponding theoretical curve. Finally, vertical lines indicate the theoretical values 0.574 and 0.704 of the optimal acceptance rate for methods with \(\nu = 3\) (MALA) and \(\nu = 5\) (fMALA and Yoshida). Clearly those values, proved to be valid in the limit \(d\rightarrow \infty\), agree very well with the simulations.

References↩︎