\section{Introduction} \label{s:intro}
Advances in methodology for the collection and storage of data have led to scientific challenges in a wide array of disciplines. This is particularly the case in Statistics as the complexity of appropriate statistical models and those we wish to fit often increases with data size. Many current state-of-the-art statistical methodologies have algorithmic cost that scale poorly with increasing volumes of data. As noted by \cite{Jordan:2013}, `many statistical procedures either have unknown runtimes or runtimes that render the procedure unusable on large-scale data' and has resulted in a proliferation in the literature of methods `\ldots{}which may provide no statistical guarantees and which in fact may have poor or even disastrous statistical properties'.\\
\\
This is particularly keenly felt in computational and Bayesian statistics, in which the standard computational tools are Markov chain Monte Carlo (MCMC), Sequential Monte Carlo (SMC) and their many variants (see for example \cite{robert2013monte}). MCMC methods are {\em exact} in the (weak) sense that they construct Markov chains which have the correct limiting distribution. Although MCMC methodology has had considerable success in being applied to a wide variety of substantive areas, they are not well-suited to this new era of `big data' for the following reasons:
\begin{enumerate}
\item[\textbullet]
Most MCMC algorithms require at least one likelihood evaluation per iteration. For example, Metropolis-Hastings algorithms require likelihood evaluations in order to perform the {\em accept/reject} step. Given a data set of size $n$, this will be at best an $\mathcal{O}(n)$ calculation. For algorithms which avoid accept/reject steps, such as Gibbs samplers, alternative $\mathcal{O}(n)$ calculations such as the evaluation of suitable sufficient statistics are required.
\item[\textbullet]
A common approach in MCMC for missing data and latent variable models is data augmentation \cite{tannerwong87}. In the context of MCMC for large data sets this approach increases the dimensionality of the state space, and in some cases this increase will be $\mathcal{O}(n)$. As we demonstrate in \secref{s:examp:2} we are able to circumvent this issue.
\end{enumerate}
\noindent Naturally, owing to the previous successes of MCMC, there has been considerable effort to modify existing methodology in order to address computational scalability. The success of such methods in reducing algorithmic cost has been mixed, and has invariably led to a compromise on exactness --- such methodologies generally construct a stochastic process with limiting distribution which is (at least hopefully) close to the desired target distribution. Broadly speaking these methods can be divided into three classes of approach: `Divide-and-conquer' methods; `Exact Subsampling' methods; and, `Approximate Subsampling' methods. Each of these approaches has its own strengths and weaknesses, and we will briefly review these in the following paragraphs.
\\ \\
Divide-and-conquer methods (for instance, \cite{Neiswanger:2013,Wang:2013,Scott:2013,Minsker:2014}) begin by splitting the data set into a large number of smaller data sets (which may or may not overlap). Inference is then conducted on these smaller data sets and resulting estimates are somehow combined. A clear advantage of such an approach is that inference on each small data set can be conducted independently, and in parallel, and so if one had access to a large cluster of computing cores then computational cost would be significantly reduced. Interaction between the cores is typically impracticable as the primary bottleneck is not computational speed, but latency. Therefore, from a parallelisation perspective this approach is perhaps the most appealing. The primary weakness of these methods is that the recombination of the separately conducted inferences is inexact. All current theory is asymptotic in the number of data points, $n$ \cite{Neiswanger:2013,Li:2016}. For these asymptotic regimes the posterior will tend to a Gaussian distribution \cite{Johnson:1970}, and it is questionable whether divide-and-conquer methods offer an advantage over simple approaches such as a Laplace approximation to the posterior \cite{Bardenet:2015}. Most results on convergence rates (e.g. \cite{Srivastava:2015}) have rates that are of the $\mathcal{O}(m^{-1/2})$, where $m$ is the number of data-points in each sub-set. As such they are no stronger than convergence rates for analysing just a single batch. One exception is in \cite{Li:2016}, where convergence rates of $\mathcal{O}(n^{-1/2})$ are obtained, albeit under strong conditions. However these results only relate to estimating marginal posterior distributions, rather than the full posterior. 
\\ \\ 
Subsampling methods are designed so that each iteration requires access to only a subset of the data. Exact approaches in this vein typically require subsets of the data of random size at each iteration. One approach is to construct unbiased estimators of point-wise evaluations of the target density using subsets of the data, which could then be embedded within the pseudo-marginal MCMC framework recently developed by \cite{Andrieu/Roberts:2009}. Unfortunately, the construction of such positive unbiased estimators is not possible in general \cite{Jacob:2013}. An approach to MCMC in which the acceptance of a proposed move is determined by using a subset of the data is proposed by \cite{Korattikara:2013}, but they choose the size of subset such thatwith high probability the correct decision is made.  However, it was shown in \cite{Bardenet:2014} that such approaches, although cheaper, require subsets of the data of size $\mathcal{O}(n)$, and thus do not even reduce the computational complexity of standard MCMC. 
\\ \\
Approximate subsampling methods have the advantage that they typically require fixed subsets of the data at every iteration, and are based on direct (although often difficult to justify) modifications of exact methods which themselves benefit from sound theoretical underpinning. In particular, one direction which has garnered considerable interest is to exploit the fact that the Langevin diffusion has known invariant distribution. If we denote our target density by $\pi$ then the Langevin diffusion 
\begin{equation} \label{eq:LangevinSDE}
\,\mathrm{d} \mathbf{X}_t=\frac{1}{2}\mathbf{\Lambda} \nabla \log \pi(\mathbf{X}_t) \,\mathrm{d} t + \mathbf{\Lambda}^{1/2}  \,\mathrm{d} \mathbf{W}_t, \quad\quad \mathbf{X}_0 = \mathbf{x} \in\mathbbm{R}^d,\,\,\,t\in[0,\infty),
\end{equation}
in which $\mathbf{W}_t$ is $d$-dimensional Brownian motion and $\mathbf{\Lambda}$ is a prespecified preconditioning matrix, has $\pi$ as its invariant distribution. So one could simply simulate this diffusion as a means to obtain Monte Carlo draws from $\pi$. Unfortunately, exact simulation from the Langevin diffusion is not possible for general $\pi$. Were it possible, even in the simplest Bayesian setting $\pi(\mathbf{x}) \propto \prod^n_{i=0} f_i(\mathbf{x})$ (where $f_0(\mathbf{x})$ is a prior and $\{f_i(\mathbf{x})\}^n_{i=1}$ are likelihood factors), evaluation of the drift of the diffusion would be an $\mathcal{O}(n)$ calculation.\\
\\ 
One natural approach to big data problems suggested by \cite{ICML:WT11} is to employ a stochastic gradient Euler-Maruyama discretisation scheme. As in a regular Euler-Maruyama discretisation scheme they suggest that the sample path is simulated iteratively over small time intervals of length $h>0$, but where the drift component of the diffusion dynamics is estimated from  a subset of the available data. So for a single iteration of their algorithm one would simulate
\begin{align}
\mathbf{X}_{t+h}|(\mathbf{X}_{t}=\mathbf{x}) \sim \text{N}\left(\mathbf{x}+h \tilde{\beta}(\mathbf{x}), h\mathbf{\Lambda}^{1/2}\right), \label{eq:stocgradprop}
\end{align}
where $\tilde{\beta}$ is an unbiased estimator of the drift of (\ref{eq:LangevinSDE}) calculated at $\mathbf{X}_t=\mathbf{x}$ using a random subset of size $m\ll n$ of the available data. One would hope that the approximation from using this sub-sampling estimate of the drift would be small if the variance of the estimator is small compared to the variance of the driving Brownian motion over a
time-step, $h\mathbf{\Lambda}^{1/2}$. The attraction of this approach is that the construction of such a $\tilde{\beta}$ is trivial, requiring only the unbiased estimation of a sum. In particular, we might draw a random subset $A$ of size $m$ without replacement from
$\{0,\ldots,n\}$,  setting 
\begin{align}
\tilde{\beta} (\mathbf{x}) = \dfrac{n+1}{2m}
\mathbf{\Lambda}\sum_{j\in A} \nabla \log f_j(\mathbf{x})\ .
\label{eq:stochgrad}
\end{align} 
Such `stochastic gradient' approximations are not exact: even if we were to use the entire data set at every iteration such an approach would, due to the discretisation error, not recover the correct target distribution \cite{Bescom}. A comprehensive theoretical study of the conditions and (fixed) step sizes required to recover a distribution which is close to the target distribution is provided by \cite{arxiv:dm15}. The most obvious way then to ensure the correct target distribution is obtained, is to embed the discretised dynamics as a proposal within a Markov chain Monte Carlo framework, which leads to the class of Metropolis-adjusted Langevin  algorithms outlined in \cite{Bescom,B:RT96}. However, due to the $\mathcal{O}(n)$ accept/reject step this MCMC embedding renders the approach impractical for large data sets.  \\
\\
Despite these theoretical concerns, using either (\ref{eq:stocgradprop})--(\ref{eq:stochgrad}) or some related approach \cite{Chen:2014,Ma:2015} has been shown to work well in practice. They can give accurate results and lead to speed-ups when compared to standard MCMC. However, for the approximation error from sub-sampling to be small we often require sub-sample sizes that are $\mathcal{O}(n)$. Furthermore, it can be hard to quantify the bias introduced by the both the subsampling and the time-discretisation of the Langevin dynamics. Some effort has been made to rectify the original na\"{i}ve stochastic gradient scheme by using an increasingly fine discretisation with iteration \cite{Teh:2014,Vollmer:2015}. Although \cite{arxiv:dm15} show that this approach recovers (asymptotically) the correct target distribution, the trade-off is a computational cost which increases nonlinearly with diffusion time, and limited model class applicability (at least with any theoretical guarantees).
\\ \\
The approach to the problem of big data proposed below is a significant departure from the current literature. Rather than building our methodology upon the stationarity of appropriately constructed Markov chains, we develop a novel approach based on the {\em quasi-limiting} distribution of suitably constructed stochastically weighted diffusion processes. A {\em quasi-stationary distribution} for a Markov process $X$ with respect to a Markov stopping time $\tau $ is the limit of the distribution of $X_t\mid \tau>t$ as $t \to \infty $ \cite{BK:QSD}, and is completely unrelated to the popular area of Quasi-Monte Carlo.  These {\em Quasi-Stationary Monte Carlo (QSMC) methods} which we develop can be used for a broad range of Bayesian problems (of a similar type to MCMC) and exhibit interesting and differing algorithmic properties. There are a number of different possible implementations of the theory which open interesting avenues for future research, in terms of  branching processes, by means of stochastic approximation methods, or (as outlined in this paper) SMC methods. One particularly interesting difference between our class of Monte Carlo algorithms and MCMC is that QSMC methods allow us to circumvent entirely the Metropolis-Hastings type accept/reject steps, while still retaining theoretical guarantees that the correct limiting target distribution is recovered. In the case of big data problems, this removes one of the fundamental $\mathcal{O}(n)$ bottlenecks in computation.
\\\\
Our Quasi-Stationary Monte Carlo methods can be applied in big data contexts by using a novel subsampling approach. 
The resulting \textit{Scalable Langevin Exact Algorithm (ScaLE)} not only has the correct limiting target distribution (and so in this sense is {\em exact}) but also has {\em sub-linear} iterative cost as a function of data size. As we show in \secref{s:scale}, our approach to subsampling decreases the computational complexity of QSMC despite requiring subsampled data sets of fixed $\mathcal{O}(1)$ size. A key component of our approach is the use of control variates to reduce the variability of sub-sampling estimators of features of the posterior, but in addition exploits properties and simulation methods for Brownian motion. Although the use of control variates for MCMC in the context of big data and sub-sampling have been previously considered \cite{Bardenet:2015,Huggins:2016}, this is the first use of control variates within an exact method. We note that closely related control variate methods are
being developed in concurrent work on MCMC (see \cite{arxiv:bfr16}) and it seems the methodology may have quite wide applicability within the field of likelihood inference for big data.
\\\\
The remainder of this paper is structured as follows. In the next section we outline the theory underpinning the methodology we develop, with Sections \ref{s:dmc} and \ref{s:dmc:trad} introducing diffusions and the theory for diffusions typically exploited in Monte Carlo methodology, while \secref{s:ergodicity} outlines the construction of our quasi-stationary process and the theory by which one can parametrise this process in order for it to recover a probability distribution of interest. \secref{s:emcfd} is focussed on developing Monte Carlo methodology based on the theory in \secref{s:dmc}. In particular, \secref{s:nsea} outlines methodology for the exact simulation of diffusions based on path-space rejection sampling and their application to the simulation of our quasi-stationary process (without normalisation), along with an importance sampling variant. \secref{s:emcfd:smc} provides a theoretical underpinning for embedding the importance sampling variant within a Sequential Monte Carlo (SMC) framework, along with algorithmic detail for our Quasi-Stationary Monte Carlo (QSMC) algorithm. In \secref{s:scale} we extend our QSMC method to deal with the particular problem of `big data' and present the Scalable Langevin Exact Algorithm (ScaLE) which is an application of these methods to QSMC. We conclude in \secref{s:examp} by providing a number of examples to which we have applied our methodology. Note that for clarity of presentation we have suppressed much of the technical and algorithmic detail, but this can be found in the appendices: \apxref{apx:polynomial} examines the conditions on the theory of our quasi-stationarity of killed Brownian motion by considering heavy tailed densities; \apxref{s:conditions} provides conditions under which the path-space rejection sampling construction of  \secref{s:nsea} can be employed; \apxref{s:pathspace layer} details from basic random number generation the mechanism by which we simulate our process; \apxrefpl{apx:prs}{apx:kbm} provide further technical and algorithmic detail on the path-space rejection sampler and killed Brownian motion algorithm respectively, which are introduced and discussed in  \secref{s:nsea}; \apxref{apx:rqsmc} provides detail of the rejection sampling variant of \secref{s:nsea} within the SMC construction of \secref{s:emcfd:smc}; in \apxref{apx:r-scale} we provide details of the sub-sampling embedding within a rejection sampling framework; and finally, \apxref{app:smc} provides more theoretical rigour to the application of SMC methodology we use.



\section{Diffusions for Monte Carlo} \label{s:dmc}
A ($d$-dimensional) diffusion $\mathbf{X}:\,\mathbbm{R}\!\to\!\mathbbm{R}^d$ over an infinite-time horizon is a Markov process with almost surely continuous sample paths, which can be defined as the solution to the following stochastic differential equation (SDE),
\begin{align}
\,\mathrm{d} \mathbf{X}_t 
& = \mathbf{\beta}(\mathbf{X}_t) \,\mathrm{d} t + \mathbf{\sigma}(\mathbf{X}_t) \,\mathrm{d} \mathbf{W}_t \quad\quad \mathbf{X}_0 = \mathbf{x} \in\mathbbm{R}^d,\,\,\,t\in[0,\infty), \label{eq:gendiffusion}
\end{align}
where $\mathbf{\beta}:\,\mathbbm{R}^d\!\to\!\mathbbm{R}^d$ denotes the (instantaneous) drift vector, $\mathbf{\sigma}:\,\mathbbm{R}^{d}\!\to\!\mathbbm{R}^{d\times d}_+$ denotes the (instantaneous) volatility matrix and $\mathbf{W}_t$ is a standard ($d$-dimensional) Brownian motion. Regularity conditions are assumed to hold to ensure the existence of a unique non-explosive weak solution (see for instance \cite{BK:SDE}).\\
\\
As discussed in \secref{s:intro}, we are interested in exploiting the known theoretical properties of diffusions in order to sample from some distribution of interest. The Fokker--Planck partial differential equation describes the law of the diffusion (\ref{eq:gendiffusion}). In particular, the evolution of the probability density for $\mathbf{X}_t$, denoted by $\nu(\mathbf{x},t)$, can be described as follows,
\begin{align}
\dfrac{\partial}{\partial t}\nu(\mathbf{x},t) 
& = -\sum_i \dfrac{\partial}{\partial \mathbf{x}_i}\left[\mathbf{\beta}_i(\mathbf{x})\nu({\mathbf{x},t})\right] +\dfrac{1}{2}\sum_{i,j}\dfrac{\partial^2}{\partial \mathbf{x}_i \partial \mathbf{x}_j}\left[\left[\mathbf{\sigma}(\mathbf{x})\mathbf{\sigma}^T(\mathbf{x})\right]_{i,j}\nu(\mathbf{x},t)\right]. \label{eq:fokker-planck}
\end{align}
We are interested in choosing the dynamics of the diffusion such that it has an invariant distribution that coincides with a distribution of interest. One can then, in principle, construct a scheme for generating samples suitable for approximating expectations with respect to that distribution. This is done by simply simulating a (finite horizon) sample trajectory from the measure induced on path-space by (\ref{eq:gendiffusion}). Mirroring the use of stationary discrete time processes in an MCMC context, such an approach is motivated by the elementary fact that if $\mathbf{X}_t\sim\nu$, where $\nu$ is the invariant distribution, then $\mathbf{X}_{t+s}|(\mathbf{X}_t=\mathbf{x}_t) \sim \nu$ for every $s>0$. As we are interested in the simulation of such a class of diffusions without any form of error, we further restrict our attention to the class of diffusions in which the volatility matrix is constant (i.e., $\mathbf{\sigma}(\mathbf{x})\equiv\mathbf{\Lambda}^{1/2}$). Setting $\nu(\mathbf{x},t)=\nu(\mathbf{x})$ in the Fokker-Planck equation it is straightforward to show that the following `Langevin' diffusion has $\nu(\mathbf{x})$ as its invariant distribution,
\begin{align}
\,\mathrm{d} \mathbf{X}_t 
& = \dfrac{1}{2}\mathbf{\Lambda} \nabla\log\nu(\mathbf{X}_t) \,\mathrm{d} t + \mathbf{\Lambda}^{1/2} \,\mathrm{d} \mathbf{W}_t, \quad \mathbf{X}_0 := \mathbf{x} \sim \nu \in\mathbbm{R}^d,\,t\in[0,\infty), \label{eq:langevin}
\end{align}
Furthermore \cite{B:RT96} show that under certain regularity conditions the Langevin diffusion will converge at an exponential rate to the correct invariant distribution from any starting point $\mathbf{X}_0:=\mathbf{x} \in\mathbbm{R}^d$.


\subsection{Traditional Exploitation of Diffusions for Monte Carlo} \label{s:dmc:trad}
Loosely speaking, in order for a computational statistician to draw from some target distribution $\pi$, one approach she could take (although not the only approach, as will become clear in the sequel) is to allow a diffusion to evolve with dynamics such that the invariant distribution $\nu$ coincides with $\pi$ and then use its occupation measure as a proxy for $\pi$. This approach requires the simulation and storage of a continuous time object over an infinite time horizon, which is not feasible. We can, however, attempt to realise and store on a computer some finite-dimensional subset of the sample path noting that, as we can exploit the strong Markov property, we can simulate a single trajectory of arbitrary length by iteratively simulating and concatenating short sub-trajectories. We can then use the realised draws of the trajectory at the end of each sub-trajectory as a proxy for the occupation measure, and in turn realisations from $\pi$. As such, we can restrict our attention to the (iterative) simulation of finite-time horizon diffusions of the following form,
\begin{align}
\,\mathrm{d} \mathbf{X}_t 
& = \mathbf{\beta}(\mathbf{X}_t) \,\mathrm{d} t + \mathbf{\Lambda}^{1/2} \,\mathrm{d} \mathbf{W}_t, \quad \mathbf{X}_0 = \mathbf{x} \in\mathbbm{R}^d,\,t\in[0,T], \label{eq:diffusion}
\end{align}
where $\mathbf{\beta}(\mathbf{x}):=\frac{1}{2}\mathbf{\Lambda}\nabla\log\nu(\mathbf{x}):\,\mathbbm{R}^d\!\to\!\mathbbm{R}^d$ denotes the (instantaneous) drift vector. In the following let $\mathbbm{Q}^{\mathbf{x}}$ denote the probability distribution induced on path-space by (\ref{eq:diffusion}), and $\mathbbm{Q}^{\mathbf{x},\mathbf{y}}$ denote $\mathbbm{Q}^{\mathbf{x}}$ conditioned on $\mathbf{X}_T=\mathbf{y}\in\mathbbm{R}$.\\
\\
The methodological approach outlined above requires us to be able to simulate from the transition density of (\ref{eq:diffusion}), which we denote by $p_{T}({\mathbf{x},\mathbf{y}}) := \mathbbm{P}_{\mathbbm{Q}^{\mathbf{x}}}(\mathbf{X}_T \in \,\mathrm{d} \mathbf{y}\,|\, \mathbf{X}_0 =\mathbf{x})/ \,\mathrm{d} \mathbf{y}$. Under conditions given in \apxref{s:conditions}, an explicit (although not tractable) 
form for the transition density can be found as detailed in \corrolref{cor:transition}. Allowing $\mathbbm{W}^{\mathbf{x}}$, $\mathbbm{W}^{\mathbf{x},\mathbf{y}}$ to denote the law of a driftless version of (\ref{eq:diffusion}) and its conditioned counterpart, defining $w_{T}({\mathbf{x},\mathbf{y}}) := \mathbbm{P}_{\mathbbm{W}^{\mathbf{x}}}(\mathbf{X}_T \in \,\mathrm{d} \mathbf{y}\,|\, \mathbf{X}_0 =\mathbf{x})/ \,\mathrm{d} \mathbf{y}$  and $\phi(\mathbf{u}) := (||\mathbf{\Lambda}^{-1/2}\mathbf{\beta}(\mathbf{u})||^2 + \text{div\,}\mathbf{\Lambda}^{-1/2}\mathbf{\beta}(\mathbf{u}))/2$, where $||\cdot||$ denotes the Euclidean norm, we can write
\begin{align}
p_{T}({\mathbf{x},\mathbf{y}})
& = w_{T}({\mathbf{x},\mathbf{y}})\cdot\mathbbm{E}_{\mathbbm{W}^{{\mathbf{x},\mathbf{y}}}}\left[\dfrac{\,\mathrm{d} \mathbbm{Q}^{\mathbf{x}}}{\,\mathrm{d} \mathbbm{W}^{\mathbf{x}}}(\mathbf{X})\right] \\
& \propto w_{T}({\mathbf{x},\mathbf{y}})
\cdot \nu^{1/2}(\mathbf{y}) \cdot \mathbbm{E}_{\mathbbm{W}^{{\mathbf{x},\mathbf{y}}}}\left[\exp\left\{- \int^T_0\phi(\mathbf{X}_s)\,\mathrm{d} s\right\}\right].\label{eq:langevin transition}
\end{align}
Furthermore, we have the classic stationarity result
\begin{align}
\lim_{T\to\infty} p_{T}({\mathbf{x}},\cdot) \rightarrow{} \nu(\cdot), \quad\quad \forall {\mathbf{x}} \in \mathbbm{R}^d, \label{eq:classconvergence}
\end{align} 
in a sense that will be made precise later, under suitable regularity conditions (see for example  the discussion in \cite{AAP:FR05}).
\\\\ 
Now, adopting the approach outlined above in which the dynamics of the Langevin diffusion are chosen such that the invariant distribution $\nu$ coincides with some target distribution $\pi$,
we require the ability to simulate from (\ref{eq:langevin transition}). One method to simulate from (\ref{eq:langevin transition}) is to deploy a numerical discretisation scheme \cite{BK:DP,BK:NSSDE,BK:NSSDEJF} (a so-called `Unadjusted Langevin Algorithm' (ULA)). However in all but a small number of trivial cases such discretisation schemes will not provide an exact draw from the transition density, and more importantly will not recover the correct invariant distribution (it may not even possess any invariant distribution \cite{Bescom}). It has been shown that the ULA can behave very badly \cite{B:RT96}.\\
\\
An alternative to a discretisation scheme to simulate from the transition density  is to deploy a rejection sampling scheme (the so called \textit{Path-Space Rejection Sampler} (PRS) 
for diffusions \cite{AAP:BR05,MCAP:BPR08,PhD:P13,MOR:CH13,B:PJR15}). Such a rejection sampling scheme would involve simulating from the alternative proposal transition density $q(\mathbf{y})\propto w_{T}({\mathbf{x},\mathbf{y}})\cdot\nu^{1/2}(\mathbf{y})$, and accepting or rejecting the proposal with probability proportional to $\mathbbm{E}_{\mathbbm{W}^{{\mathbf{x},\mathbf{y}}}}[\exp\{- \int^T_0\phi(\mathbf{X}_s)\,\mathrm{d} s\}]$. As detailed in \apxref{s:conditions} all the conditions necessary to implement a PRS are met for the class of Langevin diffusions we consider (\ref{eq:langevin}), but the implementation of such schemes would require the construction of a proposal density which incorporates the invariant distribution $\nu$. If this were readily feasible, there would be no need to use a $\nu$-invariant diffusion at all, and so the implementation of such a scheme is not practical. \\
\\
This paper provides alternative methodology to access $\pi$ by means of simulating normalised killed Brownian motion, entirely avoiding the issues outlined above associated with discretisation and rejection-based schemes.\\ 
\\
To begin briefly outlining our approach (further details can be found in \secref{s:ergodicity}), we let $\{\mathbbm{K}^{\mathbf{x}}_T,\  T\ge 0\}$ be a collection of probability measures where $\mathbbm{K}_T^{\mathbf{x}}$ describes a probability law on $C[0,T]$ such that the following \rnd exists,
\begin{align}
\dfrac{\,\mathrm{d} \mathbbm{K}^{\mathbf{x}}_T}{\,\mathrm{d} \mathbbm{W}^\mathbf{x}_T} \left( \mathbf{X} \right) 
& \propto \exp{\left\{-\int_0^T\phi(\mathbf{X}_s)\,\mathrm{d} s \right\}}. \label{eq:k killing}
\end{align}
Letting the marginal density of $\mathbbm{K}^{\mathbf{x}}_T$ at time $T$ be denoted by $k_T({\mathbf{x},\mathbf{y}})$, we establish that,
\begin{align}
k_T({\mathbf{x},\mathbf{y}}) \propto w_T({\mathbf{x},\mathbf{y}}) \cdot \underbrace{\mathbbm{E}_{\mathbbm{W}_T^{{\mathbf{x},\mathbf{y}}}}\left[\dfrac{\,\mathrm{d} \mathbbm{K}^{\mathbf{x}}_T}{\,\mathrm{d} \mathbbm{W}^{\mathbf{x}}_T} \left( \mathbf{X} \right)\right]}_{=:P(\mathbf{X})}, \label{eq:killed transition}
\end{align}
which follows by noting that if, say, $\phi$ is everywhere non-negative, then $\mathbbm{K}_T^{\mathbf{x}}$ is the law of killed Brownian motion conditioned to survive until time $T$ with instantaneous killing rate $\phi(\mathbf{X}_s)$. Intuitively, recalling (\ref{eq:langevin transition}) and  (\ref{eq:classconvergence}), we may naturally posit that
\begin{align}
\lim_{T\to\infty}  k_T({\mathbf{x},\mathbf{y}}) \propto p_{T}({\mathbf{x},\mathbf{y}}) \cdot \nu^{-1/2}(\mathbf{y}) \to \nu ^{1/2}(\mathbf{y}),
\end{align}
and so if $\nu ^{1/2}$ is itself integrable we may expect that $k_{T}({\mathbf{x}},\cdot)$ converges to a probability density proportional to $\nu^{1/2}$. Indeed, this is established in the following theorem,
\begin{theorem}[Quasi-Stationarity of Killed Brownian motion] \label{thm:qsd}
Suppose that $\nu^{1/2}\in L^1$, is bounded, and there exists some $\gamma \geq 1/2$ such that
\begin{align}
\label{eq:Q1}
\liminf_{\mathbf{x} \to \infty } \left(
{\Delta \nu ({\mathbf{x}}) \over \nu ^{\gamma + 1/2}({\mathbf{x}} ) } - { \gamma \| \nabla \nu  ({\mathbf{x}})\|^2
\over \nu ^{\gamma + 3/2}({\mathbf{x}})}
\right) >0   \hspace{2cm}
\end{align}
where $\Delta $ represents the Laplacian:
$$
\Delta f( \mathbf{x} )= \sum_{i=1}^d {\partial^2 f \over \partial x_i^2},
$$
then we have the following convergence in $L^1$:
\begin{align}
\lim_{T \to \infty }
  k_{T}({\mathbf{x}},\cdot )= \frac{ \nu^{1/2} (\cdot )}{\int _\mathbbm{R}^d
 \nu^{1/2} (\mathbf{z}) d\mathbf{z} }.
\end{align}
\proof See \secref{apx:qstheory}. \endproof
\end{theorem}
\noindent 
Consequently, a computational statistician who wishes to draw from a target distribution $\pi$ has a novel approach to explore. She can choose the dynamics of (\ref{eq:langevin}) such that the invariant distribution $\nu\equiv\pi^2$, and generate a trajectory of killed Brownian motion from $\mathbbm{K}_T^{\mathbf{x}}$. Then by \thmref{thm:qsd} we have (noting as required we have $\nu^{1/2}\in L^1$ as $\nu^{1/2}\equiv\pi$), as $T\to\infty$:
\begin{align}
k_{T}({\mathbf{x}},\cdot) \xrightarrow{L^1} \pi(\cdot), \quad\quad \forall {\mathbf{x}} \in \mathbbm{R}^d. \label{eq:l1 convergence2}
\end{align}
To generate a sample trajectory from $\mathbbm{K}_T^{\mathbf{x}}$, we need only sample Brownian motion killed instantaneously at a state-dependent rate $\phi (\mathbf{X}_s)$, and this can be done by thinning a dominating Poisson process. This provides an unexplored avenue for sampling (asymptotically) from a target distribution $\pi$ without the need to resort to an unknown approximation.\\
\\
In the remainder of this section we provide a proof and discussion of \thmref{thm:qsd}. We return in \secref{s:emcfd} to provide detail on how to exploit practically \thmref{thm:qsd} and (\ref{eq:l1 convergence2}).


\subsection{A Novel Quasi-Stationary Approach for Monte Carlo} \label{s:dmc:qsd}\label{s:ergodicity}\label{apx:qstheory}

In \secref{s:dmc:trad} we obtained an interesting limit result for the distribution of killed Brownian motion at a given time conditioned on not being killed by that time --- a so called {\em quasi-limiting} distribution. The literature on quasi-limiting distributions and quasi-stationarity is vast (see for example \cite{BK:QSD} for an entry into the area) and contains many unexpected subtleties which distinguish it from the theory of stationarity for Markov chains and processes.\\
\\
Here we present a proof of \thmref{thm:qsd}.


\begin{proof}[Proof (Theorem \ref{thm:qsd})]
The diffusion in (\ref{eq:langevin}) has  generator given by
$$
 \mathfrak{A} f ({\mathbf{x}} )= {1\over 2} \Delta f ({\mathbf{x}}) +  {1\over 2}  \nabla \log \nu ({\mathbf{x}} ) \cdot \nabla f ({\mathbf{x}} ).
$$
As $\nu $ is bounded, we assume without loss of generality that its upper bound is $1$.  Our proof shall proceed by checking the conditions of Corollary 6  of \cite{AAP:FR05}, which establishes the result: we need to check that the following are satisfied:
\begin{enumerate}
\item For all $\delta >0$, the discrete time chain $\{X_{n\delta }, n=0, 1, 2,\ldots \}$
is irreducible;
\item
All closed bounded sets are petite;
\item
We can find  a drift function $V({\mathbf{x}} ) =  \nu({\mathbf{x}} )^{-\gamma}$ for some $\gamma >0$, that satisfies the condition
\begin{align}
\label{eq:drift}
{\mathfrak{A}} V^{\eta }({\mathbf{x}} ) \le - c_\eta V({\mathbf{x}})^{\eta - \alpha} 
\end{align}
for ${\mathbf{x}} $ outside some bounded set, for each $\eta \in [\alpha ,1]$ with associated positive constant $c_{\eta }$ and where $\alpha = 1-(2\gamma )^{-1}$.
\end{enumerate}
The first condition holds for any regular diffusion since the diffusion possesses positive continuous transition densities over time intervals $t>0$; and positivity and continuity of the density also implies the second condition. For the final condition we require that
\begin{align}
\label{eq:ratiodrift}
\limsup_{|{\mathbf{x}} | \to \infty }
{{\mathfrak{A}} V^{\eta }({\mathbf{x}} )
\over
V^{\eta - \alpha } ({\mathbf{x}})}
<0\ .
\end{align}
Now by direct calculation
\begin{align}
{\mathfrak{A}} V^{\eta }({\mathbf{x}} ) = {\eta \gamma \over 2}
  \nu  ({\mathbf{x}} )^{-\gamma \eta -2} 
\left[
{\eta \gamma \|\nabla \nu ({\mathbf{x}} ) \|^2  }
-
\nu ({\mathbf{x}} ) \Delta \nu ({\mathbf{x}} )
\right]
\end{align}
 so that
\begin{align}
{{\mathfrak{A}} V^{\eta }({\mathbf{x}} )\over V({\mathbf{x}} )^{\eta - \alpha}}
=
{\eta  \gamma \nu  ({\mathbf{x}} )^{-3/2-\gamma } \over 2}
\left[
{\eta \gamma \|\nabla \nu ({\mathbf{x}} ) \|^2 }
-
\nu ({\mathbf{x}} ) \Delta \nu ({\mathbf{x}} )
\right].
\end{align}
Therefore (\ref{eq:ratiodrift}) will hold whenever (\ref{eq:Q1}) is true since we have the constraint that $\eta \le 1$ and $\|\nabla \nu ({\mathbf{x}} ) \|^2$ is clearly non-negative. As such the result holds as required.
\end{proof}
\noindent 
Note that the condition in (\ref{eq:Q1}) is essentially a condition on the tail of $\nu$. This will hold even for heavy-tailed distributions, and we show this is the case for a class of 1-dimension target densities in \apxref{apx:polynomial}.


\section{Exact Monte Carlo for Diffusions} \label{s:emcfd}

In this section we develop Monte Carlo methodology for the practical exploitation of the theory developed in \secref{s:dmc}. We begin in \secref{s:nsea} by outlining methodology for the exact simulation of diffusions based on path-space rejection sampling, with particular focus on their application to our desired quasi-stationary process. In \secref{s:nsea} we focus on an importance sampling construction which requires a number of trajectories in order to conduct normalisation. The drawback of this approach is the variance of the importance sampling weights increases with diffusion time. We address this issue in \secref{s:emcfd:smc} by embedding the methodology developed in \secref{s:nsea} within a Sequential Monte Carlo (SMC) framework, providing a theoretical underpinning for doing so. In \secref{s:emcfd:smc} we also provide algorithmic detail for the resulting {\em Quasi-Stationary Monte Carlo (QSMC)} algorithm.\\
\\
We omit at this stage any discussion of big data, returning to this context in \secref{s:scale}.


\subsection{Path-space Rejection Samplers} \label{s:emcfd:pims} \label{s:langevin}\label{s:nsea}
Path-space Rejection Samplers (PRS) is the name we give to a class of rejection samplers operating on diffusion path-space over a finite time horizon, in which appropriately weighted finite dimensional subsets of sample paths are drawn from some target measure (in our case $\mathbbm{K}_T^{\mathbf{x}}$), by means of simulating from a tractable equivalent measure with respect to which the target has bounded \rnd. As established in \thmref{thm:qsd}, one such equivalent measure is provided by Brownian motion with an appropriate pre-conditioning matrix $\mathbf{\Lambda}$, denoted by $\mathbbm{W}^{\mathbf{x}}$. In this paper for simplicity we restrict our attention to diagonal pre-conditioning matrices, $\mathbf{\Lambda}$, but note there are no theoretical or methodological limitations which would prevent us from considering other non-diagonal matrices. $\mathbbm{W}^{\mathbf{x}}$ is a convenient measure as finite dimensional subsets of sample paths can be drawn without error (by means of simulating scaled standard $d$-dimensional Brownian motion, see \cite[\textsection 2.8]{PhD:P13}). In the interest of clarity we defer, wherever possible, the technical conditions of PRS to \apxref{s:conditions}.\\
\\
We first consider the simplest path-space rejection sampling setting (introduced by \cite{AAP:BR05}). Provided that there exists $M<\infty$ such that $\frac{\,\mathrm{d} \mathbbm{K}_T^{\mathbf{x}}}{\,\mathrm{d} \mathbbm{W}_T^{\mathbf{x}}}(\mathbf{X})\leq M$, we can proceed as in standard rejection sampling and draw $\mathbf{X}\sim\mathbbm{W}^{\mathbf{x}}$ and accept the sample path $(I=1)$ with probability $P(\mathbf{X}) := \frac{1}{M}\frac{\,\mathrm{d} \mathbbm{K}_T^{\mathbf{x}}}{\,\mathrm{d} \mathbbm{W}_T^{\mathbf{x}}}(\mathbf{X})\in[0,1]$, then $(\mathbf{X}|I=1)\sim\mathbbm{K}_T^{\mathbf{x}}$. Under \conref{cond:phi} of \apxref{s:conditions}, such a bound can be found and
\begin{align}
P(\mathbf{X}) = e^{\Phi T}\cdot \exp\left\{-\int^T_0 \phi(\mathbf{X}_s)\,\mathrm{d} s\right\}\in [0,1]. \label{eq:waccprob}
\end{align}
Direct simulation of an event of probability $P(\mathbf{X})$ is not possible as evaluation of the integral in (\ref{eq:waccprob}) requires the realisation of an entire (infinite dimensional) sample path. However, as noted in \cite{AAP:BR05,B:BPR06,MCAP:BPR08,PhD:P13,MOR:CH13,B:PJR15}, an indirect and unbiased approach to simulate $P(\mathbf{X})$, which can be implemented using only a finite dimensional realisation of the sample path is possible.\\
\\
This approach to simulate $P(\mathbf{X})$ requires upper and lower bounds for $\phi(\mathbf{X}_{[0,T]})$ (to ensure that it can be computed in finite expected time). Such a condition on $\phi(\mathbf{X}_{[0,T]})$ is not restrictive as it is possible to simulate $\mathbf{X}\sim\mathbbm{W}^{\mathbf{x}}$ in conjunction with a sequence of stopping times $\{\tau_0:=0,\tau_1,\tau_2,\ldots\}$, such that for the interval between any two successive stopping times the sample path is constrained to a hypercuboid, $\mathbf{X}_{[\tau_{i-1},\tau_i]}\in[\mathbf{\ell}^{(i)}_{\mathbf{X}},\mathbf{\upsilon}^{(i)}_{\mathbf{X}}]$ where $\mathbf{\ell}^{(i)}_{\mathbf{X}},\mathbf{\upsilon}^{(i)}_{\mathbf{X}}\in\mathbbm{R}^d$ and for any $a = (a_1,\ldots,a_d),b=(b_1,\ldots,b_d) \in\mathbb{R}^d$, $[a,b] := \otimes_{i=1}^d [a_i,b_i]$. As a consequence of \resref{res:compact}, such a {\em path-space layer} provides upper and lower bounds on $\phi(\mathbf{X})$ over any finite time interval (which will vary depending on the interval being considered).\\
\\
The construction and simulation of the path-space layer is provided in detail in \apxref{s:pathspace layer}. The user chooses the extent of the state space inscribed by the hypercuboid, $[\mathbf{\ell}^{(i)}_{\mathbf{X}},\mathbf{\upsilon}^{(i)}_{\mathbf{X}}]$, in each time-interval. The algorithm then simulates the time and location, ($\tau_{i},\mathbf{X}_{\tau_i}$), at which the process first leaves this hypercuboid. Furthermore, the construction detailed in \apxref{s:pathspace layer} allows one to  simulate the sample path at any other desired intermediate point.  As such we further denote by $R^{(i)}_{\mathbf{X}}$ the path-space information and by $U^{(i)}_{\mathbf{X}}\in\mathbbm{R}$ and $L^{(i)}_{\mathbf{X}}\in\mathbbm{R}$ as the upper and lower bounds for $\phi(\mathbf{X}_{[\tau_{i-1},\tau_i]})$ (i.e. for the $i^{\text{th}}$ stopping time interval, $[\tau_{i-1},\tau_i]\subseteq[0,T]$). Conditional on the simulated layer ($R:=R(\mathbf{X})\sim\mathcal{R}$, comprised of $n_R:=\inf\{n:\tau_n\geq T\}$ segments) we can represent our acceptance probability for the draw $\mathbf{X}\sim\mathbbm{W}^{\mathbf{x}}_T$ as follows
{\small
\begin{align}
{P}(\mathbf{X}) 
& = \underbrace{\exp\left\{\Phi T -\sum^{n_R}_{i=1} L^{(i)}_{\mathbf{X}}\cdot[(\tau_i\wedge T)\!-\!\tau_{i-1}]\right\}}_{=:P^{(1)}(\mathbf{X}) \in [0,1]} \cdot \prod^{n_R}_{i=1}\Bigg[\underbrace{\exp\left\{-\!\!\int^{\tau_i\wedge T}_{\tau_{i-1}}\!\left(\phi(\mathbf{X}_s)\!-\!L^{(i)}_{\mathbf{X}}\right)\,\mathrm{d} s\right\}}_{=:P^{(2,i)}(\mathbf{X})}\Bigg] \label{eq:p1decomp} \\
& = \prod^{n_R}_{i=1} \Bigg[\underbrace{\exp\left\{(\Phi - L^{(i)}_{\mathbf{X}})\cdot[(\tau_i\wedge T)\!-\!\tau_{i-1}]\right\}}_{=:P^{(1,i)}(\mathbf{X}) \in [0,1]} \cdot \underbrace{\exp\left\{-\int^{\tau_i\wedge T}_{\tau_{i-1}}\!\left(\phi(\mathbf{X}_s)\!-\!L^{(i)}_{\mathbf{X}}\right)\,\mathrm{d} s\right\}}_{=:P^{(2,i)}(\mathbf{X})}\Bigg]. \label{eq:p1decomp2}
\end{align} }
\noindent The path-space layer ensures $\phi(\mathbf{X}_{[\tau_{i-1},\tau_i]})\in[L^{(i)}_{\mathbf{X}},U^{(i)}_{\mathbf{X}}]$, and so $\forall i$
\begin{align}
{P}^{(2,i)}(\mathbf{X})\in\left[e^{-(U^{(i)}_{\mathbf{X}}-L^{(i)}_{\mathbf{X}})T},1\right] \subseteq (0,1]. \label{eq:pois ap}
\end{align}
From (\ref{eq:p1decomp}) note that the simulation of an event of probability ${P}(\mathbf{X})$ can be conducted by first simulating an event of probability ${{P}}^{(1)}(\mathbf{X})$ (which is straight-forward to simulate as it is directly determined using the path-space layer simulated), and then $n_R$ events of probability ${{P}}^{(2,i)}(\mathbf{X})$ (which are functions of the infinite-dimensional proposal sample path). Acceptance of ${P}(\mathbf{X})$ is equivalent to accepting all events in this extended construction. We can simulate an event of probability $P^{(2,i)}(\mathbf{X})$ indirectly and using a finite dimensional realisation of the sample path by means of the following \thmref{thm:acceptprob2}.
\begin{theorem}[Segment Acceptance Probability] \label{thm:acceptprob2}
Let $\mathbbm{U}^{(i)}_R$ be the distribution of the tuple $\{\xi^{(i)}_1,\xi^{(i)}_2,\ldots{}\}\cap[\tau_{i-1},(\tau_i\wedge T)]$, where $\xi^{(i)}_j = \xi^{(i)}_{j-1}+\zeta^{(i)}_j$, $\zeta^{(i)}_j \sim \mathrm{Exp}(U^{(i)}_{\mathbf{X}}-L^{(i)}_{\mathbf{X}})$ and $\xi^{(i)}_0:=\tau_{i-1}$. Denoting by $\kappa_i$ the cardinality of this tuple, we have,
\begin{align} \label{eq:expacceptprob}
{P}^{(2,i)}(\mathbf{X})
 = \mathbbm{E}_{\mathbbm{U}^{(i)}_R}\Bigg[\prod^{\kappa_i}_{j=1}\underbrace{\Bigg(\dfrac{U^{(i)}_{\mathbf{X}}\!-\!\phi\big(\mathbf{X}_{\xi^{(i)}_j}\big)}{U^{(i)}_{\mathbf{X}}\!-\!L^{(i)}_{\mathbf{X}}}\Bigg)}_{=:P^{(2,i,j)}(\mathbf{X})}\,\Bigg|\,\mathbf{X}\Bigg]. 
\end{align}
\proof Follows directly from \cite[\textsection 3.1]{B:PJR15}, noting the law induced by the uniform distribution of a heterogenous Poisson process with intensity which changes at the time of each event so that it is, informally, $U^{(i)}_{\mathbf{X}}-L^{(i)}_{\mathbf{X}
$ on the interval $[\tau_{i-1},(\tau_i\wedge T)]$ has an equivalent distribution to the procedure above (see for instance \cite{PhD:P13,BK:PP}). \endproof
\end{theorem}
\noindent \thmref{thm:acceptprob2} allows us to simulate a Bernoulli random variable of success probability $P^{(2,i)}(\mathbf{X})$, by sampling the sample path at a finite number of exponentially distributed times, and simulating events of probability $P^{(2,i,j)}(\mathbf{X})$ (all of which need to be accepted for success). As such we now have a path-space rejection sampler for $\mathbbm{K}^{\mathbf{x}}_T$. Over a pre-specified time horizon, then (\ref{eq:p1decomp}) and \thmref{thm:acceptprob2}  together allow us to determine acceptance of $\mathbf{X}\sim\mathbbm{W}^{\mathbf{x}}_T$, but only using a finite dimensional path skeleton. For clarity of presentation we omit the technical and algorithmic details of our PRS for $\mathbbm{K}^{\mathbf{x}}_T$, which can instead be found in \apxref{apx:prs} and \algref{alg:altuea}.
\\\\
Recall that we are ultimately interested in the transition density $k_T(\mathbf{x},\mathbf{y})$ when $T\to\infty$. As such, we may more naturally be interested in identifying the time and location at which the Brownian motion sample path proposal is killed (denoted by ($\bar{\tau},\mathbf{X}_{\bar{\tau}}$)), as opposed to determining over a pre-specified time interval whether or not killing has occurred (as in the PRS for $\mathbbm{K}^{\mathbf{x}}_T$ outlined above). Note that killing will almost-surely occur in finite time (unless $\Phi=U_\mathbf{X}$). As the proposal process is strongly Markov, by appealing to (\ref{eq:p1decomp2})
the following is an informal description of a valid (although not unique) weight process in the interval $[\tau_{i-1},\tau_i]$, 
\begin{align}
\,\mathrm{d} w_t \label{eq:weightprocess}
& = \left\{ 
\begin{array}{l l}
-w_{t-}\cdot (L^{(i)}_\mathbf{X} -\Phi) \cdot \,\mathrm{d} t & \text{w.p. } (1-\,\mathrm{d} t)\cdot (U^{(i)}_\mathbf{X}-L^{(i)}_\mathbf{X}),\\
-w_{t-} \cdot \dfrac{\phi(\mathbf{X}_t)-L^{(i)}_\mathbf{X}}{U^{(i)}_\mathbf{X}-L^{(i)}_\mathbf{X}} & \text{w.p. } \,\mathrm{d} t \cdot (U^{(i)}_\mathbf{X}-L^{(i)}_\mathbf{X}).
 \end{array}\right.
\end{align}
The first line of (\ref{eq:weightprocess}) can be interpreted as a constant weight degradation in an increment of time arising from $P^{(1,i)}(\mathbf{X})$ of (\ref{eq:p1decomp2}), and the second line of (\ref{eq:weightprocess}) the weight degradation if the increment coincides with an event time $\xi^{(i)}_{\cdot}$ as in \thmref{thm:acceptprob2}.\\
\\
Consequently, identifying the pair ($\bar{\tau},\mathbf{X}_{\bar{\tau}}$) is trivial once one identifies the killing time (where $u\sim\text{U}[0,1]$),
\begin{align}
\bar{\tau} := \inf\{t\in[0,\infty): 1 - w_t(\mathbf{X}_{[0,\infty)}) \leq u\}.
\end{align}
Fortuitously the simulation techniques available to us for our   Brownian motion proposal constrained to a path-space layer (as detailed in \apxref{s:pathspace layer}), is sufficiently tractable to allow us to simulate ($\bar{\tau},\mathbf{X}_{\bar{\tau}}$) exactly. We term this algorithm \textit{Killed Brownian Motion (KBM)}, full details of which can be found in \apxref{apx:kbm} and \algref{alg:kbm}.
\\\\
The Monte Carlo algorithm we focus on in this paper is an importance sampling variant of the KBM algorithm described above, in which Brownian motion with path-space information is simulated and weighted continuously in time, with evaluations of the trajectory occurring at exponentially distributed times as given in \thmref{thm:acceptprob2}. If $L_{\mathbf{X},j}:=\{L^{(i)}_\mathbf{X}:\xi_j\in[\tau_{i-1},\tau_i]\}$ and $U_{\mathbf{X},j}:=\{U^{(i)}_\mathbf{X}:\xi_j\in[\tau_{i-1},\tau_i]\}$), then the weight of the trajectory at time $t\in[0,\infty)$ is simply
\begin{align}
w_t(\mathbf{X}) = \exp\{\Phi t\}\cdot \underbrace{\exp\bigg\{ -\!\!\!\sum_{i: t\in[0,\tau_i]}\! \!\!L^{(i)}_{\mathbf{X}}\cdot[(\tau_i\wedge t)-\tau_{i-1}]\bigg\} \cdot 
\prod_{j:\xi_j\leq t}\left(\dfrac{U_{\mathbf{X},j}-\phi(\mathbf{X}_{\xi_j})}{U_{\mathbf{X},j}-L_{\mathbf{X},j}}\right)}_{=:w^*_t(\mathbf{X})}. \label{eq:is weight2}
\end{align}
In practice we only need to calculate this weight up to a common constant of proportionality, so we calculate $w^*_t(\mathbf{X})$.  We term this approach the Importance Sampling Killed Brownian Motion (IS-KBM) algorithm, which we present in \algref{alg:is-kbm} and schematically in \figref{Fig:IS}. The advantage of this approach over KBM (presented in \apxref{apx:kbm}) is that we can recover for any 
$t\in[0,\infty)$ a weighted trajectory, but to do so we need to halt the algorithm. If we require the trajectory at time $t$, this can be done by simply  setting $\xi_j:= \xi_j\wedge t$ after \algstref{alg:is-kbm}{alg:qsmc:scalesim2}, halting upon simulation of the trajectory at time $t$, and returning $\mathbf{X}_t$ along with the (un-normalised) importance weight $w^*_t$.\\
\\
Repeating this procedure several times allows us to use self-normalised importance sampling to approximate expectations with respect to the law of the killed process conditioned to remain alive until any desired time.
\begin{figure}[h] 
\centering
\subfigure[Layers simulated for an illustrative (un-realised) Brownian motion sample path. Red horizontal lines mark the bounds of the layer for a time segment, red vertical lines mark stopping times in which there is a change in layer.]{
\includegraphics[width=0.475\textwidth]{Example_Layer.pdf}} \hspace*{0.05cm}
\subfigure[Brownian motion sample path simulated at the finite collection of event times (black dots) as given by exogenous randomness, and path-space layers (red lines), that are required within IS-KBM.]{
\includegraphics[width= 0.475\textwidth]{Example_X.pdf}} 
\subfigure[Evaluation of $\phi$ at the event times of the sample path given by (b).]{
\includegraphics[width= 0.475\textwidth]{Example_Phi.pdf}} 
\hspace*{0.05cm}
\subfigure[Un-normalised importance weight process of sample path, comprising exponential growth between event times and discrete jumps at event times.]{
\includegraphics[width= 0.475\textwidth]{Example_Weight.pdf}} 
\caption{Schematic of the Importance Sampling Killed Brownian Motion (IS-KBM) Algorithm in 1 dimension, simulated over an interval of length $T=4$. To summarise, we would begin by simulating a set of path-space layers over the desired interval (as shown in (a)), realising the sample path at a collection of event times given by exogenous randomness (as shown in (b)), evaluating $\phi$ at those event times (as shown in (c)), and then recovering the (un-normalised) weight process which has exponential growth between weight times (with exponent $-L^{(i)}_{\mathbf{X}}$ as given by the associated path-space layer for $\phi$) and discrete jumps at an event time $\xi$ of size $(\phi(X_\xi)-L^{(i)}_{\mathbf{X}}) /(U^{(i)}_{\mathbf{X}}-L^{(i)}_{\mathbf{X}})$ (as shown in (d)).} \label{Fig:IS}
\end{figure}

\begin{algorithm}[h]
	\caption{Importance Sampling Killed Brownian Motion (IS-KBM) Algorithm} \label{alg:is-kbm}
	\begin{enumerate}
	\item Initialise: Input initial value $\mathbf{X}_0$. Set $i=1$, $j=0$, $\tau_0=0$, $w^*_0=1$.
	\item $R$: Simulate layer information $R^{(i)}_{\mathbf{X}}\sim\mathcal{R}$ as per \apxref{s:pathspace layer}, obtaining $\tau_i$, $L^{(i)}_\mathbf{X}$, $U^{(i)}_\mathbf{X}$. \label{alg:is-kbm:layer}
	\item $\zeta$: Simulate $\zeta \sim \text{Exp}(U^{(i)}_\mathbf{X}-L^{(i)}_\mathbf{X})$. \label{alg:is-kbm:loop} 
	\item $\xi_j$: Set $j=j+1$ and $\xi_j = (\xi_{j-1}+\zeta)\wedge \tau_i$. \label{alg:qsmc:scalesim2}
	\item $w^*_{\xi_j}$: Set $w^*_{\xi_j}=w^*_{\xi_{j-1}}\cdot\exp\{-L^{(i)}_\mathbf{X}[ \xi_j-\xi_{j-1}]\}$.
	\item $\mathbf{X}_{\xi_j}$: Simulate $\mathbf{X}_{\xi_j}\sim \left.\text{MVN}(\mathbf{X}_{\xi_{j-1}},(\xi_j-\xi_{j-1})\mathbf{\Lambda})|R^{(i)}_\mathbf{X} \right.$.
	\item $\tau_i$: If $\xi_j=\tau_i$, set $i=i+1$, and return to \stepref{alg:is-kbm:layer}. Else set $w^*_{\xi_j}=w^*_{\xi_j}\cdot (U^{(i)}_\mathbf{X}-\phi(\mathbf{X}_{\xi_j}) )/ (U^{(i)}_\mathbf{X} -  L^{(i)}_\mathbf{X}) $ and return to \stepref{alg:is-kbm:loop} .\label{alg:is-kbm:w}
	\end{enumerate}
\end{algorithm}

\FloatBarrier

\subsection{Normalisation through Sequential Monte Carlo} \label{s:emcfd:smc} \label{s:convergence}
 \label{s:qsmc}
Although the IS-KBM framework we have presented thus far is conceptually appealing, in order to normalise we require a number of trajectories, and the variance of the importance sampling weights of these trajectories will increase with time. We present one possible strategy for addressing this difficulty in this section.
\\
\\
Sequential Monte Carlo (SMC) methods combine importance sampling and resampling techniques in order to approximate a sequence of distributions of interest. First used by \cite{smc:methodology:GSS93} in the context of target tracking, such methods apply naturally when one is interested in approximating particular time marginals of a sequence of distributions defined upon an increasing state space (see, for example, \cite[Section 24.3]{amj26:DJ11} for a general presentation). They propagate an ensemble of $N$ weighted samples over time, attaching importance weights to each sample and using resampling to replicate samples with large weights and to eliminate those with small weights according to a stochastic rule which ensures that the sample remains properly weighted for the actual target distribution. They are particularly applicable in the current setting as the resampling step allows computational effort to be focussed on trajectories which have survived. Loosely speaking, one can think of resampling as a mechanism in which samples with low survival probability are eliminated and those with a high probability are allowed to branch.
\\\\
We apply such an algorithm here in order to provide, iteratively, a population of samples (or \textit{particles}) suitable for approximating the law of the killed process conditioned to survive until the current time. This can be done for both the rejection sampling (KBM) and importance sampling (IS-KBM) versions of the algorithm presented in \secref{s:nsea}. For clarity we present details of the KBM embedding in SMC in \apxref{apx:rqsmc}, and instead consider here the importance sampling version as it more closely aligns with typical SMC constructions.\\
\\
Considering the embedding of IS-KBM (\algref{alg:is-kbm}) within SMC, 
we need to consider the times at which to conduct resampling. One could in principle devise a scheme to resample in continuous time, however for simplicity and theoretical reasons we only conduct resampling at pre-specified times $t_0:=0< t_1<\ldots{}<t_m:=T$ (we would naturally choose these such that $t_i-t_{i-1}:=T/m$). This is presented as the \textit{Quasi-Stationary Monte Carlo (QSMC)} algorithm in \algref{alg:qsmc} in which we define ${N_{\text{eff}}}$ as the effective sample size of \cite{smc:methodology:KLW94}, and $N_{\text{th}}$ as a user chosen threshold for which will conduct (on a given iteration) resampling of the particle set at the resampling time. The algorithm outputs the weighted particles at the end of each iteration.
\begin{algorithm}[h]
	\caption{Quasi-Stationary Monte Carlo Algorithm (QSMC) Algorithm.} \label{alg:qsmc}
    \begin{enumerate}
	\item \textbf{Initialisation Step ($i=0$)}
	\begin{enumerate}
	\item Input: Starting value, $\hat{\mathbf{x}}$, number of particles, $N$, and  set of $m$ times $t_{1:m}$.
	\item $\mathbf{X}^{(\cdot)}_0$: For $k$ in $1$ to $N$ set $\mathbf{X}^{(1:N)}_{t_0}=\hat{\mathbf{x}}$ and $w^{(1:N)}_{t_0}=1/N$.
	\end{enumerate}
	\item \textbf{Iterative Update Steps ($i=i+1$ while $i\leq m$)}
	\begin{enumerate}
	\item $N_{\text{eff}}$: If ${N_{\text{eff}}}\leq N_{\text{th}}$ then for $k$ in $1$ to $N$ resample $\mathbf{X}^{(k)}_{t_{i-1}} \sim \tilde{\pi}^N_{t_{i-1}}$, the empirical distribution defined by the current set of weighted particles, and set $w^{(k)}_{t_{i-1}} = 1/N$.
	\item For $k$ in $1$ to $N$,
	\begin{enumerate}
	\item $\mathbf{X}^{(\cdot)}_{t_i}$ Simulate $\mathbf{X}^{(k)}_{t_i} | \mathbf{X}^{(k)}_{t_{i-1}}$ along with un-normalised weight increment $w^*_{t_i-t_{i-1}}$ as per \algref{alg:is-kbm}.\label{alg:qsmc:scalesim1} 
	\item ${w'}^{(\cdot)}_{\!\!t_i}$: Calculate un-normalised weights, ${w'}^{(k)}_{\!\!t_i} = w^{(k)}_{t_{i-1}}\cdot w^{*}_{t_i-t_{i-1}}$. \label{alg:qsmc:accrej1}
	\end{enumerate}
	\item $w^{(\cdot)}_{t_i}$: For $k$ in $1$ to $N$ set $w^{(k)}_{t_i} = {w'}^{(k)}_{\!\!t_i}/\sum^N_{l=1} {w'}^{(l)}_{\!\!t_i}$.
	\item $\tilde{\pi}^N_{t_i}$: Set $\tilde{\pi}^N_{t_i}(\!\,\mathrm{d} \mathbf{x}) := \sum^N_{k=1} w^{(k)}_{t_i}\cdot\delta_{\mathbf{X}^{(k)}_{t_i}}(\!\,\mathrm{d} \mathbf{x})$.
	\end{enumerate}	
	\end{enumerate}
\end{algorithm}
\\\\ \noindent Given the output from \algref{alg:qsmc}, we can then estimate the target distribution $\pi$. We first choose a burn-in time, $t^*$, such that we believe we have achieved convergence to the quasi-stationary distribution
of killed Brownian motion. Our approximation to the law of the killed process is then simply the weighted occupation measures of the particle trajectories in the interval $[t^*,T]$. More precisely, using the output of the QSMC algorithm then
\begin{align}
\pi(\!\,\mathrm{d} \mathbf{x}) \approx \hat{\pi}(\!\,\mathrm{d} \mathbf{x}) 
& := \dfrac{1}{|\{i:t_i\in[t^*,T]\}|}\sum^m_{i=1:t_i\in[t^*,T]} \sum^N_{k=1}w^{(k)}_{t_i}\cdot \delta_{\mathbf{X}^{(k)}_{t_i}}(\!\,\mathrm{d} \mathbf{x}). \label{eq:occupation2}
\end{align}
In \apxref{app:smc} we use existing theory for SMC to obtain a central limit theorem for estimators based on the output QSMC. We defer the theoretical characterisation of the SMC approach used to \secref{sec:smctheory}, at which point we can study QSMC with full application to our big data context. Provided the timescale upon which resampling is conducted is not too coarse, and low variance resampling strategies are employed (for instance, \cite{DCM05}), one would not expect qualitative differences in performance between the approaches.


\FloatBarrier
\section{Sub-sampling} \label{s:scale} 
In \secref{s:emcfd} we presented two methodological implementations for simulating our desired quasi-limiting distribution -- in particular, we simulated a number of trajectories of Brownian motion which were either killed or re-weighted after an exponentially distributed amount of time (with state space dependent intensity). As we will outline, in the case of a big data setting the most na\"{i}ve implementations of the killing and re-weighting mechanisms are the $\mathcal{O}(n)$ bottlenecks that we need to overcome.\\
\\
We begin here by studying a simple property of our estimator of $P(\mathbf{X})$, the survival probability of a Brownian motion path. In \secref{s:repest} we exploit this property to overcome the $\mathcal{O}(n)$ discussed above.\\
\\ 
To recap, in \secref{s:nsea} both our rejection sampler (KBM) and our importance sampler (IS-KBM) used the estimator
\begin{align} \label{eq:Px}
  {P}(\mathbf{X})  
  & = \prod^{n_R}_{i=1} \Bigg[
  \underbrace{\exp\left\{(\Phi - L^{(i)}_{\mathbf{X}})\cdot[(\tau_i\wedge T)\!-\!\tau_{i-1}]\right\}}_{=:P^{(1)}(\mathbf{X}) \in [0,1]}  \nonumber\\
  & \quad\quad\quad\quad\quad\quad \cdot 
 \mathbbm{E}_{\mathbbm{U}^{(i)}_R}\Bigg[\prod^{\kappa_i}_{j=1}\underbrace{\Bigg(\dfrac{U^{(i)}_{\mathbf{X}}\!-\!\phi\big(\mathbf{X}_{\xi^{(i)}_j}\big)}
 {U^{(i)}_{\mathbf{X}}\!-\!L^{(i)}_{\mathbf{X}}}\Bigg)}_{=:P^{(2,i,j)}(\mathbf{X}) \in [0,1]}\,\Bigg|\,\mathbf{X}\Bigg]
  \Bigg],                                               
\end{align}
where the outer product is over the different layers, the inner product is over the  events within a layer, and the expectation is with respect to the number and position of these events (for which we have observed the proposal Brownian motion sample path). Note that the $\mathcal{O}(n)$ bottleneck within this estimator is the evaluation of $\phi(\cdot)$.\\
\\
We can replace each term within the inner product, $P^{(2,i,j)}(\mathbf{X})$, by an unbiased estimator, and this identity will still hold in expectation. That is we introduce an auxiliary random variable $A\sim\mathcal{A}$, and construct estimator $\tilde{\phi}_A(\cdot)$ such that
\begin{align} \label{eq:phiunbiased}
 \mathbbm{E}_{\mathcal{A}}\left[\tilde{\phi}_A(\cdot)\right] = \phi(\cdot).
\end{align}
Recall that to ensure in (\ref{eq:Px}) that ${P}(\mathbf{X})\in[0,1]$, we used the bounds provided by each layer to find $L^{(i)}_{\mathbf{X}},U^{(i)}_{\mathbf{X}}\in\mathbbm{R}$ such that $\phi(\mathbf{X}_{[\tau_{i-1},\tau_i]})\in[L^{(i)}_{\mathbf{X}},U^{(i)}_{\mathbf{X}}]$. To ensure ${P}(\mathbf{X})\in[0,1]$ on replacement of the terms in the inner product,  $P^{(2,i,j)}(\mathbf{X})$, we would require new bounds $\tilde{\Phi}\leq\tilde{L}^{(i)}_{\mathbf{X}}\leq L^{(i)}_{\mathbf{X}}$ and $\tilde{U}^{(i)}_{\mathbf{X}}\geq U^{(i)}_{\mathbf{X}}$ such that $\forall A\sim\mathcal{A}$ we have $\tilde{\phi}_A(\mathbf{X}_{[\tau_{i-1},\tau_i]})\in[\tilde{L}^{(i)}_{\mathbf{X}},\tilde{U}^{(i)}_{\mathbf{X}}]$. Note under this modification we would have the same number of layers, $n_R$, but the number of events in each layer would increase.  The resulting modified estimator, ${P}(\mathbf{X})$, is as follows
\begin{align} \label{eq:Px2}
  {P}(\mathbf{X})  
  & = \prod^{n_R}_{i=1} \Bigg[
  \underbrace{\exp\left\{(\tilde{\Phi} - \tilde{L}^{(i)}_{\mathbf{X}})\cdot[(\tau_i\wedge T)\!-\!\tau_{i-1}]\right\}}_{=:P^{(1)}(\mathbf{X}) \in [0,1]}  \nonumber\\
  & \quad\quad\quad\quad\quad\quad \cdot 
 \mathbbm{E}_{\tilde{\mathbbm{U}}^{(i)}_R, \mathcal{A}}\Bigg[\prod^{\tilde{\kappa}_i}_{j=1}\underbrace{\Bigg(\dfrac{\tilde{U}^{(i)}_{\mathbf{X}}\!-\!\tilde{\phi}_{A_j}\big(\mathbf{X}_{\xi^{(i)}_j}\big)}
 {\tilde{U}^{(i)}_{\mathbf{X}}\!-\!\tilde{L}^{(i)}_{\mathbf{X}}}\Bigg)}_{=:P^{(2,i,j)}(\mathbf{X}) \in [0,1]}\,\Bigg|\,\mathbf{X}\Bigg]
  \Bigg],                                               
\end{align}
The expectation in (\ref{eq:Px2}) is now also with respect to the auxiliary random variables, and we need to use independent copies of $A$ for each time we calculate the estimator $\tilde{\phi}_A$. Proving this identity is straight-forward, as we just take expectations with respect to the auxiliary random variables, and use the unbiasedness property of $\tilde{\phi}_A$, to get back to (\ref{eq:Px}) (albeit with a different choice of lower and upper bound). One should note that (\ref{eq:Px}) holds for all valid lower and upper bounds, including for those chosen in (\ref{eq:Px2}). As such, one could study the direct effect of the unbiased estimator replacement in (\ref{eq:phiunbiased}) by simply contrasting (\ref{eq:Px}) and (\ref{eq:Px2}) using the same bounds ($\tilde{L}^{(i)}_{\mathbf{X}},\tilde{U}^{(i)}_{\mathbf{X}}$).

\subsection{Constructing a scalable replacement estimator} \label{s:repest}
Consider the estimator we developed in (\ref{eq:Px2}). If we construct our auxiliary random variable $A\sim\mathcal{A}$ and estimator $\tilde{\phi}_A(\cdot)$ with $\mathcal{O}(1)$ computational complexity, 
then, if for a given layer we define $\tilde{\lambda}:=\tilde{U}^{(i)}_{\mathbf{X}}-\tilde{L}^{(i)}_{\mathbf{X}}$ (which will determine the intensity of events in the $i^{\text{th}}$ layer), 
the computational cost of the modified algorithm is simply proportional to $\tilde{\lambda}$. Note that the intensity of events will increase relative to 
the estimator in (\ref{eq:Px}) as $\lambda:=U^{(i)}_{\mathbf{X}}-L^{(i)}_{\mathbf{X}}\leq \tilde{\lambda}$. Therefore there is a trade-off between the computational loss in this regard, with that saved from employing the unbiased estimator $\tilde{\phi}_A(\cdot)$.\\
\\
Recalling the construction of KBM in \secref{s:emcfd}, 
\begin{align}
\phi(\mathbf{x}) := (||\mathbf{\Lambda}^{1/2} \nabla\log\pi(\mathbf{x})||^2 + \text{div\,}\mathbf{\Lambda}^{1/2} \nabla\log\pi(\mathbf{x}))/2-\Phi. \label{eq:targetPP}
\end{align}
Developing an unbiased estimator with $\mathcal{O}(1)$ complexity to determine acceptance is trivial. For instance considering a single component of 
(\ref{eq:targetPP}), $\mathbf{\Lambda}^{1/2} \nabla\log\pi(\mathbf{x})$, then (as outlined in \secref{s:intro} and as motivated by \cite{ICML:WT11}) the
na\"{i}ve estimator with complexity $\mathcal{O}(1)$ is simply $(n+1)\mathbf{\Lambda}^{1/2} \nabla\log f_I(\mathbf{x})$, where $I\sim \text{U}\{0, \ldots ,n\}$. Unfortunately our na\"{i}ve estimator of $\mathbf{\Lambda}^{1/2} \nabla\log\pi(\mathbf{x})$ has a bound which grows at least as $\mathcal{O}(n^{1/2})$. Therefore the computational speed-up of employing an $\mathcal{O}(1)$ unbiased subsampling estimator of (\ref{eq:targetPP}) would be counteracted by having a $\tilde{\lambda}$ (and consequently an increase in the intensity of the events required in (\ref{eq:Px2})) which grows as $\mathcal{O}(n)$. 
By comparison, the bound achieved when not employing the sub-sampling estimator, $\lambda$, will be $\mathcal{O}(n)$ times smaller than $\tilde{\lambda}$.\\
\\
The argument we develop in the remainder of this section is that by employing a control variate approach for the construction of the unbiased estimator, $\tilde{\phi}_A(\cdot)$, we can control the computational growth of the corresponding $\tilde{\lambda}$ with $n$.\\
\\
We begin by choosing a point close to a mode of the posterior distribution $\pi$, denoted by $\hat{\mathbf{x}}$. In fact for our scaling arguments to hold, we require ${\hat{\mathbf{x}}}$ to be within $\mathcal{O}(n^{-1/2})$ of the true model, and achieving this is a less demanding task than actually locating the mode. Moreover we note that this operation is only required to be done once, and not at each iteration.\\
\\ 
Letting $\mathcal{A}$ be the law of $I\sim \text{U}\{0, \ldots ,n\}$ we have
\begin{align}
\mathbbm{E}_{\mathcal{A}}\big[\underbrace{(n+1)\cdot\left[\nabla\log f_I(\mathbf{x})-\nabla\log f_I(\hat{\mathbf{x}})\right]}_{=:\tilde{\alpha}_I(\mathbf{x})}\big] = \underbrace{\nabla\log\pi(\mathbf{x})-\nabla\log\pi(\hat{\mathbf{x}})}_{=:\alpha(\mathbf{x})}. \label{eq:alphatilde}
\end{align}
As such, we can re-express $\phi(\mathbf{x})$ as
\begin{align}
\phi(\mathbf{x}) 
& = (\alpha(\mathbf{x})^{T}\mathbf{\Lambda}\left(2\nabla\log\pi(\hat{\mathbf{x}})+ \alpha(\mathbf{x})\right) + \mathbf{\Lambda}^{1/2}\text{div\,}\alpha(\mathbf{x}))/2 + C,
\end{align}
where $C := (||\mathbf{\Lambda}^{1/2} \nabla\log\pi(\hat{\mathbf{x}})||^2 +  \text{div\,}\mathbf{\Lambda}^{1/2} \nabla\log\pi(\hat{\mathbf{x}}))/2 -\Phi$ is a constant. Letting $\mathcal{A}$ now be the law of $I,J\overset{\text{iid}}{\sim} \text{U}\{0, \ldots ,n\}$ we can construct the following unbiased estimator of $\phi$,
\begin{align}
\mathbbm{E}_{\mathcal{A}}\big[\underbrace{(\tilde{\alpha}_I(\mathbf{x}))^{T}\mathbf{\Lambda}\left(2\nabla\log\pi(\hat{\mathbf{x}})+ \tilde{\alpha}_J(\mathbf{x})\right) + 
\mathbf{\Lambda}^{1/2}\text{div\,}\tilde{\alpha}_I(\mathbf{x}))/2 + C}_{=:\tilde{\phi}_A(\mathbf{x}) }\big] = \phi(\mathbf{x}). \label{eq:subsamplelambda}
\end{align}
The construction of our estimator requires evaluation of the constants $\nabla\log\pi(\hat{\mathbf{x}})$ and $\text{div\,}\nabla\log\pi(\hat{\mathbf{x}})$. Although both are be $\mathcal{O}(n)$ evaluations they only have to be computed once, and furthermore their calculation can be done entirely in parallel. For a given target density we would then proceed to compute a suitable $\tilde{\lambda}$ for our sub-sampler. We return at the end of this section to consider finding an appropriate $\tilde{\lambda}$ and its growth with $n$.\\
\\
Embedding our sub-sampling estimator described above within the QSMC algorithm of \secref{s:qsmc}, we arrive at \algref{alg:scale} which we call the \textit{Scalable Langevin Exact algorithm (ScaLE)}. A similar modification could be made to the rejection-sampling version, R-QSMC, which we discussed in \secref{s:qsmc} and detailed in \apxref{apx:rqsmc}. We term this variant \textit{Rejection Scalable Langevin Exact algorithm (R-ScaLE)} and provide full algorithmic details in \apxref{apx:r-scale}. 
\begin{algorithm}[h]
	\caption{The ScaLE Algorithm (as per \algref{alg:qsmc} unless stated otherwise).} \label{alg:subscale3} \label{alg:scale}
	\begin{enumerate}
    \item[0.] Choose $\hat{\mathbf{x}}$ and compute $\nabla\log\pi(\hat{\mathbf{x}})$, $\text{div\,}\nabla\log\pi(\hat{\mathbf{x}})$.
        \item[\ref{alg:qsmc:scalesim1}.] On calling \algref{alg:is-kbm},
    \begin{enumerate}
    \item Replace $L^{(i)}_{\mathbf{X}},U^{(i)}_{\mathbf{X}}$ in \stepref{alg:is-kbm:layer} with $\tilde{L}^{(i)}_{\mathbf{X}},\tilde{U}^{(i)}_{\mathbf{X}}$.
    \item Replace \stepref{alg:is-kbm:w} with: $\tau_i$: If $\xi_j=\tau_i$, set $i=i+1$, and return to \stepref{alg:is-kbm:layer}. Else simulate $A_j=(I_j,J_j)$, with $I_j,J_j\overset{\text{iid}}{\sim} \text{U}\{0, \ldots ,n\}$, and set $w^*_{\xi_j}=w^*_{\xi_j}\cdot (\tilde{U}^{(i)}_\mathbf{X}-\tilde{\phi}_{A_j}(\mathbf{X}_{\xi_j}) )/ (\tilde{U}^{(i)}_\mathbf{X} -  \tilde{L}^{(i)}_\mathbf{X})$ (where $\tilde{\phi}_{A_j}$ is defined as in (\ref{eq:subsamplelambda})) and return to \algstref{alg:is-kbm}{alg:is-kbm:loop}.
    \end{enumerate}
	\end{enumerate}
\end{algorithm}
\\\\ \noindent As previously motivated, given our auxiliary random variable $A\sim\mathcal{A}$ and estimator $\tilde{\phi}_A(\cdot)$ both have $\mathcal{O}(1)$ computational complexity, the iterative computational complexity and scalability of the ScaLE algorithm is determined entirely by $\tilde{\lambda}$. In particular, we will analyse the computational cost of ScaLE by studying the growth of the bound on (\ref{eq:subsamplelambda}) as a function of $n$.\\
\\
For simplicity, the argument we present assumes that we have a bound on {the contribution to} the second-derivative of the log-likelihood {for each datum}. That is we assume {$\forall I\in\{0,\ldots{},n\}$}
\begin{align}
 |\nabla\nabla^T\log f_I(\mathbf{x})|\leq K_{{n}}, \label{eq:hessiancon}
\end{align}
for some $K_{n}>0$.  Frequently this bound will hold globally. For example, considering densities with Gaussian and heavier tails, then the Hessian of the log-likelihood is typically uniformly bounded (in both data and parameter). However, as the argument we present below illustrates, we {actually} only require (\ref{eq:hessiancon}) to hold locally within a bounded region. {We return at the end of this section to discuss the dependence of $K_{n}$ on the size of the data set, $n$.}\\
\\
Recall from \secref{s:nsea} that the practitioner is able to select an appropriate diagonal pre-conditioning matrix, {which we write in the form $\mathbf{\Lambda} = n^{-1}\mathbf{\Lambda}_0$, where $\mathbf{\Lambda}_0$ approximately homogenises the scale of the components.
 We assume such a scaling of the preconditioning matrix as {the size of the movements of our killed Brownian motion} over a fixed time-interval {will then be} of $O(n^{-1/2})$, and hence contract at the same rate as the posterior, at least in the regular inference case.\\
\\
Recalling the layer construction of \secref{s:langevin} {for a single trajectory of killed Brownian motion,}  we can ensure that over any finite time interval we have $\mathbf{x} \in R_{\mathbf{X}}:=\mathbf{\Lambda} ^{1/2}  [\mathbf{\ell_{X}},\mathbf{\upsilon_{X}}]$  {(where 
$[\mathbf{\ell_{X}},\mathbf{\upsilon_{X}}]$ denotes the hypercuboid with given corners)}. We shall set $\mathbf{x}^*= \mathbf{\Lambda} ^{1/2}  (\mathbf{\ell_{X}}+\mathbf{\upsilon_{X}})/2$.  \\
\\
The practitioner has complete freedom to choose $R_{\mathbf{X}}$, and following
the guidance given in (\ref{eq:preconchoice}) of \apxref{ss:fptfinal})}, it makes sense that
$\mathbf{\ell_{X}}$ and $\mathbf{\upsilon_{X}}$ be $\mathcal{O}(1)$, which in turn ensures that $||\mathbf{x}-\mathbf{x}^*||< Cn^{-1/2}$ for some $C>0$ and for all $\mathbf{x} \in R_{\mathbf{X}}$.\\
\\
As we are interested in the {magnitude} of $\tilde{\lambda}$ for $\mathbf{x} \in R_{\mathbf{X}}$, this can be set to be twice  the maximum of $|\tilde{\phi}_A(\mathbf{x})|$ for $\mathbf{x} \in R_{\mathbf{X}}$, over all possible realisations of {$A$}.
To bound $|\tilde{\phi}_A(\mathbf{x})|$, we begin by first considering our elementary estimator in (\ref{eq:alphatilde}). {By imposing our condition in (\ref{eq:hessiancon}) we have}
\begin{align}
{{\max_{\mathbf{x} \in R_{\mathbf{X}},I\in\{0,\ldots{},n\}}}}|\tilde{\alpha}_I(\mathbf{x})| \leq (n+1)\cdot K_{{n}} \cdot {\max_{\mathbf{x} \in R_{\mathbf{X}}}}||\mathbf{x}-\hat{\mathbf{x}}||.
\end{align}
We can proceed to find a bound for our estimator in (\ref{eq:subsamplelambda}) as follows
\begin{align}
 \max_{\mathbf{x} \in R_{\mathbf{X}}{,A\in\mathcal{A}}}|\tilde{\phi}_A(\mathbf{x})|
& \leq K_{{n}}\cdot  n^{-1}\cdot (n+1)\cdot \max_{\mathbf{x} \in R_{\mathbf{X}}} ||\mathbf{x}-\hat{\mathbf{x}}||\cdot ||\mathbf{\Lambda}_0|| \nonumber\\
& \quad\quad\quad \cdot (2\nabla\log\pi(\hat{\mathbf{x}}) +(n+1)\cdot K_{{n}} \cdot \max_{\mathbf{x} \in R_{\mathbf{X}}} ||\mathbf{x}-\hat{\mathbf{x}}||)/2 \nonumber\\
&  \quad + K_{{n}}\cdot  n^{-1/2}\cdot ||\mathbf{\Lambda}_0^{1/2}||\cdot \text{div\,} ((n+1) \cdot \max_{\mathbf{x} \in R_{\mathbf{X}}}||\mathbf{x}-\hat{\mathbf{x}}||)/2 + C.
\end{align}
Now, as we have $\max_{\mathbf{x} \in R_{\mathbf{X}}}||\mathbf{x}-\hat{\mathbf{x}}|| \leq ||\mathbf{x}^*-\hat{\mathbf{x}}||+Cn^{-1/2}$, then
\begin{align}
 \tilde{\lambda} 
 & = \mathcal{O}\Big(K_{{n}} \cdot \big(||\mathbf{x}^*-\hat{\mathbf{x}}||+Cn^{-1/2}\big) \cdot \big(1+n K_{{n}} \cdot  \big(||\mathbf{x}^*-\hat{\mathbf{x}}||+Cn^{-1/2}\big)\big)  \nonumber\\
 & \quad\quad\quad\quad  + n^{1/2} \cdot K_{{n}} \cdot \big(||\mathbf{x}^*-\hat{\mathbf{x}}||+Cn^{-1/2}\big)\Big).
\end{align}
At stationarity $\mathbf{x}^*$ will be a draw from the support of the posterior. So, if we assume posterior contraction at the regular $\sqrt{n}$ rate, then $||\mathbf{x}^*-\hat{\mathbf{x}}||=\mathcal{O}_p(n^{-1/2})$, and hence $\tilde{\lambda}= \mathcal{O}_p(1)$.\\
\\
Now, returning to consider typical values of {$K_{n}$}, we first note that for smooth densities with Gaussian and heavier tails, the Hessian of the log-likelihood is typically uniformly bounded (in both data and parameter). As such it is automatic that {$K_{n}$} should be constant in this case and so the ScaLE algorithm has $\mathcal{O}(1)$ iterative computational complexity.\\
\\
For light-tailed likelihoods, the size of {$K_{n}$} will likely depend upon the most ``extreme'' data. Since maxima of heavy-tailed distributions are typically $\mathcal{O}((\log n)^{\eta })$ for some constant $\eta \le 1/2$, this can be expected to affect {$K_{n}$} making it grow logarithmically as a function of $n$. We shall omit detailed calculations here. This discussion does however suggest that non-linear transformations could potentially be used to completely eliminate this problem should it be considered important.\\
\\
{\bf Example}:\ 
As a toy example to illustrate the effect of subsampling and control variates, we carried out an analysis of the Bayesian posterior for a $t_5$ model
\begin{align}
\pi (\theta | {\bf y}) \propto \prod_{i=1}^n  \left(\dfrac{1}{5 + (y_i - \theta )^2}\right)^{3},
\end{align}
with data generated from the same model with $\theta =0$. \figref{fig:compcost} empirically illustrates the scaling properties of the approach, showing that the for this example the computational cost per iteration remains more or less constant as a function of $n$, as predicted by the above arguments, so that there is essentially no iterative cost associated with $n$.
\begin{figure}[h] 
\centering
\includegraphics[width=0.95\textwidth]{compcost.pdf}
\caption{Iterative computational cost of ScaLE (black line) when applied to increasing volumes of data generated from a one dimensional $t_5$ distribution, with algorithmic parameters chosen to give a similar overall Monte Carlo error for each data-set size. Blue line depicts linear computational cost, vertical hatching sub-linear computational cost, and horizontal hatching super-linear computational cost.}\label{fig:compcost}
\end{figure}


\subsection{Theoretical Properties} \label{sec:smctheory}
SMC algorithms in both discrete and continuous time have been studied extensively in the literature. In order to avoid a lengthy technical diversion, we restrict ourselves here to studying a slightly simplified version of the problem in order to obtain the simplest and most interpretable possible form of results. We defer the technical details of this construction to \apxref{app:smc} and give here only a qualitative description intended to guide intuition and the key result: that the resulting estimator satisfies a Gaussian central limit theorem with the usual Monte Carlo rate.\\
\\
We consider a variant of the algorithm in which (multinomial) resampling occurs at times $kh$ for $k\in \mathbb{N}$ where $h$ is a timestep specified in advance and consider the behaviour of estimates obtained at these times. Extension to resampling at a random subset of these resampling times would be possible using the approach of \cite{smc:theory:DDJ12}, considering precisely the QSMC algorithm presented in \algref{alg:qsmc} and the ScaLE algorithm in \algref{alg:scale} would require additional technical work somewhat beyond the scope of this paper; we did not observe any substantial difference in behaviour.\\
\\
In order to employ standard results for SMC algorithms it is convenient to consider a discrete time embedding of the algorithms described. We consider an abstract formalism in which between the specified resampling times the trajectory of the Brownian motion is sampled, together with such auxiliary random variables as are required in any particular variant of the algorithm. Provided the potential function employed to weight each particle prior to resampling has conditional expectation (given the path) proportional to the exact killing rate integrated over these discrete time intervals we recover a valid version of the ScaLE algorithm.\\
\\
This discrete time formalism allows for results on more standard SMC algorithms to be applied directly to the ScaLE Algorithm. We provide  in the following proposition a straightforward corollary to \cite[Chapter 9]{smc:theory:Del04}, which demonstrates that estimates obtained from a single algorithmic time slice of the ScaLE algorithm satisfy a central limit theorem. 
\begin{proposition}[Central Limit Theorem]
  In the context described, under mild regularity conditions (see references given in \apxref{app:smc}):
\begin{align*}
\lim_{N \to \infty} \sqrt{N} \left[\frac{1}{N} \sum_{i=1}^N \varphi(X^i_{hk}) -
  \mathbb{E}_{\mathbb{K}_{hk}^x}\left[\varphi(X^i_{hk})\right] \right] \Rightarrow \sigma_k(\varphi) Z
\end{align*}
where, $\varphi:\mathbb{R}^d \to\mathbb{R}$,  $Z$ is a standard normal random variable, $\Rightarrow$ denotes
convergence in distribution, and $\sigma_k(\varphi)$ depends upon the precise choice of subsampling scheme as well as the test function of interest and is specified in \apxref{app:smc}. 
  \end{proposition}


\section{Examples} \label{s:examp}
In this section we present two applications of the methodology developed in this paper. In particular, in \secref{s:examp:1} for exposition we consider the application of ScaLE to a simple example of logistic regression with a moderately large data set. In \secref{s:examp:2} we consider a problem in which we wish do parameter inference for a more complicated mixture model, motivated by big data application and with $n=2^{27}$.


\subsection{Example 1: Menarche} \label{s:examp:1}
In this subsection we consider the R \cite{software:R} dataset `Menarche', due originally to \cite{MS66}, which appears in the `MASS' package \cite{software:MASS}. Menarche is a dataset comprising the age at which $n=3918$ adolescent female children reach menarche. This dataset is frequently used as a simple textbook example to which we can fit a single feature logistic regression model. We will contrast the likelihood-based estimates provided by the R glm function for this dataset, with those provided by the ScaLE algorithm (\algref{alg:scale}) we developed in \secref{s:repest}. More precisely, we have binary observations of and fit a model of the form,
\begin{equation}
y_i = \left\lbrace \begin{array}{ll} 
1 & \quad\text{with probability } \dfrac{\exp\{\mathbf{x}^T_i \mathbf{\beta}\}}{1+\exp\{\mathbf{x}^T_i \mathbf{\beta}\}},
\\*[10pt] 0 & \quad\text{otherwise}, \label{eq:logistic} \end{array}  \right.
\end{equation}
noting that we can readily compute for each datum its contribution to the likelihood along with its derivatives. We normalise the design matrix (subtracting the sample mean and dividing by standard deviation); the glm [R] package estimates $\hat{\mathbf{\beta}} = [1.410,4.658]^T$. We show the glm estimate as the blue line in \figref{fig:logfit}.
\begin{figure}[h] 
\centering
\includegraphics[width=0.95\textwidth]{logistic_fit.pdf}
\caption{Contrast of glm (blue) and ScaLE (red) fit to the MASS Menarche dataset.}\label{fig:logfit}
\end{figure}\\
\noindent In contrast, we initialise the ScaLE algorithm at a point one standard deviation from $\hat{\beta}$, with a particle set of size $N=100$. In total we run ScaLE for diffusion time of $T=1000$, and approximate the law of the killed process as in (\ref{eq:occupation2}) using observations of the trajectory at a resolution of $t_i-t_{i-1}=0.1$. Each particle trajectory at each time $t\in[0,T]$ is associated with a sub-sample of the full dataset of size $2$. The contrast with glm of the fit to the dataset is shown as the red line in \figref{fig:logfit}. An example of a typical run can be found in \figref{fig:trace}. In \figref{fig:loglike} we compute the profile of the log-likelihood along with the estimates of glm and ScaLE, contrasting the profile with a kernel density estimate computed using the locations of the weighted trajectories at the terminal time slice (as in (\ref{eq:is weight2})).
\begin{figure}[h] 
\centering
\includegraphics[width=0.95\textwidth]{logistic_trace.pdf}
\caption{Particle trajectories of ScaLE applied to the Menarche dataset.}\label{fig:trace}
\end{figure}
\begin{figure}[h] 
\centering
\includegraphics[width=0.95\textwidth]{logistic_likelihood.pdf}
\caption{Log-likelihood surface and ScaLE estimate for the Menarche dataset.}\label{fig:loglike}
\end{figure}



\subsection{Example 2: Contaminated Mixture} \label{s:examp:2}
In this subsection we consider the application of ScaLE to conduct parameter inference for a contaminated mixture model. This is motivated by big data sets obtained from internet applications in which the dataset obtained is corrupted with noisy observations. In particular, we the case in which each datum comprises two features and we consider fitting a model in which the likelihood of an individual observation ($y_i$) is as follows,
\begin{align}
F_i := \dfrac{1-p}{\sqrt{2\pi\sigma^2}}\exp\left\{-\dfrac{1}{2\sigma^2}\left(\alpha\cdot x_{i,1} + \beta \cdot x_{i,2} -y_i\right)^2\right\} + \dfrac{p}{\sqrt{2\pi\phi^2}}\exp\left\{-\dfrac{1}{2\phi^2}y_i^2\right\}. \label{eq:contamlike}
\end{align}
In this model $p$ represents the level of corruption and $\phi$ the variance of the corruption. A common approach if we wished to conduct parameter inference for this model using MCMC is data augmentation \cite{tannerwong87}. However, for large data sets this is not feasible as the dimensional of the auxiliary variable vector will be $\mathcal{O}(n)$.
\\
\\
We generate a dataset of size $n=2^{27}$ from the model with parameters $\mathbf{\mu}=[\alpha,\beta,\sigma,\phi,p]=[2,5,1,10,0.05]$, and conduct parameter inference using ScaLE initialised at a point several deviations from that used to generate the data. For this application we use a particle set of size $N=2048$, run ScaLE for diffusion time of $T=200$, and approximate the law of the killed process as in (\ref{eq:occupation2}) using observations of the trajectory at a resolution of $t_i-t_{i-1}=0.1$. Each particle trajectory at each time $t\in[0,T]$ is associated with a sub-sample of the full dataset of size $32$. An example of a typical run can be found in \figref{fig:contamtrace}, and in \figref{fig:contammarginal} estimate the marginal density of the parameters using the occupation measure of the trajectories in the interval $t\in[100,200]$.
\begin{figure}[h] 
\centering
\includegraphics[width=0.9\textwidth]{contam_trace.png}
\caption{Particle trajectories of ScaLE applied to the contaminated mixture dataset. Red lines mark the parameter values used to initialise the algorithm, blue lines mark the mode of the marginal density.}\label{fig:contamtrace}
\end{figure}
\begin{figure}[h] 
\centering
\includegraphics[width=0.9\textwidth]{contam_marg.pdf}
\caption{Marginal densities obtained by ScaLE for the contaminated mixture dataset. Red lines mark the parameter values for which the data set was generated.}\label{fig:contammarginal}
\end{figure}

\FloatBarrier
\section{Conclusions} \label{s:conclusions}
In this paper we have introduced a new class of \textit{Quasi-Stationary Monte Carlo (QSMC)} methods which are genuinely continuous-time algorithms for simulating from complex target distributions. We have emphasised its particular effectiveness in the context of {\em big data} by developing novel sub-sampling approaches and the \textit{Scalable Langevin Exact (ScaLE)} algorithm. Unlike its immediate competitors, our subsampling approach within ScaLE is essentially computationally free and does not result in any approximation to the target distribution. Our methodology is embedded within an SMC framework, which we use to provide theoretical results. In addition we give examples to which we apply ScaLE, demonstrating its robust scaling properties for large data sets.\\
\\
We see this paper as a first step in what we hope will be a fruitful new direction for Computational Statistics. Many ideas for variations and extensions our implementation exist and will stimulate further investigation.\\
\\
Firstly, the need to simulate a quasi-stationary distribution creates particular challenges. Although quasi-stationarity is underpinned by an elegant mathematical theory, the development of numerical methods for quasi-stationarity is understudied.  We have presented an SMC methodology for this problem, but alternatives exist. For instance, \cite{blanchet2015theoretical} suggest alternative approaches.\\
\\
Even within an SMC framework for extracting the quasi-stationary distribution, there are interesting alternatives we have not explored. For example, by a modification of our re-weighting mechanism we can relate the target distribution of interest to the limiting {\em smoothing} distribution of the process, as opposed to the filtering distribution as we do here. Within the quasi-stationary literature this is often termed the type II quasi-stationary distribution. As such, the rich SMC literature offers many other variations on the procedures we adopt here.\\
\\
We have also concentrated solely on killed (or re-weighted) Brownian motion in our paper. We have shown that this strategy has robust convergence properties. However, given existing methodology for the exact simulation of diffusions in \cite{AAP:BR05,MCAP:BPR08,B:PJR15}, there is scope to develop methods which use proposal measures which much better {\em mimick} the shape of the posterior distribution.\\
\\
Our sub-sampling and control variate approaches offer dramatic computational savings. However there may be scope to extend these ideas still further. For instance, it might be possible to sub-sample dimensions and thus reduce the dimensional complexity for implementing each iteration.\\
\\
We conclude by noting that as a byproduct, the theory behind our methodology offers new insights into problems concerning the existence of quasi-stationary distributions for diffusions killed according to a state-dependent hazard rate, complementing and extending current state-of-the-art literature \cite{steinsaltz2007quasistationary}.

\section*{Acknowledgments}
MP was supported by the EPSRC [grant number EP/K014463/1]. PF was supported by [grant numbers EP/N031938/1, EP/K014463/1]. GOR was supported by [grant numbers EP/K034154/1, EP/K014463/1, EP/D002060/1].

\FloatBarrier
\bibliographystyle{imsart-number}
