Sequential Monte Carlo samplers

This chapter develops sequential Monte Carlo (SMC) samplers as a population-based approach to inference in Bayesian deep learning.^[1] Building on the sequential Monte Carlo methods of Chapter Introduction to sampling, we focus on the two extensions that make them practical for neural network posteriors: multilevel SMC and scalable parallel SMC.

Motivation¶

Chapter Introduction to sampling introduced SMC as a particle filter Doucet et al., 2001Doucet & Johansen, 2011: a latent state evolving in time, tracked by a population of particles over a sequence of growing spaces. This chapter applies the same population-based approximation to a static target—a Bayesian posterior on a fixed parameter space. Without any dynamics to generate them, the intermediate distributions are now constructed by design, and the resulting methods are known as SMC samplers Del Moral et al., 2006.

It is useful to place the neighbouring methods relative to this sampler. Tempering is the mechanism inside the sampler that builds the intermediate distributions, annealing the likelihood from the prior to the posterior Geyer, 1991Gelman & Meng, 1998Neal, 2001 (data-tempering Chopin, 2002 is a common alternative). Annealed importance sampling (AIS) is the reduced case in which resampling is omitted Neal, 2001, so it is best viewed as a simplified SMC sampler rather than a distinct method. The chapter’s two main contributions are developments of the SMC sampler, along two different axes: multilevel SMC (MLSMC) improves statistical efficiency by combining a hierarchy of coarse-to-fine approximations in a telescoping estimator Beskos et al., 2017Beskos et al., 2018, while parallel SMC improves scalability by running many samplers concurrently and weighting them by their normalizing-constant estimates, avoiding the heavy communication of island particle methods Vergé et al., 2015Whiteley et al., 2016Liang et al., 2025.

These methods matter for Bayesian deep learning because neural network posteriors are high-dimensional, multi-modal with many symmetric modes, and costly to evaluate—the regime where importance sampling degenerates and single-chain MCMC mixes slowly. An SMC sampler anneals an entire population towards the posterior, provides an unbiased estimate of the model evidence, and parallelises naturally. The remainder of the chapter develops the two extensions that make this practical at scale: MLSMC with trace-class neural network priors Sell & Singh, 2023Chada et al., 2025, and scalable parallel SMC Liang et al., 2025.

Introduction¶

Monte Carlo based Bayesian algorithms of the type described in Chapter Introduction to sampling typically require generating many samples from a target distribution which at best can be evaluated up to a normalizing constant. This can be expensive and cumbersome in high-dimensions, and so methods to improve complexity are paramount. As discussed in Chapter Introduction to sampling, the best computational complexity achievable by a Monte Carlo algorithm in general is cost $\propto1/$MSE, i.e. achieving a mean square error (MSE) of $\mathcal{O}(\varepsilon^{2})$ for $\varepsilon>0$ will asymptotically require $\mathcal{O}(\varepsilon^{-2})$ samples. This is hence referred to as the canonical rate, as it holds for any intractable integral for which we can generate i.i.d. unbiased samples of the integrand.

The sequential Monte Carlo (SMC) sampler Del Moral et al., 2006Dai et al., 2022Chopin & Papaspiliopoulos, 2020 was developed at the turn of the millennium Jarzynski, 1997Berzuini & Gilks, 2001Gilks & Berzuini, 2001Neal, 2001Chopin, 2002, and is well suited to this high-dimensional regime. Forcing the (successive) importance distributions to be close provides a generic mechanism to overcome the “curse-of-dimensionality”, i.e. complexity* constant* scaling like $e^P$ for dimension $P$ Chatterjee & Diaconis, 2018Agapiou et al., 2017. Under suitable high-dimensional stability assumptions, SMC samplers can avoid the exponential degradation of direct importance sampling and may achieve polynomial cost scaling in dimension Beskos et al., 2014. The SMC sampler also handles bad initialization or multi-modality on par with other population methods Lee et al., 2010, and facilitates adaptive tuning of the Markov kernel Buchholz et al., 2021, with comparable efficiency to MCMC approaches like the No-U-Turn Sampler (NUTS) Hoffman & Gelman, 2014. These benefits are not shared by the particle filter, which does not utilize MCMC methodology. It also delivers an* unbiased* estimator of the normalizing constant, or model evidence, which can be useful in practice Chopin & Papaspiliopoulos, 2020.

The simplest method for improving complexity is parallelization. Standard i.i.d. Monte Carlo simulation is “embarrassingly parallel”, as the samples can all be simulated independently on different machines. If the cost of each simulation dominates the $\mathcal{O}(\varepsilon^{-2})$ communication and summation operations, which should both have extremely small scalar constants which are often invisible to the algorithm ($<10^{-9}$ smaller than likelihood computations), then the constant and hence final time complexity can be reduced significantly in practice, even down to effectively $\mathcal{O}(1)$ when the merge is invisible. Monte Carlo methods for Bayesian simulation do not typically deliver i.i.d. samples, and so they are not a priori parallel. SMC samplers admit a certain degree of parallelism in their raw form, and we will show that they can also be provably “embarrassingly” parallelized beyond this Liang et al., 2025.

If the underlying problem we are simulating from additionally requires a discrete approximation of a function then the complexity compounds. For example, suppose that achieving a bias$^2$ of $\mathcal{O}(\varepsilon^2)$ demands a cost of $\mathcal{O}(\varepsilon^{-\xi})$. Then the MSE complexity compounds as $\mathcal{O}(\varepsilon^{-2-\xi})$. To overcome this issue, one remedy to reduce the computational cost is through a modified Monte Carlo estimator which is known as multilevel Monte Carlo (MLMC) Heinrich, 2001Giles, 2008Giles, 2015. MLMC introduces a hierarchy of intermediate convergent discretization levels and constructs an estimator as a telescopic sum of a coarse and cheap estimator with many samples and successive increment estimators with increasing cost and decreasing sample size, ultimately improving the complexity to $\mathcal{O}(\varepsilon^{-(2\wedge \xi)})$ instead of the product form, which can translate to huge scaling gains for small target MSE.

MLMC methodology in the Bayesian context was first developed for Markov chain Monte Carlo (MCMC) Hoang et al., 2013Dodwell et al., 2015 and sequential Monte Carlo (SMC) samplers Beskos et al., 2017. We consider a multilevel sequential Monte Carlo (MLSMC) sampler based on the application of a function-space version of the latter algorithm Beskos et al., 2018 to a strongly convergent class of neural network models known as trace-class neural networks (TNN) Sell & Singh, 2023Chada et al., 2025. These are Bayesian feedforward neural networks (BNN) which impose a* width decay* on the weights instead of a standard $L^2$ width decay prior. This simple and natural idea eliminates exchangeability and non-identifiability of the weights and delivers strong convergence with a tunable smoothness parameter. The limiting Gaussian measures are defined in function spaces with similar properties to Sobolev spaces Bogachev, 2007. The resulting method was originally developed in Chada et al. (2025), where the authors were able to demonstrate the canonical complexity is achievable. We will use function evaluations as the fundamental unit of computational cost, which is also a convenient proxy for wall-clock time which is agnostic to particulars of the implementation or hardware.

The organization of this chapter is as follows. First, we introduce the SMC sampler in Section SMC sampler. In Section MLSMC for TNN we present the MLSMC sampler method for TNN. In Section Parallel SMC we present the parallel SMC method.

SMC sampler¶

Suppose that we aim to simulate from a target distribution

One option is to use an MCMC kernel $\mathcal{M}$ such that $p \mathcal{M} = p$ to simulate $\theta^{(k)} \sim \mathcal{M}(\theta^{(k-1)}, \cdot)$, but this method produces correlated local samples and is prone to getting stuck in a given basin of attraction. If that basin of attraction is in the set of indistinguishable basins that we want then it may not particularly matter in practice. Gradient-based MCMC methods such as HMC can reduce random-walk behaviour, but they can still mix poorly across isolated posterior modes or permutation-symmetric basins.

Population methods such as importance sampling (IS) are able to mitigate this issue, but they suffer from their own difficulties. In particular, in its raw form, IS suffers from a curse of dimensionality in both parameter and data Chatterjee & Diaconis, 2018, and as such it is not suitable for high dimensional problems and big data. The SMC sampler bypasses this problem by

• interpolating between an easy-to-sample-from initial distribution $p_0$ and the target distribution of interest $p$ with several intermediate distributions $p_t$ such that $p_T = p$, and

• interleaving sequential importance sampling (SIS) with resampling and mutation by a sequence of appropriate MCMC kernels $\mathcal{M}_t$ such that $p_t \mathcal{M}_t = p_t$.

In the Bayesian context, typically the prior $p$ is a suitable initial distribution. The sequence of intermediate targets can be built to gradually introduce the data, for example either from growing subsets of the data Chopin, 2002 or with a tempering schedule $0=\lambda_0<\lambda_1 < \dots < \lambda_T=1$ Neal, 2001:

Three quantities are tracked along this path. The normalizing constant at temperature $\lambda_t$ is

which at $\lambda_T=1$ equals the model evidence $p(Y \mid X)$. Algorithm 1 maintains an unbiased estimate $Z_t^{M}$ at each step, with $Z_T^{M}$ returned as the estimate of the evidence. Successive targets are linked by the incremental importance weight that reweights particle $i$ from $\pi_{t-1}$ towards $\pi_t$,

How evenly these weights are spread is measured by the effective sample size (ESS),

Values near $M$ indicate balanced weights, whereas small values mean a few particles dominate. The ESS drives the algorithm in two ways: a low value motivates resampling, and in the adaptive variant the next temperature is chosen to preserve a target $\mathrm{ESS}=\alpha M$. Data-tempering and likelihood-tempering construction schemes for the intermediate targets can also be mixed and matched. We will consider another option in Section MLSMC for TNN. See Algorithm 1 for the version with tempering for Bayesian inference. The steps which require communication during execution are denoted with (c). It is noteworthy that if we skip the resampling step then the algorithm is called annealed importance sampling (AIS) and was introduced in Neal (2001). The normalizing constant then only needs to be computed once at the end. This version of the algorithm is convenient because it is embarrassingly parallel – communication is only required at resampling times. However, the purpose of resampling is to prevent degeneracy and improve stability, and this benefit often outweighs the added variance Chopin, 2002Chopin & Papaspiliopoulos, 2020.

MLSMC for TNN¶

In this section we will build up the MLSMC for TNN method. First we introduce the MLMC method, followed by the TNN. Then we present the MLSMC sampler and provide theoretical justification, and numerical results validating and extending the theory.

MLMC method¶

We shall begin with a short review of MLMC. MLMC is useful when accurate simulations are expensive but coarse approximations are cheap. Rather than estimating the finest-level expectation directly, one estimates a coarse expectation plus a sequence of corrections. If consecutive levels are strongly coupled, the correction variances decay rapidly, so fewer samples are needed at expensive fine levels. Let us assume that we are given a probability density $p$, on a state-space $\mathsf{U}$. It is of interest to compute expectations of $p-$integrable functions, $\varphi:\mathsf{U}\rightarrow\mathbb{R}$; $p(\varphi):=\int_{\mathsf{U}}\varphi(u)p(u)du$. Now, we assume that we need to approximate $p$ by a density $p_l$ on a state-space $\mathsf{U}_l\subseteq\mathsf{U}$ such that:

1. $\lim_{l\rightarrow\infty}p_l(\varphi)=p(\varphi)$, for any integrable $\varphi:\mathsf{U}\rightarrow\mathbb{R}$.

2. Computing with $p_l$ grows progressively more expensive as $l$ increases.

Fix a finest level $L\in\{1,2,\dots\}$. The finest-level expectation decomposes into a cheap coarse term plus a sum of level-wise corrections,

where $[p_l-p_{l-1}](\varphi)$ abbreviates $p_l(\varphi)-p_{l-1}(\varphi)$. MLMC is built on exactly this decomposition: estimating the right-hand side term by term turns out to be cheaper than attacking the left-hand side directly. The vehicle for this is a coupling of each consecutive pair: a joint density $\check{p}_l$ on $\mathsf{U}_l\times\mathsf{U}_{l-1}$, for $l\in\{1,\dots,L\}$, whose marginals recover the pair, i.e. $\int_{\mathsf{U}_{l-1}}\check{p}_l(u_l,u_{l-1})du_{l-1}=p_{l}(u_{l})$ and $\int_{\mathsf{U}_l}\check{p}_l(u_l,u_{l-1})du_l=p_{l-1}({u}_{l-1})$.

Equipped with these couplings, the estimator is assembled in two steps.

1. Draw $M_0\in\mathbb{N}$ i.i.d. samples $U_0^1,\dots,U_0^{M_0}$ from the coarsest density $p_0$.

2. For each increment level $l\in\{1,\dots,L\}$, and mutually independently of everything else, draw $M_l\in\mathbb{N}$ i.i.d. pairs
   $(U_l^1,\tilde{U}_{l-1}^1),\dots,(U_l^{M_l},\tilde{U}_{l-1}^{M_{l}})$ from the coupling $\check{p}_l$.

The resulting estimator reads

Using i.i.d. samples $U_L^1,\dots,U_L^M$ from $p_L$ to approximate $p_L(\varphi)$, we have

Neither estimator carries sampling bias, and the mean square error of any unbiased estimator decomposes into its variance plus the squared discretization bias; for (7),

where $\operatorname{Var}$ denotes the variance operator. The same decomposition holds for (8) with an identical bias term (which an optimized allocation balances against the variance), so any multilevel advantage must enter through the variance. For the multilevel estimator this is

while the single-level estimator (8) has

Whenever the couplings force the increment variances $\operatorname{Var}[\varphi(U_l^1)-\varphi(U_{l-1}^1)]$ to decay rapidly in $l$, few samples are required at the expensive fine levels, and the multilevel estimator (7) attains the same order of MSE (in the sense of (9)) as (8) at strictly lower cost. The following theorem quantifies when this occurs.

The cost (13) is asymptotically below that of the single-level estimator (8). To see this, note that (8) requires $\mathcal{O}(\varepsilon^{-2})$ samples, each at per sample cost $\mathcal{O}(\varepsilon^{-\gamma/\alpha})$, for a total of $\mathcal{O}(\varepsilon^{-2-\gamma/\alpha})$, which exceeds the MLMC complexity whenever $\beta>0$. When $\beta > \gamma$ the multilevel cost reaches $\mathcal{O}(\varepsilon^{-2})$, known as the canonical rate: the cost of i.i.d. sampling of a tractable integrand, and hence unimprovable in general.

Trace class neural networks¶

We first fix the setting. The data consist of $N\in\mathbb{N}$ pairs $\mathcal{D} = \left((x_1,y_1),\dots,(x_N,y_N)\right)$, with $x_i \in \mathsf{X}$ and $y_i\in \mathsf{Y}$ for $i\in\{1,\dots,N\}$. The goal is to learn from these a predictive model $f: \mathsf{X} \rightarrow \mathsf{Y}$, and a standard route is a parametric family $f: \mathsf{X}\times\Theta \rightarrow \mathsf{Y}$ with $\Theta\subseteq\mathbb{R}^{P}$. The inputs $x_{1:N}$ are treated as deterministic, i.e. the discriminative rather than generative view of supervised learning.

Regression. If $\mathsf{Y}=\mathbb{R}^C$, the observations are modelled, for $i\in\{1,\dots,N\}$, as

where the noise terms are independent across $i\in\{1,\dots,N\}$ (denoted $\textrm{ind}$) and $\mathcal{N}_C(\mu,\Sigma)$ is the Gaussian distribution on $\mathbb{R}^C$ with mean $\mu$ and covariance $\Sigma$. The corresponding likelihood is

with $\phi_C(y;\mu,\Sigma)$ the density of $\mathcal{N}_C(\mu,\Sigma)$ evaluated at $y$.

Classification. If $\mathsf{Y} = \{1,\dots,C\}$ for some $C\in\mathbb{N}$, for convenience of the present exposition we will define $f:\mathsf{X}\times\Theta\rightarrow\mathbb{R}^C$, with $f(x,\theta)=(f_1(x,\theta),\dots,f_C(x,\theta))$, and then separately the softmax function as

The labels are modelled as $y_i \sim S(x_i,\theta)$, independently over $i\in\{1,\dots,N\}$; here $S(x,\theta)=(S_1(x,\theta),\dots,S_C(x,\theta))$ is read as a categorical distribution over the $C$ classes associated with input $x$. In this case, the likelihood is

The predictive $S$ is Lipschitz, so conclusions about $f$ translate immediately. We conclude this part by stating the key advantages of TNNs: one is able to have a stable infinite limit related to the width, as well as high-dimensional scalability, due to the fact that one does not need to define the covariance operator structure.

Bayesian neural networks¶

Define the element-wise activation function as $\sigma$. Let $\mathsf{X} = \mathbb{R}^D$ and $\mathsf{Y}=\mathbb{R}^C$. A DNN is specified by layer dimensions $(D_0,\dots,D_L)\in\mathbb{N}^{L+1}$, where necessarily $D_0=D$ (input layer) and $D_L=C$ (output layer), together with weight matrices $A_d\in \mathbb{R}^{D_{d} \times D_{d-1}}$ and bias vectors $b_d\in \mathbb{R}^{D_d}$ for $d\in\{1,\dots,L\}$. Collecting the parameters as $\theta := \left((A_1,b_1),\dots,(A_L,b_L)\right)$, so that $\theta\in\Theta=\bigotimes_{d=1}^L\{ \mathbb{R}^{D_d\times D_{d-1}}\times\mathbb{R}^{D_d}\}$, the network output is built by the recursion

with the final-layer output $f(x,\theta)$ defining the DNN. The BNN is given by placing a prior $p$ on $\Theta$.

TNN¶

We now introduce the trace class neural network (TNN) priors, which were first proposed in Sell & Singh (2023), and differ from standard BNN priors in two fundamental ways. The first is that the prior on the weights and biases depends on the rows and columns, as opposed to standard isotropic weight decay. This is referred to as width-decay and it eliminates exchangeability and the associated non-identifiability. The second thing is that nested approximations of the well-defined TNN function-space converge strongly to a non-Gaussian process in the limit of infinite width. This stands in contrast to standard isotropic priors, under which the appropriately rescaled network converges weakly to a Gaussian process Neal, 1996G. Matthews et al., 2018, so that the infinite-width limit is a kernel method and the compositional structure of the prior is lost. This will be made precise below.

These priors were introduced to mimic Gaussian measure priors $p \sim \mathcal{N}(0,\mathcal{C})$ for inverse problems over function-space, for which Gaussian random fields are commonly simulated through the Karhunen-Loève expansion

with $(\lambda_j,\Phi_j)_{j \in \mathbb{Z}^+}$ the eigenpairs of the covariance operator $\mathcal{C}$ and $\{\iota_j\}_{j \in \mathbb{Z}^+}$ Gaussian white noise; see Lord et al. (2014) for the derivation of (19) and its use in stochastic numerics. Priors built from (19) scale poorly with input dimension, however, and this shortcoming was the original motivation for TNN priors, in which the weight and bias variances are collected in a trace-class diagonal covariance operator $\mathcal{C}$. See Sell & Singh, 2023. More precisely, the TNN prior is given by

The tuning parameter $s$ controls how much information one believes concentrates on the first nodes, and as such also controls the smoothness, expressiveness, and flexibility of the prior in terms of the functions it can represent. In the case of $s>1$, we refer to the prior as trace-class, which is the motivation for the name trace-class neural network prior.

A convenient property of these priors is strong convergence, which is made precise with the following proposition. The proof is given in Chada et al. (2025). This is the key result that enables application of the MLSMC algorithm to TNN. The advantage of TNN is that it allows one not to prespecify the covariance operator structure in the prior, and it is also capable of having an infinite-width limit.

Multilevel TNN¶

The remaining architectural freedom is the hidden-layer widths $D_{L'}$, $L'\in\{1,\dots,L-1\}$; the input and output widths are dictated by the problem, $D_0=D$ and $D_L=C$. We take the depth $L$ as fixed and give every hidden layer a common width governed by a resolution parameter $l\in\mathbb{N}$, namely $D_l=2^l$, after which the per-layer width variables $D_{L'}$ are no longer needed. Denote the corresponding vector of parameters by $\theta_l := \left((A_{1}^l,b_{1}^l),\dots,(A_{D}^l,b_{D}^l)\right) \in \Theta_l \subset \Theta$. The limiting infinite width NN output function $f(x,\theta)$ (18) approximated at finite resolution $l$ is denoted by $f_l(x,\theta_l)$, the likelihood by $p_l(Y|\theta_l,X)$, and the posterior distribution by

Both $\pi_l$ and the network $f_l(x,\theta_l)$ should be read as finite-width surrogates for their counterparts under the non-parametric limiting DNN as $l \rightarrow \infty$, whenever that limit is well defined.

MLSMC TNN method¶

This section presents the algorithm for estimating expectations under the posterior $\pi$, in particular the posterior predictive $\mathbb{E}_{\pi}[f(x,\theta)],$ together with the mathematical results that justify it. The guarantees transfer unchanged to any objective $\varphi \circ f$ with $\varphi$ Lipschitz, a class covering the usual performance metrics, negative log likelihood and accuracy among them. Everything is stated for the posterior predictive first, then extended as a corollary. Splitting off the discretization error via $\mathbb{E}_{\pi}[f(x,\theta)] = \mathbb{E}_{\pi_L}[f_L(x,\theta_L)] + (\mathbb{E}_{\pi}[f(x,\theta)]-\mathbb{E}_{\pi_L}[f_L(x,\theta_L)])$, the plan is to bring the MLMC machinery to bear on the telescoping identity

It will be shown that estimating the summands on the R.H.S. independently achieves a MSE matching $(\mathbb{E}_{\pi}[f(x,\theta)]-\mathbb{E}_{\pi_L}[f_L(x,\theta_L)])^2$ at lower optimal cost than directly approximating $\mathbb{E}_{\pi_L}[f_L(x,\theta_L)]$.

Algorithm¶

Our construction is that of Beskos et al. (2017)Beskos et al. (2018), and we adopt notation compatible with those works so that their results can be invoked directly. The input $x\in\mathsf{X}$ and the finest level $L\in\mathbb{N}$ are held fixed throughout. The parameter spaces are built recursively, $\Theta_l = \Theta_{l-1} \times \Delta_l$ for $l\in\{1,\dots,L\}$, with $\Delta_0=\Theta_0$; the role of these spaces will emerge shortly. Write $\theta_0=\delta_{0}\in\Theta_0$ and, for $l\in\{1,\dots,L\}$,

The coordinates $(\delta_1,\dots,\delta_{l})$ record the new parameters introduced as the width grows from level $l-1$ to level $l$; the algorithm to be presented is organized around these objects. Next introduce proposals: a strictly positive density $q_0(\theta_0)$ on $\Theta_0$ and, for each level, a strictly positive conditional density $q_l(\cdot|\theta_{l-1})$ on $\Delta_l$. The incremental weights are then $G_0(\theta_0) = \kappa_0(\theta_0)/q_0(\theta_0)$ together with

For each $l\in\{1,\dots,L-1\}$, take a Markov kernel $\mathcal{K}_l$ leaving $\pi_l$ invariant, and compose it with the proposal for the next increment, giving

in which $\theta_{l+1}' = (\theta_{l}',\delta_{l+1})$ and $d\delta_{l+1}$ denotes Lebesgue measure of the matching dimension. Step 4 of Algorithm 2 applies exactly this kernel.

Define $\eta_0(\theta_0)=q_0(\theta_0)$ and, for each $l\in\{1,\dots,L\}$,

The sampler below produces particle approximations of these densities, and thereby of expectations taken under them. In the subsequent exposition, given $(\theta_l^1,\dots,\theta_l^M)\in\Theta_l^M$, the so-called $M-$empirical measure will be denoted $\eta_l^M$. In other words we will have access to $M\in\mathbb{N}$ samples $(\theta_l^1,\dots,\theta_l^M)\in\Theta_l^M$ so that

Recalling (23), we will approximate the increments as follows

Note that $\eta_{l}^{M_{l}}(f_{l-1})$ will converge in probability (as $M_{l}\rightarrow\infty$) to $\pi_{l-1}(f_{l-1})$ Del Moral, 2004. Then $\eta_{l}^{M_{l}}(G_{l})$ will converge to

and $\eta_{l}^{M_{l}}(G_{l}f_l)$ converges to

Consequently, one can use the following approximation of $\pi_L(f_L)$:

Ordinarily we choose $q_l$ so that for each $\theta_l\in\Theta_l$

This means that

The main result is given in Chada et al. (2025).

Numerical experiments¶

Regression Problem¶

The first numerical experiment will be based on a well-specified Bayesian regression problem, as given in (14), i.e. the data is generated from the model. We let $\Sigma_i^2 =0.01^2 I$, $\sigma(z) = \tanh(z)$, $D=10$, $D_l=2^l$, and $N=200$. The inputs are simulated as $x_i \sim \mathcal{N}(2,0.5)$. We use a high-resolution parameter $l=9$ for the ground truth. We compare SMC and MLSMC samplers, using 100 replications to compute MSE. We report the complexity rate $\xi$ such that cost $\propto$ MSE$^{-\xi}$ for different values of $\beta=2s-1$. Recall that $C_l = \mathcal{O}(2^{2 l})$ so one expects to attain the canonical rate of convergence $\xi=1$ when $s>1.5 \rightarrow \beta>2$. We conduct our numerical experiment with levels $L \in \{3,4,\ldots,7\}$ and smoothness parameters $s\in\{3,1.1\}$, delivering canonical and sub-canonical convergence behaviour. The results are presented in Figure Figure 1, including credible sets around the MSE values, given by the thin blue and red curves. Further details on implementation of the MLSMC sampler in general can be found in Beskos et al. (2018).

The MLSMC sampler reaches any given accuracy more cheaply than its single-level counterpart. In particular: (i) the error-versus-cost slopes differ, with MLSMC approximately attaining the canonical $1/$MSE rate for $s=3$; (ii) at the smallest MSE the cost gap is roughly a factor of 10 for $s=3$, so the asymptotic advantage is already realized at finite resolution; and (iii) the advantage shrinks in the sub-canonical case $s=1.1$.

Binary MNIST Classification¶

As a second proof of concept, we consider a simplified binary MNIST LeCun et al., 2010 classification task using $D=100$ principal components (PCs) of the original $28\times28$ image dataset, with $N=400$. The results are similar to the regression case, and are presented in Figure Figure 2.

Parallel SMC¶

For this section, we will assume a fixed and finite network architecture $f$ of the standard form, as introduced in Section Bayesian neural networks. This gives rise to a computable target distribution (up to a normalizing constant) $\pi = \kappa/Z$.

A first source of parallelism is internal to a single SMC run: the MCMC mutation of each particle can proceed concurrently Lee et al., 2010Paige et al., 2014Syed et al., 2024, and it is in the mutations that the likelihood evaluations, the dominant computational cost, take place. Two hardware caveats apply. Each core must have sufficient memory available, or else the “memory wall” is hit Ivanov et al., 2021. Resampling forces communication among all particles, a potential bandwidth bottleneck, although on single-instruction-multiple-data (SIMD) hardware such as a single multi-core CPU or GPU this is usually benign. Being simple to implement, and sometimes exhibiting strong parallel scaling Lee et al., 2010, this intra-SMC parallelism is the most widely used.

Beyond this, one can run several SMC samplers at once and combine them, as studied in Vergé et al. (2015)Whiteley et al. (2016). The schemes recommended there typically make all samples communicate, which aids stability but works against scalability once a large model must be spread over many SIMD nodes whose inter-connect is slower, or absent altogether. In the island particle model Vergé et al., 2015, the total sample budget $M$ is divided into $R \leq R_{\rm max}$ SMC islands of $M/R$ samples apiece. If the SMCs never interact and are combined by naive (equal-weight) averaging, an asymptotic bias penalty of order $(R/M)^2$ appears Crisan et al., 2018Vergé et al., 2015. However,* weighting each SMC appropriately removes this penalty* Whiteley et al., 2016Dai et al., 2022. The scalable parallel sampler that results is denoted SMC$_\parallel$.

Algorithm 3 displays the SMC$_\parallel$ method, and we define the consistent (in $R$, for finite $M$ suitably large) SMC$_\parallel$ ratio estimator as

This estimator is consistent, which is not the case for a naive unweighted average. See Proposition 3 and the discussion above.

Parallel MCMC¶

MCMC can also exploit parallel resources: one can retain a single sample from each of $M$ parallel chains, each run for $b$ warm-up epochs, achieving $\text{MSE}=O(e^{-b}+1/M)$ Margossian et al., 2025, which is indistinguishable from the single-chain result for $b \gg \log M$, at a parallel cost equivalent to SMC when the integrated autocorrelation time matches the total number of mutation epochs. However, unlike SMC, the resulting estimator is not consistent for finite $b$, which can potentially spoil convergence. A systematic empirical comparison of parallel SMC and parallel MCMC for Bayesian deep learning is given in Liang et al. (2024).

Theoretical Result¶

The convergence result is made rigorous in the following proposition. Suitable assumptions and proof are provided in Liang et al. (2025).

Numerical experiments¶

As above, we will measure wall-clock time complexity in terms of function evaluations, as an implementation and hardware agnostic proxy. We will refer to a likelihood evaluation as an epoch, in analogy with stochastic gradient descent (SGD) Robbins & Monro, 1951Bottou, 2010 based optimization methods like Adam Kingma & Ba, 2014. However, note that the time to compute an epoch of mini-batches depends heavily on the size of the problem and the details of the hardware and implementation. If all data fits on a single computational unit (either a node or a core), then a serial epoch of mini-batch likelihood computations will typically require much longer wall-clock time in comparison to a single likelihood evaluation. In contrast, if only a single mini-batch fits on a node, and interconnect is slow, then even data parallel may be slower than serial mini-batches. It is also worth noting that mini-batch gradients can deliver comparable gain to full gradient steps far from convergence, which can significantly accelerate progress as measured by epochs Bottou et al., 2018.

Proof of concept¶

First, we will consider logistic regression on the Australian Credit data Quinlan, 1987, illustrating comparable results to state-of-the-art (SOTA) parallel MC algorithms and the $1/R$ convergence guaranteed by Proposition 3.

The dataset has $d=15$ covariates and $m=690$ data and Bayesian logistic regression is used. We compare SMC$_\parallel$-pCN (preconditioned Crank–Nicolson, pCN) with annealed importance sampling Neal, 2001, various MCMC$_\parallel$ methods, and recent synchronous-parallel MCMC Schwedes & Calderhead, 2021. Figure Figure 3 reports how the empirical variance of the posterior mean estimates decays with the number of samples and of parallel processes; the decay matches the $1/R$ rate of Proposition 3. Among the baselines in the figure that avoid communication scaling with $R$, SMC$_\parallel$-pCN achieves the lowest empirical variance. SMC and AIS incur a per-sample overhead in cost ($\sum_t m_t$, where $m_t$ are the mutations at tempering step $t$), but this is offset by intra-parallelism ($\sum_t m_t/M \ll 1$).

IMDb dataset¶

We now consider a large problem in semantic classification of natural language, using the IMDb dataset Maas et al., 2011 of $50{,}000$ movie reviews, split evenly into training and test sets.

Architecture. Each review is mapped to a fixed feature vector by SBERT embeddings Reimers & Gurevych, 2019 from the model all-mpnet-base-v2 Song et al., 2020^[2], whose weights are frozen throughout: they are used only to produce the 768 dimensional [CLS] embedding and are never sampled. On top of these embeddings we place (i) one hidden layer with 128 neurons, (ii) ReLU activation, (iii) a final linear layer, and (iv) a softmax output. Inference is over* all* weights and biases of this network, of dimension $P = 98690$, and not merely the final classifier head. The whole train (25000 data) and test dataset (25000 data) are considered.

Sampling setup. The likelihood is the softmax (categorical cross-entropy) likelihood over the whole training set of 25000 reviews, and inference uses an adaptive likelihood-tempering SMC sampler with HMC mutations of length $T$. We compare SMC$_\parallel$ against two baselines that share the same prior, likelihood, and frozen embeddings: HMC$_\parallel$, which averages over $M R$ independent HMC chains of length $T$, and a single serial HMC run of length $M R T$; in every method the HMC step size is adapted to a common acceptance-rate target, so the comparison is made at matched cost.

Experiment. Denote by $m$ the number of evaluations of $\kappa$ per mutation step, which we refer to as epochs, and by $T$ the number of tempering-and-mutation steps. The product $mT$ (or $\sum_{t=1}^T m_t$ for an adaptive schedule) is the total simulation cost, which we call the* irreducible serial time* $\tau_{\sf irr}$: the serial simulation the sampler needs in order to mix, which cannot be removed by adding more parallel runs. In practice $\tau_{\sf irr}$ is read off empirically as the budget at which accuracy and NLL plateau, here $\tau_{\sf irr}\approx 10^4$ epochs. When $m$ and/or $T$ are too small, so that the budget falls well short of $\tau_{\sf irr}$, the sampler has not mixed and the method* catastrophically fails*: test accuracy collapses and NLL inflates. Figure Figure 4 summarises this behaviour. The left panel shows performance as a function of $R$ at $\tau_{\sf irr}$, where SMC$_\parallel$ improves steadily with $R$ and matches the far more expensive serial HMC baseline; the right panel fixes $R=8$ and reduces the budget below $\tau_{\sf irr}$, exposing the catastrophic failure of SMC$_\parallel$ and HMC$_\parallel$ once the budget falls short of the irreducible time.