Informed Sub-Sampling MCMC: Approximate Bayesian Inference for Large Datasets

TL;DR

This paper presents ISS-MCMC, a scalable Bayesian inference method that uses data subsampling guided by summary statistics to efficiently approximate the posterior distribution in large datasets, balancing accuracy and computational cost.

Contribution

The paper introduces Informed Sub-Sampling MCMC, a novel subsampling framework that maintains simplicity and provides theoretical bias quantification for large-scale Bayesian inference.

Findings

ISS-MCMC achieves efficient approximate inference with limited computational resources.

The choice of summary statistics as MLE is theoretically justified and improves performance.

The method's bias is rigorously analyzed and shown to be acceptable in practice.

Abstract

This paper introduces a framework for speeding up Bayesian inference conducted in presence of large datasets. We design a Markov chain whose transition kernel uses an (unknown) fraction of (fixed size) of the available data that is randomly refreshed throughout the algorithm. Inspired by the Approximate Bayesian Computation (ABC) literature, the subsampling process is guided by the fidelity to the observed data, as measured by summary statistics. The resulting algorithm, Informed Sub-Sampling MCMC (ISS-MCMC), is a generic and flexible approach which, contrary to existing scalable methodologies, preserves the simplicity of the Metropolis-Hastings algorithm. Even though exactness is lost, i.e. the chain distribution approximates the posterior, we study and quantify theoretically this bias and show on a diverse set of examples that it yields excellent performances when the computational budget is limited. If available and cheap to compute, we show that setting the summary statistics as the maximum likelihood estimator is supported by theoretical arguments.

Tables

Table 1. Table 1: (Example 2 : Probit model) Comparison of the KL divergence between π 𝜋 \pi and the optimal π ~ n ∈ 𝖥 n subscript ~ 𝜋 𝑛 subscript 𝖥 𝑛 \tilde{\pi}_{n}\in\mathsf{F}_{n} ( ‖ Δ n ​ ( U ) ‖ = 3 norm subscript Δ 𝑛 𝑈 3 \|\Delta_{n}(U)\|=3 ) and other distributions in 𝖥 n subscript 𝖥 𝑛 \mathsf{F}_{n} .

$n$	$\Vert {\Delta }_n\left(U\right)\Vert$	$\mathrm{KL}\{\pi ,{\tilde{\pi }}_n\left(U\right)\}$	$B(Y,U)$
$1,000$	$3$	$0.004$	$0.04$
$1,000$	$14$	$0.11$	$0.18$
$1,000$	$23$	$0.19$	$0.29$
$100$	$33$	$0.41$	$0.54$

Table 2. Table 2: Comparison between ABC and Informed Sub-Sampling, an adaptation of ABC designed for situations where N ≫ 1 much-greater-than 𝑁 1 N\gg 1 and likelihood simulation is not possible. The exact algorithms provide samples from π 𝜋 \pi while the noisy algorithms sample from approximation of π 𝜋 \pi given in ( 21 ) and ( 22 ).

step	ABC		Informed Sub-Sampling
(i)	$\vartheta \sim p$		-
(ii)	$Z\sim f(\cdot |\vartheta )$		$Z=Y_U$, $U\sim \text{unif}\left(𝖴\right)$
(iii)	exact	noisy	exact	noisy
	if $Z=Y$,	with proba. $\propto$	if $Z=Y$	with proba. $\propto$
	set $\theta =\vartheta$	$e^{-ϵ{\Vert S\left(Y\right)-S\left(Z\right)\Vert }^2}$	draw $\theta \sim \pi \left(U\right)$	$e^{-ϵ{\Vert S\left(Y\right)-\left(N/n\right)S\left(Y_U\right)\Vert }^2}$
		set $\theta =\vartheta$		draw $\theta \sim \pi \left(U\right)$

Table 3. Table 3: (Example 4 : Logistic regression) Tradeoff Bias-Variance of the Monte Carlo estimator from Metropolis-Hastings, ISS-MCMC, Stochastic Gradient Langevin Dynamics (SGLD) Welling and Teh ( 2011 ) , the Subsampled likelihoods (SubLhd1) Bardenet et al. ( 2014 ) and the improved Confidence Sampler (SubLhd2) Bardenet et al. ( 2017 ) for a fixed computational budget of 500 seconds. Those results were replicated using 100 replications of each algorithm. Note that for SubLhd1 and SubLhd2, n 𝑛 n corresponds to the initial subset size and not to the actual subsample size that was actually used in each iteration, a parameter which is chosen by the algorithms.

algorithm	time/iter.(s)	iter. completed	RMSE	$\text{var}\left\{\widehat{\pi \left(D\right)}\right\}$
M-H	$10$	$50$	$0.1417$	$0.004$
ISS-MCMC, $n=1,000$	$0.05$	$10,000$	$0.1016$	$0.0104$
ISS-MCMC, $n=5,000$	$0.08$	$6,250$	$0.0351$	$0.0012$
ISS-MCMC, $n=10,000$	$0.13$	$3,840$	$0.0267$	$0.0007$
SGLD, $n=1,000$	$0.08$	$6,000$	$0.1370$	$0.0157$
SGLD, $n=5,000$	$0.11$	$5,250$	$0.0996$	$0.0100$
SGLD, $n=10,000$	$0.12$	$4,500$	$0.0326$	$0.0011$
SubLhd1, $n=1,000$	$1.45$	$350$	$0.0762$	$0.0042$
SubLhd1, $n=5,000$	$1.56$	$323$	$0.0680$	$0.0046$
SubLhd1, $n=10,000$	$2.24$	$223$	$0.0656$	$0.0044$
SubLhd2, $n=1,000$	$0.10$	$5,046$	$0.0304$	$0.0002$
SubLhd2, $n=5,000$	$0.14$	$3,581$	$0.0260$	$0.0006$
SubLhd2, $n=10,000$	$0.19$	$2,631$	$0.0195$	$0.0002$

Equations

$$\pi (\mathrm{d}\theta |Y_1,\dots ,Y_N)\propto f(Y_1,\dots ,Y_N|\theta )p(\mathrm{d}\theta ).$$

$$Y_i=\varphi ({\theta }_{J(i)})+{\sigma }^2{ϵ}_i,$$

$$d(t)=\sum_{j=1}^{5}\|\theta_{j}^{\ast}-\mu(\theta_{j,1:\kappa(t)})\|,\qquad\text{where}\quad\left\{\begin{array}[]{l}\forall t\in\mathbb{R},\;\kappa(t)=\max_{k\in\mathbb{N}}\{t\geq\tau_{k}\}\,,\\ \tau_{k}\;\text{is\,the\,time\,at\,the\,end\,of\,the\,}k\,\text{-th\,iteration}\,,\\ \forall k\in\mathbb{N},\;\mu(\theta_{j,1:k})=(1/k)\sum_{\ell=1}^{k}\theta_{j,\ell}\,,\end{array}\right.$$

$$f(y\,|\,\theta)=\exp\left\{g(\theta)^{T}S(y)\right\}\bigg{/}L(\theta)\,,\qquad L(\theta)=\int_{Y\in\mathsf{Y}}\exp\left\{S(y)^{T}g(\theta)\right\}\mathrm{d}y\,,$$

$$\pi(\theta\,|\,Y_{1:N})=p(\theta)\frac{\exp\left\{\sum_{k=1}^{N}S(Y_{k})^{T}g(\theta)\right\}}{L(\theta)^{N}}\bigg{/}Z(Y_{1:N})\,,$$

$$Z(Y_{1:N})=\int p(\mathrm{d}\theta )\frac{\exp \left\{{\sum }_{k=1}^{N}S(Y_k{)}^{T}g(\theta )\right\}}{L(\theta {)}^N}.$$

$${\tilde{\pi }}_n(\theta |Y_U)\propto p(\theta )f(Y_U|\theta {)}^{N/n}.$$

$$\mathrm{KL}(\pi ,\tilde{\pi })={\mathbb{E}}_{\pi }\left\{\log \frac{\pi (\theta )}{\tilde{\pi }(\theta )}\right\}.$$

$${\Delta }_n(U)=\sum _{k=1}^{N}S(Y_k)-(N/n)\sum _{k\in U}S(Y_k).$$

$$\mathrm{KL}\left\{\pi ,{\tilde{\pi }}_n(U)\right\}\le B(Y,U),$$

$$B(Y,U)=\log {\mathbb{E}}_{\pi }\exp \left\{\Vert {\mathbb{E}}_{\pi }(g(\theta ))-g(\theta )\Vert \Vert {\Delta }_n(U)\Vert \right\}$$

$$\Vert {\Delta }_n(U_1)\Vert \le \Vert {\Delta }_n(U_2)\Vert \Rightarrow B(Y,U_1)\le B(Y,U_2).$$

$${\mathsf{U}}_n^{\star }:=\left\{U\in {\mathsf{U}}_n,\frac{1}{N}\sum _{k=1}^{N}S(Y_k)=\frac{1}{n}\sum _{k\in U}S(Y_k)\right\}.$$

$$\hat{\pi }(\cdot |Y_{1:N}):=\mathcal{N}\left({\theta }^{\ast }(Y_{1:N}),I^{-1}\left({\theta }_0\right)/N\right),$$

$${\nu }_{n,ϵ}(U)\propto \exp \left\{-ϵ\Vert {\Delta }_n(U){\Vert }^2\right\},\text{for all}U\in {\mathsf{U}}_n.$$

$${\mathbb{E}}_{n,ϵ}\left\{\frac{f(Y|\theta )}{f(Y_U|\theta {)}^{N/n}}\right\}<M,$$

$$Y_{k}=\left\{\begin{array}[]{ll}1,&\text{if}\;X_{k}>0,\\ 0,&\text{otherwise}.\end{array}\right.$$

$$f(Y_k|\theta )=\alpha (\theta {)}^{Y_k}(1-\alpha (\theta ){)}^{(1-Y_k)}=\left(1-\alpha (\theta )\right){\left(\frac{\alpha (\theta )}{1-\alpha (\theta )}\right)}^{Y_k},$$

$$\pi (\theta |Y_{1:N})\propto p(\theta ){\left(1-\alpha (\theta )\right)}^N{\left(\frac{\alpha (\theta )}{1-\alpha (\theta )}\right)}^{{\sum }_{k=1}^N{Y}_k},$$

$${\tilde{\pi }}_n(\theta |Y_U)\propto p(\theta ){\left(1-\alpha (\theta )\right)}^N{\left(\frac{\alpha (\theta )}{1-\alpha (\theta )}\right)}^{(N/n){\sum }_{k\in U}{Y}_k},$$

$$\pi (\mathrm{d}\theta |Y_{1:N})\propto p(\mathrm{d}\theta )f(Y_{1:N}|\theta ).$$

$$\bar{S}(Y_U)=S(Y_U)/n,{\bar{\Delta }}_n(U)=S(Y)/N-S(Y_u)/n.$$

$$\alpha ({\theta }_i,\theta )=1\wedge a({\theta }_i,\theta ),a({\theta }_i,\theta )=\frac{\pi (\theta |Y)Q(\theta ,{\theta }_i)}{\pi ({\theta }_i|Y)Q({\theta }_i,\theta )}$$

$$\beta (U_i,U)=1\wedge b(U_i,U),b(U_i,U)=\exp \left\{ϵ\left(\Vert {\Delta }_n(U_i){\Vert }^2-\Vert {\Delta }_n(U){\Vert }^2\right)\right\}$$

$$\tilde{\alpha }({\tilde{\theta }}_i,\tilde{\theta })=1\wedge \tilde{a}({\tilde{\theta }}_i,\theta |U_{i+1}),\tilde{a}({\tilde{\theta }}_i,\tilde{\theta }|U_{i+1})=\frac{{\tilde{\pi }}_n(\tilde{\theta }|Y_{U_{i+1}})Q(\tilde{\theta },{\tilde{\theta }}_i)}{{\tilde{\pi }}_n({\tilde{\theta }}_i|Y_{U_{i+1}})Q({\tilde{\theta }}_i,\tilde{\theta })}$$

$${\hat{\pi }}_{\text{ABC}}(\mathrm{d}\theta |Y):\propto p(\mathrm{d}\theta ){\hat{f}}_{\text{ABC}}(Y|\theta )=p(\mathrm{d}\theta )\int f(\mathrm{d}Z|\theta )\exp \{-ϵ\Vert S(Z)-S(Y){\Vert }^2\},$$

$$\begin{gathered} {\hat{\pi }}_n(\mathrm{d}\theta |Y):\propto p(\mathrm{d}\theta )\hat{f}(Y|\theta )= \\ p(\mathrm{d}\theta )\sum _{U\in {\mathsf{U}}_n}{f}^{(N/n)}(Y_U|\theta )\exp \{-ϵ\Vert (N/n)S(Y_U)-S(Y){\Vert }^2\}. \end{gathered}$$

$$\Vert \pi -{\tilde{\pi }}_i\Vert =(1/2){\int }_{\Theta }|\pi (\theta )-{\tilde{\pi }}_i(\theta )|\mathrm{d}\theta .$$

$$\Vert \pi -K^i({\theta }_0,\cdot )\Vert \le C({\theta }_0){\varrho }^i.$$

$$\underset{{\theta }_0\in \Theta }{\sup }\Vert \pi -K^i({\theta }_0,\cdot )\Vert \le C{\varrho }^i.$$

Full text

1 Introduction

The development of statistical methodology that scale to large datasets represents a significant research frontier in modern statistics. This paper presents a generic and flexible approach to directly address this challenge when a Bayesian strategy is followed. Given a set of observed data $(Y_{1},\ldots,Y_{N})$, a specified prior distribution $p$ and a likelihood function $f$, estimating parameters $\theta\in\Theta$ of the model proceeds via exploration of the posterior distribution $\pi$ defined on $(\Theta,\mathcal{B}(\Theta))$ by

[TABLE]

Stochastic computation methods such as Monte Carlo methods allow one to estimate characteristics of $\pi$. In Bayesian inference, Markov chain Monte Carlo (MCMC) methods remain the most widely used strategy. Paradoxically, improvements in data acquisition technologies together with increased storage capacities, present a new challenge for these methods. Indeed, the size of the data set $N$ (along with the dimension of each observation) can become so large, that even a routine likelihood evaluation is made prohibitively computationally intensive. As a consequence, MCMC methods such as the Metropolis-Hastings algorithm (Metropolis et al., 1953) cannot be considered for reasonable runtime. This issue has recently generated a lot of research activity, see Bardenet et al. (2017) for a comprehensive review.

Most of the scalable MCMC methods proposed in the literature are based on approximations of the Metropolis-Hastings (M-H) algorithm. In the sequel, we will refer to as exact approximations, algorithms that produce samples from the target distribution when the chain is in the stationary regime, as opposed to approximate methods that do not. Central to those scalable MCMC approaches is the idea that only the calculation of the likelihood of a subset of data would be required to simulate a new state of the Markov chain. Following the development of pseudo-marginal algorithms (Andrieu and Roberts, 2009; Andrieu and Vihola, 2015), a first direction has been to replace the likelihoods in the M-H acceptance ratio by positive unbiased estimators (based on a subset of data). Although appealing since exact, this approach remains (for now) mostly theoretical because such estimators are in general not available (Jacob et al., 2015). Attempts to circumvent the positivity and unbiasedness requirements of the estimator have been studied in Quiroz et al. (2016) and Quiroz et al. (2015), respectively. In both cases, the authors resort to sophisticated control variates, which can be computationally expensive to compute.

Other authors have proposed to approximate the log-likelihood ratio by subsampling data points (Korattikara et al. (2014); Bardenet et al. (2014, 2017)), the objective being to mimic the accept/reject decision that would be achieved by the Metropolis-Hastings algorithm. Even though the resulting algorithms are not exact, the Confidence sampler proposed in Bardenet et al. (2014) and refined in Bardenet et al. (2017) is designed such that the accept/reject decision is, with an arbitrarily high probability, identical to that taken by the Metropolis-Hastings algorithm. The construction of this algorithm, based on concentration inequalities, allows to bound the L1 distance between the stationary distribution of the algorithm and $\pi$. The price to pay is that the number of likelihood evaluations is not fixed but adaptively set by the algorithm at each iteration and, as noted in Bardenet et al. (2017), it is of order $\mathcal{O}(N)$ when the chain reaches equilibrium. This number can be brought down if an accurate proxy of the log-likelihood ratio, acting as control variates, is available, as demonstrated in Bardenet et al. (2017).

More recently, a stream of research has shed light on the use of continuous time Markov processes (Zig-Zag process, Langevin diffusion) to perform Bayesian analysis of tall dataset (Bierkens et al., 2018; Pollock et al., 2016; Fearnhead et al., 2016). The computational bottleneck for this class of methods is the calculation of the gradient of the log-likelihood and it has been shown that provided that an unbiased estimate of this gradient is used, they remain exact. Here again, the use of control variates to reduce the variance of the estimator is in practice essential to reach the full potential of these methods. However, we note that those approaches represent a significant departure from the M-H algorithm and as such lose its implementational simplicity.

In this paper, we propose Informed Sub-Sampling MCMC (ISS-MCMC), a novel methodology which aims to make the best use of a computational resource available for a given computational run-time, while still preserving the celebrated simplicity of the standard M-H sampler. The state space $\Theta$ is extended with an $n$-dimensional vector of unique integers $U_{k}\subset\{1,\ldots,N\}$ identifying a subset of the data used by the Markov transition kernel at the the $k$-th iteration of the algorithm, where $n\ll N$ is set according to the available computational budget. Central to our approach is the fact that each subset is weighted according to a similarity measure with respect to the full set of data through summary statistics, in the spirit of Approximate Bayesian Computation (ABC) (see e.g.  Marin et al., 2012). The subset variable is randomly refreshed at each iteration according to the similarity measure. The Markov chain transition kernel only uses a fraction $n/N$ of the available data which is by construction –and contrary to Maclaurin and Adams (2015), Korattikara et al. (2014) and Bardenet et al. (2014)– held constant throughout the algorithm. Moreover, unlike most of the papers mentioned before, our method can be applied to virtually any model (involving i.i.d. data or not), as it does not require any assumption on the likelihood function nor on the prior distribution. Our algorithm can be cast as a noisy MCMC method since the marginal in $\theta$ of our Markov chain targets an approximation of $\pi$ that we quantify using the framework established in Alquier et al. (2016). In the special case where the data are i.i.d. realizations from an exponential model, we prove that when the summary statistics is set as the sufficient statistics, this yields an optimal approximation, in the sense of minimizing an upper bound of the Kullback-Leibler (KL) divergence between $\pi$ and the marginal target of our method. In the general case, we show that setting the summary statistics as the maximum likelihood estimator allows to bound the approximation error (in L1 distance) of our algorithm. We connect our work to a number of recent papers including Rudolf and Schweizer (2018); Huggins and Zou (2016); Dalalyan (2017) that bound approximation error of MCMC algorithms, using the Wasserstein metric.

To summarize, the main contribution of our work is to show that, under verifiable conditions, it is possible to infer $\pi$ through a scalable approximation of the M-H algorithm where the computational budget of each iteration is fixed (through the subset size $n$). To do so, it is necessary to draw the subsets according to a similarity measure with respect to the full data set and not uniformly at random, as previously explored in the literature. We show that setting the similarity measure as the squared L2 distance between the full dataset and subsample maximum likelihood estimators is supported by theoretical arguments.

Section 2 presents a striking real data example which we hope will help the reader to understand the problem we address and motivate the solution we propose, without going into further technical details at this stage. In Section 3, we provide theoretical results concerning exponential-family models, which we illustrate through a probit example. This section allows us to justify our motivations supporting the Informed Sub-Sampling general methodology which is rigorously presented in Section 4. In Section 5, we study the transition kernel of our algorithm and show that it yields a Markov chain targeting, marginally, an approximation of $\pi$. The approximation error is quantified and we provide theoretical justifications for setting up the Informed Sub-Sampling tuning parameters, including the choice of summary statistics. Finally, in Section 6, our method is used to estimate parameters of an autoregressive time series and a logistic regression model. It is also illustrated to perform a binary classification task. In the latter example, we compare the performance of our algorithm with the SubLikelihoods approach proposed in Bardenet et al. (2014).

2 An introductory example

We showcase the principles of our approach on a first real data example. The problem at hand is to infer some template shapes of handwritten digits from the MNIST database (http://yann.lecun.com/exdb/mnist/).

Example 1**.**

The data $Y_{1},Y_{2},\ldots$ are modelled by a deformable template model (Allassonnière et al., 2007). Each data $Y_{i}$ is a $15\times 15$ pixel image representing an handwritten digit whose conditional distribution given its class $J(i)\in(0,1,\ldots,9)$ is a random deformation of the template shape, parameterized by a $d=256$ dimensional vector $\theta_{J(i)}$. Assuming small deformations, the model is similar to a standard regression problem:

[TABLE]

where $Y_{i}$ is regarded as a vector $\mathbb{R}^{225}$, $\phi:\mathbb{R}^{256}\to\mathbb{R}^{225}$ is some deterministic mapping and $\sigma>0$ is the standard deviation of the additive noise $\epsilon_{i}\sim\pazocal{N}(0_{225},\mathrm{Id}_{225})$.

Given a set of $N$ labeled images $Y_{1},Y_{2},\ldots,Y_{N}$ and a prior distribution for $\theta=\{\theta_{1},\ldots,\theta_{9}\}$, one can estimate $\theta$ through its posterior distribution $\pi$, for example using the Metropolis-Hastings (M-H) algorithm (Metropolis et al., 1953). However, since the regression function $\phi$ in (2) is quite sophisticated, even a single likelihood evaluation is expensive to calculate. As a result, the M-H efficiency can be questioned as computing the $N$ likelihoods in the M-H ratio dramatically slows down each transition.

At this stage, we do not provide precise details on the Informed Sub-Sampling MCMC method but we simply provide an insight of the rationale of our approach. It designs a Markov chain whose transition kernel targets a scaled version of the posterior distribution of the parameter of interest $\theta$ given a random subset of $n$ images ($n\ll N$). More specifically, we inject in the standard M-H transition a decision about refreshing the subset of data, which, as a result, will change randomly over time. In this example, we use the knowledge of the observation labels to promote subsets of images in which the proportion of each digit is balanced.

We consider $N=10,000$ images of five digits $1,\ldots,5$, subsets of size $n=100$ and a non-informative Gaussian prior for $\theta$, as specified in Allassonnière et al. (2007). Figure 1 indicates a striking advantage of our method compared to a standard M-H using the same $N=10,000$ images. In this scenario, we allow a fixed computational budget (1 hour) for both methods and compare the estimation of the mean estimate of the two Markov chains. Qualitatively, the upper part of Figure 1 compares the estimated template shapes of the five digits at different time steps and shows that ISS-MCMC allows one to extract template shapes much quicker than the standard M-H, while still reaching an apparent similar graphical quality after one hour. This fact is confirmed quantitatively, in the lower part of Figure 1, which plots, against time and for both methods, the Euclidean distance between the Markov chain mean estimate and the maximum likelihood estimate $(\theta^{\ast}_{1},\ldots,\theta^{\ast}_{5})$ obtained using a stochastic EM (Allassonnière et al., 2007). More precisely, we compare the real valued function $\left\{d(t),\;t\in\mathbb{R}\right\}$ defined as

[TABLE]

where we have defined for $(j,k)\in\{1,\ldots,5\}\times\mathbb{N}$, $\theta_{j,k}$ as the $j$-th class parameter obtained after $k$ iterations of the Markov chains. For a vector $x\in\mathbb{R}^{n}$, $\|\cdot\|$ will refer to the usual L2 norm on $\mathbb{R}^{n}$, unless stated otherwise.

One can see that the transient phase of the Informed Sub-Sampling Markov chain is significantly shorter than that of the Metropolis-Hastings chain. More details on Example 1 can be found at Section 6. In particular, Figure 14 shows that the stationary distribution of Informed Sub-Sampling matches reasonably well $\pi$, which is a primary concern in Bayesian inference.

Our algorithm provides very encouraging results for this real data example. We motivate and formalize our method in Sections 3 and 4 and provide theoretical arguments supporting it at Section 5.

3 Approximation of the posterior distribution in exponential models: an optimality result

In this section, we consider the case of $N$ independent and identically distributed (i.i.d. ) observations from an exponential model. Sampling from the posterior distribution of such models using the Metropolis-Hastings algorithm is effortless since the information conveyed by the $N$ observations is contained in the sufficient statistics vector, which needs to be calculated only once.

The existence of sufficient statistics in this type of models allows us to establish a number of theoretical results that will be used to design and justify our Informed Sub-Sampling methodology that approximately samples from posterior distributions in general contexts, i.e non-i.i.d. observations from general likelihood models without sufficient statistics. More precisely, Propositions 1 and 2 put forward an optimal approximation of the posterior distribution $\pi$ by a distribution $\tilde{\pi}_{n}$ of the parameter of interest given only a subsample of $n$ observations. Finally, Proposition 3 justifies the introduction of a probability distribution on the set of subsamples. This is an essential element of our work as it represents a significant departure from all existing subsampling methodologies proposed in the Markov chain Monte Carlo literature, that have assumed uniform distribution on the subsamples.

3.1 Notation

Let $(Y_{1},\ldots,Y_{N})\in\mathsf{Y}^{N}$ be a set of i.i.d. observed data $(\mathsf{Y}\subseteq\mathbb{R}^{m},\,m>0)$ and define

$Y_{i:j}=(Y_{i},\ldots,Y_{j})$ if $1\leq i\leq j\leq N$ with the convention that $Y_{i:j}=\{\emptyset\}$, otherwise.

• 2.

$Y_{U}=\{Y_{k},\,k\in U\}$, where $U\subseteq\{1,\ldots,N\}$.

In this section, we assume that the likelihood model $f$ belongs to the exponential family and is fully specified by a vector of parameters $\theta\in\Theta$, ($\Theta\subseteq\mathbb{R}^{d},\,d>0$), a bounded mapping $g:\Theta\to\mathsf{S}$ and a sufficient statistic mapping $S:\mathsf{Y}\to\mathsf{S}$ ($\mathsf{S}\subseteq\mathbb{R}^{s},s>0$) such that

[TABLE]

is the density of the likelihood distribution with respect to the Lebesgue measure. The posterior distribution $\pi$ is defined on the measurable space $(\Theta,\mathcal{B}(\Theta))$ by its density function

[TABLE]

where

[TABLE]

$p$ is a prior distribution defined on $(\Theta,\mathcal{B}(\Theta))$ and with some abuse of notation, $p$ denotes also the probability density function on $\Theta$.

For all $n\leq N$, we define $\mathsf{U}_{n}$ as the set of possible combinations of $n$ different integer numbers less than or equal to $N$ and $\mathcal{U}_{n}$ as the powerset of $\mathsf{U}_{n}$. In the sequel, we set $n$ as a constant and wish to compare the posterior distribution $\pi$ (3) with any distribution from the family $\mathsf{F}_{n}=\{\tilde{\pi}_{n}(U),\,U\in\mathsf{U}_{n}\}$, where for all $U\in\mathsf{U}_{n}$, we have defined $\tilde{\pi}_{n}(U)$ as the distribution on $(\Theta,\mathcal{B}(\Theta))$ with probability density function

[TABLE]

3.2 Optimal subsets for the Kullback-Leibler divergence between $\pi$ and $\tilde{\pi}_{n}$

Recall that for two measures $\pi$ and $\tilde{\pi}$ defined on the same measurable space $(\Theta,\mathcal{B}(\Theta))$, the Kullback-Leibler (KL) divergence between $\pi$ and $\tilde{\pi}$ is defined as:

[TABLE]

Although not a proper distance between probability measures defined on the same space, $\mathrm{KL}(\pi,\tilde{\pi})$ is used as a similarity criterion between $\pi$ and $\tilde{\pi}$. It can be interpreted in information theory as a measure of the information lost when $\tilde{\pi}$ is used to approximate $\pi$, which is our primary concern here. We now state the main result of this section.

Proposition 1**.**

For any subset $U\in\mathsf{U}_{n}$, define the vector of difference of sufficient statistics between the whole dataset and the subset $Y_{U}$ as

[TABLE]

Then, the following inequality holds:

[TABLE]

where

[TABLE]

and $\|\cdot\|$ is the L2 norm.

The proof is detailed in A.1 and follows from straightforward algebra and applying Cauchy-Schwartz inequality. Note that by definition of $B$, we remark that for any two subsets $(U_{1},U_{2})\in\mathsf{U}_{n}^{2}$,

[TABLE]

The following corollary is an immediate consequence of Proposition 1.

Corollary 1**.**

Define the set:

[TABLE]

If $\mathsf{U}^{\star}_{n}$ is non-empty, then for any $U\in\mathsf{U}^{\star}_{n}$, then $\pi(\theta\,|\,Y)=\tilde{\pi}_{n}(\theta\,|\,Y_{U})$, $\pi$-almost everywhere.

A stronger result can be obtained under the assumption that a Bernstein-von Mises Theorem Van der Vaart (2000); Le Cam (1986) holds for the concentration of $\pi$ to its Normal approximation:

[TABLE]

where $\pazocal{N}$ denotes the Normal distribution, $\theta^{\ast}(Y_{1:N})=\arg\max_{\theta\in\Theta}f(Y_{1:N}\,|\,\theta)$, $\theta_{0}\in\Theta$ is some parameter and $I(\theta)$ is the Fisher information matrix given $Y_{1:N}$ at $\theta$.

Proposition 2**.**

Let $(U_{1},U_{2})\in\mathsf{U}_{n}^{2}$. Assume that for all $i\in\{1,\ldots,d\}$, $|\Delta_{n}(U_{1})^{(i)}|\leq|\Delta_{n}(U_{2})^{(i)}|$, where $|\Delta_{n}(U_{1})^{(i)}|$ refers to the $i$-th element of $\Delta_{n}(U_{1})$ (7). Then $\widehat{\mathrm{KL}}_{n}(U_{1})\leq\widehat{\mathrm{KL}}_{n}(U_{2})$, where $\widehat{\mathrm{KL}}_{n}(U)$ is the Kullback-Leibler divergence between the asymptotic approximation of the posterior $\hat{\pi}$ (11) and $\tilde{\pi}_{n}(U)$ (5).

The proof is detailed in A.2. Note that the asymptotic approximation is for $N\to\infty$ and for a fixed $n$ and is thus relevant to the context of our analysis.

3.3 Weighting the subsamples

Consider the distribution $\nu_{n,\epsilon}$ on the discrete space $\mathsf{U}_{n}$ defined for all $\epsilon\geq 0$ by:

[TABLE]

$\nu_{n,\epsilon}$ assigns a weight to any subset according to their representativeness with respect to the full dataset. When $\epsilon=0$, $\nu_{n,\epsilon}$ is uniform on $\mathsf{U}_{n}$ while when $\epsilon\to\infty$, $\nu_{n,\epsilon}$ is uniform on the set of subset(s) that minimize(s) $U\mapsto\|\Delta_{n}(U)\|$. Proposition 2 suggests that for exponential models, the optimal inference based on subsamples of size $n$ is obtained by picking the subposterior $\pi_{n}(U)$ (5) using the distribution $U\sim\nu_{n,\epsilon}$ with $\epsilon\to\infty$.

We now state Proposition 3. This result is important even though somewhat obscure at this stage. Indeed, we will show that it is a necessary condition for the method we introduce in Section 4 to converge. In fact, moving away to general models (i.e non i.i.d. and non exponential) amounts to relax the sufficient statistics existence assumption as well as the $\epsilon\to\infty$ condition. This will be achieved by constructing a class of summary statistics for the model at hand for which a similar result to Proposition 3 holds.

Proposition 3**.**

For any $\theta\in\Theta$ and $\epsilon>0$, there exists $M<\infty$ such that:

[TABLE]

where $\mathbb{E}_{n,\epsilon}$ is the expectation under $\nu_{n,\epsilon}$, as defined in (12).

The proof is postponed to A.3. Note that Proposition 3 essentially holds because $\log\nu_{n,\epsilon}$ is quadratic in $\|\Delta_{n}(U)\|$. Other weighting schemes for the subsets (e.g. uniform weights or weights $\propto\exp\{-\epsilon\|\Delta_{n}(U)\|\}$) would not necessarily allow to bound $\mathbb{E}_{n,\epsilon}\{f(Y\,|\,\theta)/f(Y_{U}\,|\,\theta)^{N/n}\}$.

3.4 Illustration with a probit model: effect of choice of sub-sample

We consider a pedagogical example, based on a probit model, to illustrate the results from the previous subsections.

Example 2**.**

A probit model is used in regression problems in which a binary variable $Y_{k}\in\{0,1\}$ is observed through the following sequence of independent random experiments, defined for all $k\in\{1,\ldots,N\}$ as:

1. (i)

Draw $X_{k}\sim\pazocal{N}(\theta^{\ast},\gamma^{2})$ 2. (ii)

Set $Y_{k}$ as follows

[TABLE]

Observing a large number of realizations $Y_{1},\ldots,Y_{N}$, we aim to estimate the posterior distribution of $\theta$. If $\gamma$ is unknown, the model is not identifiable and for simplicity we considered it as known here. The likelihood function can be expressed as

[TABLE]

where $\alpha(\theta)=\int_{0}^{\infty}(2\pi\gamma^{2})^{-1/2}\exp\{-(1/2\gamma^{2})(t-\theta)^{2}\}\mathrm{d}t$ and clearly belongs to the exponential family. The pdf of the posterior distribution $\pi$ and any distribution $\tilde{\pi}_{n}(U)\in\mathsf{F}_{n}$ writes respectively as

[TABLE]

[TABLE]

where $p$ is a prior density on $\theta$. Again, in this example, the posterior density is easy to evaluate pointwise, even when $N$ is extremely large, as it only requires to sum over all the binary variables $Y_{1},\ldots,Y_{N}$. As a consequence, samples from $\pi$ can routinely be obtained by a standard M-H algorithm and similarly for any distribution $\tilde{\pi}_{n}({U})\in\mathsf{F}_{n}$.

We simulated $N=10,000$ simulated data $Y_{1},\ldots,Y_{N}$ from (14), with true parameter $\theta^{\ast}=1$. We used the prior distribution $p=\pazocal{N}(0,10)$. In this probit model, $S$ is simply the identity function, implying that $\|\Delta_{n}(U)\|$ gives the absolute value of the difference between the scaled proportion of 1 and 0’s between the full dataset and the subset $Y_{U}$. Figure 2 reports the density functions of $\pi$ and several other distributions $\tilde{\pi}_{n}(U)\in\mathsf{F}_{n}$, for $n=100$ and $n=1,000$, with different values for the quantity $\|\Delta_{n}(U)\|$ (7). This plot, as well as the quantitative result of Table 1 are consistent with the statement of Corollary 1: when learning from a subsample of $n$ data, one should work with a subset $U$ featuring a perfect match with the full dataset, i.e $\|\Delta_{n}(U)\|=0$, or as small as possible to achieve an optimal approximation of $\pi$. Finally, Figure 3 illustrates Proposition 3: assigning the distribution $\nu_{n,\epsilon}$ (12) to the subsamples allows one to control the expectation of the likelihood ratio $f(Y\,|\,\theta)/f(Y_{U}\,|\,\theta)^{N/n}$ around 1.

4 Informed Sub-Sampling MCMC

In this section, we do not assume any particular correlation pattern for the sequence of observations, nor any specific likelihood model and simply write the posterior distribution $\pi$ as

[TABLE]

The Informed Sub-Sampling MCMC (ISS-MCMC) methodology that we describe now can be regarded as an extension of the approximation detailed in the previous section to non-exponential family models with possibly dependent observations.

4.1 Motivation of our approach

Central to our approach is the idea that all subsamples $Y_{U}$ $(U\in\mathsf{U}_{n})$ are not equally valuable for inferring $\pi$. Here, we do not assume the existence of a sufficient statistic mapping for the models under consideration. Thus, in order to discriminate between different subsamples, we introduce an artificial summary statistic mapping $S:\mathsf{Y}_{n}\to\mathsf{S}$ ($n\leq N$), where $\mathsf{S}\subseteq\mathbb{R}^{s}$. The choice of the summary statistics $S$ is problem specific and is meant to be the counterpart of the sufficient statistic mapping for general models (hence sharing, slightly abusively, the same notation). Since the question of specifying summary statistics also arises in Approximate Bayesian Computation (ABC), one can take advantage of the abundant ABC literature on this topic to find some examples of summary statistics for usual likelihood models (see e.g.  Nunes and Balding, 2010; Csilléry et al., 2010; Marin et al., 2012; Fearnhead and Prangle, 2012). More details on validation of summary statistics are discussed in Section 5.3.2.

Because the statistics used to assess the representativeness of a subsample $Y_{U}$ w.r.t. the full dataset $Y$ are only summary and not sufficient, the results of Section 3 are no longer valid. In particular, should an optimal subset $U^{\ast}$ minimising a distance between $S(Y_{U})$ and $S(Y)$ exist, inferring $\pi$ through the approximation $\tilde{\pi}_{n}(U^{\ast})$ is in no sense optimal. In fact, as shown in several examples of Section 6, this approximation is usually poor. In such a setting, it is reasonable to consider extending the set of subsamples of interest to a pool of good subsamples. This naturally suggests using the distribution $\nu_{n,\epsilon}$ (12) to discriminate between the subsamples, replacing sufficient by summary statistics and relaxing the assumption $\epsilon\to\infty$, in order to account for a collection of good subsamples. Before proceeding to the presentation of our algorithm, we define the following quantities related to a subset $U\in\mathsf{U}_{n}$:

[TABLE]

4.2 Informed Sub-Sampling MCMC: the methodology

Informed Sub-Sampling MCMC is a scalable adaptation of the Metropolis-Hastings algorithm (Metropolis et al., 1953), designed for situations when $N$ is prohibitively large to perform inference on the posterior $\pi$ in a reasonable time frame. ISS-MCMC relies on a Markov chain whose transition kernel has a bounded computational complexity, which can be controlled through the parameter $n$. We first recall how the Metropolis-Hastings algorithm produces a $\pi$-reversible Markov chain $\{\theta_{i},\,i\in\mathbb{N}\}$, for any distribution $\pi$ known up to a normalizing constant. The index $i$ is used hereafter to refer to the Markov chain iteration.

4.2.1 Metropolis-Hastings

Let $Q$ be a transition kernel on $(\Theta,\mathcal{B}(\Theta))$ and assume that the Metropolis-Hastings Markov chain is at state $\theta_{i}$. A transition $\theta_{i}\to\theta_{i+1}$ consists in the two following step:

1. (a)

propose a new parameter $\theta\sim Q(\theta_{i},\cdot)$ 2. (b)

set the next state of the Markov chain as $\theta_{i+1}=\theta$ with probability

[TABLE]

and as $\theta_{i+1}=\theta_{i}$ with probability $1-\alpha(\theta_{i},\theta)$.

Algorithm 1 details how to simulate a Metropolis-Hastings Markov chain $\{\theta_{i},\,i\in\mathbb{N}\}$.

4.2.2 Informed Sub-Sampling MCMC

To avoid any confusion, we denote by $\{\tilde{\theta}_{i},\,i\in\mathbb{N}\}$ the sequence of parameters generated by the Informed Sub-Sampling Markov chain, by contrast to the Markov chain $\{\theta_{i},\,i\in\mathbb{N}\}$ produced by the Metropolis-Hastings algorithm (Alg. 1). The pool of good subsamples used in the Informed Sub-Sampling inference is treated as a sequence of missing data $U_{1},U_{2},\ldots$ and is thus simulated by our algorithm. More precisely, ISS-MCMC produces a Markov chain $\{(\tilde{\theta}_{i},U_{i}),\,i\in\mathbb{N}\}$ on the extended space $\Theta\times\mathsf{U}_{n}$. Inspired by the analysis of Section 3, the sequence of subsamples $\{U_{i},\,i\in\mathbb{N}\}$ is randomly updated in a way that favours those subsets whose summary statistics vector is close to that of the full dataset. Let $R$ be a symmetric transition kernel on $(\mathsf{U}_{n},\mathcal{U}_{n})$, a transition $(\tilde{\theta}_{i},U_{i})\to(\tilde{\theta}_{i+1},U_{i+1})$ consists in the two following steps:

1. i-

2. (a)

propose a new subset variable $U\sim R(U_{i},\,\cdot\,)$ 2. (b)

set $U_{i+1}=U$ with probability

[TABLE]

and $U_{i+1}=U_{i}$ with probability $1-\beta(U_{i},U)$. $\Delta_{n}$ is defined at Eq. (7). 2. ii-

1. (a)

propose a new parameter $\tilde{\theta}\sim Q(\tilde{\theta}_{i},\cdot)$ 2. (b)

set $\tilde{\theta}_{i+1}=\tilde{\theta}$ with probability

[TABLE]

and as $\tilde{\theta}_{i+1}=\tilde{\theta}_{i}$ with probability $1-\tilde{\alpha}(\tilde{\theta}_{i},\tilde{\theta}\,|\,U_{i+1})$.

Algorithm 2 details how to simulate an Informed Sub-Sampling Markov chain. Note that at step 11, if $U_{i}=U_{i-1}$, the quantity $\tilde{\pi}_{n}(\theta_{i-1}\,|\,U_{i})$ has already been calculated at the previous iteration.

4.3 Connection with noisy ABC

Approximate Bayesian Computation (ABC) is a class of statistical methods, initiated in Pritchard et al. (1999), that allows one to infer $\pi$ in situations where the likelihood $f$ is intractable but forward simulation of pseudo data $Z\sim f(\,\cdot\,|\,\theta)$ is doable. More precisely, the algorithm consisting of (i) $\vartheta\sim p$, (ii) $Z\sim f(\,\cdot\,|\,\vartheta)$ and (iii) set $\theta=\vartheta$ only if $\{Z=Y\}$, does produce a sample $\theta$ whose distribution is $\pi(\,\cdot\,|\,Y)$. Regarding the situation $N\to\infty$ as a source of intractability, one could attempt to borrow from ABC to sample from $\pi$. However, since $N\to\infty$, sampling from the likelihood model is impossible and a natural idea is to replace step (ii) by drawing subsamples $Y_{U}$ ($U\in\mathsf{U}_{n}$), leading to what we refer as Informed Sub-Sampling, as opposed to Informed Sub-Sampling MCMC described in the previous Subsection. Obviously, the event $\{Y_{U}=Y\}$ is impossible except in the trivial situation where $N=n$. Overcoming situations where $\{Y=Z\}$ is impossible or very unlikely has already been addressed in the ABC literature (see Fearnhead and Prangle (2012) and Wilkinson (2013)), leading to approximate ABC algorithms. In particular, step (iii) is replaced by a step that sets $\theta=\vartheta$ with probability $\propto\exp\{-\epsilon\|S(Z)-S(Y)\|^{2}\}$ where $S$ is a vector of summary statistics and $\epsilon>0$ a tolerance parameter. We build on this analogy to propose a noisy Informed Sub-Sampling algorithm, see Table 2 for more details.

The Noisy ABC algorithm replaces direct inference of $\pi$ by the following surrogate distribution

[TABLE]

where the exact likelihood is replaced by $\hat{f}_{\text{ABC}}$. Similarly, the approximation of $\pi$ stemming from Informed Sub-Sampling is:

[TABLE]

This analogy shows that there is a connection between the ABC and the Informed Sub-Sampling in the way both approximate $\pi$, see (21) and (22). However, since sampling from $\nu_{n,\epsilon}$ and $\pi_{n}(U)$ are not feasible, this approach cannot be considered, hence motivating the use of Markov chains instead, i.e Informed Sub-Sampling MCMC. Moreover, quantifying the approximation of $\pi$ by $\hat{\pi}_{n}$ (22) is technically challenging while resorting to the Informed Sub-Samping Markov chain allows to use the Noisy MCMC framework developed in Alquier et al. (2016) to quantify this approximation. This is the purpose of the following Section.

5 Theoretical Analysis of Informed Sub-Sampling MCMC

By construction, ISS-MCMC samples a Markov chain on an extended state space $\{(\tilde{\theta}_{i},U_{i}),\,i\in\mathbb{N}\}$ but the only useful outcome of the algorithm for inferring $\pi$ is the marginal chain $\{\tilde{\theta}_{i},\,i\in\mathbb{N}\}$. In this section, we study the distribution of the marginal chain and denote by $\tilde{\pi}_{i}$ the distribution of the random variable $\tilde{\theta}_{i}$. Note that $\{\tilde{\theta}_{i},\,i\in\mathbb{N}\}$ is identical to the Metropolis-Hastings chain $\{\theta_{i},\,i\in\mathbb{N}\}$, up to replacing $\alpha$ by $\tilde{\alpha}$ in the accept/reject step. This change, from which the computational gain of our method originates, has important consequences on the stability of the Markov chain and, in particular, implies that $\pi$ is not the stationary distribution of $\{\tilde{\theta}_{i},\,i\in\mathbb{N}\}$. Interest lies in quantifying the distance between $\tilde{\pi}_{i}$ and $\pi$. In this paper, our results are expressed in total variation distance but the recent works of Rudolf and Schweizer (2018) and Johndrow and Mattingly (2017) suggest that carrying out the analysis using the Wasserstein metric may lead to more accurate bounds when $\Theta$ is not compact. We first recall the definition of the total variation distance which, for two distributions with density function $\pi$ and $\tilde{\pi}_{i}$ respectively w.r.t. the same common dominating measure, denoted $\mathrm{d}\theta$, can be expressed as

[TABLE]

5.1 Assumptions

Let $K$ denote the exact MH transition kernel, with proposal $Q$, described in Algorithm 1. $Q$ is fixed and set as a random walk kernel that achieves a reasonable acceptance rate. By construction, $K$ is $\pi$-reversible and thus $\pi$-invariant. Moreover, $K$ is assumed to be ergodic i.e $\|K(x,\,\cdot\,)-\pi\|\to 0$ at a geometric rate and the convergence is either simple (Assumption A.5.1) or uniform (Assumption A.5.1).

A 1. Geometric ergodicity

There exists a constant $\varrho\in(0,1)$ and a function $C:\Theta\to\mathbb{R}^{+}$ such that for all $(\theta_{0},i)\in\Theta\times\mathbb{N}$

[TABLE]

A 2. Uniform ergodicity

There exists two constants $C<\infty$ and $\varrho\in(0,1)$ such that for all $i\in\mathbb{N}$

[TABLE]

As observed in Remark 2 (A.6), the ISS-MCMC marginal Markov chain $\{\tilde{\theta}_{i},\,i\in\mathbb{N}\}$ is time inhomogeneous. Indeed, conditionally on $\tilde{\theta}_{i}$, the probability of the transition $\tilde{\theta}_{i}\to\tilde{\theta}_{i+1}$ depends on the iteration index $i$. This complicates the analysis of ISS-MCMC as most results on perturbation of Markov chains are established for time homogeneous Markov chains. For simplicity, we present in this section an analysis of a slight variation of ISS-MCMC that assumes independence between the different subsets $\{U_{i},\,i\in\mathbb{N}\}$ (Assumption A.5.1).

A 3. IID subsets

The subsets $U_{1},U_{2},\ldots$ are independent and identically distributed under $\nu_{n,\epsilon}$.

In practice, Assumption A.5.1 is satisfied when steps (3)-(9) of Algorithm 2 are repeated a large number of times to simulate $U_{i+1}$ given $U_{i}$. Under A.5.1, $\{\tilde{\theta}_{i},\,i\in\mathbb{N}\}$ is a time homogeneous Markov chain whose transition kernel $\tilde{K}_{n,\epsilon}$ is

[TABLE]

where for all $\theta\in\Theta$, $K(\theta,\,\cdot\,|\,u)$ is the exact MH transition kernel conditionally on some $\theta\in\Theta$ with proposal $Q$ that targets $\tilde{\pi}_{n}(\,\cdot\,|\,Y_{u})$ (5).

The results stated in this section hold under Assumption A.5.1. We nevertheless note that this assumption might be relaxed. In particular, we show how the analysis carried out under uniform ergodicity assumption can be extended even if A.5.1 does not hold, see A.6. In the geometric ergodic case, a similar extension may be doable, see e.g. (Douc et al., 2004, Theorem 8), but this is out of the scope of this paper. In general, the perturbation bounds with time inhomogeneous kernels are more obscure to interpret. Note that the numerical illustrations of ISS-MCMC presented at Section 6 were performed without satisfying Assumption A.5.1, i.e implementing Algorithm 2, and lead to satisfactory results.

Finally, we consider the following assumption for the summary statistics mapping. This assumption is motivated at two levels. First, it is necessary to have some assumptions on the summary statistics to derive theoretical results for ISS-MCMC in absence of sufficient statistics. Second, it offers a way to validate empirically the choice of summary statistics for a given model, see Section 6.

A 4. Summary Statistics

There exists a constant $\gamma_{n}<\infty$, such that for all $(\theta,U)\in\Theta\times\mathsf{U}_{n}$

[TABLE]

Assumption A.5.1 imposes a condition simultaneously on the model $f$ and the summary statistics $S$. In particular, it assumes that for any $\theta\in\Theta$, the variation of the scaled likelihood of the subsamples $Y_{U}$ $(U\in\mathsf{U}_{n})$ around $f(Y\,|\,\theta)$ is controlled by the distance between the full dataset $Y$ and the subsample $Y_{U}$, as measured through their summary statistics. This is a strong assumption which is unlikely to hold if $\Theta$ is not a compact set. It implies that even in absence of sufficient statistics, a result similar to Proposition 3 exists. One can also note that when $n\to N$, the constant $\gamma_{n}$ in Eq. (26) goes to zero.

5.2 $K$ is geometrically ergodic

Our main result is that for a sufficiently large size of subsample $n$, ISS-MCMC admits a stationary distribution. This follows from an adaptation of the work of Medina-Aguayo et al. (2016) to the context of ISS-MCMC.

Proposition 4**.**

Assume that assumptions A.5.1, A.5.1 and A.5.1 hold, then there exists an $n_{0}\leq N$ such that for all $n>n_{0}$, $\tilde{K}_{n,\epsilon}$ is also geometrically ergodic for all $\epsilon>0$.

The proof is outlined to A.4.

A direct consequence of Proposition 4, see for instance (Meyn and Tweedie, 2009, Theorem 16.0.1), is that for $n$ sufficiently large, $\tilde{K}_{n,\epsilon}$ admits a stationary distribution and that this stationary distribution converges to $\pi$ as $n\to N$. In most cases, it is difficult to obtain a rate of convergence under the assumption that $K$ is geometrically ergodic. We nevertheless note that this rate is related to rate of convergence of $\gamma_{n}$ to [math] as hinted by Medina-Aguayo et al. (2016, Theorem 4.1).

5.3 $K$ is uniformly ergodic

In addition to admitting a stationary distribution for a large enough $n$, we now show that under the assumption of uniform ergodicity it is possible to quantify the rate of convergence. Our main result follows from an adaptation of the work of Alquier et al. (2016) to the context of ISS-MCMC.

Proposition 5**.**

Define

[TABLE]

where for all $(\theta,U)\in(\Theta\times\mathsf{U}_{n})$, we have set $\phi_{U}(\theta):=f(Y_{U}\,|\,\theta)^{N/n}/f(Y\,|\,\theta)$ and

[TABLE]

Assume that assumptions A.5.1, A.5.1 and A.5.1 hold, then there exists a constant $\kappa<\infty$ such that for all $i\in\mathbb{N}$

[TABLE]

and

[TABLE]

Moreover, for a large enough subset size $n$, the marginal Markov chain produced by ISS-MCMC admits an invariant distribution $\tilde{\pi}_{n}$ that satisfies

[TABLE]

The proof of Proposition 5 is postponed to A.5. Note that an extension of this result to the case where A.5.1 does not hold is presented at A.6.

Since for any two measures $(\mu,\mu^{\prime})$, $\|\mu-\mu^{\prime}\|\leq 1$, the upper bounds of Proposition 5 are only informative if there are smaller than 1. Those bounds are a product of two expectations. We now show how those two expectations are controlled respectively through the choice of proposal kernel and the choice of summary statistics.

5.3.1 Choice of the proposal kernel

Assuming a Gaussian random walk proposal with covariance matrix $\Sigma^{T}\Sigma$, $B_{n}$ can be expressed as $B_{n}(\theta)=\sup_{U\in\mathsf{U}_{n}}D_{1}(U,\theta)$ where $D_{1}$ is defined as

[TABLE]

where $\Phi_{d}$ is the standard Gaussian distribution in dimension $d=\text{dim}(\Theta)$. When $N\gg 1$, the Bernstein-von Mises theorem states that, under conditions on the likelihood function, the posterior distribution can be approximated by a Gaussian with mean set as the maximum likelihood estimator $\theta^{\ast}$ and covariance $I(\theta_{0})^{-1}/N$ where $I$ is the Fisher information matrix and $\theta_{0}\in\Theta$ some parameter. Since ISS-MCMC aims at sampling from an approximation of $\pi$, setting $\Sigma=(1/\sqrt{N})M$ where $M^{T}M$ is an approximation of $I(\theta_{0})^{-1}$ is a reasonable choice. Proposition 6 shows that $D_{1}$ can, in this scenario, be arbitrarily brought down close to 0.

Proposition 6**.**

Under the assumption that the proposal kernel $Q$ is a Gaussian Random Walk with covariance matrix $\Sigma=(1/\sqrt{N})M$, we have

[TABLE]

where $R(x)=_{x\to 0}o(x)$ and for any square matrix $M$ of dimension $\mathbb{R}^{d}$, we have set $\|M\|_{1}:=\sum_{1\leq i,j\leq d}|M_{i,j}|$, $\|M\|_{2}:=\{\sum_{1\leq i,j\leq d}M_{i,j}^{2}\}^{1/2}$ and $|\!|\!|\,\cdot\,|\!|\!|$ is the operator norm.

The proof is postponed to A.7. Under regularity assumptions on the likelihood model, the gradient of $\log\pi$ and $\phi_{U}$ and the Hessian of $\phi_{U}$ are bounded and the upper bound of $D_{1}(U,\theta)$ can be brought down arbitrarily to [math], uniformly in $(U,\theta)$, through $M$ when $N\gg 1$.

5.3.2 Choice of the summary statistics

In Eq. (27), the likelihood of each subsample is raised at the power $N/n$ (i.e typically several orders of magnitude) and therefore subsamples unlikely under $f(\,\cdot\,|\,\theta)$ will contribute to make $A_{n}$ very large. Ideally the choice of $S$ would guarantee that subsamples $Y_{U}$ having a very small likelihood $f(Y_{U}\,|\,\theta)$ are assigned to a weight $\nu_{n,\epsilon}(U)\approx 0$ to limit their contribution. In other words, $S$ should be specified in a way that prevents $f(Y_{U}\,|\,\theta)$ to go to 0 at a rate faster than $\nu_{n,\epsilon}(U)$. This is ensured if Assumption A.5.1 holds. Indeed, in such a case

[TABLE]

where we have defined $\mathsf{E}_{n}(\theta):=\{U\in\mathsf{U}_{n},\,\sup_{\theta\in\Theta}f(Y\,|\,\theta)/f(Y_{U}\,|\,\theta)^{N/n}<1\}$ and $Z_{n}(\epsilon)=\sum_{U\in\mathsf{U}_{n}}\exp\{-\epsilon\|\Delta_{n}(U)\|^{2}\}$. Clearly, if $\epsilon$ has the same order of magnitude as $\gamma_{n}$, each term of the sum remains bounded when $\|\Delta_{n}(U)\|\to\infty$. Conversely, setting $\epsilon=0$ is equivalent to choosing $\nu_{n,\epsilon}$ as the uniform distribution on $\mathsf{U}_{n}$ and may not allow to bound $A_{n}$, see Figure 3 related to the probit example.

Potential summary statistics can be empirically validated by checking that they satisfy Assumption A.5.1. This validation can be performed graphically, by repeating the following operations for a number of parameters $\theta_{k}\sim_{\text{i.i.d}}p$:

1. (i)

draw subsets $U_{1},U_{2},\ldots$ uniformly at random in $\mathsf{U}_{n}$, 2. (ii)

plot the points with coordinates

[TABLE]

The statistics are validated if there exists $\gamma_{n}<\infty$ such that the points $(x_{k,i},y_{k,i})$ satisfy $|y_{k,i}/x_{k,i}|\leq\gamma_{n}$, as illustrated at Figure 4.

In situations where the maximum likelihood estimator $\theta^{\ast}(Y_{1:n})$ is easy and quick to evaluate numerically, we recommend setting $\bar{S}(Y_{1:n})=\theta^{\ast}(Y_{1:n})$. In the case of independent observations of a well-specified model, setting the summary statistics as the maximum likelihood estimate is justified by the following Proposition which implies that Assumption A.5.1 holds, asymptotically, up to a constant.

Proposition 7**.**

We assume that the whole dataset comprises $N=\rho n$ independent observations and there exists some $\theta_{0}\in\Theta$ such that $Y_{i}\sim f(\,\cdot\,|\,\theta_{0})$. Let $\theta^{\ast}$ be the MLE of $Y_{1},\ldots,Y_{N}$ and $\theta^{\ast}_{U}$ be the MLE of the subsample $Y_{U}$ ($U\in\mathsf{U}_{n}$). Then, there exists a constant $\beta$, a metric $\|\cdot\|_{\theta_{0}}$ on $\Theta$ and a non-decreasing subsequence $\{\sigma_{n}\}_{n\in\mathbb{N}}$, ($\sigma_{n}\in\mathbb{N}$) such that for all $U\subset\{1,2,\ldots,\rho\sigma_{n}\}$ with $|U|=\sigma_{n}$, we have for $p$-almost all $\theta$ in an neighborhood of $\theta_{0}$:

[TABLE]

where

[TABLE]

The proof is detailed in B and follows from a careful application of a Bernstein-von Mises theorem. Note that an extension of Proposition 7 to cases where the observations are not independent may exist provided that a Bernstein-von Mises theorem holds for the model at hand, which is the case for dependent observations if the likelihood model satisfies local asymptotic normality conditions Le Cam (1953, 1986).

Note that since Assumption A.5.1 is mostly used in Propositions 4 and 5 to guarantee that the log-likelihood ratio between the likelihood and the scaled likelihood of a subsample is bounded, the constant $\beta$ in Proposition 7 is not a major concern. In addition, it is straightforward to see that this constant vanishes when the subset size grows faster than the full dataset, i.e $\rho\downarrow 1$, once in the asymptotic regime of Eq. (34).

We remark that Proposition 7 is in line with the results regarding optimal summary statistics for ABC established in Fearnhead and Prangle (2012). The authors show that the quadratic error loss between the ABC estimate based on $\hat{\pi}_{\text{ABC}}$ (Eq. 21) and the true parameter is minimized when setting the summary statistics as the posterior mean, a choice which asymptotically coincides with the maximum likelihood estimator.

Finally, we note that, similarly to any approximate MCMC method, ISS-MCMC does not guarantee Law of Large number for $\pi$-integrable functionals. However, assume that the MH chain $K$ is geometrically ergodic, it is straightforward to establish that for a large enough $n$,

[TABLE]

where $\pi f:=\int f\mathrm{d}\pi$ and $\|f\|=\sup_{\theta\in\Theta}|f(\theta)|$. If in addition, $K$ is uniformly ergodic Proposition 5 can help to bound the asymptotic error that typically arises in MCMC estimation of $\pi f$.

6 Illustrations

We evaluate the efficiency of ISS-MCMC on three different applications: inferring a time series observed at $N=10^{6}$ contiguous time steps, a logistic regression with $N=10^{6}$ observations and a Gaussian binary classification problem based on $N=10^{7}$ data.

6.1 Implementation details of Informed Sub-Sampling MCMC

Before illustrating the ISS-MCMC algorithm on the different examples, we address a few technical implementation details.

On the subset size $n$: this parameter is essentially related to the computational budget available to the user. In the following examples we have used $n\propto N^{1/2}$ which achieves a substantial computational gain at a price of a negligible asymptotic bias.

• 2.

On the sufficient statistics $S$: to reduce the bias resulting from the Metropolis-Hastings approximation, $S$ should be constructed so that Assumption A.5.1 holds. If the maximum likelihood estimator $\theta^{\ast}(Y)$ is quick to compute then Proposition 7 suggests that setting $S(Y)=\theta^{\ast}(Y)$ will theoretically satisfy A.5.1. Other sufficient statistics mapping can be used, typically those arising in the Approximate Bayesian Computation literature. In any case, we recommend checking A.5.1 graphically (see Section 5).

• 3.

On the bandwidth parameter $\epsilon$: the theory shows that when $\epsilon\approx\gamma_{n}$, the asymptotic bias is controlled ($\gamma_{n}$ is the constant in A.5.1). In practice, this may prove to be too large and could potentially cause the algorithm to get stuck on a very small number of subsets. To avoid such a situation, we suggest monitoring the refresh rate of subsamples that should occur with probability of at least 1%.

• 4.

On the initial subset $U_{0}$: in theory, one would run a preliminary Markov chain $\{U_{1}^{(0)},\ldots,U_{L}^{(0)}\}$ (for some $L>0$) targeting $\nu_{n,\epsilon}$, and set $U_{0}=U_{L}^{(0)}$ in order for the results of Section 5.3.2 to hold. In practice, a more efficient approach is to use a simulating annealing Metropolis-Hastings algorithm, see Geyer and Thompson (1995). It introduces a sequence of tempered distributions $\nu_{k}:=\nu_{n,\epsilon_{k}}$, such that $\epsilon_{k}=t_{k}\epsilon$ $(k\in\{1,\ldots,L\})$ where $t_{1}=0$ and $t_{L}=1$. The transition kernel of the $k$-th iteration of the preliminary Markov chain is designed to be $\nu_{k}$ invariant. This technique facilitates sampling from a proxy of $\nu_{n,\epsilon}$ in a relative short time period as the successive tempered distributions help identifying those subsamples belonging to the high probability sets of $\nu_{n,\epsilon}$.

6.2 Inference of an AR(2) model

Example 3**.**

An autoregressive time series of order 2 AR(2) $\{Y_{k},\,k\leq N\}$ is defined recursively by:

[TABLE]

where $\theta\in\Theta\subset\mathbb{R}^{3}$. The likelihood of an observed time series for this model writes

[TABLE]

such that for all $k\geq 1$,

[TABLE]

where $x\to\Phi_{1}(x;\,m,v)$ is the pdf of the univariate Gaussian distribution with mean $m$ and variance $v$.

This model has been used in Chib and Greenberg (1995, Section 7.2) to showcase the Metropolis-Hastings (M-H) algorithm. We follow the same setup here and in particular we use the same true parameter $\theta^{\ast}=(1,-.5,1)$, same prior distribution and proposal kernel $Q$ ; see Chib and Greenberg (1995) for more details. We sampled a time series $\{Y_{k},\,k\leq N\}$ according to (36), with $N=10^{6}$ under $\theta^{\ast}$. Of course, in such a setup, M-H is prohibitively slow to be used in practice to sample from $\pi$ as it involves evaluating the likelihood of the whole time series at each iteration. We nevertheless use M-H to obtain a ground truth of $\pi$.

For simplicity, we restrict the set of subsamples to $n$ contiguous observations:

[TABLE]

This induces a set of subset $\bar{\mathsf{U}}_{n}\subset\mathsf{U}_{n}$ defined such that a subset $U\in\bar{\mathsf{U}}_{n}$ is identified with its starting index, i.e for all $i\leq|\bar{\mathsf{U}}_{n}|$, $\bar{\mathsf{U}}_{n}\ni U_{i}:=\{i,i+1,\ldots,i+n-1\}$. Indeed, using such subsamples yields a tractable likelihood (37) as otherwise, missing variables need to be integrated out, hence loosing the simplicity of our approach.

With some abuse of notation, the proposal kernel $R$ can be written as a transition kernel on the alphabet $\{0,\ldots,N-n+1\}$. It is defined in this example as:

[TABLE]

The rationale is to propose a new subset through a mixture of two distributions: the first gives higher weight to local moves and the latter allows jumps to remote sections of the time series. In this example, we have used $\omega=0.9$ and $\lambda=0.1$. We study the efficiency of ISS-MCMC in function of $n$, $\epsilon$ and $S$.

For any subsample $Y_{U}$, $U\in\bar{\mathsf{U}}_{n}$, we have set the summary statistics $\bar{S}(Y_{U})$ to the solution of the AR(2) Yule-Walker equations for the dataset $Y_{U}$. As shown in Figure 4, this choice of summary statistics satisfies (graphically) A.5.1 with $\gamma\approx 5\,10^{6}$. We therefore set $\epsilon=5.0\,10^{6}$ to make sure that $A_{n}$ (27) is bounded. Theoretically, Proposition 5 guarantee that the bias is controlled. This is illustrated graphically in Figure 5 where $\pi$ is compared to $\tilde{\pi}_{i}$ ($i=50,000$). We also report the distribution of the Informed Sub-Sampling chain when the $\epsilon=0$, i.e when the subsampling is actually uninformed and all subsamples have the same weight. In the latter case, $A_{n}$ is not bounded which explains why the bias on $\lim_{i\to\infty}\|\pi-\tilde{\pi}_{i}\|$ is not controlled. Figures 5 and 6 illustrate the distribution of $\{\tilde{\theta}_{i},\,i\in\mathbb{N}\}$ for some runs of ISS-MCMC with $\epsilon=0$ and $\epsilon=5.0\,10^{6}$. Finally, Figure 9 gives a hint at the computational efficiency of ISS-MCMC. Metropolis-Hastings was compared to ISS-MCMC with $n\in\{1,000;5,000;10,000\}$ and $\epsilon\in\{0;1;5.0\,10^{6}\}$. The performance indicator is defined as the average of the marginals Total Variation distance, i.e

[TABLE]

where $\pi^{(j)}$ and $\tilde{\pi}_{t}^{(j)}$ are respectively the true $j$-th marginal and the $j$-th marginal of the chain distribution after a runtime of $t$ seconds. The true marginals were estimated from a long Metropolis-Hastings chain, at stationarity. $\tilde{\pi}_{t}^{(j)}$ was estimated using 500 independent chains starting from the prior. The Matlab function ksdensity in default settings was applied to estimate $\|\pi^{(j)}-\tilde{\pi}_{t}^{(j)}\|_{\mathrm{TV}}$ from the chains samples, hence the variability. On the one hand, when $\epsilon=0$, there is no informed search for subsamples which makes the algorithm much faster than the other setups but yields a significant bias (larger than $0.5$). On the other hand, setting $\epsilon>0$ adds to the computational burden but allows to reduce the bias. In fact, for $n=5,000$ and a computational budget of $t=1,000$ seconds, the bias of ISS-MCMC is similar to that of Metropolis Hastings but converges 100 times as fast. Finally, note that as expected by the theory, when the unrepresentative subsamples are not penalized enough (e.g. by setting $\epsilon=1$), ISS-MCMC yields a significant bias which hardly improves on uninformed subsampling when $n=10,000$. Following Assumption A.5.1, we see that setting $\epsilon=5.0\,10^{6}$ significantly reduces the bias. Note that setting $\epsilon>5.0\,10^{6}$ could potentially reduce further the bias but may fail the algorithm: indeed, when $\epsilon$ is too large the chain $\{U_{k},\,k\in\mathbb{N}\}$ gets easily stuck on a set of the best subsamples (for this choice of summary statistics) and may considerably slow down the convergence of the algorithm.

Other choices of summary statistics can be considered. Since $Y$ is modelled as an autoregressive time series, an option would be to set the summary statistics as the empirical autocorrelation function. Figure 8 shows that it is not a recommended choice. The left panel suggests that Assumption A.5.1 does not hold for this type of summary statistics: good subsets yield a large value for $\log f(Y\,|\,\theta)-(N/n)\log f(Y_{U}\,|\,\theta)$ and conversely for bad subsets, hence generating a bias (see Eq. (30)). As a consequence, the right panel which shows a clear mismatch between $\pi^{(3)}$ and the Informed Sub-Sampling third marginal is not surprising.

6.3 Logistic Regression Example

Example 4**.**

A $d$-dimensional logistic regression model is parameterized by a vector $\theta=(\theta_{1},\ldots,\theta_{d})\in\Theta\subset\mathbb{R}^{d}$. Observations are realizations of the following model:

simulate covariates $X_{i}=(X_{i,1},\ldots,X_{i,d})\sim\pazocal{N}(0,(1/d)^{2})$

• 2.

simulate $Y_{i}$ given $\theta$ and $X_{i}$ as

[TABLE]

We have simulated $N=10^{6}$ observations $Y_{1},Y_{2},\ldots$ under the true parameter $\theta^{\ast}=(1,2,-1)$ $(d=3)$.

The summary statistics were set as the maximum likelihood estimator returned by the Matlab routine glmfit and were graphically validated, as in Figure 4. The tolerance parameter was consequently set $\epsilon=5.0\,10^{6}$. We study the influence of $n$ on the Informed Sub-Sampling chain marginal distributions in Figure 10. We note that as soon as $n\geq 5,000$, the bias vanishes and that when random subsampling is used (i.e $\epsilon=0$), the bias is much larger. Of course, Figure 10 only gives information about the marginal distributions. To complement the study, we consider estimating the probability $\pi(D)$ where $D$ is the domain defined as:

[TABLE]

in order to check that the joint distribution $\pi$ is reasonably inferred. Numerical integration using a long Metropolis-Hastings algorithm, gave the ground truth $\pi(D)=0.1$. The top panel of Figure 11 illustrates the Monte Carlo estimation of $\pi(D)$ based on $i=10,000$ iterations of ISS-MCMC implemented with $n\in\{1,000\,;\,5,000\,;\,10,000\}$ and compares it to Metropolis-Hastings. As expected ISS-MCMC has a negligible bias and the variance of the estimator decreases when $n$ increases. Indeed, when $n$ increases, the Informed Sub-Sampling process is less likely to pick irrelevant subsets, which in turns lower the variability of the chain. The Monte-Carlo estimation based on ISS-MCMC with $n=10,000$ and Metropolis-Hastings are very similar. However, when we normalize the experiment by the CPU time, Metropolis-Hastings is clearly outperformed by ISS-MCMC. The lower panel of Figure 11 assumes that only $t=500$ seconds of computation are available. All the chains are started from $\theta^{\ast}$. Table 3 reports the quantitative details of this experiment. In such a situation, one should clearly opt for the Informed Sub-Sampling approach as the Metropolis-Hastings algorithm only achieves 50 iterations for this amount of computation and as such fails to reach stationarity.

We also compare ISS-MCMC with two other algorithms that approximate the Metropolis-Hastings algorithm by using subset of data, drawn, unlike ISS-MCMC, uniformly at random. More precisely, we have implemented the Stochastic Gradient Langevin Dynamic (SGLD) from Welling and Teh (2011) and the Subsampled likelihoods MH (SubLhd1) algorithm from Bardenet et al. (2014) along with an improved version of this algorithm that makes use of control variates, referred to as the Improved Confidence sampler in Bardenet et al. (2017) but abbreviated here as SubLhd2 for simplicity. Those algorithms have been implemented in their default version, following the parameterization prescribed in their original article. All those methods are inexact and we are interested in comparing the bias/variance tradeoff per CPU time unit. Results in terms of convergence in distribution and Monte Carlo estimation are reported respectively in Figure 10 and Table 3. For this model, SGLD and SubLhd1 show a larger bias than ISS-MCMC and SubLhd2: they need larger subset size $n$ to achieve a similar precision than ISS-MCMC or SubLhd2, see Figure 10. SubLhd2 seems to outperform ISS-MCMC when $n$ is low in terms of distribution bias but the two methods perform equally good when $n\geq 5,000$. Quantitatively, the Monte Carlo estimation of $\pi(D)$ appears better with SubLhd2 than any other method for any subset size, as indicated by the RMSE reported at Table 3. However, looking at the comparative boxplot representing the distribution of the Monte Carlo estimator of $\pi(D)$ in time normalized experiments (Figure 12), one can see that when $n$ is larger than $5,000$, estimators from ISS-MCMC and SubLhd2 are quite similar confirming the qualitative impression of Figure 10 and perhaps moderating the RMSE-based assessment made at Table 3.

6.4 Binary Classification

Example 5**.**

A training dataset consisting of $N=10^{7}$ labeled observations $Y=\{Y_{k},\,k\leq N\}$ from a 2 dimensional Gaussian mixture model is simulated with

[TABLE]

where $\mu_{1}=[\theta_{1}\,,\,0]$, $\mu_{2}=[\theta_{2}\,,\,0]$, $\Gamma_{1}=\mathrm{diag}([\theta_{3}/2\,,\,\theta_{3}])$ and $\Gamma_{2}=\mathrm{diag}([\theta_{4}/2\,,\,\theta_{4}])$. We define $\theta=(\theta_{1},\theta_{2})$ with $\theta_{i}\in\mathbb{R}\times\mathbb{R}_{*}^{+}$ for each model $i\in\{1,2\}$. A prior distribution $(\theta_{1},\theta_{2})\sim_{\textit{i.i.d.}\,}p:=\pazocal{N}(0,1/2)\otimes\Gamma(1,2)$ ($\Gamma(a,b)$ is the Gamma distribution with shape $a$ and rate $b$) is assigned to $\theta$. Consider an algorithm a that simulates a Markov chain

[TABLE]

targeting the posterior distribution of $\theta$ given $Y$, perhaps approximately. We consider the real-time supervised classifier $I_{\texttt{a}}^{\ast}(t)$, driven by $\theta_{\texttt{a}}$, for the test dataset $Y^{\ast}=\{Y_{k}^{\ast},\,k\leq N_{\text{test}}\}$ ($N_{\text{test}}=10^{4}$) and defined as:

[TABLE]

where $\kappa_{\texttt{a}}(t)=\sup_{k\in\mathbb{N}}\{\tau_{k}^{\texttt{a}}\leq t\}$ and $\bar{\theta}_{i,k}^{(\texttt{a})}=(1/k)\sum_{\ell=1}^{k}\theta_{i,\ell}^{(\texttt{a})}$. We have defined $\tau_{k}^{\texttt{a}}$ as the wall clock time to generate $k$ iterations of algorithm a. We define the live classification error rate as $\epsilon_{\texttt{a}}(t)=\|I_{\texttt{a}}^{\ast}(t)-I^{\ast}\|_{1}$ where $I_{\texttt{a}}^{\ast}=(I_{\texttt{a},1}^{\ast},\ldots,I_{\texttt{a},N}^{\ast})$ and $I_{k}^{\ast}$ is the true class of $Y_{k}^{\ast}$. We compare $\epsilon_{\texttt{a}}$ for three different algorithms a: ISS-MCMC, Metropolis-Hastings and Subsampled Likelihoods (Bardenet et al., 2014).

In this simulation example, we have used the true value $\theta^{\ast}=(-1,1/2,1,1/2)$ and simulated $Y$ such that it contains the same number of observations from model 1 and model 2 i.e $N/2$. The three algorithms were implemented with the same proposal kernel, namely a single site random walk with adaptive variance that guarantees an acceptance rate between 0.40 and 0.50, see Roberts et al. (2001); Haario et al. (2001). ISS-MCMC was implemented with parameters $n=1,000$ and $\epsilon=10^{7}$. The summary statistics were taken as $S(Y_{U})=0$ if $\sum_{k\in U}\mathds{1}_{\{I_{k}=1\}}\neq\sum_{k\in U}\mathds{1}_{\{I_{k}=2\}}$ and

[TABLE]

otherwise. This choice allows to keep the right proportion of data from the two models in any subsample used for the inference. The statistics in (42) are sufficient for each model, taken separately. Subsampled Likelihoods was implemented with the default parameters prescribed in the introduction of Section 4 in Bardenet et al. (2014).

Figure 13 compares the live classification error rate achieved by the three algorithms. We also report the optimal Bayes classifier which achieves $\epsilon_{B}(t)=0.0812$ classifying $Y_{k}^{\ast}$ in class 1 if $Y_{k,1}^{\ast}<0$ and in class 2 alternatively. Unsurprisingly, Metropolis-Hastings is penalized because it evaluates the norm of a $N=10^{7}$ dimensional vector at each iteration. Subsampled Likelihoods does slightly better than M-H but suffers from the fact that close to stationary regime, the algorithm ends up drawing the quasi-entire dataset with high probability, a fact which was explained in Bardenet et al. (2014).

6.5 Additional details for the handwritten digit inference (Example 1)

In the handwritten digit example (Example 1), we have used batches of $n=100$ data. Since the initial dataset comprises $2,000$ observations per digit, the summary statistics were defined in a way that any subsample contains $20$ observations from each class. More precisely, we have set for any subsample $Y_{U}$, $S(Y_{U})=0$ if for at least one class $i\in\{1,\ldots,5\}$, $(1/n)\sum_{k\in U}\mathds{1}_{\{J(k)=i\}}\neq 20$ and

[TABLE]

otherwise. We set the bandwidth to $\epsilon=10^{5}$. The proposal kernels of the Informed Sub-Sampling chain and M-H were defined as the same Random Walk kernel. In particular, at each iteration only a bloc of the template parameter of one of the five classes is updated. The variance parameter of the Random Walk is adapted according to the past trajectory of the chain, so as to maintain an acceptance rate of $.25$.

Figure 14 reports the empirical marginal distribution of one component for each vector $\theta_{1},\ldots,\theta_{5}$ obtained from ISS-MCMC and from M-H. Those distributions are estimated from $50,000$ iterations of both algorithms, in stationary regime. This shows that the distribution of those parameters are in line with each other.

7 Conclusion

When the available computational budget is limited, inferring a statistical model based on a tall dataset in the Bayesian paradigm using the Metropolis-Hastings (M-H) algorithm is not computationally efficient. Several variants of the M-H algorithm have been proposed to address this computational issue (Bardenet et al., 2014; Banterle et al., 2015; Korattikara et al., 2014; Maclaurin and Adams, 2015). However, (i) they often lose the original simplicity of M-H, (ii) they are only applicable in situations where the data are independent and (iii) the computational cost of one iteration is stochastic which can potentially compromise any computational saving. Informed Sub-Sampling MCMC pushes the approximation one step forward: the computational cost of one iteration is deterministic and is controlled through a user specified parameter, the size of the subsamples.

The aforementioned methods rely on subsampling the whole dataset uniformly at random, at each iteration. Using such subsamples in lieu of the whole dataset in the original M-H chain leads to an algorithm whose approximation of $\pi$ comes with no guarantee when the subsamples size is fixed. Our main contribution is to show that assigning a distribution to the subsamples that reflects their fidelity to the whole dataset allows one to control the L1 approximation error, even when the subsample size is fixed, which is of computational interest.

The algorithm we propose to achieve this task, Informed Sub-Sampling MCMC, offers an alternative to situations where other scalable Metropolis-Hastings variants cannot be implemented (because the model does not satisfy the assumptions e.g. independence of the data, existence of a concentration inequality for the model or a cheap lower bound of the likelihood) or are inefficient (because the method ends up using nearly the whole dataset at each iteration). However, in scaling up Metropolis-Hastings there is no free lunch neither. In particular, our method replaces the uniform subsampling approach by a more sophisticated subsampling mechanism involving summary statistics. In this regard, even though our method is in principle widely applicable, it will only be useful in situations where a cheap summary statistics function satisfying Assumption 5.1 is available. In particular, we have shown that our method will give meaningful results when the maximum likelihood estimator is cheap to compute, which somewhat correlates with the optimality of summary statistics in ABC established in Fearnhead and Prangle (2012). We note that it is possible to construct a counterpart to ISS-MCMC in the Sequential Monte Carlo framework, as a sequence of scaled subposteriors naturally arises in our algorithm. We leave the design and evaluation of such a method for future research.

Similarly to the other noisy subsampling methods that approximate MH (e.g. Welling and Teh (2011), Korattikara et al. (2014), Bardenet et al. (2017), etc.), our theoretical results address the convergence in distribution of the ISS-MCMC chain to the posterior distribution. Questions of interest to practitioners encompass establishing Law of Large Numbers and Central limit theorems for those types of algorithms. Addressing those questions would allow a better understanding of approximate MCMC applied to large dataset contexts, from a practical perspective. In the noisy Monte Carlo literature, the closest result is perhaps Theorem 2.5 of Johndrow et al. (2015) which comes in the form of a non-asymptotic error bound in the L2 norm between $\pi f$ and its MCMC estimate. The assumptions are however quite strong which prevent applying it to ISS-MCMC and we will study those questions in a future work. Finally, recent developments in the understanding of approximate Markov chains have been carried out in Rudolf and Schweizer (2018). These authors provide explicit convergence bounds in the Wasserstein metric under the Wasserstein ergodicity assumption, a notion which is closely related to geometric ergodicity in V-norm, hence milder than uniform ergodicity. We regard this contribution as a very promising research avenue that unveils new ways of deriving quantitative error bounds for the practical problem of approximate Metropolis-Hastings algorithms that make use of subsets of data, under more realistic assumptions.

Appendix A Proofs

A.1 Proof of Proposition 1

Proof.

For notational simplicity and without loss of generality, we take here $g$ as the identity on $\Theta$. Let $n<N$ and $U$ be a subset of $\{1,\ldots,N\}$ with cardinal $n$. Consider the power likelihood:

[TABLE]

and the corresponding power posterior:

[TABLE]

where

[TABLE]

For any $\theta$ such that $p(\theta)\neq 0$, write:

[TABLE]

and the KL divergence between $\pi(\,\cdot\,|\,Y_{1:N})$ and $\tilde{\pi}(\,\cdot\,|\,Y_{U})$, denoted $\text{KL}_{n}(U)$, simply writes

[TABLE]

where $\Delta_{n}(U)=\sum_{k=1}^{N}S(Y_{k})-(N/n)\sum_{k\in U}S(Y_{k})$. Now, note that

[TABLE]

Plugging (45) into (44) yields:

[TABLE]

Finally, Cauchy-Schwartz inequality provides the following upper bound for $\text{KL}_{n}(U)$:

[TABLE]

∎

A.2 Proof of Proposition 2

Proof.

Under some weak assumptions, Bernstein-von Mises theorem states that $\pi(\,\cdot\,|\,Y_{1:N})$ is asymptotically (in $N$) a Gaussian distribution with the maximum likelihood $\theta^{\ast}$ as mean and $\Gamma_{N}=I^{-1}(\theta^{\ast})/N$ as covariance matrix, where $I(\theta)$ is the Fisher information matrix at $\theta$. Let us denote by $\Phi$ the pdf of $\pazocal{N}(\theta^{\ast},\Gamma_{N})$. Under this approximation, $\mathbb{E}_{\pi}(\theta)=\theta^{\ast}$ and from (45), we write:

[TABLE]

by integration of a multivariate Gaussian density function. Eventually, (48) yields the following approximation:

[TABLE]

∎

A.3 Proof of Proposition 3

Proof.

Let $\mathsf{U}_{n}\supset A_{n}(\theta):=\left\{U\in\mathsf{U}_{n},\;g(\theta)^{T}\Delta_{n}(U)\leq 0\right\}$ and remark that using Cauchy-Schwartz inequality, we have:

[TABLE]

Now, define $\bar{\Delta}_{n}(U):=\bar{S}(Y)-\bar{S}(Y_{U})$ where $\bar{S}$ is the normalized summary statistics vector, i.e if $U\in\mathsf{U}_{n}$, $\bar{S}(Y_{U})=S(Y_{U})/n$. Clearly, when $N\to\infty$, some terms

[TABLE]

will have a large contribution to the sum. More precisely, any mismatch between summary statistics of some subsamples $\{Y_{U},\,U\in\mathsf{U}_{n}\backslash A_{n}(\theta)\}$ with respect to the full dataset will be amplified by the factor $N$, whereby exponentially inflating the upper bound. However, assigning the distribution $\nu_{n,\epsilon}$ (12) to the subsamples $\{Y_{U},\,U\in\mathsf{U}_{n}\}$, allows to balance out this effect. Indeed, note that

[TABLE]

where $Z(\epsilon)=\sum_{U\in\mathsf{U}_{n}}\exp\{-\epsilon\|\Delta_{n}(U)\|^{2}\}$ and we have, for a fixed $n$ and when $N\to\infty$, that

[TABLE]

Since $g$ is bounded, then $\mathbb{E}\left\{{f(Y\,|\,\theta)}/{f(Y_{U}\,|\,\theta)^{N/n}}\right\}$ is bounded too. ∎

A.4 Proof of Proposition 4

We preface the proof Proposition 4 with five Lemmas, some of which are inspired from Medina-Aguayo et al. (2016). For notational simplicity, the dependence on $(n,\epsilon)$ of any ISS-MCMC related quantities is implicit. For all $(\theta,U)\in\Theta\times\mathsf{U}_{n}$, we denote by $\phi_{U}(\theta)=f(y_{U}\,|\,\theta)^{N/n}/f(y\,|\,\theta)$ and recall that $a(\theta,\theta^{\prime})$ is the (exact) MH acceptance ratio so that $\alpha(\theta,\theta^{\prime})=1\wedge a(\theta,\theta^{\prime})$. Unless stated otherwise, $\mathbb{E}$ is the expectation taken under $\nu_{n,\epsilon}$. For simplicity, $\tilde{K}_{n,\epsilon}$ is written as $\tilde{K}_{n}$.

Lemma 1**.**

For any $(\theta,\theta^{\prime})\in\Theta^{2}$, we have

[TABLE]

Proof.

This follows from a slight adaptation of Lemma 3.3 in Medina-Aguayo et al. (2016):

[TABLE]

where we have used Jensen’s inequality and the fact that the inequality $1\wedge ab\leq(1\wedge a)b$ holds for $a>0$ and $b\geq 1$. ∎

Lemma 2**.**

For any $\theta\in\Theta$ and all $\delta>0$, we have

[TABLE]

Proof.

The proof is identical to proof of Lemma 3.2 in Medina-Aguayo et al. (2016) by noting that Lemma 3.1 in the same reference holds for two random variables $\phi_{U}(\theta)$ and $\phi_{U}(\theta^{\prime})$ that are not independent, i.e for all $(\theta,\theta^{\prime})\in\Theta^{2}$ any $U\in\mathsf{U}_{n}$ and all $\delta\in(0,1)$

[TABLE]

∎

Lemma 3**.**

Assume that A.5.1 holds. Then we have

[TABLE]

Proof.

Using Cauchy-Schwartz inequality, we write that for all $(\theta,\theta^{\prime})\in\Theta^{2}$,

[TABLE]

Now for all $\theta\in\Theta$, we define the event $\mathcal{E}_{\theta}:=\{U\in\mathsf{U}_{n}\,,\;f(Y\,|\,\theta)\leq f(Y_{U}\,|\,\theta)^{N/n}\}$ so that

[TABLE]

and we note that for all $(\theta,U)\in\Theta\times\mathsf{U}_{n}$, Eq. (26) writes

[TABLE]

but also

[TABLE]

so that

[TABLE]

A similar argument gives the same upper bound for $\mathbb{E}\left\{{f(Y\,|\,\theta)}/{f(Y_{U}\,|\,\theta)^{N/n}}\right\}^{2}$ so that Eq. (50) yields

[TABLE]

The proof is completed by noting that for three numbers $a,b$ and $c$, $c>b\Rightarrow a\vee b\leq a\vee c$ and $\gamma\|\Delta_{n}(U)\|>0$. ∎

Lemma 4**.**

Assume that A.5.1 holds. Then we have for all $\theta\in\Theta$ and $\delta>0$

[TABLE]

Proof.

With the same notations as in proof of Lemma 3 and roughly with the same reasoning we have for all $\theta\in\Theta$ and all $\delta>0$

[TABLE]

where the first inequality follows by inclusion (on $\mathcal{E}_{\theta}$) of

[TABLE]

and similarly for the second term. Now, note that for all $x>0$, $\log(1+x)<-\log(1-x)$ so that

[TABLE]

where the last inequality follows from Markov inequality. ∎

We study the limiting case where $N$ is fixed and $n\to N$.

Lemma 5**.**

Assume $N$ is fixed and let $n\to N$. Then,

[TABLE]

Proof.

It follows from the fact that when $n\to N$, $\nu_{n,\epsilon}$ converges to the dirac on $U^{\dagger}=\{1,\ldots,N\}$ and therefore,

[TABLE]

∎

We can now prove Proposition 4:

Proposition

Assume that A.5.1 and A.5.1 hold. If the marginal MH chain $K$ is geometrically ergodic, i.e A.5.1 holds, then there exists an $n_{0}\leq N$ such that for all $n>n_{0}$, $\tilde{K}_{n}$ is also geometrically ergodic.

Proof.

By (Meyn and Tweedie, 2009, Theorems 14.0.1 & 15.0.1), there exists a function $V:\mathsf{X}\to[1,\infty[$, two constants $\lambda\in(0,1)$ and $b<\infty$ and a small set $S\subset\mathsf{X}$ such that $K$ satisfies a drift condition:

[TABLE]

We now show how to use the previous Lemmas to establish the geometric ergodicity of $\tilde{K}_{n}$ for some $n$ sufficiently large. This reasoning is very similar to that presented in (Medina-Aguayo et al., 2016, Theorem 3.2).

[TABLE]

Combining Eq. (51) with Eq. (52), we have that

[TABLE]

Fix $\epsilon>0$. From Lemma 5, there exists $(n_{1},n_{2})\in\mathbb{N}^{2}$ such that

[TABLE]

Combining Eqs. (52) and (54) yields that for all $n\geq n_{0}:=\max(n_{1},n_{2})$, we have

[TABLE]

Taking $\delta=\epsilon/2$ in Eq. (55) gives

[TABLE]

To show that $\tilde{K}_{n}$ (for $n>n_{0}$) satisfies a geometric drift condition, it is sufficient to take $\epsilon<(1-\lambda)/(1+\lambda)$ and to check that $S$ is also small for $\tilde{K}_{n}$. This is demonstrated exactly as in the proof of Medina-Aguayo et al. (2016, Theorem 3.2). ∎

A.5 Proof of Proposition 5

This proof borrows ideas from the perturbation analysis of uniformly ergodic Markov chains. First, note that by straightforward algebra we have that

[TABLE]

Now, under A.5.1 and using Mitrophanov (2005, Corollary 3.1) we have that for any starting point $\theta_{0}\in\Theta$,

[TABLE]

where $\lambda=\lceil\log(1/C)/\log\rho\rceil$. Combining Eqs (56) and (57) leads to Eq. (29) with $\kappa=\lambda+{C\rho^{\lambda}}/{1-\rho}$. Moreover, note that using Eq. (29) we have

[TABLE]

and taking the limit when $i\to\infty$ leads to Eq. (30). Finally, for a large enough $n$, we know from Proposition 4 that the marginal Markov chain $\{\tilde{\theta}_{i}\,,i\in\mathbb{N}\}$ produced by ISS-MCMC is geometrically ergodic and we denote by $\tilde{\pi}_{n}$ its stationary distibution. For such a $n$, we have for any $\theta_{0}\in\Theta$

[TABLE]

and taking the limit as $i\to\infty$ yields Eq. (31).

A.6 Extension of Proposition 5 beyond the time homogeneous case

We start with the two following remarks relative to the Informed Sub-Sampling Markov chain.

Remark 1**.**

Assume $U_{0}\sim\nu_{n,\epsilon}$ and $\tilde{\theta}_{0}\sim\mu$ for some initial distribution $\mu$ on $(\Theta,\vartheta)$. The distribution of $U_{i}$ given $\tilde{\theta}_{i}$ is for some $u\in\mathsf{U}_{n}$,

[TABLE]

where $\bar{K}(\theta,U;\mathrm{d}\theta^{\prime},U^{\prime}):=K(\theta,\mathrm{d}\theta^{\prime}\,|\,U)H(U,U^{\prime})$ and $H$ is the transition kernel of the Markov chain $\{U_{i},\,i\in\mathbb{N}\}$. As a consequence $\mathbb{P}(U_{i}\in\,\cdot\,|\,\tilde{\theta})$ depends on $\tilde{\theta}$ and $i$.

Remark 2**.**

The marginal Markov chain $\{\tilde{\theta}_{i},\,i\in\mathbb{N}\}$ produced by ISS-MCMC algorithm is time inhomogeneous since for all $A\in\mathcal{X}$,

[TABLE]

and $\mathbb{P}(U_{i}=u\,|\,\tilde{\theta}_{i})$ depends on $i$ (Remark 1). We thus denote by $\tilde{K}_{i}$ the marginal transition kernel $\tilde{\theta}_{i-1}\to\tilde{\theta}_{i}$. However, we observe that if the random variables $\{U_{i},\,i\in\mathbb{N}\}$ are i.i.d. with distribution $\nu_{n,\epsilon}$, $K_{i}$ becomes time homogeneous as $\mathbb{P}(U_{i}=u\,|\,\theta_{i})=\nu_{n,\epsilon}(u)$ for all $i$.

A consequence of Remark 2 is that Mitrophanov (2005, Theorem 3.1) does not hold when Assumption A.5.1 is not satisfied. Indeed, $\{\tilde{\theta}_{i},\,i\in\mathbb{N}\}$ is not a time homogeneous Markov chain in this case and we first need to generalize the result from Mitrophanov in order to apply it to our context. This is presented in Lemma 6.

Lemma 6**.**

Let $K$ be the transition kernel of an uniformly ergodic Markov chain that admits $\pi$ as stationary distribution. Let $\tilde{K}_{i}$ be the $i$-th transition kernel of the ISS-MCMC Markov chain. In particular, let $p_{i}(\,\cdot\,|\,\theta):=\mathbb{P}(U_{i}\in\,\cdot\,|\,\theta)$ be the distribution of the random variable $U_{i}$, used at iteration $i$ of the noisy Markov chain given $\theta$. We have:

[TABLE]

where $\delta_{i}:\Theta\times\Theta\to\mathbb{R}^{+}$ is a function that satisfies

[TABLE]

and the expectation is under $p_{i}(\,\cdot\,|\,\theta)$.

Proof.

In addition of the notations of Section 4, we define the following quantities for a Markov transition kernel regarded as an operator on $\mathcal{M}$, the space of signed measures on $(\Theta,\mathcal{B}(\Theta))$: $\tau(K):=\sup_{\pi\in\mathcal{M}_{0,1}}\|\pi K\|$ is the ergodicity coefficient of $K$, $\|K\|:=\sup_{\pi\in\mathcal{M}_{1}}\|\pi K\|$ is the operator norm of $K$ and $\mathcal{M}_{1}:=\{\pi\in\mathcal{M},\,\|\pi\|=1\}$ and $\mathcal{M}_{0,1}:=\{\pi\in\mathcal{M}_{1},\,\pi(\Theta)=0\}$.

Remarks 1 and 2 explain why, in general, $\{\tilde{\theta}_{i},\,i\in\mathbb{N}\}$ is a time-inhomogeneous Markov chain with transition kernel $\{\tilde{K}_{i},\,i\in\mathbb{N}\}$. For each $i\in\mathbb{N}$, define $\pi_{i}$ as the distribution of $\theta_{i}$ produced by the Metropolis-Hastings algorithm (Alg. 1) with transition kernel $K$, referred to as the exact kernel hereafter. Our proof is based on the following identity:

[TABLE]

for each $i\in\mathbb{N}$. Equation (60) will help translating the proof of Theorem 3.1 in Mitrophanov (2005) to the time-inhomogeneous setting and in particular, we have for each $i\in\mathbb{N}$:

[TABLE]

Following the proof of Theorem 3.1 in Mitrophanov (2005), we obtain

[TABLE]

where $\lambda=\left\lceil{\log_{\rho}(1/C)}\right\rceil$. Without loss of generality, we take $\pi_{0}=\tilde{\pi}_{0}$ and since $\|\pi-\tilde{\pi}_{i}\|\leq\|\pi-\pi_{i}\|+\|\pi_{i}-\tilde{\pi}_{i}\|$ we have for all $i>\lambda$ that

[TABLE]

Taking the limit as $i\to\infty$ leads to

[TABLE]

Using a similar derivation than in the proof of Corollary 2.3 in Alquier et al. (2016), we obtain

[TABLE]

where the expectation is under $p_{i}(\,U\,|\,\theta)$ and which combined with (64) leads to

[TABLE]

where the expectation is under $Q(\theta,\cdot)\otimes p_{i}(\,\cdot\,|\,\theta)$. Any upper bound $\delta_{i}(\theta,\theta^{\prime})$ of the expectation on the right hand side yields (59). ∎

By straightforward algebra, we have:

[TABLE]

where we have defined $\phi_{U}(\theta)=f(Y_{U}\,|\,\theta)^{N/n}/f(Y\,|\,\theta)$. Using Lemma 6, we have that

[TABLE]

which is the counterpart of (30) when Assumption A.5.1 does not hold. We note that the second supremum in Eq. (66) is in fact $B_{n}$ defined at Eq. (28) and, as such, can be controlled as described in Section 5.3.1. However, this is not clearly the case for the first supremum in Eq. (66) which differs from $A_{n}$ defined at Eq. (27):

[TABLE]

We now show that, under two additional assumptions (A.6 and A.6), the control based on the summary statistics also applies to the time inhomogeneous case when Assumption A.5.1 does not hold.

A 5. One-step minorization

For all $i\in\mathbb{N}$ and all $A\in\vartheta$, there exists some $\eta>0$ such that $p_{i}(A)>\eta\lambda(A)$ where $\lambda$ is the Lebesgue measure.

This assumption typically holds if $\Theta$ is compact or if the chain $\{\tilde{\theta}_{i},U_{i}\}_{i}$ admits a minorization condition. Since we assume, in this discussion, that the exact MH Markov chain is uniformly ergodic and as such satisfy a minorization condition, see e.g. Meyn and Tweedie (2009, Thm 16.2.3) and Hobert and Robert (2004). We may study conditions on which $\{\tilde{\theta}_{i}\}_{i}$ inherits this property and leave this for future work but already note that Assumption A.6 is not totally unrealistic.

**A 6. ** The marginal Markov chain $\{U_{i}\}_{i}$ has initial distribution $U_{0}\sim\nu_{n,\epsilon}$.

Even though this assumption is difficult to meet in practice as $|\mathsf{U}_{n}|$ may be very large, the discussion at the beginning of Section 6.1 indicates an approach to set the distribution of $U_{0}$ close from $\nu_{n,\epsilon}$.

Again, while the Assumptions 6 and 6 are perhaps challenging to guarantee, Proposition 8 aims at giving some level of confidence to the user that the ISS-MCMC method is useful, even when Assumption A.5.1 does not hold. In addition, it reinforces the importance of choosing summary statistics that satisfy Assumption A.5.1.

Proposition 8**.**

Assume that Assumptions A.5.1, A.5.1, A.6 and A.6 hold. Then there exists a positive number $M>0$ such that

[TABLE]

where $A_{n}$ and $\tilde{A}_{n}$ have been defined at Eq. (67).

Corollary 2**.**

Under the same Assumptions as Proposition 8, the control explained in Section 5.3.2 is also valid in the time inhomogeneous case.

Proof of Proposition 8.

From Assumption A4, there exists some $\gamma>0$ such that

[TABLE]

where $\mathrm{d}p_{i}(U\,|\,\theta)=p_{i}(U\,|\,\theta)\mathrm{d}U$ and $\mathrm{d}U$ is the counting measure. Now, the conditional probability writes:

[TABLE]

On the one hand, Lemma 1 shows that there exists a bounded function $f_{i}$ such that $\mathbb{P}(U_{i}\in\cdot\,,\,\tilde{\theta}_{i}\in\mathrm{d}\theta)\leq f_{i}(\tilde{\theta})\mathrm{d}\theta\nu_{n,\epsilon}(\,\cdot\,)$. On the other hand, Assumption 6 guarantees that there exists some $\eta>0$ such that for all $\tilde{\theta}\in\Theta$, $\mathbb{P}(\tilde{\theta}_{i}\in\mathrm{d}\theta)>\eta\mathrm{d}\theta$. Combining those two facts allows to write that

[TABLE]

Plugging Eq. (70) into Eq. (69), yields to

[TABLE]

which completes the proof, setting $M:=\sup_{\theta}\sup_{i}f_{i}(\theta)/\eta$. ∎

Lemma 1**.**

Assume that Assumptions A.5.1, A.5.1, A.6 and A.6 hold. In addition, let us assume that $U_{0}\sim\nu_{n,\epsilon}$. Then $p_{i}(\theta,U)$ is dominated by $\mathrm{d}\theta\mathrm{d}U$ where $\mathrm{d}\theta$ and $\mathrm{d}U$ implicitly refer to the Lebesgue and the counting measure, respectively. In other words there is a sequence of bounded functions $\{f_{i}:\Theta\to\mathbb{R}^{+}\}$ such that

[TABLE]

Proof.

We proceed by induction. Defining $\varrho(\tilde{\theta}\,|\,U)$ as the probability to reject a MH move for the parameter $\tilde{\theta}$ when the subset variable is $U$, we recall that $\varrho(\tilde{\theta}\,|\,U)<1$ and $\tilde{\alpha}(\tilde{\theta},\tilde{\theta}^{\prime}\,|\,U)<1$. By assumption on the proposal kernel, it satisfies $Q(\tilde{\theta},\,\mathrm{d}\tilde{\theta}^{\prime}\,)=Q(\tilde{\theta},\tilde{\theta}^{\prime})\mathrm{d}\tilde{\theta}^{\prime}$ and define the function $\overline{Q}:\theta\mapsto\sup_{\tilde{\theta}^{\prime}\in\Theta}Q(\tilde{\theta}^{\prime},\theta)$. Similarly, we define the function $\overline{\varrho}:\theta\mapsto\sup_{U\in\mathsf{U}_{n}}\varrho(\theta\,|\,U)$. Deriving the calculation separately for the continuous and the diagonal parts of the Metropolis-Hastings kernel $K(\theta,\cdot\,|\,U)$ (see Eq. (25)), we have:

[TABLE]

where the last equality follows from the $\nu_{n,\epsilon}$-stationarity of $H$. In this derivation, we have defined $\mu$ as the initial distribution of the Markov chain $\{\tilde{\theta}_{i}\}_{i}$ and $\nu$ as a shorthand notation for $\nu_{n,\epsilon}$. Now, let us assume that there is a bounded function $f_{i-1}$ such that $\mathrm{d}p_{1}(U,\tilde{\theta})\leq f_{i-1}(\tilde{\theta})\mathrm{d}\tilde{\theta}\mathrm{d}\nu(U)$. Using the notation $\mu K:=\int\mu(\mathrm{d}x)K(x,\cdot)$ for any Markov kernel $K$ and a measure $\mu$ on some measurable space $(\mathsf{X},\mathcal{X})$ and recalling that $\bar{K}$ is the transition kernel of ISS-MCMC on the extended space $\Theta\times\mathsf{U}_{n}$, we have:

[TABLE]

and $f_{i}$ is bounded. The first term in the third inequality follows from noting that

[TABLE]

∎

A.7 Proof of Proposition 6

Proof.

Note that for all $(\theta,\zeta)\in\Theta\times\mathbb{R}^{d}$, a Taylor expansion of $\pi(\theta)$ and $\phi_{U}(\theta)$ at $\theta+\Sigma\zeta$ in (32) combined to the triangle inequality leads to:

[TABLE]

where the expectation is under $\Phi_{d}$ and $R(x)=o(x)$ at [math]. Applying Cauchy-Schwartz gives:

[TABLE]

Now, we observe that:

$\mathbb{E}\{\|M\zeta\|\}=\mathbb{E}\{\sum_{i=1}^{d}(\sum_{j=1}^{d}M_{i,j}\zeta_{j})^{2}\}^{1/2}\leq\mathbb{E}\{\sum_{i=1}^{d}|\sum_{j=1}^{d}M_{i,j}\zeta_{j}|\}\leq\mathbb{E}\{\sum_{i=1}^{d}\sum_{j=1}^{d}|M_{i,j}||\zeta_{j}|\}=\sum_{i=1}^{d}\sum_{j=1}^{d}|M_{i,j}|\mathbb{E}\{|\zeta_{i}|\}=\sqrt{\frac{2}{\pi}}\|M\|_{1}$

• 2.

$\mathbb{E}\{\|M\zeta\|^{2}\}=\mathbb{E}\{\sum_{i=1}^{d}(\sum_{j=1}^{d}M_{i,j}\zeta_{j})^{2}\}=\sum_{i=1}^{d}\mathbb{E}\{(\sum_{j=1}^{d}M_{i,j}\zeta_{j})^{2}\}=\sum_{i=1}^{d}\text{var}(\sum_{j=1}^{d}M_{i,j}\zeta_{j})=\sum_{i=1}^{d}\sum_{j=1}^{d}M_{i,j}^{2}\text{var}(\zeta_{j})=\|M\|_{2}^{2}$

• 3.

considering the quadratic form associate to the operator $T(U,\theta)=M^{T}\nabla_{\theta}^{2}\phi_{U}(\theta)M$, noting that $T(U,\theta)$ is symmetric its eigenvalues $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{d}$ are real and we have

[TABLE]

so that:

[TABLE]

where for any square matrix $A$, we have defined $|\!|\!|A|\!|\!|=\sup_{x\in\mathbb{R}^{d},\|x\|=1}\|Ax\|$ as the operator norm.

∎

Appendix B Proof of Proposition 7

In this section, we are assuming that there is an infinite stream of observations $(Y_{1},Y_{2},\ldots)$ and a parameter $\theta_{0}\in\Theta$ such that $Y_{i}\sim f(\,\cdot\,|\,\theta_{0})$. Let $\rho>1$ be a constant defined as the ratio $N/n$ i.e the size of the full dataset over the size of the subsamples of interest. The full dataset is thus $Y_{1:\rho n}$. We define the set

[TABLE]

such that $Y_{U}$ ($U\in\mathsf{U}_{n}^{\rho}$) is the set of subsamples of interest. We study the asymptotics when $n\to\infty$ i.e we let the whole dataset and the size of subsamples of interest grow at the same rate.

Proposition 7**.**

Let $\theta^{\ast}_{\rho n}$ be the MLE of $Y_{1},\ldots,Y_{\rho n}$ and $\theta^{\ast}_{U}$ be the MLE of the subsample $Y_{U}$ ($U\in\mathsf{U}_{n}^{\rho}$). Assume that there exists a compact set $\kappa_{n}\subset\Theta$ such that $(\theta^{\ast}_{\rho n},\theta_{0})\in\kappa_{n}^{2}$ and for all $U$, there exists a compact set $\kappa_{U}\subset\Theta$ such that $(\theta^{\ast}_{U},\theta_{0})\in\kappa_{U}^{2}$. Then, there exists a constants $\beta$, a metric $\|\cdot\|_{\theta_{0}}$ on $\Theta$ and a non-decreasing subsequence $\{\sigma_{n}\}_{n\in\mathbb{N}}$, ($\sigma_{n}\in\mathbb{N}$) such that for all $U\in\mathsf{U}_{\sigma_{n}}^{\rho}$, we have for $p$-almost all $\theta\in\kappa_{n}\cap\kappa_{U}$

[TABLE]

where

[TABLE]

Proof.

Fix $n\in\mathbb{N}$. Consider the case where the prior distribution $p$ is uniform on $\kappa_{n}$. In this case, the posterior is

[TABLE]

and from corollary 3, we know that there exists a subsequence $\tau_{n}\subset\mathbb{N}$ such that for $p$-almost all $\theta\in\kappa_{n}$

[TABLE]

where $\theta\mapsto\Phi_{\rho\tau_{n}}(\theta)$ is the pdf of $\pazocal{N}(\theta^{\ast}_{\rho\tau_{n}},I(\theta_{0})^{-1}/\rho\tau_{n})$. Similarly, there exists another subsequence $\gamma_{n}\subset\mathbb{N}$ such that for all $U\in\mathsf{U}_{\gamma_{n}}^{\rho}$ and for $p$-almost all $\theta\in\kappa_{U}$

[TABLE]

where $\theta\mapsto\Phi_{U}(\theta)$ is the pdf of $\pazocal{N}(\theta^{\ast}_{U},I(\theta_{0})^{-1}/|U|)$. Let $\{\sigma_{n}\}_{n\in\mathbb{N}}$ be the sequence defined as $\sigma_{n}=\max\{\tau_{n},\gamma_{n}\}$. We know from (73) and (74) that for all $\varepsilon>0$ and all $\eta>0$, there exists $n_{1}\in\mathbb{N}$ such that for all $U\in\mathsf{U}_{\sigma_{n}}^{\rho}$ and for all $n\geq n_{1}$

[TABLE]

Now, by straightforward algebra, we have for any $U\in\mathsf{U}_{\sigma_{n}}^{\rho}$

[TABLE]

where we have used Lemma 3 for the first inequality and the triangle inequalities for the second. Combining (76) with (75) yields (34). ∎

Lemma 2**.**

Consider a posterior distribution $\pi_{n}$ given $n$ data $Y_{1:n}$ where $p$ is the prior distribution and its Bernstein-von Mises approximation is $\Phi_{n}=\pazocal{N}(\theta^{\ast}(Y_{1:n}),I(\theta_{0})^{-1}/n)$. There exists a subsequence $\{\tau_{n}\}_{n}\subset\mathbb{N}$ such that

[TABLE]

Proof.

This follows for the fact that convergence in $L_{1}$ implies pointwise convergence almost everywhere of a subsequence, i.e there exists a subsequence $\{\tau_{n}\}_{n\in\mathbb{N}}\subset\mathbb{N}$ such that

[TABLE]

Eq. 77 follows from combining the Bernstein-von Mises theorem and Eq. (78):

[TABLE]

∎

Corollary 3**.**

There exists a subsequence $\{\tau_{n}\}_{n\in\mathbb{N}}\in\mathbb{N}$ such that

[TABLE]

Proof.

Follows from Lemma 2, by continuity of the logarithm. ∎

Lemma 3**.**

For any $U\in\mathsf{U}_{n}$, let $\theta\mapsto\Phi_{U}(\theta)$ be the pdf of $\pazocal{N}(\theta_{U}^{\ast},I(\theta_{0})^{-1}/n)$ and $\Phi_{\rho n}$ be the pdf of $\pazocal{N}(\theta_{\rho n}^{\ast},I(\theta_{0})^{-1}/\rho n)$ be the Bernstein-von Mises approximations of respectively $\pi(\,\cdot\,|\,Y_{U})$ and $\pi(\cdot\,|\,Y_{1:\rho n})$ where $U\subset\mathsf{U}_{n}(Y_{1:\rho n})$. Then we have for all $\theta\in\Theta$

[TABLE]

where for any $d$-squared symmetric matrix $M$, we have defined by $\|\cdot\|_{M}$ the norm associated to the scalar product $\left\langle u,v\right\rangle_{M}=u^{T}Mv$.

Proof.

This follows from straightforward algebra and noting that

[TABLE]

∎

Acknowledgements

The Insight Centre for Data Analytics is supported by Science Foundation Ireland under Grant Number SFI/12/RC/2289. Nial Friel’s research was also supported by an Science Foundation Ireland grant: 12/IP/1424. Pierre Alquier’s research was funded by Labex ECODEC (ANR - 11-LABEX-0047) and by the research programme New Challenges for New Data from LCL and GENES, hosted by the Fondation du Risque.

We thank the Associate Editor and two anonymous Referees for their contribution to this work.