# Optimal design of the Barker proposal and other locally-balanced
Metropolis–Hastings algorithms

Jure<EMAIL_ADDRESS>Department of Statistics, University
of Warwick  Samuel Livingstone  Department of Statistical Science, University
College London  Giacomo Zanella  Department of Decision Sciences, BIDSA &
IGIER, Bocconi University

###### Abstract

We study the class of first-order locally-balanced Metropolis–Hastings
algorithms introduced in [9]. To choose a specific algorithm within the class
the user must select a balancing function $g:\mathbb{R}\to\mathbb{R}$
satisfying $g(t)=tg(1/t)$, and a noise distribution for the proposal
increment. Popular choices within the class are the Metropolis-adjusted
Langevin algorithm and the recently introduced Barker proposal. We first
establish a universal limiting optimal acceptance rate of 57% and scaling of
$n^{-1/3}$ as the dimension $n$ tends to infinity among all members of the
class under mild smoothness assumptions on $g$ and when the target
distribution for the algorithm is of the product form. In particular we obtain
an explicit expression for the asymptotic efficiency of an arbitrary algorithm
in the class, as measured by expected squared jumping distance. We then
consider how to optimise this expression under various constraints. We derive
an optimal choice of noise distribution for the Barker proposal, optimal
choice of balancing function under a Gaussian noise distribution, and optimal
choice of first-order locally-balanced algorithm among the entire class, which
turns out to depend on the specific target distribution. Numerical simulations
confirm our theoretical findings and in particular show that a bi-modal choice
of noise distribution in the Barker proposal gives rise to a practical
algorithm that is consistently more efficient than the original Gaussian
version.

##  1 Introduction

Markov chain Monte Carlo algorithms are the workhorse of many contemporary
statistical analyses, and an essential part of the modern data science
toolkit. Despite many advances, however, reliable inference using Markov chain
Monte Carlo can still be a cumbersome task. It is common for practitioners to
dedicate much effort to making careful algorithm design choices and adjusting
algorithmic tuning parameters to ensure that performance is adequate for a
given problem. Failure to do this can be catastrophic; examples for which a
well-designed algorithm performs adequately but a less carefully-chosen
alternative does not are ubiquitous (e.g. [19]).

Suitable guidelines on the intelligent design and implementation of Markov
chain Monte Carlo methods are therefore important. They are not always easy to
offer, however, the best choice of method can depend on the user and the
problem at hand. In some contexts, a simpler algorithm with less need for
adjustment and for which potential problems are easy to diagnose may be
preferable. In others contexts, one may be comfortable with more complex
methods, which can perform adequately on a larger class of problems if enough
fine tuning is done.

For Metropolis–Hastings algorithms, perhaps the most celebrated guidelines
concern the choice of optimal acceptance rate [15]. Rigorous theoretical
justification for certain values tend to be restricted to the case in which
dimension tends to infinity and the distribution from which samples are
desired has a particular structure (such as a product form), but empirically
the same values are known to be appropriate in many other settings [15]. The
apparent lack of dependence of these optimal choices on the target
distribution allows particularly simple recommendations to be offered to the
user for a given algorithm.

Adaptive Markov chain Monte Carlo methods have also facilitated efficient
implementation [1, 16]. Users can implement adaptive algorithms in which
algorithmic tuning parameters are automatically adjusted towards guideline
values, using ideas from stochastic optimisation and controlled Markov chains.
When combined with appropriate theory, adaptive algorithms can therefore allow
users to implement their chosen method on a given problem without the need for
cumbersome hand-tuning. Such innovations have made it possible to develop
popular tailored software packages for users of Markov chain Monte Carlo [18,
3].

Not all adaptive algorithms are created equally, however. Empirically it has
long been observed that certain approaches are more sensitive to tuning than
others [11]. In recent work [9] provided some theoretical justification for
this phenomenon, in particular showing that for popular gradient-based
approaches such as the Metropolis-adjusted Langevin algorithm and Hamiltonian
Monte Carlo spectral gaps decay exponentially quickly to zero as the tuning
parameters are perturbed from their optimal values. By contrast, spectral gaps
for the simpler random walk Metropolis decay at a polynomial rate, indicating
that the algorithm is much more robust to tuning. This has a compounding
effect if the tuning parameters are learned adaptively, as adaptive algorithms
typically learn based on past samples from the Markov chain, and if these past
samples are very poor as a result of the initial tuning parameters being sub-
optimal it can mean that the learning occurs very slowly. The moral of the
story is that algorithms can still perform poorly in practice even if an
optimally-tuned version would in theory perform well.

These findings present a conundrum, as gradient-based algorithms are
considered the state-of-the-art in Markov chain Monte Carlo to sample from
continuous and smooth distributions when properly tuned. To explore the
phenomenon in more detail [9] introduce a general class of gradient-based
algorithms, termed _first order locally-balanced_ Metropolis–Hastings, of
which the Metropolis-adjusted Langevin algorithm is a special case.
Constructing a member of the class requires a Markov kernel, which can be
thought of as the initial noise distribution for the transition, together with
a _balancing function_ , which must satisfy certain properties (see Section
2). The authors consider different choices from within the class, and in
particular construct a method called _the Barker proposal_. This algorithm has
spectral gaps that are robust to tuning as in the random walk Metropolis. The
authors also establish sufficient conditions for geometric ergodicity and some
preliminary results on scaling with dimension, suggesting that relaxation
times are $\mathcal{O}(n^{1/3})$, where $n$ is the dimension of the state.
Empirical results in the paper show that the Barker algorithm pairs extremely
well with adaptive learning of tuning parameters, and enables reliable
sampling on complex examples in which other gradient-based methods may not,
despite being remarkably simple to implement. More discussion and a
pedagogical derivation of the Barker algorithm is provided in [7].

Several unexplored questions remain regarding locally-balanced
Metropolis–Hastings algorithms. The initial noise distribution in the Barker
algorithm is simply chosen to be Gaussian in [9], but no justification besides
convenience is given for this choice. It could be that a different choice
leads to a more effective algorithm. Similarly, general guidelines on the
optimal acceptance rate for the Barker algorithm are not established. More
generally, little discussion is provided on other first-order locally-balanced
Metropolis–Hastings. It is natural to wonder whether all members of the class
will exhibit $\mathcal{O}(n^{1/3})$ relaxation times, if the
Metropolis–adjusted Langevin is the most efficient choice when optimally
tuned, and indeed whether such a direct quantitative comparison of methods is
possible in general. These questions are of both theoretical and practical
interest, as they have direct implications for the optimal design of
algorithms.

In this paper we make several new contributions. First we present universal
results on the optimal choice of acceptance rate and scaling with dimension of
any algorithm within the class of first order locally-balanced Markov
processes (under mild regularity conditions). In particular in Section 3 we
show that the 57% guideline acceptance rate for the Metropolis-adjusted
Langevin algorithm also holds for the Barker proposal and several other
methods, as does the $\mathcal{O}(n^{1/3})$ scaling with dimension as measured
by expected squared jump distance. Despite having the same optimal acceptance
rate and scaling with dimensionality, however, all such schemes have a
different asymptotic efficiency, which we explicitly characterize, enabling
for principled and generic optimization of the algorithmic design. We first
consider optimal design of the Barker proposal, in particular with respect to
the noise distribution, which is chosen to be Gaussian in [9]. We find, both
theoretically and empirically, that it is in fact beneficial to choose a bi-
modal noise distribution for each coordinate, and offer some discussion and
practical user guidelines in Section 4.1. We then consider the case in which
the noise distribution is fixed and an optimal balancing function is chosen in
Section 4.2, and the general scenario in which both the noise distribution and
balancing function are optimized over in Section 4.3. Both cases yield
surprising results, such as optimality being reached by having positive
probability of keeping some coordinates fixed at each iteration. We conduct
numerical experiments to verify the theory in Section 5, and provide a
discussion in Section 6. Our theoretical results build on the recently
introduced optimal scaling framework of [21, 24]. One powerful aspect of this
approach is the ability to analyze fairly generic schemes without requiring
overly case-specific calculations (e.g. those related to proposal
distributions with Gaussian noise, linear drift, etc.), while still obtaining
explicit expressions for the asymptotic performances that can directly be
compared among algorithms. This allows characterization of the quantitative
interplay between fine-scale properties of the target and proposal
distributions (such as moments of the noise, aggressiveness of the balancing
function and derivatives of the target) in the resulting asymptotic
efficiency, thus enabling precise methodological guidance.

##  2 Locally-balanced Markov processes

###  2.1 General framework

Consider a Markov transition kernel $Q$ defined on a Borel space
$(\mathbb{X},\mathcal{F})$. We restrict attention to
$\mathbb{X}\subset\mathbb{R}^{n}$ for some finite $n$. We say $Q$ satisfies
the detailed balance equations with respect to a probability measure $\pi$ if

$\int f(x)h(y)\pi(dx)Q(y,dx)=\int f(x)h(y)\pi(dy)Q(y,dx)$ (1)

for any $f,h\in L^{2}(\pi)$. When $Q$ does not satisfy (1), a new kernel can
be constructed using the concept of a balancing function. Let
$g:[0,\infty)\to[0,\infty)$ be such that $g(0)=0$ and for $t>0$

$g(t)=tg(1/t),$ (2)

and note that by Proposition 1 in [20] there exists a symmetric set
$\mathcal{R}\times\mathcal{R}\in\mathbb{X}\times\mathbb{X}$ such that the
Radon–Nikodym derivative

$t(x,y)=\frac{\pi(dy)Q(y,dx)}{\pi(dx)Q(x,dy)}$ (3)

is well-defined and such that $0<t(x,y)<\infty$ if $x,y\in\mathcal{R}$ and
$t(x,y)=0$ otherwise. Then the kernel

$\tilde{\mathcal{P}}(x,dy)=g\left\\{\frac{\pi(dy)Q(y,dx)}{\pi(dx)Q(x,dy)}\right\\}Q(x,dy)$
(4)

satisfies (1). However, the kernel $\tilde{\mathcal{P}}$ is not necessarily
Markov. One way of enforcing that (4) integrates to one is to restrict
attention to $g\leq 1$, ensuring that $\tilde{\mathcal{P}}(x,\mathbb{X})\leq
1$, and then combine with
$r(x,dy)=\\{1-\tilde{\mathcal{P}}(x,\mathbb{X})\\}\delta_{x}(dy)$, where
$\delta_{x}(A)=1$ if $x\in A$ and $0$ otherwise. The resulting kernel
$\tilde{\mathcal{P}}(x,dy)+r(x,dy)$ is Metropolis–Hastings (e.g. [20]).

An alternative strategy introduced in [23, 12, 9] is to instead allow any $g$
for which $\mathcal{Z}(x)=\tilde{\mathcal{P}}(x,\mathbb{X})$ is finite, and
then set

$\mathcal{P}(x,dy)=\frac{\tilde{\mathcal{P}}(x,dy)}{\mathcal{Z}(x)}.$ (5)

Note that $\mathcal{P}$ does not satisfy (1) in general, in fact $\mathcal{P}$
is invariant with respect to the measure $\mathcal{Z}(x)\pi(dx)$. A
$\pi$-invariant Markov jump process can be constructed, however, by
introducing a holding time $\mathcal{Z}(x)$ at each state $x$, and then
choosing the next state according to $\mathcal{P}$. This construction is
called a _locally-balanced Markov process_ (see [12, 7] for more detail).

###  2.2 First order locally-balanced processes

The function $\mathcal{Z}(x)$ will not be tractable in general, meaning
further work is needed to design a sampling algorithm based on a locally-
balanced Markov process. One approach is to restrict attention to symmetric
$Q$, and $\pi$ absolutely continuous with respect to the Lebesgue measure on
$\mathbb{R}^{n}$ with differentiable Lebesgue density $\pi(x)$. In this case
(3) reduces to $\pi(y)/\pi(x)$. From this point several natural first order
approximations of this ratio can be taken to construct a new more tractable
kernel. It is argued in [7, 9] that a good choice is the component-wise
approximation found by letting
$Q(x,dy)=\prod_{i}\sigma^{-1}\mu\\{(dy_{i}-x_{i})/\sigma\\}$, where $\mu$ is a
centred and symmetric distribution on $\mathbb{R}$ and $\sigma>0$, and setting

$\tilde{P}(x,dy)=\prod_{i=1}^{n}g\left(e^{(y_{i}-x_{i})\partial_{i}\log\pi(x)}\right)\mu\left(\frac{dy_{i}-x}{\sigma}\right),$

where $\partial_{i}=\partial/\partial x_{i}$ and for any event $A$ the set
$(A-x_{i})/\sigma=\\{z\in\mathbb{R}:x_{i}+\sigma z\in A\\}$, and its Markovian
counterpart

$P(x,dy)=\frac{\tilde{P}(x,dy)}{Z(x)}$ (6)

where $Z(x)=\tilde{P}(x,\mathbb{X})$. With this approximation certain choices
of $g$ and each $Q_{i}$ lead to familiar forms of $P$. Choosing $g(t)=\surd t$
and $\mu$ to be standard Gaussian, for example, leads to the unadjusted
Langevin algorithm [17]. The class of kernels obtained by (6) is much broader,
however, and is currently relatively unexplored.

###  2.3 The choice of balancing function

[9] suggest the choice of balancing function $g(t)=2t/(1+t)$, as popularised
by [2] in the context of Metropolis–Hastings. With this choice a sample from
$P$ can be easily drawn in the following manner. First sample $z_{i}\sim\mu$
for each $i$, then set $\beta_{x,i}=\partial_{i}\log\pi(x)$ and flip the sign
of each $z_{i}$ with probability $F(\beta_{x,i}z_{i})$, where
$F(x)=e^{x}/(1+e^{x})$. Finally add this to the current coordinate $x_{i}$. To
construct a $\pi$-invariant Markov chain a Metropolis–Hastings correction is
then applied to this Barker proposal. See Algorithm 1 for more detail.

Algorithm 1 Simulate from the Barker proposal. Require: current point
$x\in\mathbb{X}$
---
For $i=1$ to $n$
Draw $z_{i}\sim\mu$, and set $\beta_{x,i}\leftarrow\partial_{i}\log\pi(x)$
Set $y_{i}\leftarrow x_{i}+z_{i}$ with probability $F(\beta_{x,i}z_{i})$, and
$y_{i}\leftarrow x_{i}-z_{i}$ otherwise
Output $y=(y_{1},...,y_{d})$

It is natural to wonder how many choices of $g$ can be made. Two other simple
possibilities are $\min(1,t)$ and $\max(1,t)$, the latter being recently
studied in [4]. The below results show that in fact the family of balancing
functions is infinitely large.

######  Proposition 1.

Let $\mathcal{H}=\\{h:\mathbb{R}\to[0,\infty)~{}:~{}h(x)=h(-x)\\}$ be the
space of positive even functions. Then for every $h\in\mathcal{H}$,
$g_{h}(t)=t^{1/2}h(\log t)$ is a balancing function. Conversely, for every
balancing function $g$, the function $h_{g}(x)=e^{-x/2}g(e^{x})$ is contained
in $\mathcal{H}$.

The above provides an explicit parametrisation of $g_{h}$ in terms of a
specific $h\in\mathcal{H}$. The function $t^{1/2}$ can also be replaced with
any other balancing function to give a different bijection. The goal, of
course, is to find choices of $g$ for which tractable sampling algorithms can
be designed. In Section 4 we design new balancing functions of this nature for
specific objectives.

##  3 A universality result on the optimal acceptance rate and scaling with
dimension

###  3.1 Preliminaries

The concept of a log-Metropolis–Hastings random variable will be crucial for
our analysis of optimal scaling. We recall some key results here, for more
detail see Section 3 of [21].

######  Definition 1.

For a probability measure $\pi$, Markov kernel $Q$ on
$(\mathbb{X},\mathcal{F})$ and $\mathcal{R}$ as in Section 2.1, let $X\sim\pi$
and $Y\sim Q(X,\cdot)$. The associated log-Metropolis-Hastings-ratio random
variable is

$\rho(X,Y)=\begin{cases}~{}\log\left\\{\frac{\pi(dy)Q(y,dx)}{\pi(dx)Q(x,dy)}\right\\}&\text{if
}(X,Y)\in\mathcal{R}\times\mathcal{R}\,,\\\ ~{}0&\text{otherwise.}\end{cases}$

Let $\pi:\mathbb{R}\to[0,\infty)$ be a probability density on $\mathbb{R}$ and
for any fixed $\sigma>0$ let $Q_{\sigma}:\mathbb{R}\times\mathbb{R}\to[0,1]$
be a Markov kernel. We introduce the product measure
$\pi_{n}(dx)=\prod_{i=1}^{n}\pi(x_{i})dx_{i}$ on $\mathbb{R}^{n}$ and the
product kernel
$\mathcal{Q}_{n}(x,dy)=\prod_{i=1}^{n}Q_{\sigma_{n}}(x_{i},dy_{i})$, where
$(\sigma_{n})_{n\in\mathbb{N}}$ is a sequence of positive real numbers. The
associated log-Metropolis–Hastings random variable is

$\rho(X_{n},Y_{n})=\sum_{i=1}^{n}\rho_{n}(X_{n,i},Y_{n,i}),$

where $X_{n}=(X_{n,1},...,X_{n,n})\sim\pi_{n}$,
$Y_{n}\sim\mathcal{Q}_{n}(X_{n},\cdot)$, and $\rho_{n}$ is the log-
Metropolis–Hastings random variable associated with $f$ and $Q_{\sigma_{n}}$.
The following is established in [21].

######  Theorem 1.

Assume that there exists a positive sequence $(a_{n})_{n\in\mathbb{N}}$ with
$\lim_{n\to\infty}a_{n}=0$ such that

$\lim_{n\to\infty}\mathbb{E}[\rho_{n}^{2}1_{\rho_{n}<-a_{n}}]/\mathbb{E}[\rho_{n}^{2}]=0,$

where $1$ denotes the indicator function. If in addition
$(\sigma_{n})_{n\in\mathbb{N}}$ is chosen such that
$\lim_{n\to\infty}n\mathbb{E}[\rho_{n}^{2}]=v^{2}$ for some constant $v>0$,
then as $n\to\infty$

$\sum_{i=1}^{n}\rho_{n}(X_{n,i},Y_{n,i})\Rightarrow
N\left(-\frac{1}{2}v^{2},v^{2}\right)\,.$ (7)

######  Remark 1.

The expectation $n\mathbb{E}[-\rho_{n}]$ denotes the Kullback–Leibler
divergence between the forward and reverse Markov transition kernels,
$\pi(dx)Q(x,dy)$ and $\pi(dy)Q(y,dx)$. In fact, it is further shown in [21]
that under the above assumptions
$\lim_{n\to\infty}\mathbb{E}[\rho_{n}]/\mathbb{E}[\rho_{n}^{2}]=-1/2$, meaning
that both the mean and variance on the right-hand side of (7) can be
interpreted in terms of the Kullback–Leibler divergence in the limiting case
$n\to\infty$.

Guaranteeing the first condition, i.e. understanding how fast
$\mathbb{E}[\rho_{n}^{2}]$ decays, is key for identifying the optimal scaling
of a Metropolis-Hastings algorithm. The other condition is technical and
related to the uniform integrability of $n\rho_{n}^{2}$ and to the conditions
required in the Lindeberg’s version of Central Limit Theorem (see Theorem 4.15
of [8]). It suffices for example, to show that $\rho_{n}$ has higher moments
that vanish faster than $1/n$.

###  3.2 The asymptotic acceptance rate for locally-balanced proposals

We will establish that the above central limit theorem holds for first order
locally-balanced Metropolis–Hastings under Assumption 1 below, and then
consider optimal acceptance rates and dimension dependence in terms of the
expected squared jump distance in each coordinate. We restrict attention to
the class of target distributions
$\pi_{n}(x)=\prod_{i=1}^{n}\exp\\{\phi(x_{i})\\}$, for some
$\phi:\mathbb{R}\to\mathbb{R}$, and impose regularity conditions on $\phi$
below. Let $b(x)=\log\\{g(e^{x})\\}$, and without loss of generality set
$g(1)=1$.

######  Assumption 1.

There exist constants $H\in(0,1)$, $\gamma>0$, $\beta\geq 0$, $\epsilon>0$
such that

  1. (i)

$\phi\in\mathcal{C}^{3+H}(\mathbb{R})$ and for
$f=\phi^{\prime\prime\prime},\phi^{\prime\prime}\phi^{\prime},\phi^{\prime
3},\phi^{\prime\prime}|\phi^{\prime}|^{1+\beta}$ the integrability condition
$\int_{\mathbb{R}}f(x)^{2+\epsilon}(1+|\phi^{\prime}(x)|^{\beta})\pi(x)dx<\infty$
as well as the mixed growth-Hölder condition

$\left|f(x+\delta)-f(x)\right|\leq K(x)\max(|\delta|^{H},|\delta|^{\gamma})$

are satisfied. Function $K$ is such that
$\int_{\mathbb{R}}K(x)^{2}(1+|\phi^{\prime}(x)|^{\beta})\pi(x)dx<\infty$.

  2. (ii)

$b\in\mathcal{C}^{3}(\mathbb{R})$ and
$b^{\prime},b^{\prime\prime},b^{\prime\prime\prime}$ are all bounded above.

  3. (iii)

$\int z^{2}\mu(dz)=1$, $\int|z|^{\xi}\mu(dz)<\infty$ for
$\xi=\max(6+3\epsilon,2+2H,2+2\gamma)$, and for all $a\in\mathbb{R}$ and some
positive $C_{\mu}>0$

$\int_{\mathbb{R}}e^{b(az)}|z|^{\xi}\mu(dz)\leq
C_{\mu}(1+|a|^{\beta})\int_{\mathbb{R}}e^{b(az)}\mu(dz)<\infty.$

Part (i) of the above refers to the target distribution, part (ii) to the
balancing function and part (iii) to the interplay between them. Part (i) is
straightforwardly satisfied for many statistical models of interest, for
example likelihoods from exponential families and suitably smooth priors. Part
(ii) is satisfied by all cases explicitly studied in the paper, such as
$g(t)=\surd t$ and $g(t)=2t/(1+t)$. Part (iii) highlights the need to control
the growth of $g$ and $b$ using the tails of $\mu$. If $g$ is bounded, as in
the Barker case, then any $\mu$ with a moment generating function is
sufficient for it to be satisfied (and for many targets actually much weaker
conditions are required). When $g(t)=\surd t$, which is not bounded above,
then stronger conditions on the tails of $\mu$ are needed, such as Gaussian
tails.

Part (i) is explicitly weaker than the typical assumptions made in the optimal
scaling literature (e.g. [14]). A form of part (i) as well as $\int
z^{6}\mu(dz)<\infty$ and $g\in\mathcal{C}^{3}$ are crucial to the analysis.
Part (iii) imposes uniform control (with respect to $x$) of measures
$e^{b(\sigma_{n}z\phi^{\prime}(x))}(Z_{\sigma_{n}}(x))^{-1}\mu(z)dz$ in terms
of only the measure $\mu(z)dz$. This is required so that the normalising
constants $Z_{\sigma_{n}}$ and their second derivatives are well defined. It
may be possible to significantly relax parts (ii) or (iii), especially in
specific settings, at the expense of strengthening elsewhere. The following
Proposition identifies some simple cases in which part (iii) is satisfied.

######  Proposition 2.

Part (iii) of Assumption 1 is satisfied in the following cases.

  1. (i)

If $\mu$ has a density with compact support, for any $g$.

  2. (ii)

If $g$ is bounded, non-decreasing and $\int|z|^{\xi}\mu(dz)<\infty$ for $\xi$
as in Assumption. 1

  3. (iii)

If $g$ satisfies part (ii) of Assumption 1 and there exists
$\tilde{C}_{\mu},\tilde{\beta}>0$ such that for all $a\in\mathbb{R}$

$\int_{\mathbb{R}}e^{az}|z|^{\xi}\mu(dz)\leq\tilde{C}_{\mu}(1+|a|^{\tilde{\beta}})\int_{\mathbb{R}}e^{az}\mu(dz)<\infty.$

  4. (iv)

If $g$ satisfies part (ii) of Assumption 1, $\mu$ has a density
$\mu\in\mathcal{C}^{1}(\mathbb{R})$ such that
$\lim_{z\to\pm\infty}e^{az}\mu(z)=0$ for any $a\in\mathbb{R}$ and there exists
constants $p>1$, $A,B>0$ for which

$|z|^{p}\mu(z)\leq A\mu(z)-Bz\mu^{\prime}(z).$

In specific examples we typically verify (i), (ii) or (iv) of Proposition 2.
For instance, choices of the form $\mu(dz)\propto e^{-|z|^{p}}dz$ for $p\geq
1$ satisfy (iv). Note that a statement analogous to (ii) but for the function
$b$ is not valid. Even if $g$ is bounded, $b$ is only bounded from infinity
above, not below. In fact, since $b(x)=x+b(-x)$ holds, $b$ can never be
bounded. These conditions are required to analyze Taylor series remainder
terms for the normalising constant. It is apparent from Proposition 2 that
less conditions on $\mu$ must be assumed for the Barker proposal, for which
$g$ is bounded, compared to the Langevin choice $g(t)=\surd t$.

######  Theorem 2.

Under Assumption 1 it holds that
$\lim_{n\to\infty}\sigma_{n}^{-6}\mathbb{E}[\rho_{n}^{2}]=\theta^{2}$ for some
$\theta\in[0,\infty)$ In addition, if $\theta>0$ and $\sigma_{n}$ is chosen
such that $\lim_{n\to\infty}n^{1/6}\sigma_{n}=\ell$, then

$\sum_{i=1}^{n}\rho_{n}(X_{n,i},Y_{n,i})\Rightarrow
N\left(-\frac{1}{2}\ell^{6}\theta^{2},\ell^{6}\theta^{2}\right).$ (8)

Denoting $\mathfrak{g}=g^{\prime\prime}(1)$,
$\mu_{4}=\int_{\mathbb{R}}z^{4}\mu(dz)$,
$\mu_{6}=\int_{\mathbb{R}}z^{6}\mu(dz)$ and
$A_{\phi}=\mathbb{E}_{\pi}[(\phi^{\prime\prime\prime})^{2}]$,
$B_{\phi}=\mathbb{E}_{\pi}[(\phi^{\prime}\phi^{\prime\prime})^{2}]$,
$C_{\phi}=\mathbb{E}_{\pi}[\phi^{\prime}\phi^{\prime\prime}\phi^{\prime\prime\prime}]$
the constant $\theta^{2}$ takes the form

$\displaystyle\theta^{2}=~{}$
$\displaystyle\mu_{6}\left\\{\frac{1}{144}A_{\phi}+\left(\frac{1}{4}+\mathfrak{g}\right)^{2}B_{\phi}\right.-\left.\frac{1}{6}\left(\frac{1}{4}+\mathfrak{g}\right)C_{\phi}\right\\}$
(9)
$\displaystyle+\mu_{4}\left\\{\frac{1}{6}\left(\frac{1}{2}+\mathfrak{g}\right)C_{\phi}\right.-\left.2\left(\frac{1}{4}+\mathfrak{g}\right)\left(\frac{1}{2}+\mathfrak{g}\right)B_{\phi}\right\\}+\left(\frac{1}{2}+\mathfrak{g}\right)^{2}B_{\phi}.$

Note that the specific choice of the scaling parameter $\sigma_{n}\propto
n^{-1/6}$ in Theorem 2 is the only rate leading to a non-trivial
distributional limit for $\sum_{i=1}^{n}\rho_{n}(X_{n,i},Y_{n,i})$, despite
the fact that
$\lim_{n\to\infty}\sigma_{n}^{-6}\mathbb{E}[\rho_{n}^{2}]=\theta^{2}$ holds
for any decay rate. Note also that the expression for $\theta^{2}$ depends on
both the balancing function $g$ and the distribution $\mu$. In Section 4 we
consider optimal ways to choose $g$ and $\mu$ for certain purposes. We
consider some example choices below.

######  Example 1.

In the Langevin case $g(t)=\surd t$ and $\mu$ is standard Gaussian, so that
$g^{\prime\prime}(1)=-1/4$ and $\mu_{4}=3$, $\mu_{6}=15$. Then

$\theta^{2}=\frac{5}{48}A_{\phi}+\frac{1}{8}C_{\phi}+\frac{1}{16}B_{\phi}$

which if
$\lim_{x\to\pm\infty}e^{\phi(x)}\phi^{\prime}(x)\phi^{\prime\prime}(x)^{2}=0$
can also be written (using integration by parts)

$\theta^{2}=\frac{5}{48}\mathbb{E}\left[(\phi^{\prime\prime\prime})^{2}\right]-\frac{1}{16}\mathbb{E}\left[(\phi^{\prime\prime})^{3}\right],$

a formula that appears in [14].

######  Example 2.

For the Barker proposal $g(t)=2t/(1+t)$ and $\mu$ can be any centred and
symmetric distribution such that $\int z^{6}\mu(dz)<\infty$. With these
choices $g^{\prime\prime}(1)=-1/2$ and

$\theta^{2}=\frac{\mu_{6}}{144}\left(A_{\phi}+6C_{\phi}+9B_{\phi}\right).$
(10)

An important consequence of Theorem 2, and in particular of (8), is a simple
expression for the asymptotic acceptance rate for a first order locally-
balanced Metropolis–Hastings algorithm (see e.g. Proposition 2.4 in [13]).

######  Corollary 1.

Setting $\alpha_{n}(X,Y)=\min\\{1,\sum_{i=1}^{n}\rho_{n}(X_{i},Y_{i})\\}$,
under the conditions of Theorem 2

$\lim_{n\to\infty}\mathbb{E}[\alpha_{n}]=2\Phi(-\ell^{3}\theta/2)$

where $\Phi$ is the standard Normal cumulative distribution function.

###  3.3 Optimal acceptance rates

Given the simplified limiting expression for $\alpha_{n}$ in Corollary 1, we
can consider optimal choices of the constant $\ell$ for a fixed $\theta$,
leading to an optimal acceptance rate. We consider optimising the expected
squared jump distance here, which is well-studied and has a strong
justification motivated by diffusion limits in various settings [15].

Using the same notation as above denote by
$(\mathcal{E}^{g,\mu}_{n})_{n\in\mathbb{N}}$ the sequence of expected squared
jump distances for the first (or any other) coordinate, defined as

$\mathcal{E}^{g,\mu}_{n}=\mathbb{E}\left[(Y_{n,1}-X_{n,1})^{2}\alpha\left(X_{n},Y_{n}\right)\right],$

where $X_{n}\sim\pi_{n}$ and $Y_{n}$ is generated from $X_{n}$ using a first
order locally-balanced proposal, defined in (6), with distribution $\mu$,
balancing function $g$ and variance parameter $\sigma_{n}$. We have the
following result.

######  Theorem 3.

Let Assumption 1 and Theorem 2 be satisfied for $\phi$, $\mu$ and $g$ and
$\theta>0$. Let $(\sigma_{n})_{n\in\mathbb{N}}$ be a positive sequence with
$\lim_{n\to\infty}\sigma_{n}=0$. If either
$\lim_{n\to\infty}n^{1/6}\sigma_{n}=0$ or
$\lim_{n\to\infty}n^{1/6}\sigma_{n}=\infty$ then as $n\to\infty$

$n^{1/3}\mathcal{E}^{g,\mu}_{n}\to 0.$

If $\lim_{n\to\infty}n^{1/6}\sigma_{n}=\ell$ for some $\ell\in(0,\infty)$,
then as $n\to\infty$

$n^{1/3}\mathcal{E}^{g,\mu}_{n}\to
h(\ell)=2\ell^{2}\Phi(-\ell^{3}\theta/2)\,,$

where $\Phi$ is the standard Normal cumulative distribution function on
$\mathbb{R}$. Furthermore, there exists a unique optimal
$\ell^{*}(=\ell^{*}(g,\mu))$ that maximizes $h(\ell)$, for which
$2\Phi\\{-(\ell^{*})^{3}\theta/2\\}\approx 0.574$. The corresponding optimal
asymptotic efficiency satisfies

$h(\ell^{*})=C_{h}\theta^{-2/3}\,,$

where $C_{h}\approx 0.652$.

The above shows that any first order locally-balanced Metropolis–Hastings
algorithm will have the same asymptotic optimal acceptance rate of 0.57, and
that algorithmic efficiency as measured by expected squared jump distance will
scale as $\mathcal{O}(n^{-1/3})$ for $n\to\infty$. This includes both Barker
and Langevin proposals as well as many other possibilities. Theorem 3 also
suggests a route to both comparison and optimal design of first order locally-
balanced Metropolis–Hastings algorithms, in the former case by comparing
$\theta^{2}$ for different choices of $\mu$ and $g$, and in the latter by
choosing $\mu$ and $g$ so that $\theta^{2}$ in Theorem 2 is minimized.
According to the same theorem, under Assumption 1 the constant $\theta^{2}$
will depend on $\phi$ through $A_{\phi},B_{\phi}$ and $C_{\phi}$, on $\mu$
only through $\mu_{4}$ and $\mu_{6}$ and on $g$ only through
$\mathfrak{g}=g^{\prime\prime}(1)$. We explore optimal design under different
constraints in the next Section.

In the Langevin proposal case, the constant $h(\ell)$ was shown to correspond
to the speed measure of an overdamped Langevin diffusion limit in [14]. We
conjecture that the same is true for locally balanced proposals in general,
but do not prove explicitly diffusion limit results in this paper. Proving
diffusion limit results for general locally balanced proposals is a non-
trivial open problem, as it would require a conditional version of the central
limit theorem in (8) that is hard to obtain in such generality.

######  Example 3.

Take the Gaussian target case $\phi(x)=-x^{2}/2$. Then $\phi^{\prime}(x)=-x$,
$\phi^{\prime\prime}(x)=-1$ and $\phi^{\prime\prime\prime}(x)=0$, meaning
$A_{\phi}=C_{\phi}=0$ and $B_{\phi}=\mathbb{E}[x^{2}]=1$. For Langevin
proposals with $g(t)=\surd t$ and $\mu$ taken as Gaussian, the constant
$\theta^{2}$ in (9) becomes $\theta^{2}_{L}=1/16$, whereas for the Barker
choice $g(t)=2t/(1+t)$ and the same $\mu$ we have $\theta_{B}^{2}=\mu_{6}/16$.
The ratio of asymptotic expected squared jump distances is therefore
$(\theta_{B}/\theta_{L})^{2/3}=\mu_{6}^{1/3}$. Here $\mu_{6}=15$ meaning that
Langevin proposals are asymptotically $15^{1/3}\approx 2.47$ times more
efficient than Barker proposals with Gaussian noise when optimally tuned. This
is consistent with experiments in Section 5.2 of [9].

######  Example 4.

Consider hyperbolic targets, $\phi(x)=(\delta^{2}+x^{2})^{1/2}$, with
$\delta^{2}=0.1$ as in [9]. Then $A_{\phi}\approx 12.99$, $B_{\phi}\approx
0.22$ and $C_{\phi}\approx 1.68$. The same calculations as above imply that
Langevin proposals are 1.18 times more efficient than Barker proposals with
Gaussian noise when optimally tuned, which is also consistent with Section 5.2
of [9].

##  4 Optimal choices among the class of locally-balanced algorithms

###  4.1 Optimal choice of noise in the Barker algorithm

In this setting we fix $g(t)=2t/(1+t)$ and minimize $\theta^{2}$ with respect
to $\mu$, for a given but arbitrary choice of $\phi$. In this case
$\theta^{2}$ is given by (10), and the only influence of $\mu$ comes from the
sixth moment $\mu_{6}$. The asymptotic expected squared jump distance can
therefore be straightforwardly maximised by minimising the sixth moment of
$\mu$ subject to the constraint that $\mu_{2}=1$. Note that by Jensen’s
inequality $\mu_{6}\geq\mu_{2}^{3}=1$, and in fact the lower bound is uniquely
attained by choosing $\mu$ to be a Rademacher distribution, such that if
$W\sim\mu$ then $W=1$ with probability $1/2$ and $W=-1$ otherwise. We state
this formally below.

######  Proposition 3.

If $g(t)=2t/(1+t)$ then $\theta^{2}$ is minimized when $W\sim\mu$ is chosen to
take values $+1$ and $-1$ each with probability $1/2$.

We can compare the relative efficiency of Barker with Rademacher versus
Gaussian noise using (10) in a similar manner to Examples 3 and 4. Doing this
shows that for any $\phi$ the Rademacher version will be $\mu_{6}^{1/3}\approx
2.47$ times more efficient than the Gaussian version. It is particularly
convenient that the optimal choice of $\mu$ does not depend in any way on
$\phi$ and therefore generic methodological guidance can be provided for the
algorithm. The comparison with the Langevin proposal is instead target
dependent, as exemplified below.

######  Example 5.

When $\phi(x)=-x^{2}/2$ as in Example 3, the Barker proposal with Rademacher
noise will be exactly as efficient as the Langevin proposal. When
$\phi(x)=(\delta^{2}+x^{2})^{1/2}$ with $\delta^{2}=0.1$ as in Example 4 then
the Rademacher proposal will be 2.08 times more efficient than the Langevin
proposal.

We compare these theoretical results with empirical performances in Section 5.
The Rademacher version of the Barker proposal is not per se a practical
sampling algorithm given that the resulting algorithm will not in general
produce a $\pi$-irreducible Markov chain. One simple alternative that we
propose is therefore to set $\mu$ to be an evenly-weighted mixture of two
Normal distributions centred at $\pm\sqrt{1-\sigma^{2}}$, each with variance
$\sigma^{2}<1$. The resulting approach, termed bi-modal Barker, will satisfy
$\mu_{6}=1+12\sigma^{2}+18\sigma^{4}-16\sigma^{6}$ and be
$15^{1/3}\mu_{6}^{-1/3}$ times more efficient than the version with Gaussian
noise. For small $\sigma$ this is close to optimal whilst also being
practical. For instance, for the choice $\sigma^{2}=0.1^{2}$, which is the one
we use in simulations below, bi-modal Barker is approximately $2.37$ times
more efficient than the Gaussian version.

The result on the Radamacher optimality may seem surprising at first given
given the lack of $\pi$-irreducibility. Similar results have, however, been
uncovered previously, for example it is known that the optimum expected
squared jump distance for the random walk Metropolis when the target
distribution is spherically symmetric is found by choosing the proposal
distribution to be uniform on a hyper-sphere of fixed radius from the current
point [10]. Given the product form of $\pi$ considered in this work, the
Rademacher structure is therefore natural. For the random walk Metropolis,
however, the benefits of choosing such an optimized proposal distribution
vanish as the dimension increases [10, 22], whereas in the Barker case they do
not.

An intuitive explanation for this may be that bi-modal Barker proposal choice
makes the MCMC method less diffusive and puts more effort on moving at least a
certain distance away. This is consistent with motivation for other kinds of
development of MCMC methods, for instance Hamiltonian Monte Carlo and non-
reversible Piece-wise deterministic Markov processes [5, 6].

###  4.2 Optimising over the choice of balancing function for a fixed noise
distribution

In this section we switch attention to the optimal choice of $g$ for a fixed
choice of $\mu$. The expression (9) in this case becomes a simple quadratic in
$\mathfrak{g}$, which can be straightforwardly solved to find an optimum
choice for a given $\phi$, as given in (11) below.

######  Proposition 4.

Given $\phi$ and a fixed noise distribution $\mu$ with finite fourth and sixth
moments $\mu_{4}<\mu_{6}<\infty$, the optimum choice of $\mathfrak{g}$ is

$\mathfrak{g}^{*}=\frac{\mu_{6}\left(C_{\phi}-3B_{\phi}\right)+\mu_{4}\left(9B_{\phi}-C_{\phi}\right)-6B_{\phi}}{12B_{\phi}\left(\mu_{6}-2\mu_{4}+1\right)}.$
(11)

Any family of balancing functions for which $\mathfrak{g}=g^{\prime\prime}(1)$
can be modified to take a desired value could therefore in principle be used
to create an optimized algorithm for a particular $\mu$ and $\phi$. Consider
the family

$g_{\gamma}(t)=\frac{1}{2}\left(t^{\frac{1}{2}+\gamma}+t^{\frac{1}{2}-\gamma}\right),$
(12)

indexed by $\gamma\geq 0$, where for $\gamma=0$ we recover the Langevin case
$g(t)=\surd t$. Any choice within the family is a balancing function, and is
such that $g_{\gamma}(1)=1$ and
$\mathfrak{g}=g_{\gamma}^{\prime\prime}(1)=\gamma^{2}-\frac{1}{4}$. For a
given $\phi$, the choice of $\gamma$ can therefore be adjusted to achieve the
optimum asymptotic efficiency provided that $\mathfrak{g}^{*}$ in (11) is
larger than $-1/4$.

Given the results of the previous section it would seem natural to set $\mu$
as a Rademacher distribution, however in this case it turns out that all
choices of $g$ give equivalent algorithms. This follows straightforwardly from
the fact that (2) implies $g(t)/\\{g(t)+g(t^{-1})\\}=1/(1+t^{-1})$, which is
independent of $g$. In fact Proposition 4 does not apply to the Rademacher
case since $\mu_{4}=\mu_{6}$. Another natural option is to fix $\mu$ to be
standard Gaussian. In this case (11) implies that the maximum efficiency is
found by choosing $\mathfrak{g}=C_{\phi}/(10B_{\phi})-1/5$. This scheme can be
implemented using the family in (12), and sampling from the resulting first
order locally-balanced proposal is straightforward as it consists in a mixture
of two Gaussians, see the supplement for details. We do not implement this
scheme in the simulations, however, in favour of the more efficient
alternatives discussed in the next section.

###  4.3 Optimising over the choice of both noise distribution and balancing
function

In this section we consider optimizing over both $g$ and $\mu$ jointly. The
following proposition identifies the best possibly achievable asymptotic
efficiency with first order locally-balanced proposals for a given target.

######  Proposition 5.

A non-negative lower bound for $\theta^{2}$ that is independent of both $\mu$
and $g$ is

$\theta^{2}\geq\frac{1}{144}\left(A_{\phi}-\frac{C_{\phi}^{2}}{B_{\phi}}\right).$
(13)

Furthermore, $\theta^{2}$ can be made arbitrarily close to the lower bound by
choosing $\mu_{4}>1$ sufficiently close to one, setting $\mu_{6}=\mu_{4}^{2}$
and choosing

$\mathfrak{g}=\frac{\mu_{4}(C_{\phi}-3B_{\phi})+6B_{\phi}}{12B_{\phi}(\mu_{4}-1)}.$
(14)

###### Proof.

Given $A_{\phi}>0,B_{\phi}>0$ and $C_{\phi}\in\mathbb{R}$ we must solve the
constrained quadratic optimisation problem of minimising $\theta^{2}$ subject
to $1\leq\mu_{4}\leq\surd\mu_{6}$. The constraints on $\mu_{4}$ and $\mu_{6}$
are necessary because $1=\mu_{2}\leq\surd\mu_{4}$ by Jensen’s inequality and
$\mu_{4}\leq\surd(\mu_{2}\mu_{6})=\surd\mu_{6}$ by Cauchy’s inequality.
Moreover the Hamburger moment problem tells us these constraints are
sufficient: if they are fulfilled then there exists a symmetric proposal
distribution on $\mathbb{R}$ that satisfies them.

Defining the new variables
$m_{1}=(\mathfrak{g}+1/4)\phi^{\prime}\phi^{\prime\prime}-\phi^{\prime\prime\prime}/12$
and $m_{2}=-(\mathfrak{g}+1/2)\phi^{\prime}\phi^{\prime\prime}$, we can
rewrite $\theta^{2}$ as

$\theta^{2}=\left(\mu_{6}-\mu_{4}^{2}\right)\mathbb{E}[m_{1}^{2}]+\mathbb{E}\left[\left(\mu_{4}m_{1}+m_{2}\right)^{2}\right]\geq\mathbb{E}\left[\left(\mu_{4}m_{1}+m_{2}\right)^{2}\right].$
(15)

where the inequality follows from $\surd\mu_{6}\geq\mu_{4}$. Expressing this
lower bound in terms of $A_{\phi},B_{\phi}$ and $C_{\phi}$ gives

$\theta^{2}\geq
B_{\phi}\left\\{\left(\mu_{4}\mathfrak{g}-\mathfrak{g}+\frac{\mu_{4}}{4}-\frac{1}{2}\right)-\frac{\mu_{4}}{12}\frac{C_{\phi}}{B_{\phi}}\right\\}^{2}+\frac{\mu_{4}^{2}}{144}\left(A_{\phi}-\frac{C_{\phi}^{2}}{B_{\phi}}\right),$

which can itself be lower bounded, giving

$\theta^{2}\geq\frac{\mu_{4}^{2}}{144}\left(A_{\phi}-\frac{C_{\phi}^{2}}{B_{\phi}}\right)\geq\frac{1}{144}\left(A_{\phi}-\frac{C_{\phi}^{2}}{B_{\phi}}\right).$

We have used three inequalities. The first, in (15), is realised if and only
if $\mu_{6}=\mu_{4}^{2}$; the second simply bounds a square below by zero and
is realised if and only if $\mathfrak{g}$ is defined as in (14), which
requires $\mu_{4}>1$; the third relies on $\mu_{4}\geq 1$ and is realised if
and only if $\mu_{4}=1$. Note that the last two equalities cannot be realised
simultaneously. The final lower bound is always non-negative due to
$B_{\phi}A_{\phi}\geq C_{\phi}^{2}$ by Cauchy’s inequality. ∎

Denote by $\nu(a)$ for $a>1$ a discrete symmetric distribution taking three
possible values $-\sqrt{a},0,\sqrt{a}$, such that the probability of a non-
zero value is $1/a$, and note that this is the unique symmetric distribution
$\mu$ with moments satisfying $\mu_{2}=1$, $\mu_{4}=a$ and $\mu_{6}=a^{2}$.
Letting $\mathfrak{g}$ be defined by (14), choosing $\mu=\nu(\mu_{4})$ and
taking $\mu_{4}$ arbitrarily close to one results in $\theta^{2}$ becoming
arbitrarily close to the lower bound (13). This three point proposal results
in an algorithm that achieves close to optimal asymptotic expected squared
jump distance among the class of first order locally-balanced samplers
provided that $\mathfrak{g}$ is chosen according to (14).

######  Remark 2.

This three point proposal is in fact also the optimal choice of $\mu$ for any
fixed choice of $g$, but the amount of mass given to point zero will vary
depending on $g$. In the Barker case, for example, this point achieves no
mass, resulting in the Rademacher choice for $\mu$.

It is natural to consider taking the limit $\mu_{4}\to 1$ and expect
optimality to be reached there. When the dimension $n$ is fixed and finite,
however, this results in a Rademacher proposal, which is suboptimal. This can
be seen by noting that the lower bound (13) is always smaller than
$B_{\phi}/16+A_{\phi}/144+C_{\phi}/24$, the value attained by the Rademacher
proposal, because

$\frac{B_{\phi}}{16}+\frac{A_{\phi}}{144}+\frac{C_{\phi}}{24}=\frac{A_{\phi}}{144}+\left(\frac{B_{\phi}^{1/2}}{4}+\frac{C_{\phi}}{12B_{\phi}^{1/2}}\right)^{2}-\frac{C_{\phi}^{2}}{144B_{\phi}}\geq\frac{1}{144}\left(A_{\phi}-\frac{C_{\phi}^{2}}{B_{\phi}}\right).$

Inspecting the proof of Theorem 2 shows that $\mathfrak{g}$ must be increased
sufficiently slowly as a function of $n$ to control the remainder terms in
order for the asymptotic expression for $\theta^{2}$ to be a valid
representation of the expected squared jump distance. In other words, as
$\mu_{4}\to 1$ it takes increasingly large $n$ for the asymptotic regime to be
representative of the finite $n$ setting. For a finite $n$, it is therefore
necessary to choose $\mu_{4}>1$. We explore this phenomenon further in the
supplement. In all simulations below, we set $\mu_{4}=2$ unless stated
otherwise.

A surprising consequence of these findings is that the three point proposal
with some mass at zero outperforms a Rademacher choice that is optimum for the
Barker proposal when the freedom to choose $\mathfrak{g}$ is given. In terms
of sampling, this suggests that efficiency gains can be made by allowing some
components of the state to remain unchanged at each iteration of the algorithm
with a probability that depends on the size of the gradient in that direction.
The same family of balancing functions introduced in (12) can again be used to
create this optimum sampler.

A particular case of interest is the Gaussian setting $\phi(x)=-x^{2}/2$, in
which case $\phi^{\prime\prime\prime}(0)$ and therefore $A_{\phi}=C_{\phi}=0$.
This means that by choosing any $\mu_{4}>1$ and $\mathfrak{g}$ according to
(14) we can achieve zero asymptotic $\theta^{2}$. The result of this is a
super-efficient sampler whose efficiency will effectively decay at a slower
rate than $n^{-1/3}$. We illustrate this surprising finding numerically in
Section 5, but also stress that this property only holds when
$\phi(x)=-x^{2}/2$ to the best of our knowledge.

##  5 Simulation Study

###  5.1 Efficiency with dimension on product targets

We examine the expected squared jump distance of the first component of two
different product form target distributions as a function of dimension. This
setting is directly captured by the theoretical results of Sections 3 and 4.
The two target distributions considered are the multi-dimensional standard
Gaussian distribution and the hyperbolic distribution of Example 4. In each
case we compare the random walk Metropolis, the Metropolis-adjusted Langevin
algorithm, Barker with Gaussian noise, Barker with Rademacher noise, Barker
with bi-modal noise as described in Section 4.1 and the optimal choice over
both balancing function and noise distribution described in Section 4.3, which
will hereafter be called the three point proposal.

The results for the Gaussian target distribution are shown in Figure 1(a). It
is clear from the plots that among the Barker algorithms the Rademacher and
bi-modal choices are comparable and perform similarly to MALA, whereas the
Barker algorithm with Gaussian noise has a lower expected squared jump
distance by a factor of 2-2.5, in accordance with the theoretical value of
2.47. The three points proposal performs best and appears to exhibit a
slightly slower than $n^{-1/3}$ decay in expected squared jumping distance
when the dimension in large. This is because in the special case of Gaussian
target $\theta^{2}$ from (9) equals zero when the choices described in Section
4.3 are made.

For the hyperbolic target results are shown in Figure 1(b). The main
difference compared to the Gaussian example is that now the Barker algorithms
with Rademacher and bi-modal noise both outperform the Langevin algorithm, as
predicted by the theory described in Section 4.1. The three points proposal is
still the best performing algorithm.

Figure 1: ESJD against dimensionality. Left: Gaussian product target. Right:
Hyperbolic product target.

###  5.2 Poisson random effects model

To consider a realistic example in which the target distribution is not of the
product form, we compare algorithms on the Poisson random effects model
described in Section 6.3 of [9]. We compare the Barker algorithm with bi-modal
noise to the Barker algorithm with Gaussian noise, the Langevin algorithm and
the random walk Metropolis. The main purpose of this example is to assess
whether or not the above theoretical guidelines for the noise distribution in
the Barker algorithm lead to good choices even when the target distribution
does not have independent and identically distributed components.

The target distribution under consideration is a $51$-dimensional posterior
distribution, $p(\mu,\eta_{1},\dots,\eta_{50}|\textbf{y})$, arising from a
Poisson random effects model defined hierarchically as
$\mu\sim\hbox{N}(0,10^{2})$, $\eta_{i}|\mu\sim\hbox{N}(\mu,\sigma_{\eta}^{2})$
and $y_{ij}|\eta_{i}\sim\hbox{Poisson}(\exp(\eta_{i}))$, independently for
$i=1,\dots,50$ and $j=1,\dots,5$. In our experiment we generate the observed
data $\textbf{y}=(y_{ij})_{ij}$ from the model likelihood, i.e. sampling
$y_{ij}\sim\hbox{Poisson}(\exp(\eta^{*}_{i}))$ independently, where
$\eta^{*}_{1},\dots,\eta^{*}_{I}$ are themselves generated independently from
a $\hbox{N}(\mu^{*},\sigma_{\eta}^{2})$ distribution with $\mu^{*}=5$. Here
$\sigma_{\eta}$ is a fixed value and two scenarios are considered: in the
first we set $\sigma_{\eta}=1$, while in the second we set $\sigma_{\eta}=3$.
Effectively, $\sigma_{\eta}$ is a parameter that governs the heterogeneity
across groups $i=1,\dots,50$ in the hierarchy. Thus, larger values of
$\sigma_{\eta}$ lead to a target distribution with more heterogeneity of
scales across coordinates, which make the adaptation and sampling process more
challenging.

In each case algorithmic tuning parameters consisting of a diagonal pre-
conditioning matrix and a global scale are learned using Algorithm 4 of [1],
in the same manner as described in Section 6.3 of [9]. We measure efficiency
in terms of effective sample size for a given number of iterations since all
algorithms under comparison, apart from Random Walk Metropolis, have a roughly
equivalent cost per iteration, which is dominated by gradient computations.
Figure 2 reports the median effective sample sizes across parameters for 100
independent runs of $5\times 10^{4}$ iterations of each algorithm. All
algorithms were randomly initialized by sampling parameter values from their
prior distributions.

Figure 2: Violin plots of median effective sample sizes (median across
parameters) for 100 independent repetitions of each algorithm. Left: low
heterogeneity across coordinates. Right: high heterogeneity across
coordinates.

Both versions of the Barker algorithm appear to be more robust to different
hyperparameter values than the Langevin algorithm, which sometimes performs
well but sometimes poorly in the first scenario and always performs poorly in
the second. This is because the Langevin algorithm is very sensitive to tuning
parameter selection, and the adaptive procedure fails to converge on sensible
values for these across the time scales of the simulation. The random walk
Metropolis also performs poorly, which is largely explained by the dimension
of the problem. The Barker algorithm with bi-modal noise is approximately two
times as efficient in terms of effective sample size as the version with
Gaussian noise in this setting. More precisely, the median improvement in
estimated effective sample size is 2.08 in scenario 1 (10th and 90th quantiles
across the 100 repetitions 2.05 and 2.11 respectively) and 2.04 in scenario 2
(10th and 90th quantiles 1.98 and 2.14 respectively). Similar numbers were
obtained when looking at minimum (rather than median) effective sample sizes
across parameters. These values suggest that the asymptotic theory developed
in this paper, which quantifies bi-modal Barker to be 2.37 times more
efficient than Gaussian Barker, is highly predictive of behaviours observed in
practice also for moderate dimensionality and targets that have neither
independent nor identically distributed coordinates. More generally, in all
our simulations, we consistently observed a improvement in efficiency when
going from Gaussian to bimodal Barker with factors typically between 2 and
2.5.

###  5.3 A correlated example

Unlike the Random Walk or Langevin algorithms, the Barker and three points
schemes rely on a choice of coordinate system. This may raise the concern of
how much performance depends on specific choices of coordinate systems, and in
particular whether the $\mathcal{O}(n^{1/3})$ scaling behaviour proved above
is sensitive to the theoretical assumption that the target factorizes across
the same coordinate axes as the proposal. Here we explore these issues
numerically, performing high dimensional scaling experiments similar to
Section 5.1 but for non-product form targets with significant correlation. In
particular, we consider Gaussian distributions with non-diagonal covariance
matrix $\Sigma$ chosen in two ways. In the first case we set $\Sigma_{ii}=1$
for $i=1,\dots,n$ and $\Sigma_{ij}=\rho$ for $i\neq j$, while in the second we
take $\Sigma_{ij}=\rho^{|i-j|}$. In both cases we set $\rho=0.99$ to depart
drastically from the independence case. As in Section 5.1 we compute the
expected square jump distance per coordinate. For all algorithms under
consideration we use isotropic proposals, meaning we do not use
preconditioning to avoid aligning proposal and target axes, and we choose a
step-size that is numerically optimized to maximize performance as measured by
expected square jump distance. The results are reported in Figure 3. As
expected, all schemes perform worse than in the product case (note the
different scales on the $y$-axes between Figure 1 and Figure 3), but the
relative comparison between different schemes remains nearly unchanged and
fully coherent with the theoretical predictions obtained from Sections 3 and
4. In particular the Langevin, Barker bi-modal and Barker Radamacher schemes
perform nearly equivalently, while Barker with Gaussian noise performs around
2-2.5 times worse. Overall, the experiment suggests that the relative
performances of the Random Walk, Langevin and Barker algorithms is not
particularly sensitive to correlation and to the specific choice of coordinate
system.

Figure 3: Expected squared jump distance against dimensionality for correlated
Gaussian targets. Left: $\Sigma_{ii}=1$ and $\Sigma_{ij}=0.99$ for $i\neq j$.
Right: $\Sigma_{ij}=0.99^{|i-j|}$.

The three point proposal performs well also in these correlated examples and
actually performs surprisingly well when $\Sigma_{ij}=\rho^{|i-j|}$. Providing
better understanding of such unexpected behaviour will be the subject of
future research. Note, however, that the three point proposal implicitly uses
knowledge about the target distribution when choosing the optimal values of
the tuning parameters $\mathfrak{g}$ and $\mu_{4}$, and thus it has been given
a somewhat unfair and potentially unrealistic advantage compared to the other
schemes considered here. In particular, in this example $\mathfrak{g}$ was
chosen according to the optimal value in (14) with $B_{\phi}=1$ and
$C_{\phi}=0$ as given by product-form Gaussian targets.

##  6 Discussion

The main results of this paper rely on a product form structure of $\pi$, and
the corresponding optimal choice of locally-balanced algorithm also has a
product form. We have shown in Section 5 that this choice is still effective
when the target distribution is no longer of the product form, and therefore
recommend the use of the bi-modal Barker algorithm in practice. It is
surprising that using a non-local noise distribution of this kind results in
such a pronounced and consistent improvement in efficiency across multiple
examples. We believe that this represents a good case study of theoretical
analysis motivating new practical methodology that would be otherwise be hard
to devise. It is also worth noting that any improvement in efficiency
discussed above essentially comes for free, since all the gradient-based
schemes considered in the paper have a comparable cost per iteration, which is
typically dominated by gradient computations, and all schemes are equally
simple to implement.

The detailed quantitative analysis and comparison of algorithms within the
locally-balanced class in the high-dimensional limit is made possible by the
mathematical framework developed in Section 3 of [21]. This framework
identifies and uses only essential Taylor series expansions related to the
limiting Kullback–Leibler divergence between a locally-balanced proposal and
its time reversal. Using this we establish optimal scaling for a broad class
of algorithms including Barker and Langevin with a single unified proof, along
with significantly weaker assumptions on the smoothness and tails of the
target distribution than those in [14]. Our results are at present restricted
to limiting expected squared jump distances, rather than diffusion limits as
in [13] or [14], but we believe that it is possible to uncover a limiting
process under the current assumptions and such a line of enquiry is being
pursued at the time of writing, continuing the axiomatic approach in [21].

One intriguing finding of this work concerns the sub-optimality of the
Langevin choice $g(t)=\surd t$ with proposal input noise $\mu$ chosen to be
Gaussian. This is by far the most historically popular choice within the first
order locally-balanced class of algorithms. The results in this paper show
that according to asymptotic efficiency as measured by expected squared jump
distance not only is this combination of $\mu$ and $g$ not optimum, but in
addition that the optimum choice of $\mu$ when $g(t)=\surd t$ is not Gaussian,
and also that the optimum choice of $g$ when using Gaussian $\mu$ is not
$\surd t$.

A natural open question is whether the insights of Section 4.3 can be used to
create a novel new algorithm based on the three point proposal scheme. We have
resisted doing so here because such an algorithm would require a problem-
specific choice of balancing function and some appropriate randomisation of
the noise distribution to prevent reducibility issues. It may be possible,
however, to design an adaptive Markov chain Monte Carlo method that is able to
learn these quantities during the simulation. We look forward to designing
practical methodology based on the insights of Section 4.3 in subsequent work.

## Acknowledgement

JV was supported by a UK Engineering and Physical Sciences Research Council
grants EP/R022100/1 and EP/T004134/1. SL is supported by a UK Engineering and
Physical Sciences Research Council grant EP/V055380/1.

##  Appendix A Sampling from optimal locally-balanced kernels with Gaussian
noise

We consider here fixing $\mu$ to be standard Gaussian, as is the choice made
in default versions of the Langevin and Barker algorithms. Setting
$\mu_{6}=15$ and $\mu_{4}=3$ into (11) implies that the maximum efficiency is
found by choosing $\mathfrak{g}=C_{\phi}/(10B_{\phi})-1/5$. Such value of
$\mathfrak{g}$ can be imposed using, e.g., the family of balancing functions
defined in (12).

Sampling from a first order locally-balanced proposal using (12) and standard
Gaussian $\mu$ can then be done using Algorithm 2 applied to each coordinate.
Also, any choice of balancing function within the family (12) and $\mu$ chosen
as a mixture of Gaussians results in a proposal density with analytically
tractable normalising constant, meaning Metropolis–Hastings acceptance rates
can be evaluated.

Algorithm 2 Simulate from the locally-balanced proposal using (12). Require:
$x\in\mathbb{R}$ and $\sigma>0$
---
Set $p_{x}\leftarrow(1+e^{-\gamma\sigma^{2}\phi^{\prime}(x)^{2}})^{-1}$
Draw $u_{x}\sim\text{Bernoulli}(p_{x})$ and set $b_{x}\leftarrow 2u_{x}-1$
Draw $z\sim N(0,1)$ and set $y\leftarrow
x+\left(1/2+b_{x}\gamma\right)\sigma^{2}\phi^{\prime}(x)+\sigma z$
Output $y$

The procedure can be viewed as taking elements of both the Langevin and Barker
proposal, as the standard Langevin proposal is modified according to a random
variable that takes value $+1$ or $-1$ with some probability that is skewed in
the direction of the gradient as in the Barker proposal. An interesting aspect
of the Gaussian noise $\mu$ is that it is independent of the coordinate
system, making the Algorithm 2 depend on the coordinate system only through
the choice of the flipping directions.

######  Proposition 6.

Algorithm 2 produces a sample from a distribution with density proportional to
$\sigma^{-1}\mu((y-x)/\sigma)g_{\gamma}(e^{\phi^{\prime}(x)(y-x)})$.

###### Proof.

First note that

$p(x)=\frac{e^{-\frac{1}{2}(\gamma+\frac{1}{2})^{2}\sigma^{2}\phi^{\prime}(x)}}{e^{-\frac{1}{2}(\gamma+\frac{1}{2})^{2}\sigma^{2}\phi^{\prime}(x)}+e^{\frac{1}{2}(\gamma-\frac{1}{2})^{2}\sigma^{2}\phi^{\prime}(x)}}$

so that

$p(x)e^{-\frac{1}{2}(\gamma+\frac{1}{2})^{2}\sigma^{2}\phi^{\prime}(x)}=(1-p(x))e^{-\frac{1}{2}(\gamma-\frac{1}{2})^{2}\sigma^{2}\phi^{\prime}(x)}\,.$

The random variable $Y$ is a mixture of two Gaussians with the same variance.
Setting $\mu_{x}^{+}=x+(1/2+\gamma)\sigma^{2}\phi^{\prime}(x)$ and
$\mu_{x}^{-}=x+(1/2-\gamma)\sigma^{2}\phi^{\prime}(x)$, the density of
$\xi(y)$ of $Y$ satisfies

$\displaystyle\sigma\sqrt{2\pi}\xi(y)$
$\displaystyle=p(x)e^{-\frac{1}{2\sigma^{2}}(y-\mu_{x}^{+})^{2}}+\left(1-p(x)\right)e^{-\frac{1}{2\sigma^{2}}(y-\mu_{x}^{-})^{2}}$
$\displaystyle=e^{-\frac{(y-x)^{2}}{2\sigma^{2}}}\left(p(x)e^{-\frac{1}{2}\left(\gamma+\frac{1}{2}\right)^{2}\sigma^{2}\phi^{\prime}(x)^{2}}\right)\left(e^{\left(\frac{1}{2}+\gamma\right)(y-x)\phi^{\prime}(x)}+e^{\left(\frac{1}{2}-\gamma\right)(y-x)\phi^{\prime}(x)}\right)$
$\displaystyle\propto\frac{1}{2}\left((e^{(y-x)\phi^{\prime}(x)})^{\frac{1}{2}+\gamma}+(e^{(y-x)\phi^{\prime}(x)})^{\frac{1}{2}-\gamma}\right)\times
e^{-\frac{(y-x)^{2}}{2\sigma^{2}}}$ $\displaystyle\propto
g_{\gamma}(e^{\phi^{\prime}(x)(y-x)})\sigma^{-1}\mu((y-x)/\sigma)$

as required. ∎

##  Appendix B Illustration of different choices of $\mu_{4}$ in the optimal
choice of first order locally-balanced proposal

Here we compare the three point proposal with different choices of $\mu_{4}$
on the Hyperbolic target distribution considered in Section 5.1 of the paper.
For any value of $\mu_{4}$, the optimal value of $\mathfrak{g}$ as given by
Proposition 5 in the paper is used. Theory suggests that optimal asymptotic
performances are given by choosing $\mu_{4}$ arbitrarily close to $1$ but, as
discussed in the paper, values of $\mu_{4}$ close to $1$ may need very large
dimensionality for the actual asymptotic regime to kick in.

Figure 4: Expected squared jump distance against dimensionality for various
three point proposals on an hyperbolic product target distribution.

Figure 4 illustrates this phenomenon. For small dimensionality, the value
$\mu_{4}=2$ is best (among the ones considered) but as dimensionality
increases $\mu_{4}=1.5$ starts outperforming and then $\mu_{4}=1.25$.
Dimension $n=35000$ is not yet sufficient to for $\mu_{4}=1.1$ to outperform
$\mu_{4}=1.25$, although theory suggests that this will eventually happen.

##  Appendix C Miscellaneous proofs

###  C.1 Section 2

###### Proof of Proposition 1.

For the first part of the statement, direct calculation gives

$tg_{h}(1/t)=t\times t^{-1/2}h(-\log t)=t^{1/2}h(\log t)=g_{h}(t).$

Similarly for the second part

$h_{g}(-x)=e^{x/2}g_{h}(e^{-x})=e^{-x/2}g_{h}(e^{x})=h_{g}(x)$

using that $g_{h}(e^{-x})=e^{-x}g_{h}(e^{x})$. Noting also that
$h_{g}:\mathbb{R}\to[0,\infty)$ shows that $h_{g}\in\mathcal{H}$. ∎

###  C.2 Section 3

###### Proof of Theorem 1.

See Section 3, particularly the proof of Theorem 8, in [21]. ∎

###### Proof of Proposition 2.

(i) Easy to verify.

(ii) Since $\int_{\mathbb{R}}z^{\xi}\mu(dz)<\infty$, we have

$\int_{\mathbb{R}}e^{b(az)}|z|^{\xi}\mu(dz)=\int_{\mathbb{R}}g(e^{az})|z|^{\xi}\mu(dz)\leq\|g\|_{\infty}\int_{\mathbb{R}}|z|^{\xi}\mu(dz)\,.$

On the other hand $b$ is non decreasing and $b(0)=0$ so

$\int_{\mathbb{R}}e^{b(az)}\mu(dz)\geq\int_{0}^{\infty}\mu(dz)=\frac{1}{2}\,.$

(iii) This essentially holds because $b^{\prime}$ is bounded and hence $b$ is
at most linear. Clearly $\mu$ has finite polynomial moments. Next note that by
fundamental theorem of calculus and $b(0)=0$ we have

$b(az)=\int_{0}^{1}azb^{\prime}(azs)ds\,.$

Hence

$\displaystyle\int_{\mathbb{R}}e^{b(az)}|z|^{\xi}\mu(z)dz$
$\displaystyle=\int_{\mathbb{R}}e^{az\int_{0}^{1}b^{\prime}(saz)ds}|z|^{\xi}\mu(z)dz$
$\displaystyle\leq\tilde{C}_{\mu}\left(1+\left|a\int_{0}^{1}b^{\prime}(saz)ds\right|^{\tilde{\beta}}\right)\int_{\mathbb{R}}e^{az\int_{0}^{1}b^{\prime}(saz)ds}\mu(z)dz$
$\displaystyle\leq\tilde{C}_{\mu}\max(1,\|b^{\prime}\|_{\infty}^{\tilde{\beta}})\left(1+|a|^{\tilde{\beta}}\right)\int_{\mathbb{R}}e^{b(az)}\mu(z)dz\,.$

(iv) We will prove this implies a special instance of (iii). Due to symmetry
of $\mu$ we may assume without loss of generality that $a>0$. Firstly, $\mu$
has finite moment generating function since for any $b>a$ and $c>0$ we have

$\int_{\mathbb{R}}e^{az}\mu(z)dz\leq
2\int_{0}^{\infty}e^{bz}e^{(a-b)z}\mu(z)dz\leq 2c\left(\sup_{z\leq
c}e^{az}\mu(z)\right)+\frac{2}{b-a}e^{-c(b-a)}\left(\sup_{z>c}e^{bz}\mu(z)\right)\,.$

Secondly, we will use inequalities
$|z|^{\lambda}1_{[-1,1]^{c}}(z)\leq|z|^{\lambda+1}$ and
$e^{a}=e^{a\int_{0}^{\infty}2\mu(z)dz}\leq
2\int_{0}^{\infty}e^{az}\mu(z)dz\leq 2\int_{\mathbb{R}}e^{az}\mu(z)dz$, which
follows by Jensen’s inequality. Integration by parts implies that for any
$\lambda>p-1$ we have

$\displaystyle\int_{\mathbb{R}}e^{az}|z|^{\lambda}\mu(z)dz$
$\displaystyle\leq\int_{[-1,1]}e^{az}\mu(z)dz+A\int_{[-1,1]^{c}}e^{az}|z|^{\lambda-p}\mu(z)dz-B\int_{[-1,1]^{c}}e^{az}|z|^{\lambda-p}z\mu^{\prime}(z)dz$
$\displaystyle=\int_{[-1,1]}e^{az}\mu(z)dz+A\int_{[-1,1]^{c}}e^{az}|z|^{\lambda-p}\mu(z)dz$
$\displaystyle\qquad+B\int_{[-1,1]^{c}}e^{az}((\lambda-p+1)|z|^{\lambda-p}+a|z|^{\lambda-p}z)\mu(z)dz+B\mu(1)(e^{a}-e^{-a})$
$\displaystyle\leq\left(1+2B\mu(1)\right)\int_{\mathbb{R}}e^{az}\mu(z)dz+\left(A+B(\lambda-p+1)+|a|B\right)\int_{\mathbb{R}}e^{az}|z|^{\lambda-p+1}\mu(z)dz\,.$

Thus, we have successfully reduced the power of $z$ in the integrand by $p-1$
at the expense of producing an $|a|$. Noting that $|z|^{\xi}\leq
1+|z|^{\lceil\frac{\xi}{p-1}\rceil(p-1)}$ and recursively applying the same
argument establishes the claim. ∎

###### Proof of Theorem 2.

Denote $b(x)=\log g(e^{x})$ and note that $b(0)=0$, $b^{\prime}(0)=1/2$ and
$b^{\prime\prime}(0)=1/4+g^{\prime\prime}(1)$. First consider $n$ as fixed. By
definition for any $y\in\mathbb{R}$

$\rho_{n}(x,y)=\phi(y)-\phi(x)+b(\phi^{\prime}(y)(x-y))-b(\phi^{\prime}(x)(y-x))-\log(Z_{\sigma_{n}}(x))+\log(Z_{\sigma_{n}}(y))\,.$
(16)

We will be using the fundamental theorem of calculus

$U(1)-U(0)=\int_{0}^{1}U^{\prime}(s)ds$ (17)

for various functions $U$. Denote $w=(y-x)/\sigma_{n}$, and note that
$Z_{\sigma_{n}}(x)=\int_{\mathbb{R}}e^{b(\phi^{\prime}(x)\sigma_{n}z)}\mu(dz)$
by an analogous substitution. Then note that

$\log\left(\frac{Z_{\sigma_{n}}(x+\sigma_{n}w)}{Z_{\sigma_{n}}(x)}\right)=\sigma^{2}_{n}w\int_{0}^{1}\phi^{\prime\prime}(x+t\sigma_{n}w)\frac{\int_{\mathbb{R}}e^{b(\phi^{\prime}(x+t\sigma_{n}w)\sigma_{n}z)}b^{\prime}(\phi^{\prime}(x+t\sigma_{n}w)\sigma_{n}z)z\mu(dz)}{\int_{\mathbb{R}}e^{b(\phi^{\prime}(x+t\sigma_{n}w)\sigma_{n}z)}\mu(dz)}dt$

using (17) with the function $U_{1}(t)=\log(Z_{\sigma_{n}}(x+t\sigma_{n}w)$.
Using (17) again with the function
$U_{2}(s)=b^{\prime}(s\phi^{\prime}(x+t\sigma_{n}w)\sigma_{n}z)$ and recalling
that $b^{\prime}(0)=1/2$, this can be written as

$\displaystyle\frac{\sigma_{n}^{2}}{2}w\int_{0}^{1}\phi^{\prime\prime}(x+t\sigma
w)\frac{\int_{\mathbb{R}}e^{b(\phi^{\prime}(x+t\sigma_{n}w)\sigma_{n}z)}z\mu(dz)}{\int_{\mathbb{R}}e^{b(\phi^{\prime}(x+t\sigma_{n}w)\sigma_{n}z)}\mu(dz)}dt$
$\displaystyle\quad+\sigma_{n}^{3}w\iint_{[0,1]^{2}}\phi^{\prime\prime}(x+t\sigma_{n}w)\phi^{\prime}(x+t\sigma_{n}w)\frac{\int_{\mathbb{R}}e^{b(\phi^{\prime}(x+t\sigma_{n}w)\sigma_{n}z)}b^{\prime\prime}(s\phi^{\prime}(x+t\sigma_{n}w)\sigma_{n}z)z^{2}\mu(dz)}{\int_{\mathbb{R}}e^{b(\phi^{\prime}(x+t\sigma_{n}w)\sigma_{n}z)}\mu(dz)}dtds.$

Using (17) in the numerator of the first term with the function
$U_{3}(u)=e^{b(u\phi^{\prime}(x+t\sigma_{n}w)\sigma_{n}z)}$ and recalling that
$b(0)=0$ gives

$\displaystyle\frac{\sigma_{n}^{3}w}{2}\iint_{[0,1]^{2}}\phi^{\prime\prime}(x+t\sigma_{n}w)\phi^{\prime}(x+t\sigma_{n}w)\frac{\int_{\mathbb{R}}e^{b(u\phi^{\prime}(x+t\sigma_{n}w)\sigma_{n}z)}b^{\prime}(u\phi^{\prime}(x+t\sigma_{n}w)\sigma_{n}z)z^{2}\mu(dz)}{\int_{\mathbb{R}}e^{b(\phi^{\prime}(x+t\sigma_{n}w)\sigma_{n}z)}\mu(dz)}dtdu$
$\displaystyle\qquad+\sigma_{n}^{3}w\iint_{[0,1]^{2}}\phi^{\prime\prime}(x+t\sigma_{n}w)\phi^{\prime}(x+t\sigma_{n}w)\frac{\int_{\mathbb{R}}e^{b(\phi^{\prime}(x+t\sigma_{n}w)\sigma_{n}z)}b^{\prime\prime}(s\phi^{\prime}(x+t\sigma_{n}w)\sigma_{n}z)z^{2}\mu(dz)}{\int_{\mathbb{R}}e^{b(\phi^{\prime}(x+t\sigma_{n}w)\sigma_{n}z)}\mu(dz)}dtds.$

And further simplification leads to the expression

$\displaystyle\sigma_{n}^{3}w$
$\displaystyle\iint_{[0,1]^{2}}\phi^{\prime\prime}(x+t\sigma_{n}w)\phi^{\prime}(x+t\sigma_{n}w)$
$\displaystyle\qquad\times\frac{\int_{\mathbb{R}}e^{b(\phi^{\prime}(x+t\sigma_{n}w)\sigma_{n}z)}\left(\frac{1}{2}b^{\prime}(s\phi^{\prime}(x+t\sigma_{n}w)\sigma_{n}z)+b^{\prime\prime}(s\phi^{\prime}(x+t\sigma_{n}w)\sigma_{n}z)\right)z^{2}\mu(dz)}{\int_{\mathbb{R}}e^{b(\phi^{\prime}(x+t\sigma_{n}w)\sigma_{n}z)}\mu(dz)}dtds$
$\displaystyle=\sigma_{n}^{3}w\phi^{\prime\prime}(x)\phi^{\prime}(x)\left(\frac{1}{2}b^{\prime}(0)+b^{\prime\prime}(0)\right)+R_{1,n}(x,w)+R_{2,n}(x,w)+R_{3,n}(x,w)$
$\displaystyle=\sigma_{n}^{3}w\phi^{\prime\prime}(x)\phi^{\prime}(x)\left(\frac{1}{2}+g^{\prime\prime}(1)\right)+R_{1,n}(x,w)+R_{2,n}(x,w)+R_{3,n}(x,w),$

where

$\displaystyle R_{1,n}(x,w)$
$\displaystyle=\sigma_{n}^{3}w\left(\frac{1}{2}b^{\prime}(0)+b^{\prime\prime}(0)\right)\int_{0}^{1}\left(\phi^{\prime\prime}(x+t\sigma_{n}w)\phi^{\prime}(x+t\sigma_{n}w)-\phi^{\prime\prime}(x)\phi^{\prime}(x)\right)dt$
$\displaystyle R_{2,n}(x,w)$
$\displaystyle=\sigma_{n}^{3}w\left(\frac{1}{2}b^{\prime}(0)+b^{\prime\prime}(0)\right)\int_{0}^{1}\phi^{\prime\prime}(x+t\sigma_{n}w)\phi^{\prime}(x+t\sigma_{n}w)\frac{\int_{\mathbb{R}}e^{b(\phi^{\prime}(x+t\sigma_{n}w)\sigma_{n}z)}\left(z^{2}-1\right)\mu(dz)}{\int_{\mathbb{R}}e^{b(\phi^{\prime}(x+t\sigma_{n}w)\sigma_{n}z)}\mu(dz)}dt$
$\displaystyle R_{3,n}(x,w)$
$\displaystyle=\sigma_{n}^{3}w\int_{[0,1]^{2}}\phi^{\prime\prime}(x+t\sigma_{n}w)\phi^{\prime}(x+t\sigma_{n}w)dt$
$\displaystyle\times$
$\displaystyle\frac{\int_{\mathbb{R}}e^{b(\phi^{\prime}(x+t\sigma_{n}w)\sigma_{n}z)}\left(\frac{1}{2}b^{\prime}(s\phi^{\prime}(x+t\sigma_{n}w)\sigma_{n}z)+b^{\prime\prime}(s\phi^{\prime}(x+t\sigma_{n}w)\sigma_{n}z)-\frac{1}{2}b^{\prime}(0)-b^{\prime\prime}(0)\right)z^{2}\mu(dz)}{\int_{\mathbb{R}}e^{b(\phi^{\prime}(x+t\sigma_{n}w)\sigma_{n}z)}\mu(dz)}dtds\,.$

To expand the remaining terms of (16) we use the exact third order Taylor
expansion

$\displaystyle b(t)$
$\displaystyle=b(0)+tb^{\prime}(0)+t^{2}\frac{b^{\prime\prime}(0)}{2}+t^{3}\frac{b^{\prime\prime\prime}(0)}{6}+\frac{t^{3}}{2}\int_{0}^{1}(b^{\prime\prime\prime}(ut)-b^{\prime\prime\prime}(0))(1-u)^{2}du$
$\displaystyle=b(0)+tb^{\prime}(0)+t^{2}\frac{b^{\prime\prime}(0)}{2}+\frac{t^{3}}{2}\int_{0}^{1}b^{\prime\prime\prime}(ut)(1-u)^{2}du$

on each of the terms involving $b$ to obtain

$\displaystyle
b(-\phi^{\prime}(x+\sigma_{n}w)\sigma_{n}w)-b(\phi^{\prime}(x)\sigma_{n}w)$
$\displaystyle=-\frac{\sigma_{n}w}{2}\left(\phi^{\prime}(x)+\phi^{\prime}(x+\sigma_{n}w)\right)$
$\displaystyle\qquad+\sigma_{n}^{2}w^{2}\left(\frac{1}{8}+\frac{g^{\prime\prime}(1)}{2}\right)\left(\phi^{\prime}(x+\sigma_{n}w)^{2}-\phi^{\prime}(x)^{2}\right)$
$\displaystyle\qquad+R_{6,n}(x,w)+R_{7,n}(x,w)\,,$

where

$\displaystyle R_{6,n}(x,w)$
$\displaystyle=-\frac{1}{2}\sigma_{n}^{3}w^{3}\phi^{\prime}(x+\sigma_{n}w)^{3}\int_{0}^{1}b^{\prime\prime\prime}(-u\phi^{\prime}(x+\sigma_{n}w)\sigma_{n}w)(1-u)^{2}du\,,$
$\displaystyle R_{7,n}(x,w)$
$\displaystyle=-\frac{1}{2}\sigma_{n}^{3}w^{3}\phi^{\prime}(x)^{3}\int_{0}^{1}b^{\prime\prime\prime}(u\phi^{\prime}(x)\sigma_{n}w)(1-u)^{2}du\,.$

Taking one half of the difference of the exact Taylor expansions

$\displaystyle\phi(x+\sigma_{n}w)-\phi(x)$
$\displaystyle=\sigma_{n}w\phi^{\prime}(x)+\sigma_{n}^{2}w^{2}\int_{0}^{1}\phi^{\prime\prime}(x+u\sigma_{n}w)(1-u)du\,,$
$\displaystyle\phi(x)-\phi(x+\sigma_{n}w)$
$\displaystyle=-\sigma_{n}w\phi^{\prime}(x+\sigma_{n}w)+\sigma_{n}^{2}w^{2}\int_{0}^{1}\phi^{\prime\prime}(x+(1-u)\sigma_{n}w)(1-u)du$

and setting $v=1-u$ reveals that

$\displaystyle\phi(x+\sigma_{n}w)-\phi(x)-\frac{\sigma_{n}w}{2}\left(\phi^{\prime}(x)+\phi^{\prime}(x+\sigma_{n}w)\right)$
$\displaystyle=\frac{\sigma_{n}^{2}w^{2}}{2}\int_{0}^{1}\phi^{\prime\prime}(x+u\sigma_{n}w)(1-2u)du$
$\displaystyle=\frac{\sigma_{n}^{2}w^{2}}{2}\int_{0}^{1}\left(\phi^{\prime\prime}(x+u\sigma_{n}w)-\phi^{\prime\prime}(x)\right)(1-2u)du$
$\displaystyle=\frac{\sigma_{n}^{3}w^{3}}{2}\int_{0}^{1}u(1-2u)\int_{0}^{1}\phi^{\prime\prime\prime}(x+uv\sigma_{n}w)dvdu$
$\displaystyle=-\frac{\sigma_{n}^{3}w^{3}}{12}\phi^{\prime\prime\prime}(x)+R_{4,n}(x,w),$

where

$R_{4,n}(x,w)=\frac{\sigma_{n}^{3}w^{3}}{2}\int_{0}^{1}u(1-2u)\int_{0}^{1}\left(\phi^{\prime\prime\prime}(x+uv\sigma_{n}w)-\phi^{\prime\prime\prime}(x)\right)dvdu\,.$

Similarly,

$\displaystyle\sigma_{n}^{2}w^{2}\left(\frac{1}{8}+\frac{g^{\prime\prime}(1)}{2}\right)\left(\phi^{\prime}(x+\sigma_{n}w)^{2}-\phi^{\prime}(x)^{2}\right)$
$\displaystyle=\sigma_{n}^{3}w^{3}\left(\frac{1}{4}+g^{\prime\prime}(1)\right)\int_{0}^{1}\phi^{\prime}\phi^{\prime\prime}(x+u\sigma_{n}w)du$
$\displaystyle=\sigma_{n}^{3}w^{3}\left(\frac{1}{4}+g^{\prime\prime}(1)\right)\phi^{\prime}(x)\phi^{\prime\prime}(x)+R_{5,n}(x,w),$

where

$R_{5,n}(x,w)=\sigma_{n}^{3}w^{3}\left(\frac{1}{4}+g^{\prime\prime}(1)\right)\int_{0}^{1}\left(\phi^{\prime}\phi^{\prime\prime}(x+u\sigma_{n}w)-\phi^{\prime}\phi^{\prime\prime}(x)\right)du.$

Next, we will denote with $\tilde{\mathbb{E}}_{n}$ the expectation with
respect to the measure
$Z_{\sigma_{n}}^{-1}(x)e^{b(\phi^{\prime}(x)\sigma_{n}w)}\pi(x)dx\mu(dw)$. We
will show that $\tilde{\mathbb{E}}_{n}[R_{i,n}^{2}]$ decays faster than
$\sigma_{n}^{6}$ for $i=1,2,\dots,7$ as $n\to\infty$. Note that Assumption 1
iii) is the tool enabling the control of the Radon-Nikodym derivatives of
various proposals (indexed by $n$) with respect to $\mu$, making it possible
for us to represent these expectations with respect to the measure
$\pi(x)dx\mu(dw)$. Using $a\lesssim b$ to denote that $a\leq cb$ for some
positive finite $c$ independent of $n$, the bound for the first term $R_{1,n}$
follows by Assumption 1 i) and iii) as

$\displaystyle\tilde{\mathbb{E}}_{n}[R_{1,n}^{2}]$
$\displaystyle\lesssim\sigma_{n}^{6}\iint_{\mathbb{R}^{2}}\int_{0}^{1}\left(\phi^{\prime\prime}(x+t\sigma_{n}w)\phi^{\prime}(x+t\sigma_{n}w)-\phi^{\prime\prime}(x)\phi^{\prime}(x)\right)^{2}w^{2}\frac{e^{b(\phi^{\prime}(x)\sigma_{n}w)}}{Z_{\sigma_{n}}(x)}dt\mu(dw)\pi(x)dx$
$\displaystyle\leq\sigma_{n}^{6}\iint_{\mathbb{R}^{2}}\int_{0}^{1}K(x)^{2}w^{2}\max(|t\sigma_{n}w|^{H},|t\sigma_{n}w|^{\gamma})^{2}\frac{e^{b(\phi^{\prime}(x)\sigma_{n}w)}}{Z_{\sigma_{n}}(x)}dt\mu(dw)\pi(x)dx$
$\displaystyle\lesssim\sigma_{n}^{6+2H}\int_{\mathbb{R}}K(x)^{2}\pi(x)\int_{\mathbb{R}}|w|^{\xi}\frac{e^{b(\phi^{\prime}(x)\sigma_{n}w)}}{Z_{\sigma_{n}}(x)}\mu(dw)dx$
$\displaystyle\lesssim\sigma_{n}^{6+2H+\beta}\int_{\mathbb{R}}K(x)^{2}(1+|\phi^{\prime}(x)|^{\beta})\pi(x)dx\,.$

To bound the second term $\tilde{\mathbb{E}}_{n}R_{2,n}$, note that
$\lim_{\sigma_{n}\to 0}\sigma_{n}^{-3}R_{2,n}=0$ for every $x$ and $w$, so we
only need to provide a dominating bound for the integrand. Using Assumption 1
iii) gives that this is

$\lesssim\iint_{\mathbb{R}^{2}}\int_{0}^{1}\phi^{\prime\prime}(x+t\sigma_{n}w)^{2}\phi^{\prime}(x+t\sigma_{n}w)^{2}\left(1+\sigma_{n}^{2\beta}\phi^{\prime}(x+t\sigma_{n}w)^{2\beta}\right)\frac{e^{b(\phi^{\prime}(x)\sigma_{n}w)}}{Z_{\sigma_{n}}(x)}w^{2}dt\mu(dw)\pi(x)dx.$

To show that this is finite we use Assumption 1 i) to control the perturbation
of functions $\phi^{\prime\prime}\phi^{\prime}$ and
$\phi^{\prime\prime}|\phi^{\prime}|^{1+\beta}$ from $x$ to $x+t\sigma_{n}w$
and Assumption 1 iii) to control the Radon-Nikodym derivative
$Z_{\sigma_{n}}(x)^{-1}e^{b(\phi^{\prime}(x)\sigma_{n}w)}$ (as in the argument
for the term $\tilde{\mathbb{E}}_{n}[R_{1,n}]$):

$\displaystyle\iint_{\mathbb{R}^{2}}\int_{0}^{1}\phi^{\prime\prime}(x+t\sigma_{n}w)^{2}\phi^{\prime}(x+t\sigma_{n}w)^{2}\left(1+|\phi^{\prime}(x+t\sigma_{n}w)|^{\beta}\right)^{2}w^{2}\frac{e^{b(\phi^{\prime}(x)\sigma_{n}w)}}{Z_{\sigma_{n}}(x)}dt\mu(dw)\pi(x)dx$
$\displaystyle\quad\lesssim\iint_{\mathbb{R}^{2}}\phi^{\prime\prime}(x)^{2}\phi^{\prime}(x)^{2}\left(1+|\phi^{\prime}(x)|^{\beta}\right)^{2}w^{2}\frac{e^{b(\phi^{\prime}(x)\sigma_{n}w)}}{Z_{\sigma_{n}}(x)}\mu(dw)\pi(x)dx$
$\displaystyle\qquad+\iint_{\mathbb{R}^{2}}\int_{0}^{1}\left(\phi^{\prime\prime}(x+t\sigma_{n}w)\phi^{\prime}(x+t\sigma_{n}w)-\phi^{\prime\prime}(x)\phi^{\prime}(x)\right)^{2}dtw^{2}\frac{e^{b(\phi^{\prime}(x)\sigma_{n}w)}}{Z_{\sigma_{n}}(x)}\mu(dw)\pi(x)dx$
$\displaystyle\qquad+\iint_{\mathbb{R}^{2}}\int_{0}^{1}\left(\phi^{\prime\prime}(x+t\sigma_{n}w)|\phi^{\prime}(x+t\sigma_{n}w)|^{1+\beta}-\phi^{\prime\prime}(x)|\phi^{\prime}(x)|^{1+\beta}\right)^{2}dtw^{2}\frac{e^{b(\phi^{\prime}(x)\sigma_{n}w)}}{Z_{\sigma_{n}}(x)}\mu(dw)\pi(x)dx$
$\displaystyle\quad\lesssim\int_{\mathbb{R}}\phi^{\prime\prime}(x)^{2}\phi^{\prime}(x)^{2}(1+|\phi^{\prime}(x)|^{\beta})^{3}\pi(x)dx+\sigma_{n}^{2H+\beta}\int_{\mathbb{R}}K(x)^{2}\pi(x)\int_{\mathbb{R}}|w|^{\xi}\frac{e^{b(\phi^{\prime}(x)\sigma_{n}w)}}{Z_{\sigma_{n}}(x)}\mu(dw)dx$
$\displaystyle\quad\lesssim\int_{\mathbb{R}}\phi^{\prime\prime}(x)^{2}\phi^{\prime}(x)^{2}(1+|\phi^{\prime}(x)|^{\beta})^{3}\pi(x)dx+\int_{\mathbb{R}}K(x)^{2}(1+|\phi^{\prime}(x)|^{\beta})\pi(x)dx.$

The term $\tilde{\mathbb{E}}[R_{3,n}]$ is handled similarly with the dominated
convergence theorem as $\sigma_{n}^{-3}R_{3}$ converges point-wise to zero as
$\sigma_{n}\to 0$. The dominating bound is very similar as for $R_{2,n}$ once
we notice that $b^{\prime}$ and $b^{\prime\prime}$ are both bounded by
Assumption 1 ii).

By Assumption 1 i) the terms $R_{4,n}$ and $R_{5,n}$ are bounded absolutely by
a constant multiplier of
$K(x)\sigma_{n}^{3+H}\max(|w|^{3+H},|w|^{3+\gamma})$, hence
$\lim_{\sigma_{n}\to
0}\sigma_{n}^{-6}\tilde{\mathbb{E}}_{n}[R_{4,n}^{2}]=\lim_{\sigma_{n}\to
0}\sigma_{n}^{-6}\tilde{\mathbb{E}}_{n}[R_{5,n}^{2}]=0$ . The square of
$R_{7,n}$ can be seen to decay faster than $\sigma_{n}^{6}$ by the Dominated
convergence theorem, since $b^{\prime\prime\prime}$ is bounded,
$b^{\prime\prime\prime}(0)=0$ and integrability guaranteed by Assumption 1 i).
Similarly

$R_{6,n}=-\frac{1}{2}\sigma_{n}^{3}w^{3}\phi^{\prime}(x)^{3}\int_{0}^{1}b^{\prime\prime\prime}(-u\phi^{\prime}(x+\sigma_{n}w)\sigma_{n}w)(1-u)^{2}du\\\
-\frac{1}{2}\sigma_{n}^{3}w^{3}\left(\phi^{\prime}(x+\sigma_{n}w)^{3}-\phi^{\prime}(x)^{3}\right)\int_{0}^{1}b^{\prime\prime\prime}(-u\phi^{\prime}(x+\sigma_{n}w)\sigma_{n}w)(1-u)^{2}du\,,$

so again the square of the first part decays faster than $\sigma_{n}^{6}$ by
the Dominated convergence theorem (as for $R_{7,n}$) and the second part is
dominated by a constant multiplier of
$K(x)\sigma_{n}^{3+H}\max(|w|^{3+H},|w|^{3+\gamma)}$ due to Assumption 1 i).

Write

$T(x,w)=w^{3}\left(-\frac{1}{12}\phi^{\prime\prime\prime}(x)+\left(\frac{1}{4}+g^{\prime\prime}(1)\right)\phi^{\prime}(x)\phi^{\prime\prime}(x)\right)-w\left(\frac{1}{2}+g^{\prime\prime}(1)\right)\phi^{\prime}(x)\phi^{\prime\prime}(x)$

and set $R_{n}(x,w)=\sum_{i=1}^{7}R_{i,n}(x,w)$. Clearly

$\tilde{\mathbb{E}}_{n}[\rho_{n}^{2}]=\tilde{\mathbb{E}}_{n}\left[(\sigma^{3}T+R_{n})^{2}\right]=\tilde{\mathbb{E}}_{n}\left[\sigma^{6}T^{2}\right]+2\tilde{\mathbb{E}}_{n}\left[\sigma^{3}TR_{n}\right]+\tilde{\mathbb{E}}_{n}\left[R_{n}^{2}\right].$

Using the inequalities
$\tilde{\mathbb{E}}_{n}\left[TR_{n}\right]^{2}\leq\tilde{\mathbb{E}}_{n}[T^{2}]\tilde{\mathbb{E}}_{n}[R_{n}^{2}]$
and $\tilde{\mathbb{E}}_{n}[R_{n}^{2}]\leq
7\sum_{i=1}^{7}\tilde{\mathbb{E}}_{n}[R_{i,n}^{2}]$ we see that the last two
terms decay faster than $\sigma_{n}^{-6}$. We can now identify the limiting
$\theta^{2}$ using the Dominated convergence theorem. We show that
$\lim_{n\to\infty}\sigma_{n}^{-6}\mathbb{E}[\rho_{n}^{2}]=\lim_{n\to\infty}\tilde{\mathbb{E}}_{n}[T^{2}]=\mathbb{E}[T^{2}]=\theta^{2}$
where the last expectation is with respect to $\pi(x)dx\mu(dw)$. Indeed,

$\tilde{\mathbb{E}}_{n}[T^{2}]=\int_{\mathbb{R}}T^{2}(x,w)\frac{e^{\sigma_{n}w\phi^{\prime}(x)}}{Z_{\sigma_{n}}(x)}\mu(dw)\pi(x)dx\,.$

The integrand converges point-wise to $T^{2}(x,w)$ and is dominated by
$T^{2}(x,w)(1+|\phi^{\prime}(x)|^{\beta})$ which is integrable by Assumption1
i). Expanding the expression for $T$ and integrating with respect to
$\mu(dw)\pi(x)dx$ leads to the claimed form of $\theta^{2}$.

Finally, to finish the proof using Theorem 1 notice that by choosing
$\sigma_{n}$ such that $n\sigma_{n}^{6}\to\ell^{6}$ we get

$\lim_{n\to\infty}n\mathbb{E}[\rho_{n}^{2}]=\lim_{n\to\infty}n\sigma_{n}^{6}\sigma_{n}^{-6}\mathbb{E}[\rho_{n}^{2}]=\ell^{6}\theta^{2}\,.$

We also need to prove that
$\lim_{n\to\infty}\sigma_{n}^{-6}\tilde{\mathbb{E}}_{n}\left[\rho_{n}^{2}1_{\rho_{n}<-\sigma_{n}}\right]=0$.
Indeed,

$\displaystyle\tilde{\mathbb{E}}_{n}\left[\rho_{n}^{2}1_{\rho_{n}<-\sigma_{n}}\right]$
$\displaystyle=\sigma_{n}^{6}\tilde{\mathbb{E}}_{n}\left[T^{2}1_{\rho_{n}<-\sigma_{n}}\right]+2\sigma_{n}^{3}\tilde{\mathbb{E}}_{n}\left[TR_{n}1_{\rho_{n}<-\sigma_{n}}\right]+\tilde{\mathbb{E}}_{n}\left[R_{n}^{2}1_{\rho_{n}<-\sigma_{n}}\right]$
$\displaystyle\leq\sigma_{n}^{6}\tilde{\mathbb{E}}_{n}\left[T^{2}1_{\rho_{n}<-\sigma_{n}}\right]+2\sigma_{n}^{3}\sqrt{\tilde{\mathbb{E}}_{n}[T^{2}]\tilde{\mathbb{E}}_{n}[R_{n}^{2}]}+2\mathbb{E}[R_{n}^{2}].$

We have already established that the last two terms decay faster than
$\sigma_{n}^{-6}$. The first term can be bounded using the Hölder and Markov
inequalities by

$\displaystyle\sigma_{n}^{6}\tilde{\mathbb{E}}_{n}\left[T^{2}1_{\rho_{n}<-\sigma_{n}}\right]$
$\displaystyle\leq\sigma_{n}^{6}\tilde{\mathbb{E}}_{n}[T^{2+\epsilon}]^{\frac{2}{2+\epsilon}}\mathbb{P}[\rho_{n}<-\sigma_{n}]^{\frac{\epsilon}{2+\epsilon}}$
$\displaystyle\leq\sigma_{n}^{6}\tilde{\mathbb{E}}_{n}[T^{2+\epsilon}]^{\frac{2}{2+\epsilon}}\left(\frac{1}{\sigma_{n}^{2}}\tilde{\mathbb{E}}_{n}[\rho_{n}^{2}]\right)^{\frac{\epsilon}{2+\epsilon}}$
$\displaystyle=\sigma_{n}^{6+\frac{4\epsilon}{2+\epsilon}}\tilde{\mathbb{E}}_{n}[T^{2+\epsilon}]^{\frac{2}{2+\epsilon}}\left(\frac{1}{\sigma_{n}^{6}}\tilde{\mathbb{E}}_{n}[\rho_{n}^{2}]\right)^{\frac{\epsilon}{2+\epsilon}}.$

Analogously as before
$\tilde{\mathbb{E}}_{n}[T^{2+\epsilon}]\to\mathbb{E}[T^{2+\epsilon}]<\infty$
by the Dominated convergence theorem and Assumption 1 i) and we have already
established that
$\sigma_{n}^{-6}\tilde{\mathbb{E}}_{n}[\rho_{n}^{2}]\to\theta^{2}$. Together
these establish that

$\lim_{n\to\infty}\sigma_{n}^{-6}\tilde{\mathbb{E}}_{n}\left[\rho_{n}^{2}1_{\rho_{n}<-\sigma_{n}}\right]=0$

which completes the proof.

∎

###### Proof of Theorem 3.

We fix $g$ and $\mu$ that satisfy Assumption 1 and suppress the notation with
respect to them for the rest of the proof. First assume
$\lim_{n\to\infty}n^{1/6}\sigma_{n}=0$, which implies that

$\sigma_{n}^{-2}\mathbb{E}\left[(Y_{n,1}-X_{n,1})^{2}\right]=\sigma_{n}^{-2}\int_{\mathbb{R}}\int_{\mathbb{R}}Z_{\sigma_{n}}^{-1}(x)g(e^{\phi^{\prime}(x)w\sigma_{n}})w^{2}\mu\left(dw\right)\pi(x)dx\to
1$ (18)

as $n\to\infty$. Convergence to one follows by the Dominated convergence
theorem as $\sigma_{n}\to 0$. The dominating bound is provided by Assumption 1
iii). Hence, $\mathbb{E}\left[(Y_{n,1}-X_{n,1})^{2}\right]$ decays as
$\sigma_{n}^{2}$ which is by definition faster than $n^{-1/3}$. Because the
acceptance rate is bounded above by one, we must have that
$\lim_{n\to\infty}n^{1/3}\mathcal{E}_{n}=0$.

Next assume $\lim_{n\to\infty}n^{1/6}\sigma_{n}=\infty$. First we will
establish that

$\mathbb{E}\left[(Y_{n,1}-X_{n,1})^{2}\left(1\wedge
e\right)\left(\sum_{i=1}^{n}\rho_{n,i}\right)\right]\leq
2\mathbb{E}\left[(Y_{n,1}-X_{n,1})^{2}\right]\mathbb{E}\left[\left(1\wedge
e\right)\left(\sum_{i=2}^{n}\rho_{n,i}\right)\right]\,.$ (19)

The statement (without the two) is clear if $\rho_{n,1}$ is negative. When it
is positive then using Lemma 1 with respect to the log-Metropolis–Hastings
random variable $\rho_{n,1}$ and $h(x,y)=(y-x)^{2}$ gives

$\displaystyle\mathbb{E}\left[(Y_{n,1}-X_{n,1})^{2}\left(1\wedge
e\right)\left(\sum_{i=1}^{n}\rho_{n,i}\right)1_{\rho_{n,1}\geq 0}\right]$
$\displaystyle\leq\mathbb{E}\left[(Y_{n,1}-X_{n,1})^{2}e^{\rho_{n,1}}\cdot\left(1\wedge
e\right)\left(\sum_{i=2}^{n}\rho_{n,i}\right)1_{[0,\infty)}\left(\rho_{n,1})\right)\right]$
$\displaystyle=\mathbb{E}\left[(Y_{n,1}-X_{n,1})^{2}\left(1\wedge
e\right)\left(\sum_{i=2}^{n}\rho_{n,i}\right)1_{[0,\infty)}\left(-\rho_{n,1}\right)\right]$
$\displaystyle\leq\mathbb{E}\left[(Y_{n,1}-X_{n,1})^{2}\right]\mathbb{E}\left[\left(1\wedge
e\right)\left(\sum_{i=2}^{n}\rho_{n,i}\right)\right]\,.$

Adding terms corresponding to the sign of $\rho_{n,1}$ implies (19).

We still need to bound the acceptance rate. To do this we split the
probability space regarding the event
$\mathcal{A}_{n}:=\\{\sum_{i=1}^{n}\rho_{n,i}\leq\sum_{i=2}^{n}\mathbb{E}[\rho_{n,i}]/2\\}$.
On the set where this holds we have, using Theorem 1 and Theorem 2

$\displaystyle\mathbb{E}\left[(1\wedge
e)\left(\sum_{i=1}^{n}\rho_{n,i}\right)1_{\mathcal{A}_{n}}\right]$
$\displaystyle\leq\exp\left(\frac{1}{2}\sum_{i=1}^{n}\mathbb{E}[\rho_{n,i}]\right)$
$\displaystyle\leq\exp\left(-\frac{n}{4}\mathbb{E}[\rho^{2}_{n,1}]\left(1-\frac{\mathbb{E}[\rho_{n,1}]+\frac{1}{2}\mathbb{E}[\rho^{2}_{n,1}]}{\mathbb{E}[\rho^{2}_{n,1}]}\right)\right)$
$\displaystyle\leq\exp\left(-\frac{n}{4}\sigma_{n}^{6}\left(\sigma^{-6}_{n}\mathbb{E}[\rho^{2}_{n,1}]\right)\left(1-\frac{\mathbb{E}[\rho_{n,1}]+\frac{1}{2}\mathbb{E}[\rho^{2}_{n,1}]}{\mathbb{E}[\rho^{2}_{n,1}]}\right)\right)$
$\displaystyle\leq\exp\left(-\frac{\theta^{2}}{8}n\sigma_{n}^{6}\right)$

for all large enough $n\in\mathbb{N}$. On the complement $\mathcal{A}_{n}^{c}$
using Markov’s inequality gives

$\displaystyle\mathbb{E}\left[(1\wedge
e)\left(\sum_{i=1}^{n}\rho_{n,i}\right)1_{\mathcal{A}^{c}_{n}}\right]$
$\displaystyle\leq\mathbb{P}\left[\mathcal{A}_{n}^{c}\right]$
$\displaystyle=\mathbb{P}\left[\sum_{i=1}^{n}\left(\rho_{n,i}-\mathbb{E}[\rho_{n,i}]\right)>-\frac{1}{2}\sum_{i=1}^{n}\mathbb{E}[\rho_{n,i}]\right]$
$\displaystyle\leq\frac{4\text{Var}[\rho_{n,1}^{2}]}{n\mathbb{E}[\rho_{n,1}]^{2}}$
$\displaystyle\leq\frac{16}{n\mathbb{E}[\rho_{n,1}^{2}]\left(1-4\frac{\mathbb{E}[\rho_{n,1}]+\frac{1}{2}\mathbb{E}[\rho_{n,1}^{2}]}{\mathbb{E}[\rho_{n,1}^{2}]}\right)}$
$\displaystyle\leq\frac{32}{\theta^{2}n\sigma_{n}^{6}}$

for all large enough $n\in\mathbb{N}$.

Together with (18) and (19) these bounds imply

$\mathbb{E}\left[(Y_{n,1}-X_{n,1})^{2}\left(1\wedge
e\right)\left(\sum_{i=1}^{n}\rho_{n,i}\right)\right]\leq
2\sigma_{n}^{2}\left(e^{-\frac{\theta^{2}}{8}n\sigma_{n}^{6}}+\frac{32}{\theta^{2}n\sigma_{n}^{6}}\right)\leq
C\frac{1}{n^{1/3}}\frac{1}{n^{2/3}\sigma_{n}^{4}}$

for an appropriate constant $C>0$. Since $n\sigma_{n}^{6}\to\infty$, this
decay rate is faster than $n^{-1/3}$, meaning $\mathcal{E}_{n}$ also decays
faster than $n^{-1/3}$.

Now assume, $\lim_{n\to\infty}n^{1/6}\sigma_{n,\ell}=\ell$. Splitting the
expectation gives

$\mathbb{E}\left[(Y_{n,1}-X_{n,1})^{2}\left(1\wedge
e\right)\left(\sum_{i=1}^{n}\rho_{n,i}\right)\right]=\mathbb{E}\left[(Y_{n,1}-X_{n,1})^{2}\right]\left[\left(1\wedge
e\right)\left(\sum_{i=2}^{n}\rho_{n,i}\right)\right]\\\
+\mathbb{E}\left[(Y_{n,1}-X_{n,1})^{2}\left(\left(1\wedge
e\right)\left(\sum_{i=1}^{n}\rho_{n,i}\right)-\left(1\wedge
e\right)\left(\sum_{i=2}^{n}\rho_{n,i}\right)\right)\right]\,.$

For the first term on the right-hand side we have by (18) and Theorem 2 that

$n^{1/3}\mathbb{E}\left[(Y^{(n)}_{1}-X_{1})^{2}\right]\left[\left(1\wedge
e\right)\left(\sum_{i=2}^{n}\rho_{\sigma_{n,\ell},i}\right)\right]\to h(\ell)$

as $n\to\infty$. The second term vanishes. To see this note that by Theorem 2
and the fact that the function $t\mapsto 1\wedge e^{t}$ is $1$-Lipschitz and
bounded imply that the expression

$\left(\left(1\wedge
e\right)\left(\sum_{i=1}^{n}\rho_{\sigma_{n,\ell},i}\right)-\left(1\wedge
e\right)\left(\sum_{i=2}^{n}\rho_{\sigma_{n,\ell},i}\right)\right)$

converges to zero in the $L^{2}$ sense, and therefore also in probability. On
the other hand (18) and Assumption 1 iii) imply that the random variables
$n^{1/3}(Y_{n,1}-X_{n,1})^{2}$ have an integrable dominating bound. Their
product is therefore uniformly integrable and converges to zero in
probability, meaning it is also true that

$n^{1/3}\mathbb{E}\left[(Y_{n,1}-X_{n,1})^{2}\left(\left(1\wedge
e\right)\left(\sum_{i=1}^{n}\rho_{n,i}\right)-\left(1\wedge
e\right)\left(\sum_{i=2}^{n}\rho_{n,i}\right)\right)\right]\to 0$

as $n\to\infty$.

Finally we will optimize over the choice of $\ell$. The function
$h(\ell)=2\ell^{2}\Phi(-\ell^{3}\theta/2)$ is smooth in $\ell$ and converges
to zero both when $\ell\to 0$ and when $\ell\to\infty$. Hence, its maximum is
attained at a stationary point. Setting $s=\ell^{3}\theta/2$ we can find the
stationary points of $s\mapsto
2\left(\frac{2}{\theta}\right)^{2/3}s^{2/3}\Phi(-s)$. This corresponds to
finding the solution of the equation $2/3=s\varphi_{N}(-s)/\Phi(-s)$, where
$\varphi_{N}$ denotes the standard Gaussian probability density function.
There exists a unique solution $s^{*}$, as the function $s\mapsto
s\varphi_{N}(-s)/\Phi(-s)$ is strictly increasing, meaning $h(\ell)$ attains
its maximal value at a specific value $\ell^{*}$ satisfying
$s^{*}=(\ell^{*})^{3}\theta/2$, corresponding to an average acceptance rate of
$2\Phi(-s^{*})$, which turns out to numerically equal to $57.4\%$ to three
decimal places. This is also implies that

$h(\ell^{*})=\theta^{-2/3}\cdot 2^{5/3}(s^{*})^{2/3}\Phi(-s^{*})\approx
0.651637\times\theta^{-2/3},$

which completes the proof. ∎

##  Appendix D Technical results

######  Lemma 1.

Let $\rho$ be a log-Metropolis–Hastings random variable associated with a
probability measure $\pi$ and a Markov kernel $Q$ on
$(\mathbb{X},\mathcal{F})$. Let $f:\mathbb{X}\to\mathbb{R}$ and
$h:\mathbb{X}\times\mathbb{X}\to\mathbb{R}$ be such that $h(X,Y)=h(Y,X)$. Then
the following are true:

  1. (i)

$-\rho(X,Y)=\rho(Y,X)\,.$

  2. (ii)

If the integrals are finite, then setting $h=h(X,Y)$ and $\rho=\rho(X,Y)$

$\mathbb{E}\left[hf(\rho e^{\rho})\right]=\mathbb{E}\left[hf(-\rho)\right]\,.$

###### Proof of Lemma 1.

The proof is given Proposition 3 of [21], with the only minor difference that
we carry a symmetric function $h$ through the entire derivation. ∎

## References

  * [1] Christophe Andrieu and Johannes Thoms.  A tutorial on adaptive mcmc.  Statistics and computing, 18(4):343–373, 2008.
  * [2] Av A Barker.  Monte carlo calculations of the radial distribution functions for a proton-electron plasma.  Australian Journal of Physics, 18(2):119–134, 1965.
  * [3] Bob Carpenter, Andrew Gelman, Matthew D Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Marcus Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell.  Stan: A probabilistic programming language.  Journal of statistical software, 76(1):1–32, 2017.
  * [4] Michael CH Choi.  Metropolis–Hastings reversiblizations of non-reversible Markov chains.  Stochastic Processes and their Applications, 130(2):1041–1073, 2020\.
  * [5] Simon Duane, Anthony D Kennedy, Brian J Pendleton, and Duncan Roweth.  Hybrid monte carlo.  Physics letters B, 195(2):216–222, 1987.
  * [6] Paul Fearnhead, Joris Bierkens, Murray Pollock, and Gareth O Roberts.  Piecewise deterministic markov processes for continuous-time monte carlo.  Statistical Science, 33(3):386–412, 2018.
  * [7] Max Hird, Samuel Livingstone, and Giacomo Zanella.  A fresh take on ‘Barker dynamics’ for MCMC.  arXiv preprint arXiv:2012.09731, 2020.
  * [8] Olav Kallenberg.  Foundations of modern probability, volume 2.  Springer, 1997.
  * [9] Samuel Livingstone and Giacomo Zanella.  The Barker proposal: combining robustness and efficiency in gradient-based MCMC.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), in press, 2021.
  * [10] Peter Neal and Gareth Roberts.  Optimal scaling of random walk metropolis algorithms with non-gaussian proposals.  Methodology and Computing in Applied Probability, 13(3):583–601, 2011.
  * [11] Radford M Neal.  Slice sampling.  The annals of statistics, 31(3):705–767, 2003.
  * [12] Samuel Power and Jacob Vorstrup Goldman.  Accelerated Sampling on Discrete Spaces with Non-Reversible Markov Processes.  arXiv preprint arXiv:1912.04681, 2019.
  * [13] Gareth O Roberts, Andrew Gelman, and Walter R Gilks.  Weak convergence and optimal scaling of random walk Metropolis algorithms.  The annals of applied probability, 7(1):110–120, 1997.
  * [14] Gareth O Roberts and Jeffrey S Rosenthal.  Optimal scaling of discrete approximations to Langevin diffusions.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 60(1):255–268, 1998.
  * [15] Gareth O Roberts and Jeffrey S Rosenthal.  Optimal scaling for various Metropolis-Hastings algorithms.  Statistical science, 16(4):351–367, 2001.
  * [16] Gareth O Roberts and Jeffrey S Rosenthal.  Examples of adaptive mcmc.  Journal of computational and graphical statistics, 18(2):349–367, 2009.
  * [17] Gareth O Roberts and Richard L Tweedie.  Exponential convergence of Langevin distributions and their discrete approximations.  Bernoulli, pages 341–363, 1996.
  * [18] Jeffrey S Rosenthal.  Amcmc: An r interface for adaptive mcmc.  Computational Statistics & Data Analysis, 51(12):5467–5470, 2007\.
  * [19] Chris Sherlock, Paul Fearnhead, and Gareth O Roberts.  The random walk metropolis: linking theory and practice through a case study.  Statistical Science, 25(2):172–190, 2010.
  * [20] Luke Tierney.  A note on Metropolis-Hastings kernels for general state spaces.  Annals of applied probability, pages 1–9, 1998.
  * [21] Jure Vogrinc and Wilfrid S Kendall.  Counterexamples for optimal scaling of Metropolis–Hastings chains with rough target densities.  The Annals of Applied Probability, 31(2):972–1019, 2021.
  * [22] Ziheng Yang and Carlos E Rodríguez.  Searching for efficient markov chain monte carlo proposal kernels.  Proceedings of the National Academy of Sciences, 110(48):19307–19312, 2013.
  * [23] Giacomo Zanella.  Informed proposals for local MCMC in discrete spaces.  Journal of the American Statistical Association, 115(530):852–865, 2020.
  * [24] Giacomo Zanella, Mylene Bédard, and Wilfrid S Kendall.  A Dirichlet form approach to MCMC optimal scaling.  Stochastic Processes and their Applications, 127(12):4053–4082, 2017.