# Metropolis–Hastings with Scalable Subsampling

Estevão Prado, Christopher Nemeth and Chris Sherlock

School of Mathematical Sciences, Lancaster University, UK  The author was
supported by EPSRC grants EP/V022636/1 and EP/Y028783/1.

###### Abstract

The Metropolis–Hastings (MH) algorithm is one of the most widely used Markov
Chain Monte Carlo schemes for generating samples from Bayesian posterior
distributions. The algorithm is asymptotically exact, flexible and easy to
implement. However, in the context of Bayesian inference for large datasets,
evaluating the likelihood on the full data for thousands of iterations until
convergence can be prohibitively expensive. This paper introduces a new
subsample MH algorithm that satisfies detailed balance with respect to the
target posterior and utilises control variates to enable exact, efficient
Bayesian inference on datasets with large numbers of observations. Through
theoretical results, simulation experiments and real-world applications on
certain generalised linear models, we demonstrate that our method requires
substantially smaller subsamples and is computationally more efficient than
the standard MH algorithm and other exact subsample MH algorithms.

Keywords: Markov Chain Monte Carlo, scalable Metropolis–Hastings algorithm,
control variates, subsampling, big data.

##  1 Introduction

In the Bayesian framework, the posterior distribution of the model parameters
contains all information needed for inference, but even with simple models,
the posterior is generally known only up to a constant of proportionality.
Markov Chain Monte Carlo (MCMC; Robert et al.,, 1999; Brooks et al.,, 2011) is
a class of simulation methods commonly used in Bayesian statistics to sample
from posterior distributions. At the core of MCMC methods is the
Metropolis–Hastings (MH) algorithm, which proposes at each iteration a new
parameter value and then decides whether or not to accept it. The acceptance
decision requires the evaluation of the full (unnormalised) posterior density
evaluated at the proposed parameter and this in turn requires the evaluation
of every single term in the likelihood. In this big-data era, where datasets
consist of millions or even billions of data points, each likelihood
evaluation is computationally expensive and it becomes infeasible to repeat
this at every iteration, as required by the standard MH algorithm.

Optimisation-based algorithms, such as variational and Laplace approximations
(Raj et al.,, 2014; Yogatama et al.,, 2014; Andersen et al.,, 2018; Winter et
al.,, 2024), are still feasible, but these are inexact. In this context, a new
class of MCMC methods has gained attention more recently: scalable MCMC
algorithms; see Fearnhead et al., (2024) for a recent review. In an attempt to
alleviate the computational burden of performing traditional MCMC sampling on
large datasets, some of the works in the scalable MCMC literature propose
splitting the data into partitions and then, in parallel, performing posterior
sampling on each partition across multiple CPUs. This approach is called
‘divide-and-conquer’, as the data are divided into multiple disjoint
subsamples, and the final posterior is obtained by combining the subposteriors
from each subsample. One of the challenges related to this approach is to
efficiently combine the subposteriors in cases where the target distribution
is non-Gaussian; see, for example, Neiswanger et al., (2013); Scott et al.,
(2016); Nemeth and Sherlock, (2018); Vyner et al., (2023); Trojan et al.,
(2024).

Another avenue to reduce the computational cost of running MCMC algorithms on
datasets with large number of observations utilises data subsampling
(Maclaurin and Adams,, 2014; Korattikara et al.,, 2014; Bardenet et al.,,
2014; Liu et al.,, 2015; Bardenet et al.,, 2017; Yuan and Wang,, 2024). In
contrast to the divide-and-conquer approach, subsampling-based MCMC methods
speed up posterior sampling by running a single MCMC algorithm and using a
random subsample of the data at each iteration. For instance, Quiroz et al.,
2018a  propose a Metropolis-based MCMC algorithm that builds upon the pseudo-
marginal approach (Beaumont,, 2003; Andrieu and Roberts,, 2009) and control
variates (Ripley,, 2009). The MH acceptance ratio is replaced with an
approximation based on subsamples. To control the error introduced by this
approximation, the MH ratio needs to be bias-corrected. Consequently, the
algorithm is no longer exact, as the invariant distribution of the Markov
chain is a perturbed version of the true posterior.

There are also subsample MCMC methods that are constructed in such a way that
the invariant distribution of the resulting Markov chain is the true posterior
distribution. All of these methods require that the full likelihood is a
product of terms (typically, one term for each data point) and require some
assumptions on the individual terms in the product. Subject to these
constraints, only a random subsample of the terms needs to be evaluated at any
given iteration of the MCMC algorithm. For instance, Maclaurin and Adams,
(2014) assume that for each observation there exists a strictly positive lower
bound on the log-likelihood that accumulates tractably when summed over many
observations; it then augments the parameter space using Bernoulli latent
variables. Cornish et al., (2019) and Zhang et al., (2020) are less
restrictive and do not augment the state space. Both of these algorithms
replace the Metropolis–Hastings acceptance ratio with an alternative that uses
a small number of subsamples and, hence, has a high variance. Cornish et al.,
(2019) uses control variates to reduce this variance, whilst Zhang et al.,
(2020) incorporates an additional parameter to mitigate the variance’s
deleterious effect on the acceptance rate.

In this paper, we propose a new exact subsample MH algorithm that makes use of
control variates and subsampling and produces a Markov chain that satisfies
detailed balance with respect to the true posterior. The criterion (3) looks
at the difference in log-likelihoods at the current and proposed values. It
requires that the discrepancy between this difference and the estimate of the
difference suggested by the control variates can be bounded by a product of a
function of the parameters and a function of the data value. We establish
general bounds that are typically much tighter than those considered in
previous works, especially when the dimension, $d$, of the parameter vector is
moderate or large. Specialising to regression models, we derive specific,
tight bounds for three generalised linear models (GLMs) including the logistic
regression model considered in previous works. More broadly, we establish that
our algorithm is optimal within a particular class of subsample algorithms.
All of this leads to a substantial speed-up in the posterior sampling. We also
conduct an asymptotic analysis comparing our proposed bounds with previous
approaches, demonstrating theoretically and empirically that the computational
cost of using our bounds is at least a factor of $d^{1/2}$ better than
previous bounds.

This paper is organised as follows. Section 2 introduces our algorithm and
illustrates why it targets the correct posterior. Section 3 discusses related
works and issues with the existing subsample-based MH algorithms. Section 4
presents the general forms of the much tighter bounds we propose for the log-
likelihood differences, as well as theoretical results on the exactness and
optimality of our proposed algorithm and guidance on its tuning. The proofs of
the results in this section are deferred to Appendix A. Sections 5 and 6
present simulation experiments and real-world applications using the logistic,
probit and Poisson regression models, where our proposed method is compared
against the equivalent Metroplis–Hastings algorithm and other subsample-based
MH algorithms. Finally, Section 7 concludes with a discussion. Additional
simulation experiments, asymptotic analyses comparing the new and existing
bounds and results related to the acceptance rate of the proposed algorithms
are presented in Appendices B, C and D.

##  2 Metropolis–Hastings with Scalable Subsampling

We are interested in sampling from a posterior distribution with density
$\pi(\theta|y)$, where $\theta\in\mathcal{X}\subseteq\mathbb{R}^{d}$ denotes
the parameter vector and $y=(y_{1},\ldots,y_{n})$ are the observed data.
Assuming a prior $p(\theta)$, and data that are conditionally independent with
log likelihood terms $\ell_{i}(\theta):=\log p(y_{i}|\theta)$, the posterior
density is

$\pi(\theta|y)\propto p(\theta)\prod_{i=1}^{n}p(y_{i}|\theta)\equiv
p(\theta)p(y_{1:n}|\theta)=p(\theta)\exp\left[\sum_{i=1}^{n}\ell_{i}(\theta)\right].$

Given a current value, $\theta$, and a proposal density,
$q(\theta^{\prime}|\theta)$, the Metropolis–Hastings acceptance probability is
$1\wedge[q(\theta|\theta^{\prime})\pi(\theta^{\prime}|y)/\\{q(\theta^{\prime}|\theta)\pi(\theta|y)\\}]$.
The contribution of the computationally expensive likelihood to this ratio is

$\exp\left[\sum_{i=1}^{n}\\{\ell_{i}(\theta^{\prime})-\ell_{i}(\theta)\\}\right].$

To allow us to subsample terms from this sum whilst controlling the variance
of the subsequent estimator, we create a deterministic approximation, called a
control variate, for each difference,
$\ell_{i}(\theta^{\prime})-\ell_{i}(\theta)$. Let
$\widehat{\theta}\in\mathcal{X}$ be (ideally, but for strict correctness of
the algorithm, not necessarily) close to the posterior mode. We choose control
variates to be the approximations to
$\ell_{i}(\theta^{\prime})-\ell_{i}(\theta)$ that arise from first- and
second-order Taylor expansion of $\ell_{i}$ about $\widehat{\theta}$:

$\displaystyle r_{i}^{(1)}(\theta,\theta^{\prime};\widehat{\theta})$
$\displaystyle:=(\theta^{\prime}-\theta)^{\top}g_{i}(\widehat{\theta}),$ (1)
$\displaystyle r_{i}^{(2)}(\theta,\theta^{\prime};\widehat{\theta})$
$\displaystyle:=(\theta^{\prime}-\theta)^{\top}g_{i}(\widehat{\theta})+(\theta^{\prime}-\theta)^{\top}H_{i}(\widehat{\theta})\left\\{\frac{1}{2}(\theta^{\prime}+\theta)-\widehat{\theta}\right\\},$
(2)

where $g_{i}(\theta)=\nabla\ell_{i}(\theta)$ is the gradient of the log-
likelihood and $H_{i}(\theta)=\nabla\nabla^{\top}\ell_{i}(\theta)$ is the
Hessian. Whichever control variate, $k\in\\{1,2\\}$, is used, our algorithm
requires the existence of bounds $M^{(k)}(\theta,\theta^{\prime})\geq 0$ and
$c^{(k)}_{1},\dots,c^{(k)}_{n}>0$ such that
$M^{(k)}(\theta,\theta^{\prime})=M^{(k)}(\theta^{\prime},\theta)$ for all
$\theta,\theta^{\prime}\in\mathcal{X}$ and

$|\ell_{i}(\theta^{\prime})-\ell_{i}(\theta)-r_{i}^{(k)}(\theta,\theta^{\prime};\widehat{\theta})|\leq
c^{(k)}_{i}M^{(k)}(\theta,\theta^{\prime}).$ (3)

In the following, for notational convenience, we drop the superscript $(k)$
from $r^{(k)}_{i}(\theta,\theta^{\prime};\widehat{\theta})$, $c^{(k)}_{i}$ and
$M^{(k)}(\theta,\theta^{\prime})$ and use the notation
$r_{i}(\theta,\theta^{\prime};\widehat{\theta})$, $c_{i}$,
$M(\theta,\theta^{\prime})$, respectively. Such bounds for logistic regression
and robust linear regression models when $k=1,2$ are given in Cornish et al.,
(2019). In Section 4.2, we provide much-improved bounds for logistic, probit
and a form of Poisson regression. We define the difference between the control
variate and the truth to be

$\Delta_{i}:=r_{i}(\theta,\theta^{\prime};\widehat{\theta})-\\{\ell_{i}(\theta^{\prime})-\ell_{i}(\theta)\\},$
(4)

and for $\gamma\in[0,1]$ and we introduce the functions $\phi_{i}$ and
$\phi_{i}^{\prime}$ for $i=1,\ldots,n$ which are chosen to satisfy the
following properties for all $\theta,\theta^{\prime}\in\mathcal{X}$:

  * F1:

$0\leq\phi_{i}\leq c_{i}M(\theta,\theta^{\prime})$.

  * F2:

$0\leq\phi^{\prime}_{i}\leq c_{i}M(\theta,\theta^{\prime})$.

  * F3:

$\phi_{i}-\phi_{i}^{\prime}=\Delta_{i}$.

The assumptions on $\phi_{i}$ and $\phi_{i}^{\prime}$ ensure that these are
bounded, non-negative functions whose difference is exactly the quantity of
interest and which can be used as the expectations of Poisson random variables
that can be sampled via Poisson thinning. In this paper, we choose

$\displaystyle\phi_{i}$
$\displaystyle=\gamma\max[0,\Delta_{i}]+(1-\gamma)\left\\{c_{i}M(\theta,\theta^{\prime})+\min[0,\Delta_{i}]\right\\}.$
(5)

We also define the equivalent quantity when
$\theta\leftrightarrow\theta^{\prime}$ are interchanged:

$\displaystyle\phi^{\prime}_{i}$
$\displaystyle=\gamma\max[0,\Delta^{\prime}_{i}]+(1-\gamma)\left\\{c_{i}M(\theta^{\prime},\theta)+\min[0,\Delta^{\prime}_{i}]\right\\},$
(6)

where
$\Delta_{i}^{\prime}=r_{i}(\theta^{\prime},\theta;\widehat{\theta})-\\{\ell_{i}(\theta)-\ell_{i}(\theta^{\prime})\\}=-\Delta_{i}$,
since
$r_{i}(\theta^{\prime},\theta;\widehat{\theta})=-r_{i}(\theta,\theta^{\prime};\widehat{\theta})$
for both $r^{(1)}$ and $r^{(2)}$. From (3) and (4), both the multiplicand of
$\gamma$ and that of $1-\gamma$ in (5) are between $0$ and
$c_{i}M(\theta,\theta^{\prime})$, proving F1, and, by an analogous argument
since $M(\theta,\theta^{\prime})=M(\theta^{\prime},\theta)$, F2. Since
$\max[0,-\Delta_{i}]=-\min[0,\Delta_{i}]$ and
$\min[0,-\Delta_{i}]=-\max[0,\Delta_{i}]$, F3 holds because
$\max[0,\Delta_{i}]+\min[0,\Delta_{i}]=0+\Delta_{i}=\Delta_{i}$. Whilst our
algorithm is valid for any $\gamma\in[0,1]$, we show in Section 4.1 that it is
always most efficient when $\gamma=0$.

The non-negativity in F1 and F2 permits us to complete the set-up for our
algorithm through the independent auxiliary variables

$S_{i}\sim\mathsf{Pois}\left(\phi_{i}\right),~{}~{}~{}i=1,\dots,n.$ (7)

$S_{i}$ will be the number of times the $i$th likelihood term appears in our
acceptance ratio. When $\phi_{i}$ is sufficiently small, this is very likely
to be $0$. From F3,

$\frac{p(s_{i}|\theta^{\prime},\theta)}{p(s_{i}|\theta,\theta^{\prime})}=\frac{\exp(-\phi_{i}^{\prime}){\phi_{i}^{\prime}}^{s_{i}}/s_{i}!}{\exp(-\phi_{i}){\phi_{i}}^{s_{i}}/s_{i}!}=\exp(\Delta_{i})\frac{{\phi_{i}^{\prime}}^{s_{i}}}{{\phi_{i}}^{s_{i}}}.$

Let $s:=(s_{1},\ldots,s_{n})$ and write
$p(s|\theta,\theta^{\prime}):=\prod_{i=1}^{n}p(s_{i}|\theta,\theta^{\prime})$.
If the Metropolis–Hastings proposal density is $q(\theta^{\prime}|\theta)$,
detailed balance holds with respect to $\pi(\theta)$ if the probability of
accepting $\theta^{\prime}$ satisfies

$\pi(\theta)q(\theta^{\prime}|\theta)p(s|\theta,\theta^{\prime})\alpha(\theta,\theta^{\prime};s)=\pi(\theta^{\prime})q(\theta|\theta^{\prime})p(s|\theta^{\prime},\theta)\alpha(\theta^{\prime},\theta;s).$
(8)

This holds, for example, when the acceptance probability is

$\displaystyle\alpha(\theta,\theta^{\prime};s)$
$\displaystyle=1\wedge\frac{\pi(\theta^{\prime})q(\theta|\theta^{\prime})p(s|\theta^{\prime},\theta)}{\pi(\theta)q(\theta^{\prime}|\theta)p(s|\theta,\theta^{\prime})}$
$\displaystyle=1\wedge\frac{q(\theta^{\prime}|\theta)p(\theta^{\prime})}{q(\theta|\theta^{\prime})p(\theta^{\prime})}\exp\left[\sum_{i=1}^{n}\ell_{i}(\theta^{\prime})-\ell_{i}(\theta)\right]\exp\left[\sum_{i=1}^{n}\Delta_{i}\right]\prod_{i=1}^{n}\left\\{\frac{\phi_{i}^{\prime}}{\phi_{i}}\right\\}^{s_{i}}$
$\displaystyle=1\wedge\frac{q(\theta^{\prime}|\theta)p(\theta^{\prime})}{q(\theta|\theta^{\prime})p(\theta^{\prime})}\exp\left[\sum_{i=1}^{n}r_{i}(\theta,\theta^{\prime};\widehat{\theta})\right]\prod_{i=1}^{n}\left\\{\frac{\phi_{i}^{\prime}}{\phi_{i}}\right\\}^{s_{i}}$
(9)
$\displaystyle=1\wedge\frac{q(\theta^{\prime}|\theta)p(\theta^{\prime})}{q(\theta|\theta^{\prime})p(\theta)}\exp\left[\sum_{i=1}^{n}r_{i}(\theta,\theta^{\prime};\widehat{\theta})\right]\exp\left[\sum_{i=1}^{n}s_{i}\left\\{\log\left(\phi_{i}^{\prime}\right)-\log\left(\phi_{i}\right)\right\\}\right].$

Sampling random Poisson variables $S_{1},\dots,S_{n}$ independently would be
an $\mathcal{O}(n)$ operation that is performed at every Metropolis-Hastings
iteration, and therefore not practical. To avoid such costly simulations, we
use Poisson thinning (Lewis and Shedler,, 1979; Bouchard-Côté et al.,, 2018;
Bierkens et al.,, 2019; Cornish et al.,, 2019; Zhang et al.,, 2020) to
efficiently sample the $S_{i}$ at each iteration. Defining
$C:=\sum_{i=1}^{n}c_{i}$, we proceed as follows:

  1. 1.

Sample $B\sim\mathsf{Pois}(CM(\theta,\theta^{\prime}))$.

  2. 2.

Conditional on $B$, sample
$(T_{1},\dots,T_{B})\sim\mathsf{Multinomial}(B;c_{1}/C,\dots,c_{i}/C,\dots,c_{n}/C)$.
The set $\\{T_{j}\\}_{j=1}^{B}$ contains the indices (possibly repeated) of
data points that could be included in the likelihood evaluation.

  3. 3.

For each $j\in 1,\dots,B$, include $T_{j}$ in $\mathcal{I}$ with a probability
$\phi_{T_{j}}/\\{c_{T_{j}}M(\theta,\theta^{\prime})\\}$. For each
$i=1,\dots,n$, $S_{i}$ is the number of times the index $i$ occurs in
$\mathcal{I}$.

Marginalising Steps 1 and 2 over $B$, the number of times each index $i$
appears in the set $\\{T_{j}\\}_{j=1}^{B}$ has a
$\mathsf{Pois}(c_{i}M(\theta,\theta^{\prime}))$ distribution, independent of
the frequencies of all other indices. Step 3 then thins down to the required
distribution. Calculating $C$ has a single $\mathcal{O}(n)$ cost at the start
of the algorithm. Similarly, after an $\mathcal{O}(n)$ set up cost, Step 2 has
a per-iteration cost of $\mathcal{O}(B)$ since it assigns an index to each of
the $B$ values of $T_{j}$ that are non-zero. As we will see in Appendix B, $B$
is, at most $\mathcal{O}(1)$. Lastly, Step 3 also has a per-iteration cost of
$\mathcal{O}(B)$. The final set up cost of the algorithm is in finding a
suitable $\widehat{\theta}$, which is also $\mathcal{O}(n)$.

Algorithm 1 presents our algorithm, Metropolis–Hasting with Scalable
Subsampling (MH-SS), using a delayed-acceptance (DA) formulation (Christen and
Fox,, 2005; Cui et al.,, 2011; Golightly et al.,, 2015; Sherlock et al.,,
2017; Quiroz et al., 2018b, ) where the MH acceptance ratio is now
$\alpha(\theta,\theta^{\prime};s)=\alpha_{1}(\theta,\theta^{\prime})\alpha_{2}(\theta,\theta^{\prime};s)$.
In essence, this uses the approximation from the Taylor expansion about
$\widehat{\theta}$ to pre-screen each proposal, reducing the computational
cost still further, with no additional effort. Here
$\alpha(\theta,\theta^{\prime};s)$ in (9) is replaced with
$\alpha_{1}(\theta,\theta^{\prime})\alpha_{2}(\theta,\theta^{\prime};s)$,
where

$\displaystyle\alpha_{1}\left(\theta,\theta^{\prime}\right)$
$\displaystyle=1\wedge\frac{p(\theta^{\prime})q(\theta|\theta^{\prime})}{p(\theta)q(\theta^{\prime}|\theta)}\exp\left[\sum_{i=1}^{n}r_{i}(\theta,\theta^{\prime};\widehat{\theta})\right],$
(10) $\displaystyle\alpha_{2,MHSS}\left(\theta,\theta^{\prime}\right)$
$\displaystyle=1\wedge\exp\left[\sum_{i\in\mathcal{I}}s_{i}\left\\{\log\left(\phi_{i}^{\prime}\right)-\log\left(\phi_{i}\right)\right\\}\right].$
(11)

The detailed balance condition (8) is still satisfed and so the algorithm
continues to target $\pi(\theta)$.

If the proposal is rejected at the pre-screening stage, then the algorithm
proceeds to the next iteration and a new proposal. Though
$\alpha_{1}(\theta,\theta^{\prime})$ involves all $n$ observations, it takes
negligible CPU time to calculate because the components in
$r_{i}(\theta,\theta^{\prime};\widehat{\theta})$ that require most of the
computational effort (i.e., $\sum_{i=1}^{n}g_{i}(\widehat{\theta})$ and, for
second-order control variates, also $\sum_{i=1}^{n}H_{i}(\widehat{\theta})$)
are pre-computed before running the algorithm. This brings down the
computational cost per iteration in evaluating the control variates to
$\mathcal{O}(d)$ for first-order and $\mathcal{O}(d^{2})$ for second-order.
When the pre-screening step has been passed successfully and
$\mathbb{E}(B)=CM(\theta,\theta^{\prime})\geq n$, then subsampling is not
used, and the remainder of a step of the (delayed-acceptance)
Metropolis–Hastings is performed on the full data. Otherwise, on passing the
pre-screening, the proposal $\theta^{\prime}$ is accepted with probability
$\alpha_{2,MHSS}(\theta,\theta^{\prime};s)$, which is calculated taking into
account at most $B$ observations, those with non-zero $s_{i}$.

Algorithm 1 Metropolis–Hastings with Scalable Subsampling (MH-SS)

1:Initialise: $\theta$.

2:for ($t=1$ to $T$) do

3: Sample $\theta^{\prime}\sim q(.|\theta)$.

4: Sample $u\sim\mathsf{Uniform}\left(0,1\right)$ and compute the DA-MH Stage
1 acceptance probability $\alpha_{1}(\theta,\theta^{\prime})$ in (10).

5: if $u\leq\alpha_{1}\left(\theta,\theta^{\prime}\right)$ then

6: if $CM(\theta,\theta^{\prime})\geq n$ then $\triangleright$
Metropolis–Hastings step on the full data

7: Compute the DA Metropolis–Hastings Stage 2 acceptance probability:
$\alpha_{2,MH}(\theta,\theta^{\prime})=1\wedge\frac{p(y_{1:n}|\theta^{\prime})}{p(y_{1:n}|\theta)}\exp\left[-\sum_{i=1}^{n}r_{i}(\theta,\theta^{\prime};\widehat{\theta})\right].$

8: Sample $u\sim\mathsf{Uniform}\left(0,1\right)$: if
$u\leq\alpha_{2,MH}\left(\theta,\theta^{\prime}\right)$, set
$\theta=\theta^{\prime}$.

9: else$\triangleright$ Metropolis–Hastings step on subsamples

10: Sample $B\sim\mbox{Poisson}(CM(\theta,\theta^{\prime}))$

11: $\mathcal{I}\leftarrow\emptyset$.

12: for ($b=1$ to $B$) do $\triangleright$ Subsample creation step

13: Sample $i_{b}$ such that $P(i_{b}=i)=c_{i}/C$ for $i=1,\ldots,n$.

14: Add $i_{b}$ to $\mathcal{I}$ with probability
$\frac{\phi_{i_{b}}(\theta)}{c_{i_{b}}M(\theta,\theta^{\prime})}$.

15: end for

16: Compute the DA MH-SS Stage 2 acceptance probability
$\alpha_{2,MHSS}\left(\theta,\theta^{\prime}\right)$ in (11).

17: Sample $u\sim\mathsf{Uniform}\left(0,1\right)$: if
$u\leq\alpha_{2,MHSS}\left(\theta,\theta^{\prime}\right)$, set
$\theta=\theta^{\prime}$.

18: end if

19: end if

20:end for

21:Output: Samples from the posterior distribution $\pi(\theta)$.

For random walk-based proposals, $q$, such as
$\theta^{\prime}\sim\mathsf{N}(\theta,V)$ for some variance matrix $V$, or for
the Crank-Nicholson proposal with pre-conditioning based on $\widehat{\pi}$
that is suggested in Cornish et al., (2019), Algorithm 1 can be applied as
written. For proposals, such as that of the Metropolis-adjusted Langevin
algorithm (MALA), that use $\nabla\log\pi(\theta)$, naively, the costs of
proposing $\theta^{\prime}$ and of evaluating both $q(\theta^{\prime}|\theta)$
and $q(\theta|\theta^{\prime})$ are all $\mathcal{O}(n)$. In this case a
cheap, unbiased low-variance estimator of $\nabla\log\pi$ at both the proposed
and current parameter values can be obtained by subsampling and using control
variates; for example:

$\widehat{\nabla}\log\pi(\theta)=\sum_{i=1}^{n}g_{i}(\widehat{\theta})+\frac{n}{m}\sum_{j=1}^{m}\\{g_{i_{j}}(\theta)-g_{i_{j}}(\widehat{\theta})\\},$

where $i_{1},\dots,i_{m}$ are the indices of the random subsample chosen at
this iteration for estimating the gradients and
$\sum_{i=1}^{n}g_{i}(\widehat{\theta})$ is pre-computed.

##  3 Issues with existing subsampling MH algorithms

Maclaurin and Adams, (2014) propose the Firefly Monte Carlo (FlyMC) algorithm,
which augments the parameter space by introducing latent Bernoulli random
variables. The Bernoulli success probabilities are defined so that the
$\theta$ marginal of the invariant distribution of the Markov chain is exactly
$\pi$. Crucially, the evaluation of the joint posterior only requires
likelihood terms corresponding to the non-zero Bernoulli variables. The
biggest limitation of this algorithm is the requirement of a global
exponential-family lower bound that is needed on each likelihood term
$B_{i}(\theta)\leq\ell_{i}(\theta)$. A secondary consideration is the trade-
off between the MCMC mixing of the Bernoulli variables and the computational
cost of updating them.

Similar to our proposed algorithm, scalable MH (SMH; Cornish et al.,, 2019)
also assumes that the log-likelihood difference is bounded as in (3); however,
the delayed-acceptance idea is applied repeatedly to give an overall
acceptance probability of

$\alpha_{1}(\theta,\theta^{\prime})\prod_{i=1}^{n}1\wedge\exp\left(\Delta_{i}\right).$
(12)

We explain the remainder of the algorithm using a formulation that is
equivalent to that in Cornish et al., (2019), but using quantities more
familiar from our algorithm. $B$ and $T_{1},\dots,T_{B}$ are (effectively)
simulated via the same mechanism as in our algorithm so that $N_{i}$, the
number of repeats of the index $i$ in the set $\\{T_{j}\\}_{j=1}^{B}$ has a
$\mathsf{Pois}(c_{i}M(\theta,\theta^{\prime}))$ distribution. Then, for each
$i$, $A_{i}\sim\mathsf{Binom}(N_{i},p_{i})$ is created, where
$p_{i}=-\log[1\wedge\exp(\Delta_{i})]/[c_{i}M(\theta,\theta^{\prime})]$. Since
$A_{i}\sim\mathsf{Pois}(-\log[1\wedge\exp(\Delta_{i})])$,
$\mathbb{P}(\sum_{i=1}^{n}A_{i}=0)$ is exactly the product term in (12). As
$(1\wedge a)(1\wedge b)\leq 1\wedge ab$, the acceptance probability (12) is
always lower than ours. Experiments with first-order control variates later in
this section show that to avoid the acceptance rate dropping unacceptably low,
the scaling parameter in the proposal $q(\theta^{\prime}|\theta)$ must be
reduced to well below the optimal value for the equivalent MH algorithm. Our
bounds $c_{i}M(\theta,\theta^{\prime})$ are also typically tighter than those
in Cornish et al., (2019), especially in moderate to high dimensions.

The TunaMH algorithm of Zhang et al., (2020) again uses (3), but without
control variates. It creates quantities $\phi_{i}$ and $\phi_{i}^{\prime}$,
but with
$\phi_{i}=[\ell_{i}(\theta)-\ell_{i}(\theta^{\prime})+c_{i}M(\theta,\theta^{\prime})]/2$.
This form of $\phi_{i}$ is a special case of (5), with $\gamma=\frac{1}{2}$.
Section 4.1 shows that the optimal choice of $\gamma$ is always $0$.

In TunaMH, where there are no control variates,
$c_{i}M(\theta,\theta^{\prime})$ is much larger (because the quantity it is
bounding is much larger); however, it is also, typically, a much tighter bound
on $\phi_{i}$ and $\phi_{i}^{\prime}$. As a consequence, with $\phi$ and
$\phi^{\prime}$ varying across the range from $0$ to
$c_{i}M(\theta,\theta^{\prime})$, the acceptance ratio in (7) has a high
variance and is often very low. TunaMH adds an extra term to each Poisson
expectation, $\phi_{i}$. This additional expectation lower bounds both the
numerator and denominator of every term in the product in (7) and, if it is
sufficiently large compared with $\phi_{i}$ and $\phi_{i}^{\prime}$, it
stabilises the acceptance rate. The size of the term is _tuna_ ble via an
additional hyperparameter, giving the algorithm its name. However, this
additional expectation comes at a cost of increasing
$\mathbb{E}\left[{B}\right]$ and hence the computational cost, leading to an
apparent trade-off between the acceptance rate and the subsample size.

In fact, it is a three-way trade-off. Decreasing the scaling parameter,
$\lambda$, for the jumps in the MH proposal, $q(\theta^{\prime}|\theta)$,
decreases $\|\theta^{\prime}-\theta\|$ and hence
$\ell_{i}(\theta^{\prime})-\ell_{i}(\theta)$ and, thus,
$M(\theta,\theta^{\prime})$. In this way, each $\phi_{i}$ and
$\phi_{i}^{\prime}$ becomes smaller compared with the additional Tuna
expectation; the smaller $\lambda$ helps stabilise the acceptance probability.
In reducing $\|\theta^{\prime}-\theta\|$, it also reduces the computational
cost via a smaller $\mathbb{E}\left[{B}\right]$ (the size of the final
subsample) through smaller $\phi_{i}$. The advice in Zhang et al., (2020) is
to add only a small amount to each expectation so that no more than one extra
observation is used on average in each subsample, and to aim for an empirical
acceptance rate of around 60%. As a consequence, the Metropolis–Hastings
scaling $\lambda$ is reduced substantially, the Markov chain mixes more slowly
and although the algorithm is more efficient than MH with _the same_ $\lambda$
but no subsampling, it brings no substantial improvement over an _optimally
scaled_ MH algorithm (which we illustrate in Figure 1).

To be more specific, consider the random-walk proposal (19) for a fixed
dimension, $d$, and examine the effect of varying $n$. When the acceptance
rate is reasonable, the mixing efficiency of the RWM is approximately
proportional to $\lambda^{2}/d$, which (with $d$ fixed) is
$\mathcal{O}(\lambda^{2})$, and does not depend on $n$. However, the
computational cost of each iteration is proportional to $n$, giving an overall
efficiency of $\mathcal{O}(\lambda^{2}/n)$. By contrast, for Tuna,
$M(\theta,\theta^{\prime})\propto\|\theta^{\prime}-\theta\|=\mathcal{O}(\lambda
n^{-1/2})$, by the Bernstein-von Mises Theorem, so the computational cost is
$\mathcal{O}(\lambda n^{1/2})$. As with the RWM, for a reasonable acceptance
rate, the mixing efficiency is proportional to $\lambda^{2}$, so the overall
efficiency is $\mathcal{O}(\lambda n^{-1/2})$. For the RWM, it is possible to
choose $\lambda=\mathcal{O}(1)$ and maintain a good acceptance rate; in
contrast, in practice for Tuna (see, for example, Table 7), to follow the
advice in Zhang et al., (2020) and keep the acceptance rate around $60\%$, we
find that $\lambda$ must be decreased in proportion to $n^{-1/2}$, leading to
the same order of efficiency as the RWM; i.e., there is no free lunch.

Figure 1 shows the effective sample size (ESS) per second, average batch sizes
and acceptance rates as dimension $d$ increases for a logistic regression
model with $n=30,000$ observations. It compares three algorithms: the RWM and
versions of the same algorithm using TunaMH and SMH with first-order control
variates, which we refer to as SMH-1. Panel (b) shows that Tuna uses a
fraction of the observations that the full-batch RWM uses, but it is clear
from Panel (a) that TunaMH is not substantially more efficient than the RWM,
especially as $d$ increases. For example, for $d=5$ and $d=60$, Tuna requires
on average per MCMC iteration subsamples of sizes $10^{1.38}\approx 24$ and
$10^{2.47}\approx 301$, respectively. As recommended in Zhang et al., (2020),
we set the Tuna parameter so that the expected number of additional data
points was typically just below $1$ and tuned the scaling parameter of the
proposal distribution to return an acceptance rate of $60\%$, as shown in
panel (c). The fact that TunaMH uses a fraction of $n$ makes it significantly
faster than the RWM, but it does not make it substantially more efficient
because of the small scaling parameter of the proposal distribution that is
needed to achieve the target acceptance rate.

(a)

(b)

(c)

Figure 1: Acceptance rates and ESS per second for SMH-1, Tuna and RWM. The
results are based on synthetic datasets generated from a logistic regression
model with $n=10^{4.5}\approx 30,000$ observations. The y-axis of panels (a)
and (b) are presented in the logarithm base 10.

All of the algorithms in Figure 1 use
$q(\theta^{\prime}|\theta)=\mbox{N}(\theta,\frac{\lambda^{2}}{d}V_{d})$ as the
proposal distribution, where $V_{d}=-[H(\hat{\theta})]^{-1}$ denotes minus the
inverse of the Hessian at an estimate $\widehat{\theta}$ near the mode and
$\lambda>0$ is a scaling parameter that requires tuning. Following the
literature on optimal scaling for the random-walk Metropolis algorithm (Gelman
et al.,, 1997; Roberts and Rosenthal,, 2001; Sherlock and Roberts,, 2009), we
set $\lambda=2.38$. This choice leads, in the limit as $d\to\infty$, to an
acceptance rate of approximately $23.4\%$, as shown in panel (c) of Figure 1.
In sharp contrast, for Tuna, $\lambda$ is set from $0.03$ ($d=5$) to $0.1$
($d=60$) to maintain an acceptance rate of around $60\%$; see Table 6 in
Appendix C.3 for all values.

The bounds on the log-likelihood differences used in the SMH algorithm are
quite loose (see Section 5), even for moderate $d$, especially when using
first-order control variates. In practice, looser bounds imply that larger
subsamples need to be evaluated in the MH ratio, leading to increased
computational time and lower efficiency. Panel (a) of Figure 1 shows that
SMH-1’s ESS per second deteriorates more quickly with dimension than that of
Tuna and the RWM. The explanation for this is two-fold. Firstly, Panel (b)
shows the steady increase in SMH’s average batch size as the number of
variables in the linear predictor of the logistic regression model increases.
The greater the batch size, the slower the algorithm. For example, for $d=5$
and $d=50$, SMH’s average batch size jumps from $10^{1.94}\approx 88$ to
$10^{4.10}\approx 12,751$. For $d=60$, almost all $n$ observations are
required. Consequently, in order to keep a constant acceptance rate at which
SMH-1 is most efficient as dimension $d$ changes, it is necessary to reduce
the scaling parameter in the proposal distribution. For the results presented
in Figure 1, $\lambda$ varied from $2$ ($d=5$) to $0.65$ ($d=60$).

In Sections 5 and 6, we empirically show through simulation experiments and
applications to real-world datasets that our proposed method, MH-SS, is
substantially more efficient than Tuna and the RWM and that it requires
substantially fewer observations than SMH for both first- and second-order
control variates. The improvements over TunaMH are a direct consequence of the
the use of control variates with a good bound on $M(\theta,\theta^{\prime})$
and of the optimal choice of $\gamma$. The improvement over SMH is due to the
control variate bound being tighter and scaling better with dimension than
that in SMH and because the MH-SS acceptance probability is a product of two
terms, rather than of many terms, all bounded above by $1$. Our tighter bounds
on the log-likelihood differences are introduced in Section 4.2. As we will
show in Section 4.1, the optimal choice of $\gamma$ in (5) leads to larger
acceptance rates, which helps maximise the computational efficiency of our
algorithms. In addition, in Appendix B (see also Section 4.3), through an
asymptotic analysis, we compare our proposed bounds with those of SMH and show
that our bounds are at least a factor of $d^{1/2}$ tighter. Finally, Section
4.4 studies both the scaling and acceptance rate that maximise the efficiency
for our proposed methods in the case of a RWM proposal. In contrast to the
standard RWM algorithm, our results suggest that the optimal acceptance rate
for the MH-SS algorithm is approximately $45\%$.

##  4 Theoretical results

We now present theoretical results to support our proposed method. First, we
demonstrate that the acceptance rate of MH-SS is maximised when $\gamma$ is
set to zero in (5). Then, we present general bounds for the log-likelihood
differences considering first- and second-order control variates. The proofs
of all theoretical results in this section are deferred to Appendix A.

###  4.1 Optimality with respect to $\phi$

######  Theorem 1.

Let $\alpha_{2}^{(\gamma)}(\theta,\theta^{\prime})$ be the acceptance
probability $\alpha_{2,MHSS}(\theta,\theta^{\prime})$ from (11) when
$\phi_{i}$ and $\phi^{\prime}_{i}$ are defined as in (5) and (6). Then for any
fixed $\theta$ and $\theta^{\prime}$,
$\alpha_{2}^{(\gamma)}(\theta,\theta^{\prime})$ is a decreasing function of
$\gamma$, so is maximised when $\gamma=0$.

_Proof_ : See Section A.1 in the Appendix. $\square$

Since $\alpha_{1}(\theta,\theta^{\prime})$ in (10) does not depend on
$\gamma$, this implies that the overall acceptance probability is maximised
when $\gamma=0$. Furthermore, the number of samples used in the calculation
depends only on the bounds $c_{i}M(\theta,\theta^{\prime})$ and so the
computational cost is not affected by $\gamma$. Since the proposal,
$q(\theta^{\prime}|\theta)$, is the same, whatever the value of $\gamma$,
Peskun, (1973) and Tierney, (1998) lead to the following result:

######  Corollary 1.

The efficiency of Algorithm 1 in terms of effective samples per second is
maximised when $\gamma=0$.

To illustrate the practical implications of Theorem 1, we perform simulation
experiments using our method with first-order control variates on a logistic
regression model. For each combination of $n=(10^{3},10^{4},10^{5})$ and
$d=(5,10,30,50,100)$, we set $\gamma$ to $0$, $0.5$ and $1$. Setting
$\gamma=0.5$ corresponds to employing the $\phi$ function in (5) used by Zhang
et al., (2020) (albeit with control variates and without the additional Tuna
term), thereby illustrating the acceptance rate of our method if we were to
use this sub-optimal choice of $\gamma$. Figure 2 presents acceptance rates of
our proposed method with first-order control variates obtained from 10 Monte
Carlo simulations. Unsurprisingly, the acceptance rates of MH-SS are maximised
when $\gamma=0$, regardless of the values of $n$ and $d$, which is consistent
with Theorem 1. Conversely, $\gamma=0.5$ consistently yields lower acceptance
rates compared to $\gamma=0$, while $\gamma=1$ exhibits the poorest results.
Since the algorithms differ only in terms of the value of $\gamma$, Corollary
1 states that the one with the highest acceptance rate inevitably exhibits the
largest effective sample size. We also investigated the effect of $\gamma$ on
acceptance rates of our method with second-order control variates. Whilst the
same ordering with $\gamma$ was evident, the differences were less pronounced,
and we omit them for the sake of conciseness.

Figure 2: Acceptance rates of the proposed method with first-order control
variates. For each combination of $\gamma$, $n$ and $d$, $10$ synthetic
datasets are generated from a logistic regression model.

###  4.2 Bounds on the remainder terms

Our method assumes that
$|\ell_{i}(\theta^{\prime})-\ell_{i}(\theta)-r_{i}^{(k)}(\theta,\theta^{\prime})|\leq
c_{i}^{(k)}M^{(k)}(\theta,\theta^{\prime})$, for $k=1,2$. Bounds on these
differences are provided in Cornish et al., (2019); however, these bounds are
loose as dimension $d$ increases. In this section, we provide tighter generic
bounds and then refine these further in the case of regression models.

####  4.2.1 General bounds on the remainders

We first prove general bounds which will be useful in all of our examples in
Sections 5 and 6 and many others besides. These general bounds then lead to
specific forms in each example; see the end of this subsection for cases of
logistic regression, probit regression and a form of Poisson regression.

######  Theorem 2.

Let
$\theta,\theta^{\prime},\widehat{\theta}\in\mathcal{X}\subseteq\mathbb{R}^{d}$
and let $\mathcal{A}$ be the triangle with vertices at
$\theta,\theta^{\prime},\widehat{\theta}$. Let
$\ell:\mathbb{R}^{d}\to\mathbb{R}$ and let $g(\theta)=\nabla\log\pi(\theta)$
and $H(\theta)=\nabla\nabla^{\top}\log\pi(\theta)$. Let $\|v\|$ indicate the
$L_{2}$ norm of $v$. Set

$\displaystyle r^{(1)}(\theta,\theta^{\prime};\widehat{\theta})$
$\displaystyle:=(\theta^{\prime}-\theta)^{\top}g(\widehat{\theta}),$ (13)
$\displaystyle r^{(2)}(\theta,\theta^{\prime};\widehat{\theta})$
$\displaystyle:=(\theta^{\prime}-\theta)^{\top}g(\widehat{\theta})+(\theta^{\prime}-\theta)^{\top}H(\widehat{\theta})\left\\{\frac{1}{2}(\theta^{\prime}+\theta)-\widehat{\theta}\right\\}.$
(14)

  * •

Suppose that there exists $x\in\mathbb{R}^{d}$ and $M_{1}\geq 0$ such that for
any $\theta_{1}\in\mathcal{A}$ and $u,v\in\mathbb{R}^{d}$,

$|u^{\top}H(\theta_{1})v|\leq|x^{\top}u|~{}|x^{\top}v|~{}M_{1}.$ (15)

Then

$\left|\ell(\theta^{\prime})-\ell(\theta)-r^{(1)}(\theta,\theta^{\prime};\widehat{\theta})\right|\leq|x^{\top}(\theta^{\prime}-\theta)|\max\left(|x^{\top}(\theta-\widehat{\theta})|,|x^{\top}(\theta^{\prime}-\widehat{\theta})|\right)M_{1}.$
(16)

If, instead, for any $\theta_{1}\in\mathcal{A}$ and $u,v\in\mathbb{R}^{d}$,
there exists $M_{2}\geq 0$ such that

$|u^{\top}H(\theta_{1})v|\leq\|u\|~{}\|v\|~{}M_{2}$

then

$\left|\ell(\theta^{\prime})-\ell(\theta)-r^{(1)}(\theta,\theta^{\prime};\widehat{\theta})\right|\leq\|\theta^{\prime}-\theta\|\max\left(\|\theta-\widehat{\theta}\|,\|\theta^{\prime}-\widehat{\theta}\|\right)M_{2}.$

  * •

Alternatively, suppose that there exists $x\in\mathbb{R}^{d}$ and $L_{1}\geq
0$ such that for any $\theta_{1},\theta_{2}\in\mathcal{A}$ and
$u,v\in\mathbb{R}^{d}$

$\left|u^{\top}\left\\{H(\theta_{2})-H(\theta_{1})\right\\}v\right|\leq|x^{\top}u|~{}|x^{\top}v|~{}|x^{\top}(\theta_{2}-\theta_{1})|~{}L_{1}.$
(17)

Then

$\displaystyle\left|\ell(\theta^{\prime})-\ell(\theta)-r^{(2)}(\theta,\theta^{\prime};\widehat{\theta})\right|$
$\displaystyle\leq\frac{1}{12}|x^{\top}(\theta^{\prime}-\theta)|^{3}~{}L_{1}$
$\displaystyle+\frac{1}{2}|x^{\top}(\theta^{\prime}-\theta)|~{}|x^{\top}(\theta-\widehat{\theta})|^{2}~{}L_{1}$
$\displaystyle+\frac{1}{2}|x^{\top}(\theta^{\prime}-\theta)|~{}|x^{\top}(\theta^{\prime}-\widehat{\theta})|^{2}~{}L_{1}.$
(18)

If, instead, for any $\theta_{1},\theta_{2}\in\mathcal{A}$ and
$u,v\in\mathbb{R}^{d}$, there exists $L_{2}\geq 0$ such that

$\left|u^{\top}\left\\{H(\theta_{2})-H(\theta_{1})\right\\}v\right|\leq\|u\|~{}\|v\|~{}\|\theta_{2}-\theta_{1}\|~{}L_{2},$

then

$\displaystyle\left|\ell(\theta^{\prime})-\ell(\theta)-r^{(2)}(\theta,\theta^{\prime};\widehat{\theta})\right|$
$\displaystyle\leq\frac{1}{12}\|\theta^{\prime}-\theta\|^{3}~{}L_{2}$
$\displaystyle+\frac{1}{2}\|\theta^{\prime}-\theta\|~{}\|\theta-\widehat{\theta}\|^{2}~{}L_{2}$
$\displaystyle+\frac{1}{2}\|\theta^{\prime}-\theta\|~{}\|\theta^{\prime}-\widehat{\theta}\|^{2}~{}L_{2}.$

_Proof_ : See Section A.2 in the Appendix. $\square$

######  Remark 1.

Let $\kappa$ be the magnitude of the largest eigenvalue of a symmetric matrix
$W$, where $W=H(\theta_{1})$ or $W=H(\theta_{2})-H(\theta_{1})$. In general,
$a^{\top}Wb\leq\kappa\|a\|~{}\|b\|$, motivating the two unnumbered conditions.
In the special case of a single observation from a generalised linear model,
however, it will turn out that we may obtain a more subtle bound by taking $x$
to be a covariate vector and using the forms in (16) and (18).

####  4.2.2 Regression models: Further improvement on bounds

We now specialise to the case where each observation is a realisation from a
regression model with a covariate vector, $x$, regression coefficients,
$\beta$, and observation, $y$, such that the likelihood for the observation is
$\ell(\beta)=\mathfrak{h}(\eta;y)$, where $\eta=x^{\top}\beta$.

In the following, we abbreviate $\mathfrak{h}(\eta;y)$ to
$\mathfrak{h}(\eta)$. Firstly,

$g(\beta)=\nabla_{\beta}\ell=(\nabla_{\beta}~{}\eta)~{}\mathfrak{h}^{\prime}(\eta)=x\mathfrak{h}(\eta)~{}~{}~{}\mbox{and}~{}~{}~{}H(\beta)=\nabla_{\beta}\nabla_{\beta}^{\top}\ell=xx^{\top}\mathfrak{h}^{\prime\prime}(\eta).$

Thus,
$|a^{\top}H(\beta)b|=|x^{\top}a|~{}|x^{\top}b|~{}|\mathfrak{h}^{\prime\prime}(\eta)|$;
also, with $\eta_{1}=x^{\top}\beta_{1}$ and $\eta_{2}=x^{\top}\beta_{2}$,

$\displaystyle|a^{\top}\\{H(\beta_{2})-H(\beta_{1})\\}b|$
$\displaystyle=|a^{\top}x\\{\mathfrak{h}^{\prime\prime}(\eta_{2})-\mathfrak{h}^{\prime\prime}(\eta_{1})\\}x^{\top}b|=|x^{\top}a|~{}|x^{\top}b|~{}|\eta_{2}-\eta_{1}|~{}|\mathfrak{h}^{\prime\prime\prime}(\eta_{*})|$
$\displaystyle=|x^{\top}a|~{}|x^{\top}b|~{}|x^{\top}(\beta_{2}-\beta_{1})|~{}|\mathfrak{h}^{\prime\prime\prime}(\eta_{*})|,$

for some $\eta_{*}\in[\eta_{1},\eta_{2}]$. So we can find appropriate bounds
for (15) and (17) by bounding $|\mathfrak{h}^{\prime\prime}(\eta;y)|$ and
$|\mathfrak{h}^{\prime\prime\prime}(\eta;y)|$.

Before doing this, however, we tighten our general bounds further in the case
of a regression model, improving over the obvious bound that arises from
repeated application of the Cauchy-Schwarz inequality. The improvement helps
in practice because in high-dimensions,“most” vectors are almost perpendicular
to each other so the inner product $\omega$ that we will define is often close
to $0$.

######  Lemma 1.

Let $u,v,x\in\mathbb{R}^{d}$ and let $\omega:=u\cdot v/(\|u\|\|v\|)$, then for
$k>0$, $|u\cdot x|^{k}|v\cdot x|\leq D_{k}(\omega)\|u\|^{k}\|v\|\|x\|^{k+1}$,
where

$D_{k}(\omega)=\frac{1}{c_{k}(\omega)(k+1)^{(k+1)/2}}\left\\{k+|\omega|c_{k}(\omega)\right\\}^{(k+1)/2},$

and
$c_{k}(\omega)=\sqrt{k+\frac{1}{4}(k-1)^{2}\omega^{2}}-\frac{1}{2}(k-1)|\omega|$.

_Proof_ : See Section A.3 in the Appendix. $\square$

######  Remark 2.

We make three points.

  1. 1.

When $k=1$ and $c_{k}=1$, $D_{k}(\omega)$ simplifies to
$D_{1}(\omega)=\\{1+|\omega|\\}/2$. Giving a typical improvement of nearly a
factor of $2$ over the obvious application of the Cauchy-Schwarz inequality.

  2. 2.

When $u$ and $v$ are aligned, then $\omega=1$ and $c_{k}=1$ implies
$D_{k}(\omega)=1,$ which retrieves the bound from the Cauchy-Schwarz
inequality.

  3. 3.

When $u$ and $v$ are orthogonal, as is usually approximately the case in high
dimensions, $\omega=0$, $c_{k}=\sqrt{k}$ and

$D_{k}=\frac{1}{\sqrt{k}}\left(1-\frac{1}{k+1}\right)^{(k+1)/2}\uparrow\frac{1}{\sqrt{k}}\exp(-1/2).$

So in typical high-dimensional cases, we gain a factor of $e^{-1/2}/\sqrt{k}$
over the Cauchy-Schwarz inequality.

Combining Lemma 1 with Theorem 2 then gives the following, which we apply in
practice for each of our regression models.

######  Corollary 2.

Consider a likelihood term $\ell(\theta)=\mathfrak{h}(x^{\top}\theta;y)$ and
define

$\omega:=\frac{(\theta-\widehat{\theta})\cdot(\theta^{\prime}-\theta)}{\|\theta-\widehat{\theta}\|\|\theta^{\prime}-\theta\|}~{}~{}~{}\mbox{and}~{}~{}~{}\omega^{\prime}:=\frac{(\theta^{\prime}-\widehat{\theta})\cdot(\theta^{\prime}-\theta)}{\|\theta^{\prime}-\widehat{\theta}\|\|\theta^{\prime}-\theta\|}.$

If $\|\mathfrak{h}^{\prime\prime}(x^{\top}\theta;y)\|\leq M(y)$ then,

$\left|\ell(\theta^{\prime})-\ell(\theta)-r^{(1)}(\theta,\theta^{\prime};\widehat{\theta})\right|\leq\|x\|^{2}\|\theta^{\prime}-\theta\|\max\left[\|\theta-\widehat{\theta}\|D_{1}(\omega),\|\theta^{\prime}-\widehat{\theta}\|D_{1}(\omega^{\prime})\right]~{}M(y).$

Alternatively, if $\|\mathfrak{h}^{\prime\prime\prime}(\theta;y)\|\leq L(y)$
then

$\displaystyle\left|\ell(\theta^{\prime})-\ell(\theta)-r^{(2)}(\theta,\theta^{\prime};\widehat{\theta})\right|\leq\frac{\|x\|^{3}\|\theta^{\prime}-\theta\|}{2}\left\\{\frac{1}{6}\|\theta^{\prime}-\theta\|^{2}+\|\theta-\widehat{\theta}\|^{2}D_{2}(\omega)+\|\theta^{\prime}-\widehat{\theta}\|^{2}D_{2}(\omega^{\prime})\right\\}L(y).$

In moderate to high dimensions, most of the time, $\theta^{\prime}-\theta$ is
close to perpendicular to both $\theta-\widehat{\theta}$ and
$\theta^{\prime}-\widehat{\theta}$, so Corollary 2 gives a reduction of
(typically) close to a factor of $2$ for the first-order bound and around
$2.5$ for the second-order bound. In moderate to high dimensions, most $x$ are
also close to perpendicular to the three vectors of interest, so the true
remainder is typically much smaller than even this tighter bound.

For the regression models in Sections 5 and 6, we use the following values
(justified in Section A.4):

  * •

Logistic regression: $M(y)=1/4$ and $L(y)=\sqrt{3}/18$;

  * •

Probit regression: $M(y)=1$ and $L(y)=0.3$;

  * •

Poisson regression with an expectation of $\log[1+\lambda\exp(\eta)]$ for any
fixed $\lambda$: $M(y)=0.25+0.168y$ and $L(y)=\sqrt{3}/18+0.061y$.

For the Tuna algorithm with logistic regression, we may use the bound
$|\ell_{i}(\theta^{\prime})-\ell_{i}(\theta)|\leq\|x_{i}\|\|\theta^{\prime}-\theta\|$
(e.g., Cornish et al.,, 2019; Zhang et al.,, 2020). For Poisson regression
with an expectation of $\log[1+\exp(\eta)]$, the bound is
$|\ell_{i}(\theta^{\prime})-\ell_{i}(\theta)|\leq\|x_{i}\|\max(1,y_{i})\|\theta^{\prime}-\theta\|$.

###  4.3 Computational Cost

In Appendix B, we consider a sequence of posterior distributions in the
multivariate-normal Bernstein-von Mises limit and examine the random walk
Metropolis proposal using a suitable scaling and variance matrix. We show that
the computational cost per iteration of MH-SS is $\Theta(d^{3/2})$ when using
$r_{i}^{(1)}$ and $\Theta(d^{3}/n^{1/2})$ when using $r_{i}^{(2)}$. We also
show that whichever control variate is used, the corresponding costs when
using the bounds in Cornish et al., (2019) are at least a factor of $d^{1/2}$
larger.

###  4.4 Tuning the MH-SS proposal

In this section, we provide guidelines for the choice of the scaling parameter
$\lambda$ in the special case of a random-walk Metropolis proposal:

$\theta^{\prime}\sim\mathsf{N}\left(\theta,\frac{\lambda^{2}}{d}V_{d}\right),$
(19)

where $d$ is the dimension of $\theta$ and $V_{d}$ is a $d$-dimensional
positive definite matrix. In our analysis, we set
$V_{d}=-[H(\hat{\theta})]^{-1}$, where $[H(\hat{\theta})]^{-1}$ represents the
inverse of the Hessian at an estimate $\widehat{\theta}$ near the mode. In the
large-$n$ Bernstein-von Mises limit, this is equivalent to exploring a
$d$-dimensional standard Gaussian posterior for
$\psi=\sqrt{V_{d}}^{-1}\theta$, with a proposal of
$\psi^{\prime}\sim\mathsf{N}\left(\psi,(\lambda^{2}/d)I_{d}\right)$ and, as
such, fits squarely within the framework of previous analyses of the high-
dimensional limit Gelman et al., (1997); Roberts and Rosenthal, (2001);
Sherlock and Roberts, (2009); Sherlock et al., (2021).

Sherlock et al., (2021) investigate the tuning of delayed-acceptance random
walk Metropolis algorithms. However, while the MH-SS algorithms do have a
delayed-acceptance component, they have two specific properties that,
respectively aid and necessitate a new treatment: the extreme
accuracy/inaccuracy of the approximate target and an approximately linear
relationship between the proposal scaling and the computational cost.

The computational cost per iteration is approximately proportional to the
number of data points used and the expectation of this is proportional to
$M(\theta,\theta^{\prime})$. With the proposal (19),
$\|\psi^{\prime}-\psi\|_{2}=\Theta(1)$, whereas both
$\|\psi-\widehat{\psi}\|_{2}$ and $\|\psi^{\prime}-\widehat{\psi}\|_{2}$ are
$\Theta(d^{1/2})$, where $\widehat{\psi}=\sqrt{V_{d}}^{-1}\widehat{\theta}$.
Provided $V_{d}$ is well-conditioned, $\|\theta^{\prime}-\theta\|_{2}$ is,
therefore, also a factor of $d^{1/2}$ smaller than either
$\|\theta-\widehat{\theta}\|_{2}$ or
$\|\theta^{\prime}-\widehat{\theta}\|_{2}$; see also the analysis in Appendix
B. Theorem 2 then shows that for both first- and second-order control
variates, $M(\theta,\theta^{\prime})\propto\lambda$. This justifies the more
formal:

######  Assumption 1.

The computational cost per iteration of both MH-SS-1 and MH-SS-2 is
proportional to $\lambda$.

In Appendix D, we argue for the following assumption, which we have observed
in all of our empirical studies (see, for example, Figure 3):

######  Assumption 2.

For first-order control variates,
$\mathbb{E}\left[{\alpha_{1}(\theta,\theta^{\prime})}\right]\approx 1$ and for
values of $\theta^{\prime}$ that pass the screening stage,
$\mathbb{E}\left[{\alpha_{2,MHSS}(\theta,\theta^{\prime})}\right]\approx\mathbb{E}\left[{\alpha_{RWM}(\theta,\theta^{\prime})}\right]$.
For second-order control variates,
$\mathbb{E}\left[{\alpha_{1}(\theta,\theta^{\prime})}\right]\approx\mathbb{E}\left[{\alpha_{RWM}(\theta,\theta^{\prime})}\right]$
and for values of $\theta^{\prime}$ that pass the screening stage,
$\mathbb{E}\left[{\alpha_{2,MHSS}(\theta,\theta^{\prime})}\right]\approx 1$.

(a)

(b)

Figure 3: Acceptance rates: $\alpha_{1}$ in (10) on the left and $\alpha_{2}$
in (11) on the right for a logistic regression target in dimension $d=100$
with $n=30,000$ observations, and with covariates and true coefficients
simulated as described in Section 5.

Thus, whichever form of control variate is used, we anticipate that the
expected acceptance rate for a given scaling will be approximately the same as
the expected acceptance rate for the RWM algorithm with that scaling. High-
dimensional asymptotics for the Metropolis–Hastings algorithm with a RWM
proposal (19) (e.g., Gelman et al.,, 1997; Roberts and Rosenthal,, 2001;
Sherlock and Roberts,, 2009) provide a limiting expected acceptance rate of
$2\Phi(-J\lambda/2)$ for some constant, $J$, that depends on the roughness of
the posterior; with our choice of $V_{d}$, $J=1$. The expected squared jumping
distance of the $i$th component of $\theta$ is

$\mathsf{ESJD}(\lambda)=\mathbb{E}\left[{\left(\theta^{(j)}_{i}-\theta^{(j+1)}_{i}\right)^{2}}\right],$

where $\theta^{(j)}\sim\pi$ is the $j$th value in the MCMC chain and
$\theta^{(j+1)}$ is the $j+1$th. This is a natural measure of the efficiency
of the algorithm and, in the case of a limiting diffusion, is exactly the
right measure of the efficiency in the high-dimensional limit as it is
proportional to the speed of the limiting diffusion. In this limit, the ESJD
is proportional to the product of the expected acceptance probability and the
square of the scaling parameter, $\mathsf{ESJD}\propto
2\lambda^{2}\Phi(-\lambda/2)$. Maximising this quantity leads to the well-
known optimal acceptance rate of $0.234$ of the RWM algorithm.

Given Assumption 1, however, the overall efficiency for MH-SS algorithms is
proportional to

$\frac{\mathsf{ESJD}(\lambda)}{\lambda}\propto 2\lambda\Phi(-\lambda/2).$

This is maximised at $\lambda_{opt}\approx 1.50$ which gives an acceptance
rate of $2\Phi(-\lambda_{opt}/2)\approx 0.452$.

We found this advice to hold true in all of our empirical studies, an example
of which is given in Figure 4.

(a)

(b)

Figure 4: Optimal scaling for the MH-SS algorithms based on synthetic data
from a logistic regression model with $n=30,000$ observations and $d=100$
covariate parameters, including an intercept. The efficiency metric
(MSJD/$\mathbb{E}(B)$) is plotted as function of (a) the scaling parameter
($\lambda$) and (b) the empirical acceptance rate.

The same argument applies to any algorithm that uses our control variates,
provided that for any reasonable scaling, $\lambda$, the overall acceptance
rate is close to that for the RWM. The SMH bounds used in Cornish et al.,
(2019) are reproduced in Appendix B.2 and have only a weak relative dependence
on $\lambda$. Consequently, we might expect SMH-2 to be optimally scaled when
the acceptance rate is around $0.234$. Figure 10 in Appendix C verifies both
of these points for a particular logistic regression scenario. It also shows
that the behaviour of the SMH-1 acceptance rate on $\lambda$ is very different
to that of the RWM. We find that the scaling must decrease more quickly with
dimension than $\lambda/\sqrt{d}$ to maintain a reasonable acceptance rate as
dimension increases. Although there is no obviously optimal acceptance rate
that we might expect to hold across different scenarios, we observed
empirically through simulation experiments with various numbers of
observations $n$ (from $30,000$ to $100,000$) and dimensions $d$ (from $10$ to
$100$) that the optimal scaling parameter for SMH-1 is in the range between
$0.5$ to $1.5$. Furthermore, for each combination of $n$ and $d$, a scaling of
$\lambda=1$ achieved at least $80\%$ of the optimal efficiency for that
combination. Consequently, in Sections 5 (simulation experiments) and 6 (real-
world applications), we set $\lambda=1$ for SMH-1.

##  5 Simulation experiments

In this section, we compare the MH-SS algorithm against the Tuna, SMH and MH
algorithms, all using random-walk proposals. We apply them to a logistic
regression model and a form of the Poisson regression model. For MH-SS and
SMH, we consider first- and second-order control variates, which we denote by
MH-SS-1 and MH-SS-2; with the same convention for SMH. We obtain an estimate
of the posterior mode $\hat{\theta}$ through the stochastic gradient descent
algorithm (Robbins and Monro,, 1951). The Tuna algorithm increases
$\mathbb{E}\left[{B}\right]$ by $\chi C^{2}M(\theta,\theta^{\prime})^{2}$,
where $\chi$ is the additional tuning parameter. We follow the heuristic
proposed by Zhang et al., (2020), choosing $\chi$ so as to increase the
expected number of simulations by less than $1$ across most combinations of
$(\theta,\theta^{\prime})$. To aid reproducibility, we record the value of the
Tuna tuning parameter, $\chi$, used each time in Table 7 in Appendix C.3. For
all methods, the proposal distribution is of the form (19) with $V_{d}$ set to
$-[H(\hat{\theta})]^{-1}$, where $H(\hat{\theta})$ denotes the Hessian at our
approximation to the posterior mode and $\lambda>0$ is a scaling parameter.
For each method, $\lambda$ was chosen based on the maximisation of the mean
squared jump distance (MSJD) divided by the average subsample size,
$\mathbb{E}(B)$. Following Section 4.4, we set $\lambda=1.5$ for MH-SS-1 and
MH-SS-2, $\lambda=1$ for SMH-1, $\lambda=2$ for SMH-2 and $\lambda=2.38$ for
the RWM. In all simulation experiments carried out in this section, the ESSs
are at least $200$ for all algorithms. To ensure the computational times are
comparable, we implemented all methods in Python $3.9.6$ (Van Rossum and
Drake,, 2009).

###  5.1 Logistic regression

To empirically demonstrate the benefits of the proposed method, we initially
test it on synthetic datasets generated from a succession of logistic
regression models with varying data size, $n$, and parameter size, $d$. We
sample each covariate (except the intercept) independently as
$x_{i,j}\sim\mathsf{N}(0,1/d)$ and simulate independent
$\beta_{j}\sim\mathsf{N}(0,1)$. These choices ensure the linear predictor,
$x_{i}^{\top}\beta$, remains $\mathcal{O}(1)$ as $d$ increases and the
probabilities obtained from the logistic function applied to the linear
predictor are well distributed around $0.5$.

For each combination of $n$ and $d$, we run once each of the $6$ algorithms
under consideration and calculate the effective sample size (ESS; Green and
Han,, 1992; Liu and Chen,, 1995) per second as a metric of efficiency.
Additionally, we examine other metrics such as the average expected batch
size, which we denote $\mathbb{E}(B)$; Appendix C presents the results in
terms of the alternative efficiency metric of ESS divided by $\mathbb{E}(B)$.
To highlight any power-law relationships and enable all algorithms to be
jointly displayed, results involving the number of observations, $n$, average
expected batch size, $\mathbb{E}(B)$, and ESS per second are presented on a
logarithmic scale.

Figure 5 shows how the average batch size varies by $n$ and $d$. As argued in
Appendix B, the algorithms based on the second-order control variates use
fewer data values than their first-order counterparts and require even fewer
as $n$ increases. Since the bounds we propose are tighter than those of SMH,
the MH-SS algorithm uses less data. When $d=10$, SMH-2 sometimes uses fewer
observations than MH-SS-1, but MH-SS-2 always uses the fewest; in all other
dimensions considered, the MH-SS algorithms used fewer observations than the
SMH algorithms. For example, when $d=50$ and $n<10^{4}$, SMH requires all $n$
observations, whereas MH-SS requires fewer than $n$ observations, even for
$n=10^{3}$. In all scenarios, Tuna and MH-SS-1 require a similar subsample
size, the only exception being $d=100$ where Tuna uses fewer observations.

Figure 5: Average batch size for MH-SS, SMH and RWM for the logistic
regression model. For RWM, the average batch size is $n$. Both axes are
presented in the logarithm base 10.

Figure 6 compares the ESS per second across all methods. As expected, the ESS
per second of MH-SS with first- and second-order control variates is greater
than that of either SMH algorithm, except when $d=10$, where SMH-2 outperforms
MH-SS-1. For $d=10$, MH-SS-2 and SMH-2 exhibit similar efficiency for
$n=10^{5}=100,000$, even though the former uses fewer observations. To
understand this, we recall the results in the left-hand upper panel of Figure
5, where it is possible to see that $\mathbb{E}(B)\approx 10^{0}=1$ for MH-
SS-2, while $\mathbb{E}(B)\approx 10^{1}=10$ for SMH-2. In practice, we
observe that the run-time differences are tiny in cases where the average
batch sizes are small, helping to explain the similarity in terms of
efficiency. However, as $d$ increases, the efficiency gap becomes more
pronounced because the SMH bounds are at least a factor of $d^{1/2}$ worse
than the equivalent MH-SS bound. Even for $d=30$, SMH-2 would be a competitive
option against RWM only for $n>10^{3.5}$, while SMH-1 demonstrates relatively
poor performance for all considered dimensions.

Figure 6: ESS per second of MH-SS, SMH and RWM for the logistic regression
model. Both axes are presented in the logarithm base 10. Some ESSs are omitted
because $\mathbb{E}(B)\geq n$, which implies the use and efficiency of the RWM
algorithm.

Though Tuna and MH-SS-1 require a similar subsample size for all scenarios
considered in Figure 5, except for $d=100$, MH-SS with first- and second-order
control variates outperform Tuna in all dimensions in terms of efficiency. As
stated in Section 3, Tuna tends to favour smaller subsample sizes in detriment
of computational efficiency and, as argued in Section 3, its performance is
comparable to the RWM algorithm. To recap, the heuristic proposed by Zhang et
al., (2020) aims for an acceptance rate of $60\%$ and minimises the impact of
the Tuna tuning parameter on the average batch size. As result, in practice
Tuna takes tiny steps when exploring the posterior distribution. This means it
uses small subsample sizes and thus has low per-iteration run-time; however,
such an approach is inefficient from the effective sample size viewpoint as
small step sizes lead to high auto-correlation in the Markov chain.

###  5.2 Decomposing the improvements over SMH

Tuna uses no more information than the RWM, so it is unsurprising that their
efficiencies are similar and that MH-SS’s careful use of control variates
leads to a substantial increase in efficiency. We now investigate the
algorithms that use control variates, and separate the drivers of the superior
performance of MH-SS over SMH.

Both SMH and our proposed method bound
$|\ell_{i}(\theta^{\prime})-\ell_{i}(\theta)-r_{i}^{(k)}(\theta,\theta^{\prime};\widehat{\theta})|\leq
c_{i}^{(k)}M^{(k)}(\theta,\theta^{\prime})$, for $k=1,2$, though with
different $c_{i}^{(k)}$ and $M^{(k)}(\theta,\theta^{\prime})$, so we created a
hybrid algorithm named SMH-NB. This hybrid approach combines the structure of
SMH, which is detailed in Algorithm 1 of Cornish et al., (2019), with our
bounds, $c_{i}^{(k)}$ and $M^{(k)}(\theta,\theta^{\prime})$, from Section
4.2.2. We compare these hybrid algorithms with the vanilla SMH and MH-SS to
separate the benefit brought by these bounds from the other benefits of MH-SS.
Unlike the vanilla SMH, the scaling parameter in the proposal distribution of
the hybrid SMH with new bounds is set to $\lambda=0.5$ for SMH-1-NB and
$\lambda=1.5$ for SMH-2-NB. These values were chosen as they empirically
maximise MSJD/$\mathbb{E}(B)$.

Panel (a) of Figure 7 demonstrates the substantial reduction in the subsample
size that SMH can achieve when using our bounds. It also shows that, on
average, SMH-NB with the first-order control variates requires smaller
subsamples than MH-SS with the same approximation. This difference arises
exclusively from the smaller scaling parameter in the proposal distribution of
SMH-1-NB (i.e., $\lambda=0.5$ for SMH-1-NB versus $\lambda=1.5$ for MH-SS). In
contrast, MH-SS-2 and SMH-2-NB have approximately the same $\mathbb{E}(B)$
since they share the same $\lambda$. A smaller $\lambda$ implies smaller jumps
from $\theta$ to $\theta^{\prime}$, negatively impacting efficiency metrics
such as ESS since (when the acceptance rate is reasonable) a small $\lambda$
tends to increase the correlation between successive values in the chain,
especially in high-dimensional settings. For example, in panel (b) of Figure
7, both SMH-NB variants exhibit greater efficiency per second than their
vanilla counterparts. However, each is consistently outperformed by the
corresponding MH-SS algorithm, even though it shares the same bound,
demonstrating superior performance of the algorithm itself, as well as the
bounds.

(a)

(b)

Figure 7: Average batch size and ESS per second for MH-SS, vanilla SMH and SMH
with new bounds for the logistic regression model ($d=30$). Both axes are
presented in the logarithmic scale.

In summary, the results presented in Figure 7 suggest that SMH requires less
data per iteration and has greater efficiency if it adopts our tighter bounds.
However, this is not sufficient to make it at least as efficient as our
proposed methods. Though SMH and MH-SS share similarities, such as the use of
Poisson auxiliary variables and the assumption that the log-likelihood
differences are bounded, the final algorithms are different. SMH, and
especially SMH-1, suffers from its acceptance rate formulation as a product of
individual acceptance probabilities. Moreover, SMH requires sampling up to
$B-1$ additional Bernoulli random variables per iteration.

###  5.3 Poisson regression

We now compare the MH-SS, Tuna, SMH and RWM algorithms using a set of
synthetic data generated from a Poisson regression model with expectation
$\log(1+e^{x^{\top}\beta})$. We simulate synthetic datasets with $d=30$,
$n=(10^{4.5},10^{5})=(31,622,100,000)$ and sample the predictors and
coefficients as for the logistic regression scenarios. For Tuna, we set (in
the notation of Zhang et al., (2020)), $\chi=2\times 10^{-5}$ in both cases
and $\lambda=0.09$ for $n=31,622$ and $\lambda=0.05$ for $n=100,000$.
Following the suggestions of Zhang et al., (2020), these combinations lead to
acceptance rates of approximately $60\%$ and minimise the impact of $\chi$ on
the average batch size.

Table 1 shows that MH-SS-1 and MH-SS-2 have the largest ESS per second. Once
again, the empirical results for $\mathbb{E}(B)$ show that the SMH bounds
scale poorly in moderate to high dimensions and that our proposed bounds are
tighter. For instance, MH-SS-1 requires $14$ times fewer observations than
SMH-1, for $n=31,622$, and it is slightly more efficient than SMH-2 in both
scenarios. We note that the computational efficiency of SMH-1 is one of the
lowest among the compared algorithms because, as pointed out in Section 3, it
degrades quickly as $d$ increases. Though Tuna requires on average the
smallest subsample sizes per MCMC iteration, its efficiency is similar to that
of the RWM. As already discussed, this low efficiency is due to the scaling
parameter, $\lambda$, of the proposal distribution for Tuna being smaller than
for any other algorithm, generating highly correlated posterior samples.

Table 1: Acceptance rate, average batch size, ESS per second and ESS/$\mathbb{E}(B)$ for the Poisson regression model applied to synthetic data in $d=30$. The values in bold denote the two largest efficiency metrics. Model | Acc. rate | $\mathbb{E}(B)$ | ESS per second | ESS/$\mathbb{E}(B)$
---|---|---|---|---
$n=31,622$
Tuna | $0.6051$ | $100$ | $0.04$ | $0.02$
MH-SS-1 | $0.4245$ | $203$ | 4.31 | $5.07$
MH-SS-2 | $0.4505$ | $19$ | 60.91 | 56.42
SMH-1 | $0.1670$ | $2,889$ | $0.05$ | $0.08$
SMH-2 | $0.3184$ | $256$ | $4.29$ | 5.11
RWM | $0.2368$ | $31,622$ | $0.04$ | $0.04$
$n=100,000$
Tuna | $0.6042$ | $111$ | $0.02$ | $0.01$
MH-SS-1 | $0.4229$ | $195$ | 6.57 | $5.43$
MH-SS-2 | $0.4568$ | $10$ | 92.40 | 116.90
SMH-1 | $0.1630$ | $2,773$ | $0.05$ | $0.08$
SMH-2 | $0.3195$ | $128$ | $6.52$ | 10.06
RWM | $0.2384$ | $100,000$ | $0.01$ | $0.01$

##  6 Results on real datasets

In the simulation section, we showed that MH-SS outperforms the SMH, Tuna and
RWM algorithms for different combinations of $n$ and $d$ using synthetic data.
In this section, we extend our comparison using real-world datasets and three
regression models: logistic, probit and Poisson with an expectation of
$\log(1+e^{\eta})$, where $\eta$ is the linear predictor. The objective is to
illustrate that the MH-SS algorithm can be applied to the probit model, as
well as in more challenging settings with larger $n$ and correlated predictors
which are not continuous and/or normally distributed. Results for the Tuna
algorithm are not presented for the probit model as it is not possible to
bound the first derivative of the log-likelihood of this model. We assess the
algorithms in terms of computational efficiency, considering metrics including
ESS/$\mathbb{E}(B)$ and ESS per second. We also present the empirical
acceptance rates and average batch sizes. In all results presented in this
section, the ESSs are at least $150$ for all algorithms.

We first compare the methods on three binary datasets using the logistic and
probit regression models only. Our fourth and last comparison uses a Poisson
regression model on a dataset with $n=298,290$ and $d=28$. The first three
datasets differ in many aspects, such as the numbers of observations and
predictors as well the correlations between predictors. The later point is of
particular interest as the simulation experiments only considered an intercept
and continuous i.i.d. Gaussian predictors. Proposals are of the same form as
in the Section 5, with $V_{d}$ and $\lambda$ set in the same way. For Tuna, we
provide the values of $\lambda$ and $\chi$ separately for each dataset.

We now point out aspects of the results common to all four analyses. Whether
in terms of ESS per second or ESS/$\mathbb{E}\left[{B}\right]$, MH-SS-2 is
always the most efficient algorithm. In terms of ESS per second, it is always
at least an order of magnitude more efficient than any of the other
algorithms. The second- and third-most efficient algorithms are always MH-SS-1
and SMH-2, though the order varies. MH-SS-2 always uses the fewest
observations, but in all of the logistic regression models it is Tuna that
uses the second fewest observations, and for the Poisson regression, it uses
the third fewest observations. Despite the relatively low numbers of
observations, Tuna always performs worse than either MH-SS-1 or SMH-2 (and
often much worse than either) because the proposal scaling parameter is so
small; see the discussion in Section 3.

###  6.1 U.S. Current Population Survey

We first illustrate the applicability of our method on the 2018 United States
Current Population Survey (CPS; Ruggles et al.,, 2024). The CPS is a monthly
household survey carried out by the U.S. Census Bureau and the U.S. Bureau of
Labor Statistics that gathers information on the labour force for the
population of the U.S. The data contain variables such as income, education,
occupation, participation in welfare programs and health insurance. Using a
sample of $n=500,000$ survey participants, we wish to model whether or not an
individual has a total pre-tax personal income above $\$25,000$ for the
previous year, based on $10$ predictors. We fit the logistic and probit models
with an intercept term and normalise the covariates to improve numerical
stability, ensuring each has a mean of zero and variance of one. As some of
the variables are categorical, the linear predictor contains $d=31$ parameters
(including the intercept). For the Tuna algorithm, we set $\chi=10^{-5}$ and
$\lambda=0.019$.

Table 2 presents the results for acceptance rate, average batch size, ESS per
second and ESS divided by $\mathbb{E}(B)$ for the Tuna, MH-SS, SMH and RWM
algorithms. Whether for the logistic or probit models, SMH-1 is not a viable
option compared to the RWM. Although SMH-1 uses approximately $10\%$ of the
data, its efficiency per second is less than that of the RWM. This is due to
its relatively low scaling parameter, $\lambda$, and the additional
computational overheads compared with the vanilla RWM, due to it being a more
complex algorithm, including sampling extra random variables and performing
vector/matrix multiplications. We stress that in all of our implementations,
operations that can be pre-computed to speed up both SMH and MH-SS are
performed accordingly.

Table 2: Acceptance rate, average batch size, ESS per second and ESS/$\mathbb{E}(B)$ for the logistic and probit regression models applied to a sample of $n=500,000$ participants of the U.S. CPS. The values in bold denote the two largest efficiency metrics. Model | Algorithm | Acc. rate | $\mathbb{E}(B)$ | ESS per second | ESS/$\mathbb{E}(B)$
---|---|---|---|---|---
Logistic | Tuna | $0.5935$ | $284$ | $0.003$ | $0.001$
MH-SS-1 | $0.4405$ | $989$ | 1.35 | 0.92
MH-SS-2 | $0.4561$ | $88$ | 27.37 | 10.62
SMH-1 | $0.1577$ | $49,060$ | $0.002$ | $0.003$
SMH-2 | $0.3191$ | $8,288$ | $0.13$ | $0.13$
RWM | $0.2391$ | $500,000$ | $0.005$ | $0.002$
Probit | MH-SS-1 | $0.4038$ | $1,543$ | 0.82 | 0.66
MH-SS-2 | $0.4108$ | $71$ | 31.32 | 14.58
SMH-1 | $0.1397$ | $57,846$ | $0.0006$ | $0.0005$
SMH-2 | $0.2758$ | $5,076$ | $0.05$ | $0.028$
RWM | $0.1933$ | $500,000$ | $0.0005$ | $0.002$

###  6.2 Detection of gas mixtures

We now apply the various MCMC sampling algorithms to a dataset available in
the UCI machine learning repository (Kelly et al.,, 2024; Fonollosa et al.,,
2015) which gathers recordings of multiple chemical sensors used to detect gas
mixtures at varying concentrations. Our goal is to predict for a sample of
size $n=250,000$ whether the concentration of Ethylene in the air, measured in
parts per million, is above zero. Though we consider $d=7$ continuous
predictors only (including an intercept), each of which corresponds to a
different model sensor, some of them are highly correlated. In practice, this
correlation is expected to increase $\mathbb{E}(B)$ for SMH and MH-SS, thereby
potentially reducing their computational efficiency relative to RWM. For the
Tuna algorithm, we set $\chi=7\times 10^{-6}$ and $\lambda=0.025$.

Table 3 presents the efficiency metrics for all algorithms applied to the
logistic and probit models. Notably, for logistic regression, Tuna is not
competitive with any of the algorithms, except for the RWM, SMH-2 uses, on
average, $4$ times more observations per MCMC iteration than MH-SS-2, while
the ratio of the average batch sizes of SMH-1 and MH-SS-1 is only just over 2.
We remind the reader that as $\mathbb{E}(B)$ is a function of the bound on the
log-likelihood differences, the looser SMH bounds result in larger average
batch sizes. For a reasonable acceptance rate, SMH-1 requires a smaller
scaling than MH-SS-1; this mitigates against the difference in batch size but
reduces the mixing of the MCMC chain. The patterns are similar for the probit
model, except that the MH-SS algorithms are more than a factor of $30$ times
more efficient than their SMH counterparts.

Table 3: Acceptance rate, average batch size, ESS per second and ESS/$\mathbb{E}(B)$ for the logistic and probit regression models applied to the gas sensor dataset. The values in bold denote the two largest efficiency metrics. Model | Algorithm | Acc. rate | $\mathbb{E}(B)$ | ESS per second | ESS/$\mathbb{E}(B)$
---|---|---|---|---|---
Logistic | Tuna | $0.6091$ | $321$ | $0.02$ | $0.006$
MH-SS-1 | $0.4778$ | $945$ | $6.44$ | $4.21$
MH-SS-2 | $0.4774$ | $165$ | 91.91 | 23.77
SMH-1 | $0.4373$ | $2,260$ | $0.48$ | $0.72$
SMH-2 | $0.3491$ | $680$ | 8.90 | 6.98
RWM | $0.2738$ | $250,000$ | $0.03$ | $0.02$
Probit | MH-SS-1 | $0.4153$ | $1,338$ | 3.53 | $2.98$
MH-SS-2 | $0.4142$ | $103$ | 100.74 | 38.77
SMH-1 | $0.3792$ | $2,715$ | $0.05$ | $0.66$
SMH-2 | $0.2869$ | $382$ | $3.14$ | 12.15
RWM | $0.2308$ | $250,000$ | $0.002$ | $0.02$

###  6.3 High-energy particle physics

We apply the scalable MCMC algorithms to the Hepmass dataset (Whiteson,, 2016)
from the UCI machine learning repository. The dataset contains information
about signatures of exotic particles obtained from a high-energy physics
experiment. The binary response variable indicates whether a new particle of
unknown mass is observed. The dataset has originally $26$ continuous
predictors and is split into a training set of $7$ million observations and a
test set of $3.5$ million. To illustrate the applicability of the proposed
methods, we compare the algorithms on a subsample of $n=1,000,000$
observations from the training set. For Tuna, we set $\chi$ and $\lambda$ to
$10^{-5}$ and $0.013$, respectively.

Table 4 summarises the results of acceptance rate, $\mathbb{E}(B)$, ESS per
second and ESS/$\mathbb{E}(B)$ for the logistic and probit regression models
applied to the Hepmass dataset. The efficiencies of the RWM, Tuna and SMH are
all around two (or more) orders of magnitude smaller than those of SMH-2 and
MH-SS-1, which themselves are at least an order of magnitude smaller than that
of MH-SS-2.

Table 4: Acceptance rate, average batch size, ESS per second and ESS/$\mathbb{E}(B)$ for the logistic and probit regression models applied to the Hepmass dataset. The values in bold denote the two largest efficiency metrics. Model | Algorithm | Acc. rate | $\mathbb{E}(B)$ | ESS per second | ESS/$\mathbb{E}(B)$
---|---|---|---|---|---
Logistic | Tuna | $0.5914$ | $295$ | $0.002$ | $0.001$
MH-SS-1 | $0.4465$ | $766$ | $0.76$ | $1.25$
MH-SS-2 | $0.4557$ | $31$ | 46.07 | 31.10
SMH-1 | $0.1685$ | $11,844$ | $0.01$ | $0.01$
SMH-2 | $0.3221$ | $493$ | 2.45 | 2.43
RWM | $0.2433$ | $1,000,000$ | $0.003$ | $0.001$
Probit | MH-SS-1 | $0.4125$ | $1,031$ | $0.71$ | $1.01$
MH-SS-2 | $0.4215$ | $18$ | 48.67 | 55.93
SMH-1 | $0.1529$ | $15,374$ | $0.001$ | $0.012$
SMH-2 | $0.2852$ | $261$ | 0.84 | 4.66
RWM | $0.1992$ | $1,000,000$ | $0.0002$ | $0.001$

###  6.4 Reported road casualties in the UK

The UK Department for Transport publishes annual road safety statistics as
part of the Statistics and Registration Service Act 2007. The data include the
accident’s geographical coordinates, severity, speed limit of the road where
the accident took place, details about the vehicles involved, weather
conditions, road conditions, as well as time and date (Lovelace et al.,,
2019). We compare the performance of MH-SS, SMH and Tuna on the consolidated
UK road casualties data from $2020$ to $2022$. We aim to model at the
individual level of $n=298,290$ accidents the number of casualties based on 8
predictors, of which 2 are continuous. In total, the linear predictor in our
Poisson model with mean $\log(1+e^{\eta})$ has $d=28$ parameters (including
the intercept). For the Tuna algorithm, we set $\chi$ and $\lambda$ to
$10^{-7}$ and $0.19$, respectively.

Table 5 shows results of acceptance rate, average batch size, ESS per second
and ESS/$\mathbb{E}(B)$ for MH-SS, SMH and Tuna. As per the analyses presented
above considering datasets with a binary response variable, it is possible to
see that MH-SS-1 and MH-SS-2 are the algorithms that on average require the
smallest number of data points per MCMC iteration. Unsurprisingly, both
versions of MH-SS present the best results in terms of computational
efficiency. Though Tuna and MH-SS-1 have similar average batch sizes, again
the latter is significantly more efficient as shown by its ESS per second and
ESS/$\mathbb{E}(B)$.

Table 5: Acceptance rate, average batch size, ESS per second and ESS/$\mathbb{E}(B)$ for a Poisson regression model applied to the UK road casualties dataset. The values in bold denote the two largest efficiency metrics. Model | Algorithm | Acc. rate | $\mathbb{E}(B)$ | ESS per second | ESS/$\mathbb{E}(B)$
---|---|---|---|---|---
Poisson | Tuna | $0.5964$ | $2,680$ | $0.015$ | $0.002$
MH-SS-1 | $0.4542$ | $1,927$ | 0.79 | 0.52
MH-SS-2 | $0.4572$ | $391$ | 8.61 | 2.59
SMH-1 | $0.1893$ | $44,262$ | $0.003$ | $0.004$
SMH-2 | $0.3219$ | $20,320$ | $0.05$ | $0.06$
RWM | $0.2414$ | $298,290$ | $0.011$ | $0.004$

##  7 Discussion

Performing posterior sampling via MCMC on datasets with a large number of
observations poses a significant challenge, particularly to traditional
methods like the Metropolis–Hastings (MH) algorithm, which requires the
evaluation of all likelihood terms at each iteration. This challenge has led
Bayesian researchers to turn away from MCMC algorithms and to instead utilise
approximate and biased inference schemes, such as variational inference, in
order to approximate the posterior distribution. Current approaches in the
scalable MCMC literature aim to alleviate the computational burden of applying
MCMC on large data sets by either partitioning the data into batches and
running independent MCMC on each partition, or employing data subsampling
techniques. In this paper, we have presented a novel exact
Metropolis–Hastings-based algorithm designed to scale posterior sampling for
large datasets through data subsampling. Our approach assumes that the log-
likelihood differences can be tightly bounded and introduces an auxiliary
Poisson variable which allows us to extend the Metropolis–Hastings ratio so
that not all observations are needed, thereby reducing the computational
effort per iteration.

Several previous approaches that utilise subsampling techniques to scale MCMC
to large datasets rely on either the log-likelihood or log-likelihood
differences being bounded. These methods typically introduce auxiliary
variables to simplify the posterior sampling. However, there have been two
main drawbacks: i) the auxiliary variables can be sub-optimally designed,
which hinders the efficiency of the algorithm, and ii) the bounds on the log-
likelihood are not sufficiently tight, especially as the number of predictors,
$d$, increases. In our work, we have overcome these limitations and enhanced
the scalability of MCMC methods to large datasets.

Our auxiliary variable approach shares features with that in Zhang et al.,
(2020), using Poisson random variables whose expectation depends on both the
current and proposed value. However, as seen in the experiments and real-world
applications of Sections 5 and 6, the Tuna algorithm (Zhang et al.,, 2020)
typically performs relatively poorly because it requires sub-optimal scaling
of the Metropolis–Hastings proposal in order to control the acceptance rate.
Our more general formulation, combined with the use of control variates, leads
to orders of magnitude greater efficiency. The scalable Metropolis–Hastings
algorithm (Cornish et al.,, 2019) uses the same generic form for the control
variates as we do, but our bounds on the discrepancy between the control
variate and the truth scale better with increasing dimension and this,
together with the removal of the large product in the acceptance probability
and another novel tightening, has a substantial impact on the overall
efficiency. We also have provided guidelines on choosing the proposal scaling,
$\lambda$, for MH-SS in the case where the Metropolis–Hastings proposal is
that of a random-walk Metropolis algorithm: choose $\lambda$ to target an
empirical acceptance of around $0.452$. This guidance has been verified
empirically and holds whether using first- or second-order control variates.

The simulation experiments presented in Section 5 utilising synthetic data
from logistic and Poisson regression models demonstrate that our proposed
algorithms require much smaller average batch sizes than those in Cornish et
al., (2019) and exhibit consistently greater computational efficiency per
second across various combinations of number of data points, $n$, and number
of predictors, $d$. The real-world applications in Section 6, featuring
datasets with observations ranging from $n=250,000$ to $n=1,000,000$ and three
models (logistic, probit and Poisson regressions), show that i) our algorithms
are highly efficient even in settings characterised by large $n$ and
correlated and/or categorical predictors and that ii) the RWM, SMH and TunaMH
algorithms do not scale as well as MH-SS, and MH-SS-2 in particular.

## References

  * Andersen et al., (2018) Andersen, M., Winther, O., Hansen, L. K., Poldrack, R., and Koyejo, O. (2018).  Bayesian structure learning for dynamic brain connectivity.  In International Conference on Artificial Intelligence and Statistics, pages 1436–1446. PMLR.
  * Andrieu and Roberts, (2009) Andrieu, C. and Roberts, G. O. (2009).  The pseudo-marginal approach for efficient Monte Carlo computations.  The Annals of Statistics, 37(2):697 – 725.
  * Bardenet et al., (2014) Bardenet, R., Doucet, A., and Holmes, C. (2014).  Towards scaling up Markov chain Monte Carlo: an adaptive subsampling approach.  In International Conference on Machine Learning, pages 405–413. PMLR.
  * Bardenet et al., (2017) Bardenet, R., Doucet, A., and Holmes, C. (2017).  On Markov chain Monte Carlo methods for tall data.  Journal of Machine Learning Research, 18(47).
  * Beaumont, (2003) Beaumont, M. A. (2003).  Estimation of population growth or decline in genetically monitored populations.  Genetics, 164(3):1139–1160.
  * Bierkens et al., (2019) Bierkens, J., Fearnhead, P., and Roberts, G. (2019).  The Zig-Zag process and super-efficient sampling for Bayesian analysis of big data.  The Annals of Statistics, 47(3):1288 – 1320.
  * Böhning, (1999) Böhning, D. (1999).  The lower bound method in probit regression.  Computational Statistics & Data Analysis, 30(1):13–17.
  * Bouchard-Côté et al., (2018) Bouchard-Côté, A., Vollmer, S. J., and Doucet, A. (2018).  The bouncy particle sampler: A nonreversible rejection-free Markov chain Monte Carlo method.  Journal of the American Statistical Association, 113(522):855–867.
  * Brooks et al., (2011) Brooks, S., Gelman, A., Jones, G., and Meng, X.-L. (2011).  Handbook of Markov Chain Monte Carlo.  CRC press.
  * Christen and Fox, (2005) Christen, J. A. and Fox, C. (2005).  Markov chain Monte Carlo using an approximation.  Journal of Computational and Graphical Statistics, 14(4):795–810.
  * Cornish et al., (2019) Cornish, R., Vanetti, P., Bouchard-Côté, A., Deligiannidis, G., and Doucet, A. (2019).  Scalable Metropolis-Hastings for exact Bayesian inference with large datasets.  In International Conference on Machine Learning, pages 1351–1360. PMLR.
  * Cui et al., (2011) Cui, T., Fox, C., and O’sullivan, M. (2011).  Bayesian calibration of a large-scale geothermal reservoir model by a new adaptive delayed-acceptance Metropolis-Hastings algorithm.  Water Resources Research, 47(10).
  * Fearnhead et al., (2024) Fearnhead, P., Nemeth, C., Oates, C. J., and Sherlock, C. (2024).  Scalable Monte Carlo for Bayesian learning.  arXiv preprint arXiv:2407.12751.
  * Fonollosa et al., (2015) Fonollosa, J., Sheik, S., Huerta, R., and Marco, S. (2015).  Reservoir computing compensates slow response of chemosensor arrays exposed to fast varying gas concentrations in continuous monitoring.  Sensors and Actuators B: Chemical, 215:618–629.
  * Gelman et al., (1997) Gelman, A., Gilks, W. R., and Roberts, G. O. (1997).  Weak convergence and optimal scaling of random-walk Metropolis algorithms.  The Annals of Applied Probability, 7(1):110–120.
  * Golightly et al., (2015) Golightly, A., Henderson, D. A., and Sherlock, C. (2015).  Delayed-acceptance particle MCMC for exact inference in stochastic kinetic models.  Statistics and Computing, 25:1039–1055.
  * Green and Han, (1992) Green, P. J. and Han, X.-l. (1992).  Metropolis methods, Gaussian proposals and antithetic variables.  In Stochastic Models, Statistical Methods, and Algorithms in Image Analysis: Proceedings of the Special Year on Image Analysis, held in Rome, Italy, 1990, pages 142–164. Springer.
  * Kelly et al., (2024) Kelly, M., Longjohn, R., and Nottingham, K. (2024).  UCI machine learning repository.
  * Korattikara et al., (2014) Korattikara, A., Chen, Y., and Welling, M. (2014).  Austerity in MCMC land: Cutting the Metropolis-Hastings budget.  In International Conference on Machine Learning, pages 181–189. PMLR.
  * Lewis and Shedler, (1979) Lewis, P. W. and Shedler, G. S. (1979).  Simulation of nonhomogeneous Poisson processes by thinning.  Naval Research Logistics Quarterly, 26(3):403–413.
  * Liu and Chen, (1995) Liu, J. S. and Chen, R. (1995).  Blind deconvolution via sequential imputations.  Journal of the American Statistical Association, 90(430):567–576.
  * Liu et al., (2015) Liu, S., Mingas, G., and Bouganis, C.-S. (2015).  An exact MCMC accelerator under custom precision regimes.  In 2015 International Conference on Field Programmable Technology (FPT), pages 120–127. IEEE.
  * Lovelace et al., (2019) Lovelace, R., Morgan, M., Hama, L., Padgham, M., Ranzolin, D., and Sparks, A. (2019).  stats 19: A package for working with open road crash data.  The Journal of Open Source Software, 4(33):1181.
  * Maclaurin and Adams, (2014) Maclaurin, D. and Adams, R. P. (2014).  Firefly Monte Carlo: Exact MCMC with subsets of data.  In Proceedings of the 30th Conference on Uncertainty in Artificial Intelligence, Quebec City, Canada, pages 543–552.
  * Neiswanger et al., (2013) Neiswanger, W., Wang, C., and Xing, E. (2013).  Asymptotically exact, embarrassingly parallel MCMC.  arXiv preprint arXiv:1311.4780.
  * Nemeth and Sherlock, (2018) Nemeth, C. and Sherlock, C. (2018).  Merging MCMC subposteriors through Gaussian Process approximations.  Bayesian Analysis, 13(2):507 – 530.
  * Peskun, (1973) Peskun, P. H. (1973).  Optimum Monte Carlo sampling using Markov chains.  Biometrika, 60(3):607–612.
  * (28) Quiroz, M., Kohn, R., Villani, M., and Tran, M.-N. (2018a).  Speeding up MCMC by efficient data subsampling.  Journal of the American Statistical Association.
  * (29) Quiroz, M., Tran, M.-N., Villani, M., and Kohn, R. (2018b).  Speeding up MCMC by delayed-acceptance and data subsampling.  Journal of Computational and Graphical Statistics, 27(1):12–22.
  * Raj et al., (2014) Raj, A., Stephens, M., and Pritchard, J. K. (2014).  fastSTRUCTURE: variational inference of population structure in large SNP data sets.  Genetics, 197(2):573–589.
  * Ripley, (2009) Ripley, B. D. (2009).  Stochastic simulation.  John Wiley & Sons.
  * Robbins and Monro, (1951) Robbins, H. and Monro, S. (1951).  A stochastic approximation method.  The Annals of Mathematical Statistics, pages 400–407.
  * Robert et al., (1999) Robert, C. P., Casella, G., and Casella, G. (1999).  Monte Carlo statistical methods, volume 2.  Springer.
  * Roberts and Rosenthal, (2001) Roberts, G. O. and Rosenthal, J. S. (2001).  Optimal scaling for various Metropolis-Hastings algorithms.  Statistical Science, 16(4):351–367.
  * Ruggles et al., (2024) Ruggles, S., Flood, S., Sobek, M., Backman, D., Chen, A., Cooper, G., Richards, S., Rogers, R., and Schouweiler, M. (2024).  IPUMS USA: Version 15.0 [dataset].  IPUMS.
  * Sampford, (1953) Sampford, M. R. (1953).  Some inequalities on Mill’s ratio and related functions.  The Annals of Mathematical Statistics, 24(1):130–132.
  * Scott et al., (2016) Scott, S. L., Blocker, A. W., Bonassi, F. V., Chipman, H. A., George, E. I., and McCulloch, R. E. (2016).  Bayes and big data: The consensus Monte Carlo algorithm.  International Journal of Management Science and Engineering Management, 11:78–88.
  * Sherlock et al., (2017) Sherlock, C., Golightly, A., and Henderson, D. A. (2017).  Adaptive, delayed-acceptance MCMC for targets with expensive likelihoods.  Journal of Computational and Graphical Statistics, 26(2):434–444.
  * Sherlock and Roberts, (2009) Sherlock, C. and Roberts, G. (2009).  Optimal scaling of the random-walk Metropolis on elliptically symmetric unimodal targets.  Bernoulli, 15(3):774 – 798.
  * Sherlock et al., (2021) Sherlock, C., Thiery, A. H., and Golightly, A. (2021).  Efficiency of delayed-acceptance random walk Metropolis algorithms.  The Annals of Statistics, 49(5):2972 – 2990.
  * Tierney, (1998) Tierney, L. (1998).  A note on Metropolis-Hastings kernels for general state spaces.  Annals of Applied Probability, pages 1–9.
  * Trojan et al., (2024) Trojan, C., Fearnhead, P., and Nemeth, C. (2024).  Diffusion generative modelling for divide-and-conquer MCMC.  arXiv preprint arXiv:2406.11664.
  * Van Rossum and Drake, (2009) Van Rossum, G. and Drake, F. L. (2009).  Python 3 Reference Manual.  CreateSpace, Scotts Valley, CA.
  * Vyner et al., (2023) Vyner, C., Nemeth, C., and Sherlock, C. (2023).  Swiss: A scalable Markov chain Monte Carlo divide-and-conquer strategy.  Stat, 12(1):e523.
  * Whiteson, (2016) Whiteson, D. (2016).  HEPMASS.  UCI Machine Learning Repository.
  * Winter et al., (2024) Winter, S., Campbell, T., Lin, L., Srivastava, S., and Dunson, D. B. (2024).  Emerging directions in Bayesian computation.  Statistical Science, 39(1):62–89.
  * Yogatama et al., (2014) Yogatama, D., Wang, C., Routledge, B. R., Smith, N. A., and Xing, E. P. (2014).  Dynamic language models for streaming text.  Transactions of the Association for Computational Linguistics, 2:181–192.
  * Yuan and Wang, (2024) Yuan, W. and Wang, G. (2024).  Markov chain Monte Carlo without evaluating the target: an auxiliary variable approach.  arXiv preprint arXiv:2406.05242.
  * Zhang et al., (2020) Zhang, R., Cooper, A. F., and De Sa, C. M. (2020).  Asymptotically optimal exact minibatch Metropolis-Hastings.  Advances in Neural Information Processing Systems, 33:19500–19510.

##  Appendix A Proof of theoretical results for MH-SS

###  A.1 Proof of Theorem 1

_Proof:_ Theorem 1 relies on a preliminary result:

######  Proposition 1.

For any $c>0$, $\Delta\in\mathbb{R}$ and $\psi>0$, let
$E(\psi):=\mathbb{E}\left[{1\wedge
c\left(\frac{\psi-\Delta}{\psi}\right)^{S}}\right]$, where
$S\sim\mathsf{Pois}(\psi)$. Then $\frac{\mathsf{d}E}{\mathsf{d}\psi}\geq 0$.

_Proof of Proposition 1:_

$E(\psi)=\exp(-\psi)\sum_{s=0}^{\infty}\frac{\psi^{s}}{s!}\left\\{1\wedge
c\left(\frac{\psi^{\prime}}{\psi}\right)^{s}\right\\}=\exp(-\psi)\sum_{s=0}^{\infty}\frac{1}{s!}\\{\psi^{s}\wedge
c(\psi-\Delta)^{s}\\}.$

Now,

$a:=\frac{\mathsf{d}}{\mathsf{d}\psi}\left\\{\psi^{s}\wedge
c(\psi-\Delta)^{s}\right\\}=s\psi^{s-1}1(\psi^{s}\leq
c(\psi-\Delta)^{s})+sc(\psi-\Delta)^{s-1}1(\psi^{s}>c(\psi-\Delta)^{s}).$

The quantity $a$ is either $s\psi^{s-1}$ or $sc(\psi-\Delta)^{s-1}$ so is at
least as large as the minimum of the two. Hence, for $s>0$,

$\frac{\mathsf{d}}{\mathsf{d}\psi}\\{\psi^{s}\wedge c(\psi-\Delta)^{s}\\}\geq
s\\{\psi^{s-1}\wedge c(\psi-\Delta)^{s-1}\\}.$

Thus

$\frac{\mathsf{d}E}{\mathsf{d}\psi}\geq\exp(-\psi)\left[\sum_{s=1}^{\infty}\frac{s}{s!}\\{\psi^{s-1}\wedge
c(\psi-\Delta)^{s-1}\\}-\sum_{s=0}^{\infty}\frac{1}{s!}\\{\psi^{s}\wedge
c(\psi-\Delta)^{s}\\}\right]=0.$

proving Proposition 1. $\square$

The maximal value of $\psi$, therefore leads to the maximal value of
$E(\psi)$.

Next, decompose the expectation over $S=(S_{1},\dots,S_{n})$ of the acceptance
rate

$\mathbb{E}\left[{\alpha_{2}(\theta,\theta^{\prime};S_{1},\dots,S_{n})}\right]=\mathbb{E}\left[{\mathbb{E}\left[{\mathbb{E}\left[{\dots\mathbb{E}\left[{\alpha_{2}(\theta,\theta^{\prime};S_{1},\dots,S_{n})|S_{2},\dots,S_{n}}\right]\dots|S_{n-1},S_{n}}\right]|S_{n}}\right]}\right].$

Applying Proposition 1 recursively with $\psi=\phi_{i}$, $i=n,\dots,1$, shows
that
$\mathbb{E}\left[{\alpha_{2}(\theta,\theta^{\prime};S_{1},\dots,S_{n})}\right]$
is maximised when $\psi$ is maximised since each of the intermediate
expectations is maximised when $\psi$ is maximised.

Finally, defining $\Delta_{i}$ as in (4),

$\frac{\mathsf{d}}{\mathsf{d}\gamma}\phi_{i}=\max[0,\Delta_{i}]-\min[0,\Delta_{i}]-c_{i}M(\theta,\theta^{\prime})=|\Delta_{i}|-c_{i}M(\theta,\theta^{\prime})\leq
0,$

so $\mathbb{E}\left[{\alpha_{2}(\theta,\theta^{\prime},S)}\right]$ is
maximised when $\phi_{i}$ is maximised, which is at $\gamma=0$. $\square$

###  A.2 Proof of Theorem 2

####  A.2.1 Preliminaries

######  Lemma 2 (Mean-Value Theorem for $a:\mathbb{R}^{d}\to\mathbb{R}$).

Let $a:\mathbb{R}^{d}\to\mathbb{R}$ be once differentiable and let
$x_{A}\in\mathbb{R}^{d}$ and $x_{B}\in\mathbb{R}^{d}$. Then
$a(x_{B})-a(x_{A})=(x_{B}-x_{A})^{\top}\nabla a(x^{\prime})$ for some
$x^{\prime}$ on the line between $x_{A}$ and $x_{B}$.

_Proof_ : Let $x=(1-t)x_{A}+tx_{B}$ and set $b(t)=a(x(t))$. The Mean-Value
Theorem for the univariate function $b$ gives
$a(x_{B})-a(x_{A})=b(1)-b(0)=b^{\prime}(t^{\prime})$ for some
$t^{\prime}\in[0,1]$. But

$\frac{\mathsf{d}b}{\mathsf{d}t}=\frac{\mathsf{d}}{\mathsf{d}t}a\left((1-t)x_{A}+tx_{B}\right)=\sum_{i=1}^{d}\left.\frac{\mathsf{d}x_{i}}{\mathsf{d}t}\right|_{t=t^{\prime}}~{}\left.\frac{\mathsf{d}a}{\mathsf{d}x_{i}}\right|_{t=t^{\prime}}=(x_{B}-x_{A})^{\top}\nabla
a(x^{\prime}),$

where $x^{\prime}=(1-t^{\prime})x_{A}+t^{\prime}x_{B}$. $\square$

Next, recall the _trapezium rule_ :
$\int_{a}^{b}f(x)\mathsf{d}x=\frac{1}{2}\\{f(a)+f(b)\\}+\frac{1}{12}(b-a)^{3}f^{\prime\prime}(\xi)$,
where $\xi\in[a,b]$. We will use a slightly different version.

######  Proposition 2.

If $f^{\prime}(x)$ is Lipschitz with a Lipschitz constant of $L$, then

$\left|\int_{a}^{b}f(x)\mathsf{d}x-\frac{1}{2}(b-a)\\{f(a)+f(b)\\}\right|\leq\frac{1}{12}(b-a)^{3}L.$

_Proof_ : Let $I=\int_{a}^{b}(x-c)f^{\prime}(x)\mathsf{d}x$, where
$c=(a+b)/2$. Then, integrating by parts,

$I=[(x-c)f(x)]_{a}^{b}-\int_{a}^{b}f(x)\mathsf{d}x=\frac{1}{2}(b-a)\\{f(a)+f(b)\\}-\int_{a}^{b}f(x)\mathsf{d}x.$

However,

$\left|\int_{a}^{b}(x-c)f^{\prime}(x)\mathsf{d}x\right|=\left|\int_{a}^{b}(x-c)\\{f^{\prime}(x)-f^{\prime}(c)\\}\mathsf{d}x\right|\leq
L\int_{a}^{b}(x-c)^{2}\mathsf{d}x,$

which evaluates to the required expression. $\square$

####  A.2.2 Proof of the theorem itself

_Proof_ : We prove the implications (15)$\Longrightarrow$(16) and
(17)$\Longrightarrow$(18); the other two results follow analogously.
Proof of (15)$\Longrightarrow$(16): Lemma 2 gives
$\ell(\theta^{\prime})-\ell(\theta)=(\theta^{\prime}-\theta)^{\top}\nabla\ell(\widetilde{\theta})$,
for some $\widetilde{\theta}$ on the line between $\theta$ and
$\theta^{\prime}$. So
$\ell(\theta^{\prime})-\ell(\theta)-r_{1}(\theta,\theta^{\prime};\widehat{\theta})=r^{\prime}$,
where
$r^{\prime}:=(\theta^{\prime}-\theta)^{\top}\left\\{\nabla\ell(\widetilde{\theta})-\nabla\ell(\widehat{\theta})\right\\}$.

Let $a(x):=\\{\nabla\ell(x)\\}^{\top}(\theta^{\prime}-\theta)$ and apply Lemma
2 again, to see that

$r^{\prime}=a(\widetilde{\theta})-a(\widehat{\theta})=(\widetilde{\theta}-\widehat{\theta})^{\top}\nabla
a(\theta^{*})=(\widetilde{\theta}-\widehat{\theta})^{\top}H(\theta^{*})(\theta^{\prime}-\widehat{\theta}).$

for some $\theta^{*}$ on the line between $\widehat{\theta}$ and
$\widetilde{\theta}$. From (15), therefore,
$r^{\prime}\leq|x^{\top}(\theta^{\prime}-\theta)|~{}|x^{\top}(\widetilde{\theta}-\widehat{\theta})|~{}M$.
Now,
$\widetilde{\theta}-\widehat{\theta}=\gamma(\theta-\widehat{\theta})+(1-\gamma)(\theta^{\prime}-\widehat{\theta})$
for some $\gamma\in[0,1]$, so

$|x^{\top}(\widetilde{\theta}-\widehat{\theta})|=\gamma|x^{\top}(\theta-\widehat{\theta})|+(1-\gamma)|x^{\top}(\theta^{\prime}-\widehat{\theta})|\leq\max\left(x^{\top}(\theta-\widehat{\theta}),x^{\top}(\theta^{\prime}-\widehat{\theta})\right),$

which leads to (16).

Proof of (17)$\Longrightarrow$(18): Let $a=\theta$, $b=\theta^{\prime}$ and
$x\in[0,1]$. Then
$f(z)=(\theta^{\prime}-\theta)^{\top}g(\theta+z(\theta^{\prime}-\theta))$ is
the component in the $\theta^{\prime}-\theta$ direction of the derivative of
$\ell$ along the line joining $\theta$ to $\theta^{\prime}$ and where $z$
gives the fraction of the total distance along the line. We also have

$f^{\prime}(z)=(\theta^{\prime}-\theta)^{\top}H\left(\theta+z\\{\theta^{\prime}-\theta\\}\right)(\theta^{\prime}-\theta),$

and from (17)

$|f^{\prime}(z_{2})-f^{\prime}(z_{1})|\leq|x^{\top}(\theta^{\prime}-\theta)|~{}|x^{\top}(\theta^{\prime}-\theta)|~{}|x^{\top}(\theta^{\prime}-\theta)|~{}L~{}|z_{2}-z_{1}|.$

Thus $f^{\prime}$ is Lipschitz. Proposition 2 applied with $a=0$ and $b=1$
then gives

$\left|\ell(\theta^{\prime})-\ell(\theta)-\frac{1}{2}(\theta^{\prime}-\theta)\left\\{g(\theta)+g(\theta^{\prime})\right\\}\right|\leq\frac{L}{12}|x^{\top}(\theta^{\prime}-\theta)|^{3}.$
(20)

Now, first order Taylor expansion of the scalar
$(\theta^{\prime}-\theta)^{\top}g(u)$ around $u=\widehat{\theta}$ gives

$\displaystyle(\theta^{\prime}-\theta)^{\top}g(\theta^{\prime})$
$\displaystyle=(\theta^{\prime}-\theta)^{\top}g(\widehat{\theta})+(\theta^{\prime}-\theta)^{\top}H(\widehat{\theta})(\theta^{\prime}-\widehat{\theta})+(\theta^{\prime}-\theta)^{\top}\\{H(\eta^{\prime})-H(\widehat{\theta})\\}(\theta^{\prime}-\widehat{\theta}),$
$\displaystyle(\theta^{\prime}-\theta)^{\top}g(\theta)$
$\displaystyle=(\theta^{\prime}-\theta)^{\top}g(\widehat{\theta})+(\theta^{\prime}-\theta)^{\top}H(\widehat{\theta})(\theta-\widehat{\theta})+(\theta^{\prime}-\theta)^{\top}\\{H(\eta)-H(\widehat{\theta})\\}(\theta-\widehat{\theta}),$

where $\eta=t\theta+(1-t)\widehat{\theta}$ and
$\eta^{\prime}=t^{\prime}\theta^{\prime}+(1-t^{\prime})\widehat{\theta}$ for
$0\leq t,t^{\prime}\leq 1$. Thus

$\frac{1}{2}(\theta^{\prime}-\theta)^{\top}\\{g(\theta)+g(\theta^{\prime})\\}-r_{2}(\theta,\theta^{\prime})=r^{\prime\prime},$

where

$r^{\prime\prime}:=\frac{1}{2}(\theta^{\prime}-\theta)\\{H(\eta^{\prime})-H(\widehat{\theta})\\}(\theta^{\prime}-\widehat{\theta})+\frac{1}{2}(\theta^{\prime}-\theta)\\{H(\eta)-H(\widehat{\theta})\\}(\theta-\widehat{\theta}).$

The triangle inequality combined with (17) then gives

$|r^{\prime\prime}|\leq\frac{L}{2}|x^{\top}(\theta^{\prime}-\theta)|+|x^{\top}(\theta-\widehat{\theta})|~{}|x^{\top}(\eta-\widehat{\theta})|+\frac{L}{2}|x^{\top}(\theta^{\prime}-\theta)|+|x^{\top}(\theta^{\prime}-\widehat{\theta})|~{}|x^{\top}(\eta^{\prime}-\widehat{\theta})|.$

However, $\eta$ is on the line between $\theta$ and $\widehat{\theta}$, so
$|x^{\top}(\eta-\widehat{\theta})|\leq|x^{\top}(\theta-\widehat{\theta})|$,
and similarly,
$|x^{\top}(\eta^{\prime}-\widehat{\theta})|\leq|x^{\top}(\theta^{\prime}-\widehat{\theta})|$.
Combining with (20) via the triangle inequality gives (18). $\square$

###  A.3 Proof of Lemma 1

_Proof_ : Let $c>0$. If $u\cdot x<0$, set $u^{\prime}\leftarrow-u$ else
$u^{\prime}\leftarrow u$. If $v\cdot x<0$ set $v^{\prime}\leftarrow-
vc||u||/||v||$ else $v^{\prime}\leftarrow cv||u||/||v||$. Then, since the
geometric mean is bounded by the arithmetic mean,

$\\{(u^{\prime}\cdot x)^{k}(v^{\prime}\cdot
x)\\}^{1/(k+1)}\leq\frac{1}{k+1}\\{ku^{\prime}\cdot x+v^{\prime}\cdot
x\\}=\frac{1}{k+1}(ku^{\prime}+v^{\prime})\cdot
x\leq\frac{1}{k+1}||ku^{\prime}+v^{\prime}||~{}||x||,$

by the Cauchy-Schwarz inequality. For any vector $z$, define
$\widehat{z}:=z/\|z\|$. Thus,

$\displaystyle\left\\{(u^{\prime}\cdot x)^{k}(v^{\prime}\cdot
x)\right\\}^{2/(k+1)}$
$\displaystyle\leq\frac{1}{(k+1)^{2}}\left\\{k^{2}||u^{\prime}||^{2}+||v^{\prime}||^{2}+2ku^{\prime}\cdot
v^{\prime}\right\\}||x||^{2}$
$\displaystyle=\frac{\|u\|^{2}}{(k+1)^{2}}\left\\{k^{2}+c^{2}+2kc\widehat{u}^{\prime}\cdot\widehat{v}^{\prime}\right\\}$
$\displaystyle\leq\frac{\|u\|^{2}}{(k+1)^{2}}\left\\{k^{2}+c^{2}+2kc|\widehat{u}^{\prime}\cdot\widehat{v}^{\prime}|\right\\},$

where the penultimate line follows as $||v^{\prime}||=c||u^{\prime}||$.
However,

$(u^{\prime}\cdot x)^{k}(v^{\prime}\cdot x)=|u^{\prime}\cdot
x|^{k}|v^{\prime}\cdot x|=c\frac{\|u\|}{\|v\|}|u\cdot x|^{k}|v\cdot x|$

and
$|\widehat{u}^{\prime}\cdot\widehat{v}^{\prime}|=|\widehat{u}\cdot\widehat{v}|=|\omega|$,
so

$\left\\{|u\cdot x|^{k}|v\cdot
x|\right\\}^{2/(k+1)}\leq\|u\|^{2k/(k+1)}\|v\|^{2/(k+1)}\|x\|^{2}\frac{h_{k}(c)}{(k+1)^{2}},$

where

$h_{k}(c):=\frac{k^{2}+c^{2}+2kc|\omega|}{c^{2/(k+1)}}.$

This is true for all $c$, and so is true at the minimiser of $h_{k}(c)$, which
solves $c^{2}+(k-1)c|\omega|=k$ and is $c=c_{k}(\omega)$ as defined in the
statement of the lemma. At $c=c_{k}(\omega)$,

$k^{2}+c^{2}+2kc|\omega|=k^{2}+(k+1)c_{k}(\omega)|\omega|+k=(k+1)\left\\{k+c_{k}(\omega)|\omega|\right\\}.$

Hence, the right hand side of the inequality is

$\|u\|^{2k/(k+1)}\|v\|^{2/(k+1)}\|x\|^{2}\frac{h_{k}(c_{k}(\omega))}{(k+1)^{2}}=\|u\|^{2k/(k+1)}\|v\|^{2/(k+1)}\|x\|^{2}\frac{k+c_{k}(\omega)|\omega|}{(k+1)c_{k}(\omega)^{2/(k+1)}}.$

Raising to the power of $(k+1)/2$ gives the required result. $\square$

###  A.4 Specifics for regression models

######  Theorem 3.

With a single observation $y$, for logistic regression
$|\mathfrak{h}^{\prime\prime}(\eta)|\leq M_{logit}(y)=1/4$ and
$|\mathfrak{h}^{\prime\prime\prime}(\eta)|\leq L_{logit}(y)=\sqrt{3}/18$, for
probit regression $|\mathfrak{h}^{\prime\prime}(\eta)|\leq M_{probit}(y)=1$,
and for Poisson regression with an expectation of
$\log\left(1+\lambda\exp(\eta)\right)$,
$|\mathfrak{h}^{\prime\prime}(\eta)|\leq M_{Pois}(y)=0.25+0.168y$.

######  Remark 3.

Plotting $\mathfrak{h}^{\prime\prime\prime}(\eta)$ for probit regression and
the relevant portion of $\mathfrak{h}^{\prime\prime\prime}(\eta)$ for Poisson
regression over a fine grid strongly suggests that for probit regression,
$|\mathfrak{h}^{\prime\prime\prime}(\eta)|<L_{probit}(y)=0.30$ and for Poisson
regression,
$|\mathfrak{h}^{\prime\prime\prime}(\eta)|<L_{Pois}(y)={\sqrt{3}}/{18}+0.061y$;
we do not formally prove these bounds.

_Proof_ : Logistic regression,
$\mathfrak{h}(\eta)=y\log[p/(1-p)]+\log(1-p)=y\eta-\log(1+\exp(\eta))$, so
$\mathfrak{h}^{\prime}(\eta)=y-\frac{\exp(\eta)}{1+\exp(\eta)}$ and

$-\mathfrak{h}^{\prime\prime}(\eta)=\frac{\exp(\eta)}{\\{1+\exp(\eta)\\}^{2}}=p(1-p).$

Hence, $0\leq-\mathfrak{h}^{\prime\prime}(\eta)\leq 1/4$. Furthermore

$\mathfrak{h}^{\prime\prime\prime}(\eta)=\frac{\exp(\eta)-\exp(2\eta)}{\\{1+\exp(\eta)\\}^{3}}.$

This is maximised at $\exp(\eta)=2-\sqrt{3}$ and minimised at
$\exp(\eta)=2+\sqrt{3}$. Substitution gives
$-\sqrt{3}/18\leq\mathfrak{h}^{\prime\prime\prime}(\eta)\leq\sqrt{3}/18$.
$\square$

_Aside_ : Notice that $|\mathfrak{h}^{\prime}(\eta)|=|y-p|\leq 1$. This leads
to the Tuna bound (no control variates).

Probit regression, $\mathfrak{h}(\eta)=y\log\Phi(\eta)+(1-y)\log\Phi(-\eta)$,

$\displaystyle\mathfrak{h}^{\prime}(\eta)$
$\displaystyle=y\frac{\phi(\eta)}{\Phi(\eta)}-(1-y)\frac{\phi(\eta)}{\Phi(-\eta)},$
$\displaystyle-\mathfrak{h}^{\prime\prime}(\eta)$
$\displaystyle=y\phi(\eta)\frac{\phi(\eta)+\eta\Phi(\eta)}{\Phi(\eta)^{2}}+(1-y)\phi(\eta)\frac{\phi(\eta)-\eta\Phi(-\eta)}{\Phi(-\eta)^{2}}.$

For $\mathfrak{h}^{\prime\prime}(\eta)$, $(y,\eta)\leftrightarrow(1-y,-\eta)$
switches the first and second terms, and only one of these terms is ever non-
zero, so it is sufficient to bound the first term.

$T(\eta)=\phi(\eta)\frac{\phi(\eta)+\eta\Phi(\eta)}{\Phi(\eta)^{2}}.$

We now show that $0\leq T(\eta)\leq 1$ by showing that it is decreasing and
finding the limits at $-\infty$ and $\infty$; from this,
$0\leq-\mathfrak{h}^{\prime\prime}(\eta)\leq 1$.

Firstly, set $g(\eta)=\phi(\eta)/\Phi(\eta)\geq 0$. Calculus gives:

$T^{\prime}=\frac{\phi(\eta)}{\Phi(\eta)^{3}}\left[1-(\phi(\eta)+\eta\Phi(\eta))(2\phi(\eta)+\eta\Phi(\eta))\right]=g(\eta)\left[1-(g(\eta)+\eta)(2g(\eta)+\eta)\right].$

Equation (4) and Section 3 of Sampford, (1953) shows that
$g(-\eta)[(g(-\eta)-\eta)(2g(-\eta)-\eta)-1]>0$ for all $\eta$; hence
$T^{\prime}(\eta)<0$ for all $\eta$.

Next, we write $T(\eta)=T_{1}(\eta)T_{2}(\eta)$, where
$T_{1}(\eta)=\frac{\phi(\eta)/\eta}{\Phi(\eta)}$ and
$T_{2}(\eta)=\frac{\eta\phi(\eta)+\eta^{2}\Phi(\eta)}{\Phi(\eta)}$. Applying
l’Hopital’s rule to $T_{1}$ gives

$\lim_{\eta\to-\infty}T_{1}(\eta)=\lim_{\eta\to-\infty}\frac{-\phi(\eta)-\phi(\eta)/\eta^{2}}{\phi(\eta)}=-1.$

Applying l’Hopital to $T_{2}$ gives

$\lim_{\eta\to-\infty}T_{2}(\eta)=\lim_{\eta\to-\infty}\frac{-\eta^{2}\phi(\eta)+\phi(\eta)+\eta^{2}\Phi(\eta)+2\eta\Phi(\eta)}{\phi(\eta)}=1+2\lim_{\eta\to-\infty}\frac{1}{T_{1}(\eta)}=-1.$

So $\lim_{\eta\to-\infty}T(\eta)=1$.

Finally, as $\eta\to\infty$, $\Phi(\eta)\to 1$ and both $\phi(\eta)$ and
$\eta\phi(\eta)\to 0$, so $T(\eta)\to 0$. $\square$

_Note_ : Böhning, (1999) gives the same bound, but does not prove the
monotonicity between the limits as $\eta\to\pm\infty$.

Poisson regression: Setting $\eta^{\prime}=\eta+\log\lambda$,
$\mathsf{d}\mathfrak{h}/\mathsf{d}\eta=\mathsf{d}\mathfrak{h}/\mathsf{d}\eta^{\prime}\times\mathsf{d}\eta^{\prime}/\mathsf{d}\eta=\mathsf{d}\mathfrak{h}/\mathsf{d}\eta^{\prime}$,
so without loss of generality we may set $\lambda=1$. Now,
$\mathfrak{h}=y\log\log[1+\exp(\eta)]-\log(1+\exp(\eta)]-\log(y!)$,

$\displaystyle\mathfrak{h}^{\prime}(\eta)$
$\displaystyle=y\frac{1}{\log[1+\exp(\eta)]}\frac{\exp(\eta)}{1+\exp(\eta)}-\frac{\exp(\eta)}{1+\exp(\eta)},$
(21) $\displaystyle-\mathfrak{h}^{\prime\prime}(\eta)$
$\displaystyle=\frac{\exp(\eta)}{\\{1+\exp(\eta)\\}^{2}}\left[y\left\\{\frac{\exp(\eta)}{\\{\log(1+\exp(\eta))\\}^{2}}-\frac{1}{\log(1+\exp(\eta))}\right\\}+1\right].$
(22)

Set

$0\leq S(\eta):=\frac{\exp(\eta)}{\\{1+\exp(\eta)\\}\log(1+\exp(\eta))}\leq
1,$

as $\log(1+x)\geq x/(x+1)$ for all $x>-1$. The coefficient of $y$ in
$-\mathfrak{h}^{\prime\prime}(\eta)$ can be written as $S(\eta)R(\eta)$, where

$R(\eta)=S(\eta)-1/(1+\exp(\eta))\geq 0,$

since $\log(1+x)\leq x$.

Firstly, by Taylor expansion,
$\log(1+\exp(x))\geq\exp(x)-\frac{1}{2}\exp(2x)$, so, for $\eta<0$,

$R(\eta)\leq\frac{1}{\\{1+\exp(\eta)\\}\\{1-\frac{1}{2}\exp(\eta)\\}}-\frac{1}{1+\exp(\eta)}=\frac{\exp(\eta)}{2+\exp(\eta)-\exp(2\eta)}\leq\frac{1}{2}\exp(\eta).$

Thus, for $\eta\leq-\log 3$, $S(\eta)R(\eta)\leq R(\eta)\leq R(-\log 3)\leq
1/6$.

Secondly, $\log(1+x)>\log(x)$, so

$R(\eta)\leq\frac{\exp(\eta)}{\eta\\{1+\exp(\eta)\\}}-\frac{1}{1+\exp(\eta)}=\frac{1}{\eta\\{1+\exp(-\eta)\\}}-\frac{1}{1+\exp(\eta)}\\\
\leq\frac{1}{\eta}.$

Hence, for $\eta\geq 6$, $S(\eta)R(\eta)\leq R(\eta)\leq R(6)\leq 1/6$.

Next, we numerically maximise $S(\eta)R(\eta)$ on the interval $[-\log 3,6]$
to obtain $SR<0.168$.

Finally, the term that does not involve $y$ is positive, with a magnitude of
at most $1/4$, which leads to the bound.

_Aside_ : To obtain the bounds for the Tuna algorithm (Zhang et al.,, 2020),
it is needed to bound $|\mathfrak{h}^{\prime}(\eta)|$. Applying l’Hopital’s
rule to the coefficient of $y_{i}$ in $\mathfrak{h}^{\prime}(\eta)$ (21),
$S(\eta)$, gives

$\lim_{\eta\to\infty}S(\eta)=0~{}~{}~{}\mbox{and}~{}~{}~{}\lim_{\eta\to-\infty}S(\eta)=1,$

and the limits of the second term in $\mathfrak{h}^{\prime}(\eta)$ are $0$ and
$-1$ as $\eta$ goes to $-\infty$ and $\infty$, respectively. So,

$\lim_{\eta\to\infty}\mathfrak{h}^{\prime}(\eta)=-1~{}~{}~{}\mbox{and}~{}~{}~{}\lim_{\eta\to-\infty}\mathfrak{h}^{\prime}(\eta)=y.$

Using $\log(1+\exp(\eta))<\exp(\eta)$ for all $\eta$, the coefficient of
$y_{i}$ in (22) is always positive. Thus,
$\mathfrak{h}^{\prime\prime}(\eta)<0$ for all $\eta$ and
$\mathfrak{h}^{\prime}(\eta)$ is decreasing and hence lies between its limits,
$-1$ and $y$. So, $|\mathfrak{h}^{\prime}(\eta)|\leq\mbox{max}(1,y)$.
$\square$

##  Appendix B Asymptotic analysis of bounds

We first analyse the tightness of the bounds from Theorem 2 and then
incorporate the computational cost of evaluating the bounds under the
assumption that the cost of evaluating the likelihood $\ell_{i}(\theta)$ is
$\mathcal{O}(d)$. We compare the costs to those of the SMH algorithm (Cornish
et al.,, 2019).

In regression examples with each $\beta_{j}\sim\mathsf{N}(0,1)$,
$j=1,\dots,d$, independently, we sampled covariates

$x_{i,j}\sim\mathsf{N}(0,1/d),~{}i=1,\dots,n,~{}j=2,\dots,d,$ (23)

independently, with $x_{i,1}=1$. The variance was chosen so that the linear
predictor, $x_{i}^{\top}\beta$ remains $\mathcal{O}(1)$ as $d$ increases.

We find that the overall costs of evaluating our algorithm’s accept-reject
probability using $r^{(1)}$ is $\mathcal{O}(d^{3/2})$, which is independent of
$n$, whereas for $r^{(2)}$ the cost is $\mathcal{O}(d^{3}n^{-1/2})$. The cost
of the SMH algorithm depends on the distribution of the covariates used. If
(as in our simulations) these are Gaussian, then each iteration of the SMH-1
algorithm costs $\mathcal{O}\left(d^{2}\log d\right)$, and for SMH-2 the cost
is $\mathcal{O}\left(d^{7/2}\\{\log d\\}^{3/2}n^{-1/2}\right)$; the terms in
$\log d$ vanish if the covariates are bounded and the terms become larger if
the covariates have heavier-than-Gaussian tails. For reference, the RWM costs
$\mathcal{O}(dn)$: reductions in the cost per-data point come at the expense
of increasing cost with dimension.

Assumptions: We assume the following:

  * •

The Bernstein-von Mises limit applies and so

$\theta\sim\mathsf{N}\left(\widetilde{\theta},\frac{1}{n}\mathbb{I}_{\mathsf{E}}^{d}(\theta_{*})^{-1}\right),$
(24)

where $\theta_{*}$ is the true parameter value and $\widetilde{\theta}$ is the
posterior mode. In regression models, for example, if the covariate vector
$X_{i}\sim\nu$, for some distribution, $\nu$, then
$\mathbb{I}_{\mathsf{E}}^{d}(\theta_{*})=\mathbb{E}_{{X\sim\nu}}\left[{\mathbb{I}_{\mathsf{E}}^{d}(\theta_{*};X)}\right]$,
the expectation over covariate vectors of the expected Fisher information
conditional on the covariates.

  * •

We consider the RWM, with a proposal of
$\theta^{\prime}\sim\mathsf{N}(\theta,\lambda^{2}_{n,d}V^{d})$ for some fixed
(as $n$ increases) positive-definite matrix, $V^{d}$.

  * •

We assume that $\Sigma^{d}:=\mathbb{I}_{\mathsf{E}}^{d}(\theta_{*})^{-1}$ is
well behaved as $d\to\infty$, satisfying the conditions in (25) and that the
user ensures that $V^{d}$ is similarly well-behaved, satisfying (26):

$\displaystyle\frac{1}{d}\mathsf{trace}(\Sigma^{d})\to t_{\pi}<\infty,$
$\displaystyle~{}\frac{1}{d}\sum_{i=1}^{d}\left(\Sigma^{d}_{i,i}\right)^{1/2}\to
t^{*}_{\pi}>0~{}~{}~{}\mbox{and}~{}~{}~{}\frac{1}{d^{2}}\sum_{i=1,j=1}^{d,d}\left(\Sigma^{d}_{i,j}\right)^{2}\to
0,$ (25) $\displaystyle\frac{1}{d}\mathsf{trace}(V^{d})\to t_{v}<\infty,$
$\displaystyle~{}\frac{1}{d}\sum_{i=1}^{d}\left(V^{d}_{i,i}\right)^{1/2}\to
t^{*}_{v}>0~{}~{}~{}\mbox{and}~{}~{}~{}\frac{1}{d^{2}}\sum_{i=1,j=1}^{d,d}(V^{d}_{i,j})^{2}\to
0.$ (26)

  * •

To keep the acceptance rate bounded as $n$ and $d$ become large, we choose

$\lambda_{n,d}=\frac{\ell}{n^{1/2}d^{1/2}}.$ (27)

The function of $n$ follows directly from (24), whereas the function of $d$
follows from Gelman et al., (1997) and the many subsequent articles on optimal
scaling of the random walk Metropolis algorithm.

  * •

We assume that $\widehat{\theta}$, the approximation to the mode, found via a
stochastic-gradient algorithm, satisfies:

$||\widehat{\theta}-\widetilde{\theta}||_{2}\leq c~{}\frac{d^{1/2}}{n^{1/2}},$
(28)

for some $c<t_{\pi}$ and not depending on $n$ or $d$. This assumption is weak:
in practice, for a reasonable approximation to the mode, we would expect
$c<<t_{\pi}$.

######  Proposition 3.

Subject to assumptions (24), (25), (26), (28) and (27):

$\displaystyle\|\theta^{\prime}-\theta\|_{2}=\Theta(n^{-1/2}),~{}~{}~{}\|\theta-\widehat{\theta}\|_{2}=\Theta(d^{1/2}n^{-1/2})$
$\displaystyle~{}~{}~{}\mbox{and}~{}~{}~{}\|\theta^{\prime}-\widehat{\theta}\|_{2}=\Theta(d^{1/2}n^{-1/2}),$
$\displaystyle\mathbb{E}\left[{\|\theta^{\prime}-\theta\|_{1}}\right]=\Theta(d^{1/2}n^{-1/2}),~{}~{}~{}\mathbb{E}\left[{\|\theta-\widehat{\theta}\|_{1}}\right]=\Theta(dn^{-1/2})$
$\displaystyle~{}~{}~{}\mbox{and}~{}~{}~{}\mathbb{E}\left[{\|\theta^{\prime}-\widehat{\theta}\|_{1}}\right]=\Theta(dn^{-1/2}),$
$\displaystyle\|\theta^{\prime}-\theta\|_{1}=\mathcal{O}(d^{1/2}n^{-1/2}),~{}~{}~{}\|\theta-\widehat{\theta}\|_{1}=\mathcal{O}(dn^{-1/2})$
$\displaystyle~{}~{}~{}\mbox{and}~{}~{}~{}\|\theta^{\prime}-\widehat{\theta}\|_{1}=\mathcal{O}(dn^{-1/2}).$

_Proof_ : We will use the following repeatedly: for any vector
$A\sim\mathsf{N}(0,W)$ for some positive definite matrix $W$:

$\mathbb{E}\left[{\sum_{i=1}^{d}A_{i}^{2}}\right]=\mathsf{trace}(W),~{}\mathsf{Var}\left[{\sum_{i=1}^{d}A_{i}^{2}}\right]=2\sum_{i,j=1}^{d,d}W_{i,j}^{2},~{}\mathbb{E}\left[{\sum_{i=1}^{d}|A_{i}|}\right]=\sqrt{\frac{2}{\pi}}\sum_{i=1}^{d}\sqrt{W_{i,i}}.$

From (26) and (27), as $d\to\infty$, for any fixed $n$,

$\displaystyle\mathbb{E}\left[{n||\theta^{\prime}-\theta||_{2}^{2}}\right]$
$\displaystyle=n\lambda^{2}_{n,d}\mathsf{trace}(V_{d})\to\ell^{2}t_{v},$
$\displaystyle\mathsf{Var}\left[{n||\theta^{\prime}-\theta||_{2}^{2}}\right]$
$\displaystyle\to 0,$

so $n||\theta^{\prime}-\theta||_{2}^{2}\to\ell^{2}t_{v}$ almost surely, giving
the bound for $\|\theta^{\prime}-\theta\|_{2}$. Analogously,
$n||\theta-\widetilde{\theta}||_{2}^{2}/d\to t_{\pi}$ almost surely, so, by
the triangle inequality,

$\frac{n^{1/2}}{d^{1/2}}\|\widetilde{\theta}-\widehat{\theta}\|_{2}-\frac{n^{1/2}}{d^{1/2}}\|\theta-\widetilde{\theta}\|_{2}\leq\frac{n^{1/2}}{d^{1/2}}\|\theta-\widehat{\theta}\|_{2}\leq\frac{n^{1/2}}{d^{1/2}}\|\widetilde{\theta}-\widehat{\theta}\|_{2}+\frac{n^{1/2}}{d^{1/2}}\|\theta-\widetilde{\theta}\|_{2}.$

In the limit as $d\to\infty$, and subject to (25) and (28), we obtain

$t_{\pi}-c\leq\frac{n^{1/2}}{d^{1/2}}\|\theta-\widehat{\theta}\|_{2}\leq
t_{\pi}+c,$

giving the bounds for $\|\theta-\widehat{\theta}\|_{2}$. The triangle
inequality then gives

$\|\theta-\widehat{\theta}\|_{2}-\|\theta^{\prime}-\theta\|_{2}\leq\|\theta^{\prime}-\widehat{\theta}\|_{2}\leq\|\theta-\widehat{\theta}\|_{2}+\|\theta^{\prime}-\theta\|_{2},$

giving the bounds for $\|\theta^{\prime}-\widehat{\theta}\|_{2}$.

For the $L_{1}$ norms, by Jensen’s inequality,

$\|A\|_{1}=\sum_{i=1}^{d}|A_{i}|=d\frac{1}{d}\sum_{i=1}^{d}|A_{i}|\leq
d\sqrt{\frac{1}{d}\sum_{i=1}^{d}A_{i}^{2}}=d^{1/2}\|A\|_{2}.$

We already know the limiting $L_{2}$ norms of $\theta^{\prime}-\theta$,
$\theta-\widehat{\theta}$ and $\theta^{\prime}-\widehat{\theta}$ so we know
that the limiting $L_{1}$ norms can be no more than a factor $d^{1/2}$ larger.
However,

$\frac{n^{1/2}}{d^{1/2}}\mathbb{E}\left[{\|\theta^{\prime}-\theta\|_{1}}\right]=\frac{n^{1/2}}{d^{1/2}}\sqrt{\frac{2}{\pi}}\frac{\ell}{n^{1/2}d^{1/2}}\sum_{i=1}^{d}\sqrt{V_{ii}}=\sqrt{\frac{2}{\pi}}\frac{\ell}{d}\sum_{i=1}^{d}\sqrt{V_{ii}}\to\ell\sqrt{\frac{2}{\pi}}t^{*}_{v}.$

Similarly

$\frac{n^{1/2}}{d}\mathbb{E}\left[{\|\theta-\widehat{\theta}\|_{1}}\right]=\frac{n^{1/2}}{d}\sqrt{\frac{2}{\pi}}\frac{1}{\sqrt{n}}\sum_{i=1}^{d}\sqrt{\Sigma_{ii}}=\sqrt{\frac{2}{\pi}}\frac{1}{d}\sum_{i=1}^{d}\sqrt{V_{ii}}\to\sqrt{\frac{2}{\pi}}t^{*}_{\pi}.$

The triangle inequality shows that
$\mathbb{E}\left[{\|\theta^{\prime}-\widehat{\theta}\|_{1}}\right]$ behaves
similarly to $\mathbb{E}\left[{\|\theta-\widehat{\theta}\|_{1}}\right]$.
$\square$

###  B.1 Our bounds

For each observation, we now apply Proposition 3 to the generic bounds from
Theorem 2, (16) and (18) passed through Lemma 1. From (23), we also have that
$\|x_{i}\|_{2}=o(1)$. Given
$\Theta(\|\theta^{\prime}-\theta\|_{2})<\Theta(\|\theta-\widehat{\theta}\|_{2})=\Theta(\|\theta^{\prime}-\widehat{\theta}\|_{2})$,

$\displaystyle\left|\ell(\theta^{\prime})-\ell(\theta)-r^{(1)}(\theta,\theta^{\prime};\widehat{\theta})\right|$
$\displaystyle\leq\Theta(\|\theta^{\prime}-\theta\|_{2})\Theta(\|\theta-\widehat{\theta}\|_{2})=\Theta\left(n^{-1/2}\right)~{}\Theta\left(d^{1/2}n^{-1/2}\right)=\Theta\left(d^{1/2}n^{-1}\right),$
$\displaystyle\left|\ell(\theta^{\prime})-\ell(\theta)-r^{(2)}(\theta,\theta^{\prime};\widehat{\theta})\right|$
$\displaystyle\leq\Theta(\|\theta^{\prime}-\theta\|_{2})\Theta(\|\theta-\widehat{\theta}\|_{2}^{2})=\Theta\left(n^{-1/2}\right)\Theta\left(d^{1/2}n^{-1/2}\right)^{2}=\Theta\left(dn^{-3/2}\right).$

The above bounds are proportional to the probability that a given observation
will be used in the accept reject stage. There are $n$ observations, and
evaluation of each $r_{i}^{(1)}$ is $\Theta(d)$, so the overall costs of
$r^{(1)}$ is
$n\times\Theta(d^{1/2}n^{-1})\times\Theta(d)=\mathcal{O}(d^{3/2})$. In
general, evaluation of each $r_{i}^{(2)}$ has a cost of $\Theta(d^{2})$. Hence
the cost for $r^{(2)}$
$n\times\Theta(dn^{-3/2})\times\Theta(d^{2})=\Theta(d^{3}/n^{1/2})$.

###  B.2 SMH bounds

For the SMH algorithm, Equation (13) of Cornish et al., (2019) has

$\displaystyle\left|\ell_{i}(\theta^{\prime})-\ell_{i}(\theta)-r_{i}^{(1)}(\theta,\theta^{\prime};\widehat{\theta})\right|$
$\displaystyle\leq\frac{1}{2}\left(||\theta-\widehat{\theta}||_{1}^{2}+||\theta^{\prime}-\widehat{\theta}||_{1}^{2}\right)\overline{U}_{2,i},$
$\displaystyle\left|\ell_{i}(\theta^{\prime})-\ell_{i}(\theta)-r_{i}^{(2)}(\theta,\theta^{\prime};\widehat{\theta})\right|$
$\displaystyle\leq\frac{1}{6}\left(||\theta-\widehat{\theta}||_{1}^{3}+||\theta^{\prime}-\widehat{\theta}||_{1}^{3}\right)\overline{U}_{3,i},$

where

$\overline{U}_{2,i}=\max_{j,k=1,\dots,d}\sup_{\theta\in\Theta}\left|\frac{\partial^{2}U_{i}}{\partial\theta_{j}\partial\theta_{k}}\right|~{}~{}~{}\mbox{and}~{}~{}~{}\overline{U}_{3,i}=\max_{j,k,l=1,\dots,d}\sup_{\theta\in\Theta}\left|\frac{\partial^{3}U_{i}}{\partial\theta_{j}\partial\theta_{k}\theta_{l}}\right|$

Repeating the total cost analysis, including the $\mathcal{O}(d)$ cost of
evaluating $r^{(1)}$ and the $\mathcal{O}(d^{2})$ cost for $r^{(2)}$, for the
scalable MH bounds, we obtain costs of
$\mathcal{O}(d^{3}n^{-1})\sum_{i=1}^{n}\overline{U}_{2,i}$ and
$\mathcal{O}(d^{5}n^{-3/2})\sum_{i=1}^{n}\overline{U}_{3,i}$ for first and
second-order control variates, respectively.

For the logistic regression model, Appendix G.1 of the SMH paper (Cornish et
al.,, 2019) derives $\overline{U}_{2,i}=\frac{1}{4}\max_{1\leq j\leq
d}|x_{i,j}|^{2}$ and $\overline{U}_{3,i}=\frac{1}{6\sqrt{3}}\max_{1\leq j\leq
d}|x_{i,j}|^{3}$. Thus, if $d^{1/2}x_{i,j}$ is bounded, so are
$d~{}\overline{U}_{2,i}$ and $d^{3/2}\overline{U}_{3,i}$. If, on the other
hand, each $x_{i,j}$ is Gaussian then
$\overline{U}_{2,i}=\mathcal{O}(d^{-1}\log d)$ and
$\overline{U}_{3,i}=\mathcal{O}\left(d^{-3/2}\\{\log d\\}^{3/2}\right)$. In
this case, we have total computational costs of $\mathcal{O}\left(d^{2}\log
d\right)$ for SMH-1 and $\mathcal{O}\left(d^{7/2}\\{\log
d\\}^{3/2}n^{-1/2}\right)$ for SMH-2.

##  Appendix C Additional simulation experiments

In this section, we present additional comparisons using ESS/$\mathbb{E}(B)$
as an efficiency metric. Figure 8 shows that the MH-SS algorithms consistently
outperform SMH, Tuna and the RWM in terms of ESS/$\mathbb{E}(B)$. Notably, the
efficiency gap between MH-SS-1 and SMH-1 is more pronounced compared to MH-
SS-2 and SMH-2, especially as $d$ increases. This discrepancy can be
attributed to two main factors. Firstly, the bounds on the log-likelihood
differences in MH-SS-1 are significantly tighter than those in SMH-1, as
demonstrated in Section 4 and Appendix B. Tighter bounds result in a smaller
$\mathbb{E}(B)$. Secondly, the acceptance rates of SMH-1 deteriorate
drastically as $d$ increases, as discussed in Section 3, whereas they tend to
remain stable for SMH-2. As efficiency metrics like ESS take into account the
autocorrelation of the posterior samples, lower acceptance rates lead to
smaller ESS. In summary, the results for SMH-1 are influenced by a combination
of loose bounds and poor acceptance rates, especially in high-dimensional
settings.

Figure 8: ESS divided by $\mathbb{E}(B)$ for MH-SS, SMH and RWM for the
logistic regression model. Both axes are presented in the logarithm base 10.
Some ESSs are omitted because $\mathbb{E}(B)\geq n$, which implies the use of
the RWM algorithm.

###  C.1 Comparison of implementations

The results presented throughout the paper consider an implementation of SMH
that evaluates the likelihood at most $B\leq n$ times in each MCMC iteration.
If the proposal $\theta^{\prime}$ is accepted, it indicates that the
likelihood was exactly evaluated $B$ times. Conversely, if $\theta^{\prime}$
is rejected, it suggests that the likelihood was most likely evaluated fewer
than $B$ times because SMH is designed to reject a proposal based on
sequential point-wise likelihood evaluations; for more details, see Algorithm
1 in Cornish et al., (2019). In contrast, both MH-SS and Tuna are originally
designed to evaluate the likelihood exactly $B$ times before deciding whether
to accept or reject the proposal.

The sequential point-wise rejection process employed by SMH implies that
likelihood evaluations cannot be computed in a vectorised fashion in
programming languages like Python and R. Consequently, while loops are
required to handle these evaluations. We point out that, for the sake of fair
run-time comparisons, the results presented in the paper for RWM, SMH, MH-SS
and Tuna were all obtained from implementations based on for-loops. In
practice, however, we find that vectorised operations in Python are often
faster compared to those relying on for-loops.

In an effort to enhance the run-time efficiency of SMH, we implemented a
vectorised approach for evaluating the likelihood. For fairness in the
comparisons, when vectorising the loops in SMH-1 and SMH-2, we also vectorised
the computations in Tuna, MH-SS-1, MH-SS-2 and in the RWM. In the ‘vectorised
SMH’ algorithm, the likelihood is evaluated exactly $B$ times regardless of
whether $\theta^{\prime}$ is accepted or rejected. We note that the impact of
SMH vectorization is exclusively observed in metrics related to run-times,
such as ESS per second. Specifically, only the denominator of the ESS per
second changes, as the ESS remains the same. Metrics like $\mathbb{E}(B)$ and
ESS/$\mathbb{E}(B)$ also remain the same. Though this modification resulted in
a significant computational speed-up, particularly for SMH-1 as it tends to
evaluate more likelihood terms, it was not enough to make SMH more
computationally efficient than MH-SS.

Figure 9 shows the results of ESS per second obtained from simulation
experiments with a logistic regression model. Panel (a) exhibits results based
on foor-loop implementations, while panel (b) shows the vectorised versions.
We stress that within each panel, the RWM, SMH, MH-SS and Tuna implementations
are either all based on for-loops or all vectorised. Overall, comparing panels
(a) and (b), we can see that the vectorised versions of SMH-1 and SMH-2 are
significantly faster than their for-loop counterparts. For instance, for
$d=50$ and $n=10^{5}$, the efficiency of SMH-1 and SMH-2 shows a 15- and
39-fold improvement, respectively, as their ESS per second jumped from
$10^{-2}$ and $10^{-0.3}$ to $10^{-0.8}$ and $10^{1.3}$, respectively.
Unsurprisingly, the gains in terms of run-times for the vectorised MH-SS are
smaller, especially for MH-SS-2 with large $n$, as it requires smaller
subsamples.

(a) For-loop implementations

(b) Vectorised implementations

Figure 9: ESS per second for SMH and MH-SS on a logistic regression model.
Except for the RWM, all algorithms have either for-loop-based implementations
(panel(a)) or are all vectorised (panel (b)). Both axes are presented in the
logarithm base 10. Some ESSs are omitted because $\mathbb{E}(B)\geq n$, which
implies the use of the RWM algorithm.

###  C.2 Empirical results of the optimal scaling for SMH

In addition to the optimal scaling results for MH-SS presented in Section 4.4,
we also performed simulation experiments to find the optimal acceptance rate
for the scalable Metropolis–Hastings algorithm (Cornish et al.,, 2019) in
order to guarantee a fair comparison against MH-SS, Tuna and RWM in the
Sections 5 and 6. This section presents the efficiency results for a set of
values for the scaling parameter $\lambda$ ranging from $0.1$ to $4$ with
increments of $0.2$. Our empirical analysis focuses primarily on the results
of the vanilla SMH-1 and SMH-2, though we also show SMH’s efficiency when
equipped with the bounds introduced in Section 4.2 (i.e., SMH-NB).

As stated in Section 4.4 and now displayed in panel (e) of Figure 10, SMH-2
has an optimal scaling around that of the RWM, which in turn leads to an
acceptance rate close to $0.234$ (panel (c)). Panel (e) also displayes that
all algorithms, except SMH-1 and SMH-1-NB, have similar acceptance rates as a
function of the scaling parameter, $\lambda$. For these algorithms, we find
that the scaling must decrease more quickly with dimension than
$\ell/\sqrt{d}$ to maintain a reasonable acceptance rate as dimension
increases. For SMH-1, however, we found that the optimal acceptance rate is
less obvious. For example, for synthetic data generated from a logistic
regression with $n=100,000$ and $d=100$, panel (a) of Figure 10 suggests that
acceptance rates from $8\%$ to $20\%$ maximise SMH-1’s efficiency. However, we