层次贝叶斯深度学习方法论的理论与实践进展

# Theoretical and Practical Advances in Hierarchical Bayesian Deep Learning Methodology
Source: [https://arxiv.org/html/2610.19074](https://arxiv.org/html/2610.19074)
Neil K\. Chada††thanks:Department of Mathematics, City University of Hong Kong, Hong Kong SAR\.neilchada123@gmail\.comSimon Cotter††thanks:School of Mathematics, University of Manchester, UK\.simon\.cotter@manchester\.ac\.ukKody J\. H\. Law††thanks:School of Mathematics, University of Manchester, UK\.kody\.law@manchester\.ac\.ukXinzhu Liang††thanks:School of Mathematics, University of Manchester, UK\.xinzhu\.liang@postgrad\.manchester\.ac\.uk

###### Abstract

Reliable uncertainty quantification remains an important challenge in AI\. Hierarchical Bayesian methods provide a principled framework for this by representing unknown quantities through posterior distributions\. Joint inference of neural\-network parameters and prior hyperparameters is challenging for conventional Markov chain Monte Carlo \(MCMC\) methods due to strong dependencies\. In this work, we introduce the multilevel particle marginal Metropolis\-Hastings \(MLPMMH\) method, which combines pseudo\-marginal inference for the hyperparameters with coupled sequential Monte Carlo \(SMC\) across hierarchy of network\-width truncations for conditional parameter inference\. The hierarchy is induced by a trace\-class neural network \(TNN\) prior with width\-dependent decay, motivated empirically by the compressible structure observed in trained neural\-network parameters\. Under standard multilevel coupling assumptions, we establish that MLPMMH achieves mean\-square error \(MSE\)𝒪⁡\(ε2\)\\mathcal\{O\}\(\\varepsilon^\{2\}\)at expected computational cost𝒪⁡\(ε−2\)\\mathcal\{O\}\(\\varepsilon^\{\-2\}\), and verify the predicted rate on a Bayesian regression problem\. To address larger\-scale applications, we further introduce S\-MLPMMH, a practical approximation that replaces the exact sampling components with a PCA\-based hierarchy and scalable Bayesian Monte Carlo\. Experiments on MNIST demonstrate competitive predictive performance and uncertainty quantification while exploiting the parallel structure of the approximation\. Software has been made publicly available at[https://github\.com/submission\-code\-anon/mlpmmh\-anon](https://github.com/submission-code-anon/mlpmmh-anon)\.

## 1Introduction

The focus of this work is on the design of Bayesian statistical algorithms for efficiently performing hierarchical Bayesian deep learning\([Papamarkou et al\., 2022](https://arxiv.org/html/2610.19074#bib.bib26);[Papamarkou et al\., 2024](https://arxiv.org/html/2610.19074#bib.bib3);[BDL\-Consortium, 2026](https://arxiv.org/html/2610.19074#bib.bib2)\)\. In this context, one is interested in estimating the posterior distribution of the joint parameterθ×μ∈Θ×𝖲\\theta\\times\\mu\\in\\Theta\\times\\mathsf\{S\}, where we have access to data𝒟∈𝔻\\mathcal\{D\}\\in\\mathbb\{D\}\. Then the joint posterior, via Bayes’ Theorem, is expressed as

π⁡\(θ,μ\|𝒟\)=π⁡\(𝒟\|θ,μ\)​π¯​\(θ\|μ\)​π0​\(μ\)Z⁡\(𝒟\)∝π⁡\(𝒟,θ\|μ\)​π0​\(μ\),\\pi\(\\theta,\\mu\|\\mathcal\{D\}\)=\\frac\{\\pi\(\\mathcal\{D\}\|\\theta,\\mu\)\\bar\{\\pi\}\(\\theta\|\\mu\)\\pi\_\{0\}\(\\mu\)\}\{Z\(\\mathcal\{D\}\)\}\\propto\\pi\(\\mathcal\{D\},\\theta\|\\mu\)\\pi\_\{0\}\(\\mu\)\\,,whereZ⁡\(𝒟\)Z\(\\mathcal\{D\}\)is a normalizing constant,π⁡\(𝒟\|θ,μ\)\\pi\(\\mathcal\{D\}\|\\theta,\\mu\)is the data\-likelihood,π¯​\(θ\)\\bar\{\\pi\}\(\\theta\)is the prior distribution for parameter andπ0​\(μ\)\\pi\_\{0\}\(\\mu\)is the prior distribution for hyperparameter\. Prior distributions, in general, allow us to represent an initial belief of the posterior, providing some level of flexibility\. However, this hierarchical structure can hinder generic Markov chain Monte Carlo \(MCMC\) inference, even for powerful Bayesian methodologies such as HMC\([Neal and others, 2011](https://arxiv.org/html/2610.19074#bib.bib35);[Betancourt and Girolami, 2015](https://arxiv.org/html/2610.19074#bib.bib36)\)\. Strong correlations between hyperparameters and parameters often make estimating the joint posterior difficult, leading to pathological behaviour\. The pseudo\-marginal approach is a general class of MCMC algorithms that target the marginal posterior \(∫Θπ⁡\(𝒟,𝑑θ\|μ\)​π0​\(μ\)\\int\_\{\\Theta\}\\pi\(\\mathcal\{D\},d\\theta\|\\mu\)\\pi\_\{0\}\(\\mu\)\) by replacing the intractable marginal likelihood \(∫Θπ⁡\(𝒟,𝑑θ\|μ\)\\int\_\{\\Theta\}\\pi\(\\mathcal\{D\},d\\theta\|\\mu\)\) of the hyperparameters with an unbiased stochastic estimator\([Andrieu and Roberts, 2009](https://arxiv.org/html/2610.19074#bib.bib30)\)\. An important sub\-class is known as particle MCMC\([Andrieu et al\., 2010](https://arxiv.org/html/2610.19074#bib.bib29);[Golightley and Wilkinson, 2011](https://arxiv.org/html/2610.19074#bib.bib10)\), which uses particle filters to build the unbiased distribution estimators inside a pseudo\-marginal scheme\. We will leverage the particle marginal Metropolis\-Hasting \(PMMH\) method here\.

Although Sequential Monte Carlo \(SMC\) methods are well suited for sampling in high dimensional settings, they still suffer from the curse of dimensionality, which limits their efficiency\. To mitigate this, Multilevel Monte Carlo \(MLMC\)\([Heinrich, 2001](https://arxiv.org/html/2610.19074#bib.bib28);[Giles, 2008](https://arxiv.org/html/2610.19074#bib.bib22);[Giles, 2015](https://arxiv.org/html/2610.19074#bib.bib6)\)has been introduced to achieve canonical computational complexity,1/1/MSE\. Recent work has applied MLMC in the context of statistical inference[Beskos et al\. \(2017\)](https://arxiv.org/html/2610.19074#bib.bib8);[Beskos et al\. \(2018\)](https://arxiv.org/html/2610.19074#bib.bib12);[Chada et al\. \(2021\)](https://arxiv.org/html/2610.19074#bib.bib37);[Cotter et al\. \(2022\)](https://arxiv.org/html/2610.19074#bib.bib27)and machine learning\([Chada et al\., 2025](https://arxiv.org/html/2610.19074#bib.bib16);[Fujisawa and Sato, 2021](https://arxiv.org/html/2610.19074#bib.bib21);[Gerstner et al\., 2021](https://arxiv.org/html/2610.19074#bib.bib18);[O\. Lye et al\., 2021](https://arxiv.org/html/2610.19074#bib.bib19);[Shi and Cornish, 2021](https://arxiv.org/html/2610.19074#bib.bib23)\)\. In this paper, we leverage MLMC techniques in the context of trace\-class neural networks \(TNN\) for prior modeling\([Sell and Singh, 2020](https://arxiv.org/html/2610.19074#bib.bib15)\)following a similar approach as in[Chada et al\. \(2025\)](https://arxiv.org/html/2610.19074#bib.bib16)\. The seminal work by[Jasra et al\. \(2018\)](https://arxiv.org/html/2610.19074#bib.bib9)applied the particle marginal Metropolis–Hastings scheme to parameter estimation for particle filter algorithms for parameter estimation in discretely observed continuous time state space models\. The present work applies a similar principle to the case of joint parameter and hyper\-parameter estimation in static BNN models with TNN prior architectures by simulating from an approximate coupling of the parameter conditional distributions at coarse and fine scales and correcting the increments a posteriori with importance weights\.

### 1\.1Contributions

The contributions of this work are summarized as follows:

- •Algorithmic contribution\.We introduce a*multilevel particle marginal Metropolis\-Hastings*\(MLPMMH\) scheme for BNN which builds increment estimators from approximate couplings of parameter conditional targets at successive resolutions simulated with SMC samplers, within the PMMH method for joint estimation of hyperparameters\. These increment estimators are combined to build an MLMC estimator\.
- •Theoretical guarantees\.Under standard multilevel coupling assumptions and the trace\-class prior, we prove a convergence theorem \(Theorem[4\.1](https://arxiv.org/html/2610.19074#S4.Thmtheorem1)\): MSE=𝒪⁡\(COST−1\)\.\\text\{MSE\}=\\mathcal\{O\}\\\!\\bigl\(\\text\{COST\}^\{\-1\}\\bigr\)\.
- •Practical Approximation algorithm\.A practically applicable approximation is introduced, where each increment is estimated from an ensemble of SBMC estimators[Liang et al\. \(2025\)](https://arxiv.org/html/2610.19074#bib.bib38), each conditional on a hyper\-parameter simulated from a Laplace approximation of the marginal\. This is referred to as S\-MLPMMH\.
- •Empirical validation\.\(i\) A toy Bayesian regression problem shows the expected MSE–COST convergence rate and illustrates increment variance and mean decay across levels\. \(ii\) Strong performance of S\-MLPMMH is demonstrated on the MNIST image classification problem\. \(iii\) Power\-law decay with rank of weights and biases from MAP training of Bayesian neural networks empirically supports the use of the TNN prior\. \(iv\) Several MLSMC variants for the inner parameter inference are compared in the appendix, including their computational cost\-accuracy behavior\.

### 1\.2Outline

The outline of this document is as follows\. In Section[2](https://arxiv.org/html/2610.19074#S2), we introduce the model formulation, trace\-class neural network prior, and associated multilevel structure\. We then move to Section[3](https://arxiv.org/html/2610.19074#S3), introducing our MLPMMH methodology for joint parameter and hyperparameter estimation and the practical S\-MLPMMH approximation\. Related theory, including the MSE–computational\-cost convergence result for MLPMMH, is provided in Section[4](https://arxiv.org/html/2610.19074#S4)\. In Section[5](https://arxiv.org/html/2610.19074#S5), we numerically test MLPMMH to verify our theory, motivate and validate the TNN width\-decay structure, and apply the S\-MLPMMH approximation on MNIST classification\. Finally, we summarize our findings in Section[6](https://arxiv.org/html/2610.19074#S6)\. Additional background on SMC and several existing and novel MLSMC methods for parameter estimation is provided in Appendix[B](https://arxiv.org/html/2610.19074#A2), together with supporting theory, numerical experiments, and proofs\.

![Refer to caption](https://arxiv.org/html/2610.19074v1/mlpmmh.png)Figure 1:Diagram of the MLPMMH program\. A technical diagram of the hierarchical TNN model itself is provided in Figure[4](https://arxiv.org/html/2610.19074#A1.F4)in the Appendix\.

## 2Model Formulation

Suppose we have data𝒟=\(\(x1,y1\),…,\(xN,yN\)\)\\mathcal\{D\}=\(\(x\_\{1\},y\_\{1\}\),\.\.\.,\(x\_\{N\},y\_\{N\}\)\),N∈ℕN\\in\\mathbb\{N\}, wherexi∈𝖷⊆ℝnx\_\{i\}\\in\\mathsf\{X\}\\subseteq\\mathbb\{R\}^\{n\}are inputs associated to output labelsyi∈𝖸⊆ℝmy\_\{i\}\\in\\mathsf\{Y\}\\subseteq\\mathbb\{R\}^\{m\}\. We want to infer the joint state of parameterθ∈Θ⊆ℝdθ\\theta\\in\\Theta\\subseteq\\mathbb\{R\}^\{d\_\{\\theta\}\}and the hyperparameterμ∈𝖲⊆ℝdμ\\mu\\in\\mathsf\{S\}\\subseteq\\mathbb\{R\}^\{d\_\{\\mu\}\}based on the data𝒟\\mathcal\{D\}within a deep neural network \(DNN\) modelg:𝖷×Θ×𝖲→ℝmg:\\mathsf\{X\}\\times\\Theta\\times\\mathsf\{S\}\\rightarrow\\mathbb\{R\}^\{m\}111Although the hyperparameterμ\\mudoes not enter the forward map explicitly onceθ\\thetais fixed, we retain the notationg⁡\(x,θ,μ\)g\(x,\\theta,\\mu\)to emphasize thatθ\\thetais distributed according to a prior parameterized byμ\\mu\. Thus, conditionally onθ\\theta,g⁡\(x,θ,μ\)=g⁡\(x,θ\)g\(x,\\theta,\\mu\)=g\(x,\\theta\), while the distribution of the network output depends onμ\\muthrough the prior onθ\\theta\.whereg⁡\(x,θ,μ\)=\(g1​\(x,θ,μ\),…,gm​\(x,θ,μ\)\)\.g\(x,\\theta,\\mu\)=\(g\_\{1\}\(x,\\theta,\\mu\),\.\.\.,g\_\{m\}\(x,\\theta,\\mu\)\)\.

### 2\.1Trace\-class neural network priors

These TNN priors have also been introduced and well\-established in the[Chada et al\. \(2025\)](https://arxiv.org/html/2610.19074#bib.bib16), and here we followed a similar procedure: ford=1,…,Dd=1,\.\.\.,D, the associated weights and bias \(θ\\theta\) are defined as

Ai​j,d\\displaystyle A\_\{ij,d\}∼𝒩⁡\(0,σA2​\(i​j\)−s\),\\displaystyle\\sim\\mathcal\{N\}\(0,\\sigma\_\{A\}^\{2\}\(ij\)^\{\-s\}\),\(2\.1\)Ai​j,1\\displaystyle\\quad A\_\{ij,1\}∼𝒩⁡\(0,σA2​i−s\),\\displaystyle\\sim\\mathcal\{N\}\(0,\\sigma\_\{A\}^\{2\}i^\{\-s\}\),Ai​j,D\\displaystyle\\quad A\_\{ij,D\}∼𝒩⁡\(0,σA2​j−s\),\\displaystyle\\sim\\mathcal\{N\}\(0,\\sigma\_\{A\}^\{2\}j^\{\-s\}\),bi,d\\displaystyle b\_\{i,d\}∼𝒩⁡\(0,σb2​i−s\)\.\\displaystyle\\sim\\mathcal\{N\}\(0,\\sigma\_\{b\}^\{2\}i^\{\-s\}\)\.
The model, presented in Figure[4](https://arxiv.org/html/2610.19074#A1.F4), is constructed based on the pushforward of a countably infinite Gaussian random variable that has a specifiedweight decayin the width direction \(which can similarly be applied to depth\)\. This instils a clear inductive bias—analogous to choosing a kernel to specify the smoothness of a Gaussian process \(GP\) while remaining fundamentally different\. These priorsconverge strongly, while GPNN typically converge only weakly, so that it allows us to apply the MLMC methodology\. Since the hyperparameterssis a tricky one, whose value determines the strong convergence rate in the MLMC methodology, we will set it to be a constant value here\. We then can let the hyperparameter beμ=\(σA,σb\)\\mu=\(\\sigma\_\{A\},\\sigma\_\{b\}\), and a suitable prior could be the 2\-dimensional log\-normal distribution\.

### 2\.2Multilevel Structure

Given the levell∈ℕl\\in\\mathbb\{N\}, for simplicity, we assumeDDis fixed and let the NN width is the same for all layer other than the input and the output, i\.e\.nd=nl=2ln\_\{d\}=n\_\{l\}=2^\{l\}ford∈\{1,…,D−1\}d\\in\\\{1,\.\.\.,D\-1\\\},n0=nn\_\{0\}=nandnD=mn\_\{D\}=m\. The corresponding vector of parameters is denoted byθl:=\(\(A1l,b1l\),…,\(ADl,bDl\)\)\\theta\_\{l\}:=\(\(A\_\{1\}^\{l\},b\_\{1\}^\{l\}\),\.\.\.,\(A\_\{D\}^\{l\},b\_\{D\}^\{l\}\)\)which is in the parameter spaceΘl=\{ℝnl×n×ℝnl\}⊗\{ℝnl×nl×ℝnl\}⊗D−2⊗\{ℝm×nl×ℝm\}\\Theta\_\{l\}=\\\{\\mathbb\{R\}^\{n\_\{l\}\\times n\}\\times\\mathbb\{R\}^\{n\_\{l\}\}\\\}\\otimes\\\{\\mathbb\{R\}^\{n\_\{l\}\\times n\_\{l\}\}\\times\\mathbb\{R\}^\{n\_\{l\}\}\\\}^\{\\otimes D\-2\}\\otimes\\\{\\mathbb\{R\}^\{m\\times n\_\{l\}\}\\times\\mathbb\{R\}^\{m\}\\\}and the prior isπ¯l​\(θl\)\\bar\{\\pi\}\_\{l\}\(\\theta\_\{l\}\)\.

Here, we introduce the canonical truncation mapTl:Θl→Θl−1T\_\{l\}:\\Theta\_\{l\}\\to\\Theta\_\{l\-1\}\. Given a fine\-level parameterθl∈Θl\\theta\_\{l\}\\in\\Theta\_\{l\}, we writeθl−1=Tl​\(θl\)\\theta\_\{l\-1\}=T\_\{l\}\(\\theta\_\{l\}\)for the corresponding coarse\-level parameter obtained by setting:Ai​j,dl=Ai​j,dl−1A^\{l\}\_\{ij,d\}=A^\{l\-1\}\_\{ij,d\}andbi,dl=bi,dl−1b^\{l\}\_\{i,d\}=b^\{l\-1\}\_\{i,d\}fori,j∈1,…,nl−1i,j\\in\{1,\\ldots,n\_\{l\-1\}\}andd∈\{2,…,D\}d\\in\\\{2,\\dots,D\\\}; a similar assignment is to be adopted for\(A1l,b1l\)\(A^\{l\}\_\{1\},b^\{l\}\_\{1\}\)and\(ADl,bDl\)\(A^\{l\}\_\{D\},b^\{l\}\_\{D\}\)\.

Then, we denote the DNN output bygl​\(x,θl,μ\)g\_\{l\}\(x,\\theta\_\{l\},\\mu\)\. Using the same output mapΨ\\Psiat every level, we define the quantity of interest byφl​\(x,θl,μ\):=\(Ψ∘gl\)​\(x,θl,μ\)=Ψ⁡\(gl​\(x,θl,μ\)\)\\varphi\_\{l\}\(x,\\theta\_\{l\},\\mu\):=\(\\Psi\\circ g\_\{l\}\)\(x,\\theta\_\{l\},\\mu\)=\\Psi\(g\_\{l\}\(x,\\theta\_\{l\},\\mu\)\)\. The posterior distribution of the parameterθl\\theta\_\{l\}conditioned onμ\\muis

πl​\(θl\|μ\)∝ℒl​\(θl\|μ\)​π¯l​\(θl\)\\pi\_\{l\}\(\\theta\_\{l\}\|\\mu\)\\propto\\mathcal\{L\}\_\{l\}\(\\theta\_\{l\}\|\\mu\)\\bar\{\\pi\}\_\{l\}\(\\theta\_\{l\}\)
πl​\(θ\|μ\)=fl​\(θ\|μ\)Zl\\pi\_\{l\}\(\\theta\|\\mu\)=\\frac\{f\_\{l\}\(\\theta\|\\mu\)\}\{Z\_\{l\}\}, wherefl​\(θl\|μ\)=ℒl​\(θl\|μ\)​π¯l​\(θl\)f\_\{l\}\(\\theta\_\{l\}\|\\mu\)=\\mathcal\{L\}\_\{l\}\(\\theta\_\{l\}\|\\mu\)\\bar\{\\pi\}\_\{l\}\(\\theta\_\{l\}\)andZl=∫Θlfl​𝑑θZ\_\{l\}=\\int\_\{\\Theta\_\{l\}\}f\_\{l\}d\\theta\. Then we construct the approximate coupling,Πl\\Pi\_\{l\}, as follows

Πl​\(θl,θl−1\)\\displaystyle\\Pi\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\)=Fl​\(θl,θl−1\|μ\)Fl​\(1\|μ\),\\displaystyle=\\frac\{F\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\}\{F\_\{l\}\(1\|\\mu\)\},Fl​\(θl,θl−1\|μ\)\\displaystyle F\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)=𝐋l​\(θl,θl−1\|μ\)​Π¯l​\(θl,θl−1\),\\displaystyle=\{\\bf L\}\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\\bar\{\\Pi\}\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\),where we choose the joint likelihood,𝐋l​\(θl,θl−1\|μ\)\{\\bf L\}\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\), as the maximum coupling

𝐋l​\(θl,θl−1\|μ\)=max⁡\{ℒl​\(θl\|μ\),ℒl−1​\(θl−1\|μ\)\},\{\\bf L\}\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)=\\max\\\{\\mathcal\{L\}\_\{l\}\(\\theta\_\{l\}\|\\mu\),\\mathcal\{L\}\_\{l\-1\}\(\\theta\_\{l\-1\}\|\\mu\)\\\},and the coupling of the prior distributionΠ¯l\\bar\{\\Pi\}\_\{l\}is

Π¯l​\(d​θl,d​θl−1\)=π¯l​\(d​θl\)​δTl​\(θl\)​\(d​θl−1\)\.\\bar\{\\Pi\}\_\{l\}\(d\\theta\_\{l\},d\\theta\_\{l\-1\}\)=\\bar\{\\pi\}\_\{l\}\(d\\theta\_\{l\}\)\\,\\delta\_\{T\_\{l\}\(\\theta\_\{l\}\)\}\(d\\theta\_\{l\-1\}\)\.

## 3Hierarchical Bayesian Estimation of Neural Networks with Trace\-Class Priors

In this section, we provide and present our primary contribution, which is our methodology for hierarchical Bayesian estimation of the hyperparameter\. This will be based on particle Markov chain Monte Carlo, specifically the particle marginal Metropolis\-Hastings method\([Chada et al\., 2021](https://arxiv.org/html/2610.19074#bib.bib37);[Jasra et al\., 2018](https://arxiv.org/html/2610.19074#bib.bib9)\)\. We review this method\([Jasra et al\., 2018](https://arxiv.org/html/2610.19074#bib.bib9)\)and design it in such a way related to our framework and trace\-class neural network priors\. Further details of related algorithms are provided in the supplementary material\.

### 3\.1MLPMMH

Now, we embed the MLSMC sampler in Algorithm[16](https://arxiv.org/html/2610.19074#alg3.l16)inside the PMMH method to estimate the marginal likelihood, which gives the multilevel pseudo\-marginal Metropolis\-Hastings \(MLPMMH\) method\. In the context of Bayesian neural networks, we infer the prior standard\-deviation hyperparameters governing the weights and biases of the TNNs model with log\-normal priors\. Due to the multilevel setup, the joint posterior distribution at each levelllis given by

πl​\(θl,μ\)∝ℒl​\(θl\|μ\)​π¯l​\(θl\)⏟fl​\(θ\|μ\)​π0​\(μ\)=:fl​\(θl,μ\)\.\\pi\_\{l\}\(\\theta\_\{l\},\\mu\)\\propto\\underbrace\{\\mathcal\{L\}\_\{l\}\(\\theta\_\{l\}\|\\mu\)\\bar\{\\pi\}\_\{l\}\(\\theta\_\{l\}\)\}\_\{f\_\{l\}\(\\theta\|\\mu\)\}\\pi\_\{0\}\(\\mu\)=:f\_\{l\}\(\\theta\_\{l\},\\mu\)\.
One can run a MLSMC sampler in Algorithm[16](https://arxiv.org/html/2610.19074#alg3.l16)to approximate the marginal distribution

Πl​\(θl,θl−1\|μ\)\\displaystyle\\Pi\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)∝𝐋l​\(θl,θl−1\|μ\)​Π¯l​\(θl,θl−1\)\\displaystyle\\propto\{\\bf L\}\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\\bar\{\\Pi\}\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\)=:Fl​\(θl,θl−1\|μ\)\.\\displaystyle=:F\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\.
Then the joint coupling posterior distribution at each level can be approximate using Algorithm[12](https://arxiv.org/html/2610.19074#alg1.l12)as

ΠlN,Ml​\(d⁡\(θ,μ\)\):=1Ml​∑k=1Mlδ\(θk,l,μk,l\)​\(d⁡\(θ,μ\)\),\\Pi\_\{l\}^\{N,M\_\{l\}\}\(d\(\\theta,\\mu\)\):=\\frac\{1\}\{M\_\{l\}\}\\sum\_\{k=1\}^\{M\_\{l\}\}\\delta\_\{\(\\theta\_\{k,l\},\\mu\_\{k,l\}\)\}\(d\(\\theta,\\mu\)\),
Give a level parameterL=\{Lmin,Lmax\}L=\\\{L\_\{\\text\{min\}\},L\_\{\\text\{max\}\}\\\}, the MLPMMH estimator is given by

φ^LMLPMMH:=∑l=LminLmax\(1Ml​∑k=1Mlφl​\(θk,l,μk,l\)​fl​\(θk,l\|μk,l\)Fl​\(θk,l,θk,l−1\|μk,l\)1Ml​∑k=1Mlfl​\(θk,l\|μk,l\)Fl​\(θk,l,θk,l−1\|μk,l\)CLOSEOPEN−1Ml​∑k=1Mlφl−1​\(θk,l−1,μk,l\)​fl−1​\(θk,l−1\|μk,l\)Fl​\(θk,l,θk,l−1\|μk,l\)1Ml​∑k=1Mlfl−1​\(θk,l−1\|μk,l\)Fl​\(θk,l,θk,l−1\|μk,l\)\),\\begin\{gathered\}\\hat\{\\varphi\}\_\{L\}^\{\\text\{MLPMMH\}\}:=\\sum\_\{l=L\_\{\\text\{min\}\}\}^\{L\_\{\\text\{max\}\}\}\\bigg\(\\frac\{\\frac\{1\}\{M\_\{l\}\}\\sum\_\{k=1\}^\{M\_\{l\}\}\\varphi\_\{l\}\(\\theta\_\{k,l\},\\mu\_\{k,l\}\)\\frac\{f\_\{l\}\(\\theta\_\{k,l\}\|\\mu\_\{k,l\}\)\}\{F\_\{l\}\(\\theta\_\{k,l\},\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)\}\}\{\\frac\{1\}\{M\_\{l\}\}\\sum\_\{k=1\}^\{M\_\{l\}\}\\frac\{f\_\{l\}\(\\theta\_\{k,l\}\|\\mu\_\{k,l\}\)\}\{F\_\{l\}\(\\theta\_\{k,l\},\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)\}\}\\\\ \-\\frac\{\\frac\{1\}\{M\_\{l\}\}\\sum\_\{k=1\}^\{M\_\{l\}\}\\varphi\_\{l\-1\}\(\\theta\_\{k,l\-1\},\\mu\_\{k,l\}\)\\frac\{f\_\{l\-1\}\(\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)\}\{F\_\{l\}\(\\theta\_\{k,l\},\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)\}\}\{\\frac\{1\}\{M\_\{l\}\}\\sum\_\{k=1\}^\{M\_\{l\}\}\\frac\{f\_\{l\-1\}\(\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)\}\{F\_\{l\}\(\\theta\_\{k,l\},\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)\}\}\\bigg\),\\end\{gathered\}\(3\.1\)whereφLmin−1:=0\\varphi\_\{L\_\{\\text\{min\}\}\-1\}:=0\.

##### Parallization\.

The MLPMMH estimator also admits several levels of parallelization\. First, the level\-wise PMMH estimators can be run independently for different \(ll\)\. Second, within each PMMH iteration, the SMC estimator used to estimate the level\-dependent normalizing constant can exploit particle\-level parallelism\. Thus the main available parallelism is across levels and across SMC particles\. The PMMH chain itself remains sequential in the iteration index \(k\), although independent chains or independent repetitions may be run in parallel\.

1:Simulate

μ0∼π0\\mu\_\{0\}\\sim\\pi\_\{0\}and run Algorithm[16](https://arxiv.org/html/2610.19074#alg3.l16)for level

llwith

NNand

μ0\\mu\_\{0\}, it outputs

\(F0,l,θ^0,l\)\(F\_\{0,l\},\\hat\{\\theta\}\_\{0,l\}\)\.

2:Define

p0,lN​\(μ0\)=1N​∑i=1NF0,l​\(θ^0,li,θ^0,l−1i\|μ0\)p\_\{0,l\}^\{N\}\(\\mu\_\{0\}\)=\\frac\{1\}\{N\}\\sum\_\{i=1\}^\{N\}F\_\{0,l\}\(\\hat\{\\theta\}\_\{0,l\}^\{i\},\\hat\{\\theta\}\_\{0,l\-1\}^\{i\}\|\\mu\_\{0\}\)and

w0,li=F0,l​\(θ^0,li,θ^0,l−1i\|μ0\)∑i=1NF0,l​\(θ^0,li,θ^0,l−1i\|μ0\)w\_\{0,l\}^\{i\}=\\frac\{F\_\{0,l\}\(\\hat\{\\theta\}^\{i\}\_\{0,l\},\\hat\{\\theta\}\_\{0,l\-1\}^\{i\}\|\\mu\_\{0\}\)\}\{\\sum\_\{i=1\}^\{N\}F\_\{0,l\}\(\\hat\{\\theta\}^\{i\}\_\{0,l\},\\hat\{\\theta\}\_\{0,l\-1\}^\{i\}\|\\mu\_\{0\}\)\}\.

3:Resample: Select

I0,l1∼\{w0,l1,…,w0,lN\}I\_\{0,l\}^\{1\}\\sim\\\{w\_\{0,l\}^\{1\},\.\.\.,w\_\{0,l\}^\{N\}\\\}, and let

θ0,l=θ^0,lI0,l1\\theta\_\{0,l\}=\\hat\{\\theta\}\_\{0,l\}^\{I\_\{0,l\}^\{1\}\}\.

4:for

k=1,…,Mlk=1,\.\.\.,M\_\{l\}do

5:Propose

μ^k∼Q⁡\(⋅,μk−1\)\\hat\{\\mu\}\_\{k\}\\sim Q\(\\cdot,\\mu\_\{k\-1\}\)\. Run Algorithm[16](https://arxiv.org/html/2610.19074#alg3.l16)for level

llwith

NNand

μ^k\\hat\{\\mu\}\_\{k\}, and output

\{θ^k,li,θ^k,l−1i,fk,li​\(θ^k,li\),fk,l−1i​\(θ^k,l−1i\),Fk,li\}i=1N\\left\\\{\\hat\{\\theta\}\_\{k,l\}^\{i\},\\hat\{\\theta\}\_\{k,l\-1\}^\{i\},f\_\{k,l\}^\{i\}\(\\hat\{\\theta\}\_\{k,l\}^\{i\}\),f\_\{k,l\-1\}^\{i\}\(\\hat\{\\theta\}\_\{k,l\-1\}^\{i\}\),F\_\{k,l\}^\{i\}\\right\\\}\_\{i=1\}^\{N\}\.

6:Define

pk,lN​\(μ^k\)=1N​∑i=1NFk,l​\(θ^k,li,θ^k,l−1i\|μ^k\)p\_\{k,l\}^\{N\}\(\\hat\{\\mu\}\_\{k\}\)=\\frac\{1\}\{N\}\\sum\_\{i=1\}^\{N\}F\_\{k,l\}\(\\hat\{\\theta\}\_\{k,l\}^\{i\},\\hat\{\\theta\}\_\{k,l\-1\}^\{i\}\|\\hat\{\\mu\}\_\{k\}\)and

wk,li=Fk,l​\(θ^k,li,θ^k,l−1i\|μ^k\)∑i=1NFk,l​\(θ^k,li,θ^k,l−1i\|μ^k\)\.w\_\{k,l\}^\{i\}=\\frac\{F\_\{k,l\}\(\\hat\{\\theta\}^\{i\}\_\{k,l\},\\hat\{\\theta\}\_\{k,l\-1\}^\{i\}\|\\hat\{\\mu\}\_\{k\}\)\}\{\\sum\_\{i=1\}^\{N\}F\_\{k,l\}\(\\hat\{\\theta\}^\{i\}\_\{k,l\},\\hat\{\\theta\}\_\{k,l\-1\}^\{i\}\|\\hat\{\\mu\}\_\{k\}\)\}\.
7:Resample: Select

Ik,l1∼\{wk,l1,…,wk,lN\}I\_\{k,l\}^\{1\}\\sim\\\{w\_\{k,l\}^\{1\},\.\.\.,w\_\{k,l\}^\{N\}\\\}, and let

θk,l=θ^k,lIk,l1\\theta\_\{k,l\}=\\hat\{\\theta\}\_\{k,l\}^\{I\_\{k,l\}^\{1\}\},

θk,l−1=θ^k,l−1Ik,l1\\theta\_\{k,l\-1\}=\\hat\{\\theta\}\_\{k,l\-1\}^\{I\_\{k,l\}^\{1\}\}\.

8:Simulate

u∼U⁡\[\(0,1\)\]u\\sim U\[\(0,1\)\]and compute

α=min⁡\{π0​\(μ^k\)​Q​\(μk−1\|μ^k\)​pk,lN​\(μ^k\)π0​\(μk−1\)​Q​\(μ^k\|μk−1\)​pk,lN​\(μk−1\),1\}\.\\alpha=\\min\\Bigg\\\{\\frac\{\\pi\_\{0\}\(\\hat\{\\mu\}\_\{k\}\)Q\(\\mu\_\{k\-1\}\|\\hat\{\\mu\}\_\{k\}\)p\_\{k,l\}^\{N\}\(\\hat\{\\mu\}\_\{k\}\)\}\{\\pi\_\{0\}\(\\mu\_\{k\-1\}\)Q\(\\hat\{\\mu\}\_\{k\}\|\\mu\_\{k\-1\}\)p\_\{k,l\}^\{N\}\(\\mu\_\{k\-1\}\)\},\\ 1\\Bigg\\\}\.
9:If

u≤αu\\leq\\alpha, then

\(μk,θk,l,pk,lN​\(μk\)\)=\(μ^k,θ^k,l,pk,lN​\(μ^k\)\)\(\\mu\_\{k\},\\theta\_\{k,l\},p\_\{k,l\}^\{N\}\(\\mu\_\{k\}\)\)=\(\\hat\{\\mu\}\_\{k\},\\hat\{\\theta\}\_\{k,l\},p\_\{k,l\}^\{N\}\(\\hat\{\\mu\}\_\{k\}\)\); otherwise,

10:

\(μk,θk,l,pk,lN​\(μk\)\)=\(μk−1,θk−1,l,pk,lN​\(μk−1\)\)\(\\mu\_\{k\},\\theta\_\{k,l\},p\_\{k,l\}^\{N\}\(\\mu\_\{k\}\)\)=\(\\mu\_\{k\-1\},\\theta\_\{k\-1,l\},p\_\{k,l\}^\{N\}\(\\mu\_\{k\-1\}\)\)\.

11:endfor

12:Define

μk,l=μk\\mu\_\{k,l\}=\\mu\_\{k\}\. Store

θk,l\\theta\_\{k,l\},

θk,l−1\\theta\_\{k,l\-1\},

fl−1​\(θk,l−1\|μk,l\)=π¯l−1​\(θk,l−1\)​ℒl−1​\(θk,l−1\|μk,l\)f\_\{l\-1\}\(\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)=\\bar\{\\pi\}\_\{l\-1\}\(\\theta\_\{k,l\-1\}\)\\mathcal\{L\}\_\{l\-1\}\(\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\),

fl​\(θk,l\|μk,l\)=π¯l​\(θk,l\)​ℒl​\(θk,l\|μk,l\)f\_\{l\}\(\\theta\_\{k,l\}\|\\mu\_\{k,l\}\)=\\bar\{\\pi\}\_\{l\}\(\\theta\_\{k,l\}\)\\mathcal\{L\}\_\{l\}\(\\theta\_\{k,l\}\|\\mu\_\{k,l\}\)and

Fl​\(θk,l,θk,l−1\|μk,l\)F\_\{l\}\(\\theta\_\{k,l\},\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)\.

Algorithm 1MLPMMH for levelll\([Jasra et al\., 2018](https://arxiv.org/html/2610.19074#bib.bib9)\)

### 3\.2S\-MLPMMH

Inspired by the SBMC framework in[Liang et al\. \(2025\)](https://arxiv.org/html/2610.19074#bib.bib38), we introduce a practical SBMC\-based approximation of MLPMMH designed to reduce its serial computational cost, called S\-MLPMMH\. The original MLPMMH method uses an outer serial PMMH chain to sample the hyperparameter, and each outer proposal requires running an inner MLSMC estimator\. To obtain a scalable practical version, we replace the serial outer PMMH chain with a Laplace\-based SBMC approximation to the hyperparameter posterior, centered at an empirical\-Bayes MAP estimate\. We also replace the inner MLSMC sampler by a MAP\-centered SBMC sampler\.

Specifically, at each levelll, we consider the transformed hyperparameterη=log⁡μ\\eta=\\log\\muand construct a local Laplace approximation centered at a level\-dependent MAP estimateηl,𝖬𝖠𝖯\\eta\_\{l,\\sf MAP\}\. The outer hyperparameter samples are independently drawn as

ηk,l∼𝒩⁡\(ηl,𝖬𝖠𝖯,Vl,η\),μk,l=exp⁡\(ηk,l\),\\eta\_\{k,l\}\\sim\\mathcal\{N\}\\left\(\\eta\_\{l,\\sf MAP\},V\_\{l,\\eta\}\\right\),\\qquad\\mu\_\{k,l\}=\\exp\(\\eta\_\{k,l\}\),fork=1,…,Mlk=1,\\ldots,M\_\{l\}, whereVl,ηV\_\{l,\\eta\}is obtained from the local curvature of the empirical\-Bayes objective\.

For each sampled hyperparameterμk,l\\mu\_\{k,l\}, we then construct an inner MAP\-centered SBMC sampler\. In particular, we compute the corresponding MAP estimateθl,𝖬𝖠𝖯\\theta\_\{l,\\sf MAP\}and construct a new inner prior centered at this estimate,

π¯^l=𝒩⁡\(θl,𝖬𝖠𝖯,sθ​μk,l2​𝖨\)\.\\hat\{\\bar\{\\pi\}\}\_\{l\}=\\mathcal\{N\}\\left\(\\theta\_\{l,\\sf MAP\},s\_\{\\theta\}\\mu\_\{k,l\}^\{2\}\\sf I\\right\)\.The MAP estimate is recomputed for each hyperparameter sampleμk,l\\mu\_\{k,l\}\.

This practical version preserves the multilevel correction structure of MLPMMH while reducing both the serial outer cost and the cost of inner sampling\. Moreover, due to the parallelizability of the SBMC framework, the outer hyperparameter samples and their corresponding inner SBMC computations can be processed independently and therefore in parallel\. Note, however, that S\-MLPMMH does not constitute an exact pseudo\-marginal sampler for the original joint posterior, and the theoretical results for MLPMMH in Section[4](https://arxiv.org/html/2610.19074#S4)do not directly apply to this practical approximation\. Further details of the empirical\-Bayes, Laplace, and MAP\-centered SBMC construction are provided in Appendix[D](https://arxiv.org/html/2610.19074#A4)\.

## 4Theoretical results

The short section is devoted to our main theoretical contribution, which is an MSE bound with our MLPMMH result\. It is requires several important additional results, which are provided in Appendix[E](https://arxiv.org/html/2610.19074#A5)and proven in Appendix[H](https://arxiv.org/html/2610.19074#A8)\.

### 4\.1Main results

###### Theorem 4\.1\.

Let us assume Assumptions[E\.1](https://arxiv.org/html/2610.19074#A5.Thmass1),[E\.3](https://arxiv.org/html/2610.19074#A5.Thmass3)and[E\.2](https://arxiv.org/html/2610.19074#A5.Thmass2)withβ\>γ\\beta\>\\gamma\. Then, for suitableφ:Θ×𝖲→ℝ\\varphi:\\Theta\\times\\mathsf\{S\}\\rightarrow\\mathbb\{R\}, there exists aC<\+∞C<\+\\inftysuch that for anyε\>0\\varepsilon\>0one can choose\(L=\{Lmin,Lmax\},\{Ml\}l=LminLmax\)\(L=\\\{L\_\{\\text\{min\}\},L\_\{\\text\{max\}\}\\\},\\\{M\_\{l\}\\\}\_\{l=L\_\{\\text\{min\}\}\}^\{L\_\{\\text\{max\}\}\}\)such that

𝔼⁡\[\(φ^LMLPMMH−π⁡\(φ\)\)2\]≤C​ε2,\\mathbb\{E\}\[\(\\hat\{\\varphi\}\_\{L\}^\{\\mathrm\{MLPMMH\}\}\-\\pi\(\\varphi\)\)^\{2\}\]\\leq C\\varepsilon^\{2\},and expected COST\(φ^LMLPMMH\)≤C​ε−2\(\\hat\{\\varphi\}\_\{L\}^\{\\mathrm\{MLPMMH\}\}\)\\leq C\\varepsilon^\{\-2\}, i\.e\. the canonical rate\. The estimator is defined in equation \([3\.1](https://arxiv.org/html/2610.19074#S3.E1)\)\.

###### Proof\.

The proof sees Appendix[H\.1](https://arxiv.org/html/2610.19074#A8.SS1)\. ∎

Whereas, the computational cost of PMMH to achieve the MSE=O⁡\(ε2\)=O\(\\varepsilon^\{2\}\)is𝒪⁡\(ε−2−γ/α\)\\mathcal\{O\}\(\\varepsilon^\{\-2\-\\gamma/\\alpha\}\)\. It implies that PMMH will often suffer from the curse of dimensionality since the ratioγα\\frac\{\\gamma\}\{\\alpha\}will increase as the dimension becomes high\.

## 5Numerical Experiments

In this section, we provide numerical experiments through our MLPMMH scheme, using TNN neural networks\. We will not only verify our complexity rate \(through Theorem[4\.1](https://arxiv.org/html/2610.19074#S4.Thmtheorem1)\) but also justify the choice of the TNN prior\. That will be done by checking its weight decay in the width limit\. Experiments will be conducted on a Bayesian regression and MNIST classification problem while providing a comparison with other methods\.

### 5\.1Numerical justification for MLPMMH

#### 5\.1\.1Synthetic regression dataset

GenerateNt​r=100N\_\{tr\}=100data points,yiy\_\{i\}, followed[F\.1](https://arxiv.org/html/2610.19074#A6.E1)at levell=7l=7usingxi=2\+0\.5⋅νi∈ℝ10x\_\{i\}=2\+\\sqrt\{0\.5\}\\cdot\\nu\_\{i\}\\in\\mathbb\{R\}^\{10\}, whereνi∼𝒩⁡\(0,I10\)\\nu\_\{i\}\\sim\\mathcal\{N\}\(0,I\_\{10\}\)fori=1,…,100i=1,\\dots,100\. Evaluation is then applied on one fixed test pointx=\[2,2,…,2\]∈ℝ10x=\[2,2,\\dots,2\]\\in\\mathbb\{R\}^\{10\}\.

In the section, we numerically verified the rate decay in Assumption[E\.2](https://arxiv.org/html/2610.19074#A5.Thmass2)and convergence rate in Theorem[4\.1](https://arxiv.org/html/2610.19074#S4.Thmtheorem1)on the synthetic regression example using the Bayesian TNN model\. The number of neural layers isD=3D=3, the input neural network width is then set ton=10n=10, and the output neural network width is set tom=1m=1\. The TNN prior uses a rate decay ofs=2s=2, the MCMC process in MLPMMH takes5050burn\-in samples and adaptive step size, and the HMC kernel is applied with1010mutation steps, one leapfrog step and adaptive step size in the MLSMC\. The hyperparameter prior is the LogNormal distribution with the location parameter as−0\.125\-0\.125and scale parameter as0\.50\.5\.

A regularity test, similar to MLSMC, is conducted to estimate the decay of bias and variance for normalized estimator differences between levels under the posterior distribution\. The results are shown in Figure[7](https://arxiv.org/html/2610.19074#A6.F7), supporting the decay rate assumed in Assumption[E\.2](https://arxiv.org/html/2610.19074#A5.Thmass2)for the convergence analysis\. Again, starting from levell=4l=4, the decay of biases and weights aligns with the expected rate, where the bias decay \(weak convergence\) is bounded byα=β/2\\alpha=\\beta/2and the strong convergence rate isβ=2​s−1=3\\beta=2s\-1=3withs=2s=2\.

TakingLmin=4L\_\{\\text\{min\}\}=4as the starting level, the convergence test for the MLPMMH is presented in Figure[2](https://arxiv.org/html/2610.19074#S5.F2), i\.e\. MSE=𝒪⁡\(COST−1\)=\\mathcal\{O\}\(\\text\{COST\}^\{\-1\}\)\. Given a target accuracy levelε\\varepsilon, we can follow a similar procedure using the mean and variance decay to compute the maximum levelLmaxL\_\{\\max\}and the number of samples per levelMlM\_\{l\}in the MLPMMH as in Section[G\.1](https://arxiv.org/html/2610.19074#A7.SS1), and set theN=100N=100samples in the MLSMC for estimating the marginal distribution\. The reference solution is computed at accuracyε=8×10−4\\varepsilon=8\\times 10^\{\-4\}using MLPMMH\.

![Refer to caption](https://arxiv.org/html/2610.19074v1/mse_vs_cost_loglog.png)Figure 2:Complexity test for MLPMMH on Bayesian regression model: We use1010realizations to compute the MSE for each experiment and plot the rate of convergence w\.r\.t computation costs\. The regression coefficient is−0\.9\-0\.9\.

### 5\.2Empirical verification of TNN prior structure

To numerically validate the TNN prior, we conduct a series of controlled MAP training experiments on the MNIST dataset\.

A key issue in interpreting width\-dependent decay is the exchangeability of hidden units\. For NoReg and weight\-decay, the ordering of hidden units is arbitrary: the same network function can often be represented after permuting the hidden\-unit indices\. Therefore, a trained network may contain a low\-effective\-width or compressible structure even if this structure is not aligned with the original width ordering\. For this reason we use a rank\-wise diagnostic, which sorts learned weight magnitudes and tests whether compressibility is present at all\. This informs the width\-decay that can be non\-intrusively imposed when we break the exchangeability\.

We sort the absolute weights in each layer in descending order and plot\|w\|\|w\|against rank on log–log axes, as shown in Figure[3](https://arxiv.org/html/2610.19074#S5.F3)\. Since this sorting removes the original width indices, the figure should be interpreted as a compressibility diagnostic rather than as direct evidence of index\-wise TNN structure\. The NoReg and weight\-decay models already show clear decay in sorted weight magnitudes, suggesting that trained neural networks naturally concentrate much of their mass in a relatively small number of effective directions\. The width\-decay model preserves a comparable leading\-rank structure while inducing a stronger decay in the tail, especially in the hidden layers\. This suggests that the TNN prior is not imposing an unrelated structure from scratch; rather, it makes an already present compressibility structure explicit and ordered\.

![Refer to caption](https://arxiv.org/html/2610.19074v1/weights_decay_s2p0_sigmaA1p0_sigmab1p0_reproduced.png)Figure 3:Rank\-wise decay of learned weight magnitudes fors=2\.0s=2\.0andσA=σb=1\.0\\sigma\_\{A\}=\\sigma\_\{b\}=1\.0\. For each layer, the absolute values of the learned weights are sorted in descending order and plotted against rank on log\-log axes\.Finally, Table[1](https://arxiv.org/html/2610.19074#S5.T1)shows that width\-decay achieves performance comparable to NoReg and weight\-decay for moderate choices of the decay parameterss\. In particular, fors≤4s\\leq 4, width\-decay maintains similar accuracy, NLL, Brier score, and ECE relative to the two baselines\. This indicates that introducing the TNN ordering does not substantially compromise predictive performance within a reasonable range of decay strengths\. However, very large values such ass=10\.0s=10\.0lead to numerical instability and poor performance, suggesting that overly aggressive width decay should be avoided\. Since the appropriate value ofssis problem\-dependent, we discuss its practical choice further in the MLPMMH experiments below\.

MethodNoRegWeight\-decayWidth\-decayss0\.51\.02\.04\.010\.0Acc \(%\)97\.92 ± 0\.2397\.73 ± 0\.0697\.82 ± 0\.1497\.82 ± 0\.2897\.65 ± 0\.1597\.07 ± 0\.1389\.28 ± 1\.05NLL0\.106 ± 0\.0150\.110 ± 0\.0050\.108 ± 0\.0050\.112 ± 0\.0120\.105 ± 0\.0110\.105 ± 0\.0050\.374 ± 0\.042Brier0\.035 ± 0\.0040\.038 ± 0\.0010\.036 ± 0\.0020\.036 ± 0\.0040\.038 ± 0\.0020\.045 ± 0\.0020\.161 ± 0\.016ECE0\.014 ± 0\.0020\.015 ± 0\.0000\.014 ± 0\.0010\.015 ± 0\.0020\.015 ± 0\.0010\.008 ± 0\.0010\.027 ± 0\.004Table 1:Metrics for three methods evaluated on the MNIST dataset: MLE, standard weight decay regularization, and TNN width\-decay regularization are used for training on the MNIST training dataset with PyTorch default initialization and testing on the MNIST test dataset\. TNN width\-decay regularization with fixedσ=1\\sigma=1and varying decay parameterss\. Results are averaged over 5 replicates, with error bar given by one standard deviation\.LmaxL\_\{\\max\}MlM\_\{l\}TtotalT\_\{\\mathrm\{total\}\}\(s\)TparallelT\_\{\\mathrm\{parallel\}\}\(s\)FLOPs \(TFLOP\)Acc\. \(%\)NLLBrierECE22\[7\]\[7\]635\.089\.070\.3283\.51±0\.050\.470±0\.0000\.239±0\.0000\.023±0\.00133\[13,5\]\[13,5\]1633\.6710\.261\.2991\.71±0\.070\.568±0\.0050\.126±0\.0000\.020±0\.00044\[30,13,5\]\[30,13,5\]4240\.7210\.474\.5395\.34±0\.100\.376±0\.0020\.079±0\.0010\.024±0\.00155\[77,33,14,6\]\[77,33,14,6\]11259\.929\.4114\.6597\.06±0\.050\.241±0\.0010\.057±0\.0000\.035±0\.00166\[200,87,37,16,7\]\[200,87,37,16,7\]29716\.969\.4344\.3297\.33±0\.040\.204±0\.0030\.051±0\.0000\.035±0\.00077\[527,229,99,43,18,8\]\[527,229,99,43,18,8\]79113\.259\.46129\.1297\.54±0\.020\.204±0\.0030\.049±0\.0000\.036±0\.000

Table 2:Metrics for S\-MLPMMH on the MNIST dataset: The starting level isK=2K=2\. The reported metrics include accuracy, NLL, Brier score, ECE, total wall\-clock runtimeTtotalT\_\{\\mathrm\{total\}\}, pseudo\-parallel runtimeTparallelT\_\{\\mathrm\{parallel\}\}, and estimated FLOPs for different choices ofLmaxL\_\{\\max\}\. The QOI is the predictive probability\. Results are averaged over 3 replicates, with error bar gives by one standard deviation\.Overall, these experiments suggest that the TNN width\-decay structure is empirically compatible with effective neural network training\. The rank\-wise diagnostic shows that trained networks already exhibit substantial compressibility, while the index\-wise diagnostic shows that width\-decay turns this compressibility into an ordered, truncation\-friendly structure\. This provides empirical support for using the TNN prior as a structural prior in the MLPMMH construction\.

### 5\.3Practical application of MLPMMH

The original MLPMMH algorithm targets the joint posterior of the hyperparameters and model parameters by combining an outer pseudo\-marginal Markov chain with an inner multilevel SMC estimator\. Although this provides a principled approach to inference under strong parameter\-hyperparameter dependence, it is computationally expensive: each outer proposal requires a complete inner MLSMC run, while the outer Markov chain remains sequential\. These features make the method difficult to scale to large datasets and high\-dimensional feature spaces\. We therefore introduce practical modifications and evaluate the resulting approximation on the full MNIST dataset, using6000060000training and1000010000test images\.

#### 5\.3\.1Feature matrix from the pre\-trained LeNet NN

We first pre\-train the full LeNet architecture on the firstN𝗍𝗋N\_\{\\sf tr\}MNIST images using plain supervised learning\. The network comprises twoConv→ReLU→MaxPool\\text\{Conv\}\\to\\text\{ReLU\}\\to\\text\{MaxPool\}blocks followed by three fully–connected layers, the latter regularized with50%50\\%dropout\. All parameters are optimized for100100epochs with cross–entropy loss and the Adam optimiser \(initial learning rate10−310^\{\-3\}, halved every1010epochs\)\. Once training is complete, the parameters in the LeNet NN are fixed, i\.e\. a frozen LeNet NN\. The output of the second MaxPool layer is then flattened into a256256\-dimensional feature vector, and these vectors are collected for allN𝗍𝗋N\_\{\\sf tr\}training images\. The resultingN𝗍𝗋×256N\_\{\\sf tr\}\\times 256matrix of features is the feature matrix used for the subsequent Bayesian analysis\.

#### 5\.3\.2PCA hierarchy and multilevel coupling\.

To obtain the TNN\-style \(isotropic\) weight–decay and the nested sequence of models required for the multilevel scheme, we apply Principal Component Analysis \(PCA\) to theN𝗍𝗋×256N\_\{\\sf tr\}\\times 256training\-feature matrix and express both training and test features in this basis\. At levelll, the input layer width is set ton0=2ln\_\{0\}=2^\{l\}, i\.e\. we keep only the first2l2^\{l\}principal components \(PCs\)\. Followed the coupling construction introduced above, levell−1l\-1uses the first2l−12^\{l\-1\}PCs from the levelll, i\.e\. the fine and coarse likelihoods share their leading coordinates, yielding the strong correlation that drives multilevel variance reduction\. The prior,π¯l\\bar\{\\pi\}\_\{l\}, for the parameter in single Bayesian logistic\-regression layer is the Gaussian prior now, i\.e\.N⁡\(0,σA2​Im×n0\)N\(0,\\sigma\_\{A\}^\{2\}I\_\{m\\times n\_\{0\}\}\)andN⁡\(0,σb2​Im×1\)N\(0,\\sigma\_\{b\}^\{2\}I\_\{m\\times 1\}\)whereμ=\(σA,σb\)=\(σ,σ\)\\mu=\(\\sigma\_\{A\},\\sigma\_\{b\}\)=\(\\sigma,\\sigma\)\. The hyperparameter prior,π0\\pi\_\{0\}, is the Log\-Normal distribution\. In our implementation, we set𝔼⁡\[μ\]=128/2\\mathbb\{E\}\[\\mu\]=\\frac\{1\}\{\\sqrt\{2^\{8\}/2\}\}222Here28=2562^\{8\}=256is the full feature dimension\. The choice𝔼⁡\[μ\]=2/256\\mathbb\{E\}\[\\mu\]=\\sqrt\{2/256\}uses the He/Kaiming standard\-deviation scale as a convenient dimension\-dependent reference value\.,Var⁡\(μ\)=0\.25​𝔼​\[μ\]2\\operatorname\{Var\}\(\\mu\)=0\.25\\,\\mathbb\{E\}\[\\mu\]^\{2\}\. So, this gives the location parameter as−2\.54\-2\.54and scale parameter as0\.470\.47\.

#### 5\.3\.3Results

We then apply the practical S\-MLPMMH approximation introduced in Section[3\.2](https://arxiv.org/html/2610.19074#S3.SS2), using the PCA hierarchy described above to construct the multilevel estimators\.

Table[2](https://arxiv.org/html/2610.19074#S5.T2)reports the results for the practical MLPMMH method as the computational budget increases\. In these experiments,LmaxL\_\{\\max\}is increased gradually, together with the corresponding number of outer SBMC samplesMlM\_\{l\}at each levelll\. The runtime is reported in two forms: the total wall\-clock runtimeTtotalT\_\{\\mathrm\{total\}\}and the pseudo\-parallel runtimeTparallelT\_\{\\mathrm\{parallel\}\}\. In our implementation, the independent outer SBMC samples could be run in parallel, and the inner SBMC samples within each outer SBMC sample could also be parallelized\. Therefore, we estimate the parallel runtime asTparallel=maxl⁡TlN​MlT\_\{\\mathrm\{parallel\}\}=\\max\_\{l\}\\frac\{T\_\{l\}\}\{NM\_\{l\}\}, whereTlT\_\{l\}denotes the total runtime at levelll,NNis the number of inner SBMC samples, andMlM\_\{l\}is the number of outer SBMC samples at levelll\. We also report the estimated floating\-point cost in TFLOPs\.

Withdl=2ld\_\{l\}=2^\{l\}forl=2,…,7l=2,\\ldots,7, the level\-dependent Bayesian models are relatively small, and the idealized parallel runtime is therefore nearly constant across levels\. The estimated FLOPs, however, increase withLmaxL\_\{\\max\}, as expected from the theoretical sample allocation\. Consequently, this experiment is not well suited to demonstrating level\-dependent wall\-clock savings, although it illustrates the predictive behaviour of the practical multilevel approximation as the maximum level increases\. The ECE remains low forLmax≤4L\_\{\\max\}\\leq 4but increases fromLmax=5L\_\{\\max\}=5onward\. To investigate this behaviour, we also compute the signed calibration gap, defined asGap=1Ntest​∑i=1Ntest\(acci−confi\)\.\\operatorname\{Gap\}=\\frac\{1\}\{N\_\{\\mathrm\{test\}\}\}\\sum\_\{i=1\}^\{N\_\{\\mathrm\{test\}\}\}\\left\(\\operatorname\{acc\}\_\{i\}\-\\operatorname\{conf\}\_\{i\}\\right\)\.The gap is consistently positive across the considered values ofLmaxL\_\{\\max\}, indicating that the model is, on average, under\-confident\. This suggests conservative uncertainty estimates, which is a positive sign for reliability and robustness\.

## 6Conclusion

This paper develops and applies a Bayesian method for joint parameter and hyperparameter estimation for neural networks\. The methodology builds on multilevel Monte Carlo particle Metropolis methods\. and a particular neural network prior\. The prior structure is based on a trace\-class neural network, which imposes decay in the width direction and enjoys strong convergence\. We provide a theoretical complexity result relating MSE to computational cost, along with numerical verification\. Numerical experiments were conducted on both a Bayesian linear regression problem and on the MNIST data set for classification\. Our results demonstrate the effectiveness of our proposed method, both in terms the model prior and the algorithm used for posterior inference\. For the former, we compared our network architecture to that of standard training, and training where there is no regularization\. For moderate values of the decay parameterss, width\-decay regularization achieves predictive accuracy and NLL comparable to the unregularized and standard weight\-decay baselines, whereas excessively strong decay degrades predictive performance\. Separately, the synthetic regression experiments support the predicted multilevel regularity and the canonical MSE\-cost rate for the proposed MLPMMH method\. We also presented an SBMC\-based multilevel approximation using frozen LeNet features and a PCA hierarchy\. Since this practical variant replaces both the outer PMMH chain and the inner MLSMC sampler, it should be viewed as a scalable approximation rather than a direct verification of the exact MLPMMH\.

## Acknowledgements

NKC is supported by an EPSRC\-UKRI AI for Net Zero Grant: “Enabling CO2 Capture And Storage Projects Using AI”, \(grant EP/Y006143/1\)\. NKC is also supported by a City University of Hong Kong Start\-up Grant, project number 7200809, and a RGC\-ECS Grant, project number 21300826\. KJHL and XL gratefully acknowledge the support of IBM and EPSRC in the form of an Industrial Case Doctoral Studentship Award\.

## References

- D\. Agrawal, T\. Papamarkou, and J\. HinkleWide neural networks with bottlenecks are deep gaussian processes\.The Journal of Machine Learning Research21\(1\),pp\. 7056–7121\.Cited by:[§A\.2](https://arxiv.org/html/2610.19074#A1.SS2.p3.1)\.
- Andrieuet al\.\(2010\)C\. Andrieu, A\. Doucet, and R\. HolensteinParticle Markov chain Monte Carlo methods\.Journal of the Royal Statistical Society: Series B \(Statistical Methodology\)72\(3\),pp\. 269–342\.Cited by:[§1](https://arxiv.org/html/2610.19074#S1.p1.2),[Algorithm 5](https://arxiv.org/html/2610.19074#alg5)\.
- Andrieu and Roberts \(2009\)C\. Andrieu and G\. O\. RobertsThe pseudo\-marginal approach for efficient Monte Carlo computations\.Cited by:[§1](https://arxiv.org/html/2610.19074#S1.p1.2),[Algorithm 5](https://arxiv.org/html/2610.19074#alg5)\.
- BDL\-Consortium \(2026\)BDL\-ConsortiumHandbook of bayesian deep learning\.External Links:[Link](https://dx.doi.org/10.5281/zenodo.22086449)Cited by:[§1](https://arxiv.org/html/2610.19074#S1.p1.1)\.
- Beskoset al\.\(2018\)A\. Beskos, A\. Jasra, K\. J\. H\. Law, Y\. Marzouk, and Y\. ZhouMultilevel sequential Monte Carlo with dimension\-independent likelihood\-informed proposals\.SIAM/ASA Journal on Uncertainty Quantification6\(2\),pp\. 762–786\.Cited by:[§B\.1](https://arxiv.org/html/2610.19074#A2.SS1.p1.1),[§B\.2\.4](https://arxiv.org/html/2610.19074#A2.SS2.SSS4.p1.1),[§E\.2](https://arxiv.org/html/2610.19074#A5.SS2.p5.1.1),[§1](https://arxiv.org/html/2610.19074#S1.p2.1)\.
- Beskoset al\.\(2017\)A\. Beskos, A\. Jasra, K\. J\. H\. Law, R\. Tempone, and Y\. ZhouMultilevel sequential Monte Carlo samplers\.Stochastic Processes and their Applications127\(5\),pp\. 1417–1440\.Cited by:[§B\.2\.4](https://arxiv.org/html/2610.19074#A2.SS2.SSS4.p1.1),[§E\.2](https://arxiv.org/html/2610.19074#A5.SS2.p5.1.1),[§1](https://arxiv.org/html/2610.19074#S1.p2.1)\.
- Betancourt and Girolami \(2015\)M\. Betancourt and M\. GirolamiHamiltonian monte carlo for hierarchical models\.Current trends in Bayesian methodology with applications79\(30\),pp\. 2–4\.Cited by:[§1](https://arxiv.org/html/2610.19074#S1.p1.2)\.
- Chadaet al\.\(2021\)N\. K\. Chada, J\. Franks, A\. Jasra, K\. J\. Law, and M\. ViholaUnbiased inference for discretely observed hidden markov model diffusions\.SIAM/ASA Journal on Uncertainty Quantification9\(2\),pp\. 763–787\.Cited by:[Remark G\.2](https://arxiv.org/html/2610.19074#A7.Thmrem2),[§1](https://arxiv.org/html/2610.19074#S1.p2.1),[§3](https://arxiv.org/html/2610.19074#S3.p1.1)\.
- Chadaet al\.\(2025\)N\. K\. Chada, A\. Jasra, K\. J\. H\. Law, and S\. S\. SinghBayesian deep learning with multilevel trace\-class neural networks\.SIAM Journal on Mathematics of Data Science7\(3\),pp\. 1210–1240\.Cited by:[§A\.1](https://arxiv.org/html/2610.19074#A1.SS1.p2.1),[§B\.2\.4](https://arxiv.org/html/2610.19074#A2.SS2.SSS4.p1.1),[§B\.2](https://arxiv.org/html/2610.19074#A2.SS2.p1.1),[§E\.3](https://arxiv.org/html/2610.19074#A5.SS3.p1.1),[§E\.3](https://arxiv.org/html/2610.19074#A5.SS3.p2.1.1),[§1](https://arxiv.org/html/2610.19074#S1.p2.1),[§2\.1](https://arxiv.org/html/2610.19074#S2.SS1.p1.1)\.
- \[10\]N\. K\. Chada, B\. Leimkuhler, D\. Paulin, and P\. A\. WhalleyUnbiased kinetic langevin monte carlo with inexact gradients\.arXiv preprint arXiv:2311\.05025\.Cited by:[4th item](https://arxiv.org/html/2610.19074#A9.I1.i4.p1.1)\.
- Cho and Saul \(2009\)Y\. Cho and L\. SaulKernel methods for deep learning\.Advances in neural information processing systems22\.Cited by:[§A\.2](https://arxiv.org/html/2610.19074#A1.SS2.p3.1)\.
- Cotteret al\.\(2022\)S\. L\. Cotter, K\. J\. Law, X\. Liang, and S\. YangA randomized Multi\-index sequential Monte Carlo method\.arXiv preprint arXiv:2210\.15390\.Cited by:[§B\.2](https://arxiv.org/html/2610.19074#A2.SS2.p1.1),[§E\.2](https://arxiv.org/html/2610.19074#A5.SS2.p4.1.1),[§1](https://arxiv.org/html/2610.19074#S1.p2.1)\.
- Cotteret al\.\(2013\)S\. L\. Cotter, G\. O\. Roberts, A\. M\. Stuart, and D\. WhiteMCMC methods for functions: modifying old algorithms to make them faster\.Cited by:[§B\.1](https://arxiv.org/html/2610.19074#A2.SS1.p3.1)\.
- Del Moralet al\.\(2006\)P\. Del Moral, A\. Doucet, and A\. JasraSequential Monte Carlo samplers\.Journal of the Royal Statistical Society: Series B \(Statistical Methodology\)68\(3\),pp\. 411–436\.Cited by:[§B\.1](https://arxiv.org/html/2610.19074#A2.SS1.p1.1)\.
- Fujisawa and Sato \(2021\)M\. Fujisawa and I\. SatoMultilevel Monte Carlo variational inference\.J\. Mach\. Learn\. Res\.22,pp\. 1–44\.Cited by:[§1](https://arxiv.org/html/2610.19074#S1.p2.1)\.
- Gerstneret al\.\(2021\)T\. Gerstner, B\. Harrach, D\. Roth, and M\. SimonMultilevel Monte Carlo learning\.Preprint\.Cited by:[§1](https://arxiv.org/html/2610.19074#S1.p2.1)\.
- Geyer \(1992\)C\. J\. GeyerPractical markov chain monte carlo\.Statistical science,pp\. 473–483\.Cited by:[§B\.1](https://arxiv.org/html/2610.19074#A2.SS1.p3.1)\.
- Giles \(2008\)M\. B\. GilesMultilevel Monte Carlo path simulation\.Op\. Res\.56,pp\. 607–617\.Cited by:[§1](https://arxiv.org/html/2610.19074#S1.p2.1),[Algorithm 3](https://arxiv.org/html/2610.19074#alg3)\.
- Giles \(2015\)M\. B\. GilesMultilevel Monte Carlo methods\.Acta Numerica24,pp\. 259\.Cited by:[§1](https://arxiv.org/html/2610.19074#S1.p2.1),[Algorithm 3](https://arxiv.org/html/2610.19074#alg3)\.
- Golightley and Wilkinson \(2011\)A\. Golightley and D\. WilkinsonBayesian parameter inference for stochastic biochemical network models using particle Markov chain Monte Carlo\.Interface Focus1,pp\. 807–820\.Cited by:[§1](https://arxiv.org/html/2610.19074#S1.p1.2)\.
- Heinrich \(2001\)S\. HeinrichMultilevel Monte Carlo methods\.InInternational Conference on Large\-Scale Scientific Computing,pp\. 58–67\.Cited by:[§1](https://arxiv.org/html/2610.19074#S1.p2.1)\.
- Jasraet al\.\(2018\)A\. Jasra, K\. Kamatani, K\. J\. H\. Law, and Y\. ZhouBayesian static parameter estimation for partially observed diffusions via multilevel Monte Carlo\.SIAM Journal on Scientific Computing40\(2\),pp\. A887–A902\.Cited by:[§H\.4](https://arxiv.org/html/2610.19074#A8.SS4.p2.1.1),[§1](https://arxiv.org/html/2610.19074#S1.p2.1),[§3](https://arxiv.org/html/2610.19074#S3.p1.1),[Algorithm 1](https://arxiv.org/html/2610.19074#alg1)\.
- Jasraet al\.\(2017\)A\. Jasra, K\. Kamatani, K\. J\. Law, and Y\. ZhouMultilevel particle filters\.SIAM Journal on Numerical Analysis55\(6\),pp\. 3068–3096\.Cited by:[§E\.4](https://arxiv.org/html/2610.19074#A5.SS4.p2.1.1)\.
- Lawet al\.\(2022\)K\. J\. H\. Law, N\. Walton, S\. Yang, and A\. JasraMulti\-index Sequential Monte Carlo ratio estimators for bayesian inverse problems\.arXiv preprint arXiv:2203\.05351\.Cited by:[§B\.2](https://arxiv.org/html/2610.19074#A2.SS2.p1.1),[§E\.2](https://arxiv.org/html/2610.19074#A5.SS2.p3.1.1),[§G\.2](https://arxiv.org/html/2610.19074#A7.SS2.p1.1)\.
- Law \(2014\)K\. J\. LawProposals which speed up function\-space mcmc\.Journal of Computational and Applied Mathematics262,pp\. 127–138\.Cited by:[§B\.1](https://arxiv.org/html/2610.19074#A2.SS1.p3.1)\.
- Leeet al\.\(2018\)J\. Lee, Y\. Bahri, R\. Novak,et al\.Deep neural networks as gaussian processes\.International Conference on Learning Representations\.Cited by:[§A\.2](https://arxiv.org/html/2610.19074#A1.SS2.p3.1)\.
- Lianget al\.\(2025\)X\. Liang, J\. M\. Lukens, S\. Lohani, B\. T\. Kirby, T\. A\. Searles, X\. Qiu, and K\. J\. LawScalable bayesian monte carlo: fast uncertainty estimation beyond deep ensembles\.arXiv preprint arXiv:2505\.13585\.Cited by:[3rd item](https://arxiv.org/html/2610.19074#S1.I1.i3.p1.1),[§3\.2](https://arxiv.org/html/2610.19074#S3.SS2.p1.1)\.
- Lordet al\.\(2014\)G\. Lord, C\. E\. Powell, and T\. ShardlowAn introduction to computational stochastic pdes\.Cambridge Texts in Applied Mathematics\.Cited by:[§A\.2](https://arxiv.org/html/2610.19074#A1.SS2.p2.1)\.
- Moralet al\.\(2017\)P\. D\. Moral, A\. Jasra, K\. J\. H\. Law, and Y\. ZhouMultilevel sequential Monte Carlo samplers for normalizing constants\.ACM Transactions on Modeling and Computer Simulation \(TOMACS\)27\(3\),pp\. 1–22\.Cited by:[§E\.2](https://arxiv.org/html/2610.19074#A5.SS2.p2.1.1)\.
- Nealet al\.\(2011\)R\. M\. Nealet al\.MCMC using hamiltonian dynamics\.Handbook of markov chain monte carlo2\(11\),pp\. 2\.Cited by:[§1](https://arxiv.org/html/2610.19074#S1.p1.2)\.
- Neal \(2012\)R\. M\. NealBayesian learning for neural networks\.Springer Science & Business Media\.Cited by:[§A\.1](https://arxiv.org/html/2610.19074#A1.SS1.p2.1)\.
- O\. Lyeet al\.\(2021\)K\. O\. Lye, S\. Mishra, and R\. MolinaroA multi\-level procedure for enhancing accuracy of machine learning algorithms\.Euro\. Journal of Applied Mathematics32\(3\),pp\. 436–46\.Cited by:[§1](https://arxiv.org/html/2610.19074#S1.p2.1)\.
- Papamarkouet al\.\(2022\)T\. Papamarkou, J\. Hinkle, T\. Young, and D\. WombleChallenges in Markov chain Monte Carlo for bayesian neural networks\.Statistical Science37\(3\),pp\. 435–442\.Cited by:[§1](https://arxiv.org/html/2610.19074#S1.p1.1)\.
- Papamarkouet al\.\(2024\)T\. Papamarkou, M\. Skoularidou, K\. Palla, L\. Aitchison, J\. Arbel, D\. Dunson, M\. Filippone, V\. Fortuin, P\. Hennig, J\. M\. Hernández\-Lobato, A\. Hubin, A\. Immer, T\. Karaletsos, M\. E\. Khan, A\. Kristiadi, Y\. Li, S\. Mandt, C\. Nemeth, M\. A\. Osborne, T\. G\. J\. Rudner, D\. Rügamer, Y\. W\. Teh, M\. Welling, A\. G\. Wilson, and R\. ZhangPosition: bayesian deep learning is needed in the age of large\-scale AI\.InForty\-first International Conference on Machine Learning,External Links:[Link](https://openreview.net/forum?id=PrmxFWI1Fr)Cited by:[§1](https://arxiv.org/html/2610.19074#S1.p1.1)\.
- Paulinet al\.\(2025\)D\. Paulin, P\. A\. Whalley, N\. K\. Chada, and B\. LeimkuhlerSampling from bayesian neural network posteriors with symmetric minibatch splitting langevin dynamics\.Proceedings of The 28th International Conference on Artificial Intelligence and Statistics, PMLR258,pp\. 5014–5022\.Cited by:[4th item](https://arxiv.org/html/2610.19074#A9.I1.i4.p1.1)\.
- Robertet al\.\(1999\)C\. P\. Robert, G\. Casella, and G\. CasellaMonte Carlo statistical methods\.Vol\.2,Springer\.Cited by:[§B\.1](https://arxiv.org/html/2610.19074#A2.SS1.p3.1)\.
- Sell and Singh \(2020\)T\. Sell and S\. S\. SinghDimension\-robust function space mcmc with neural network priors\.Preprint\.Cited by:[§A\.2](https://arxiv.org/html/2610.19074#A1.SS2.p1.1),[§A\.2](https://arxiv.org/html/2610.19074#A1.SS2.p2.1),[2nd item](https://arxiv.org/html/2610.19074#A9.I1.i2.p1.1),[§1](https://arxiv.org/html/2610.19074#S1.p2.1)\.
- Shi and Cornish \(2021\)Y\. Shi and R\. CornishOn multilevel Monte Carlo unbiased gradient estimation for deep latent variable models\.Proceedings of The 24th International Conference on Artificial Intelligence and Statistics, PMLR130,pp\. 3925–3933\.Cited by:[§1](https://arxiv.org/html/2610.19074#S1.p2.1)\.

## Appendix AModel Formulation

In this section, we present and review our model setup for the paper\. This includes a review of our Bayesian posterior distribution of interest, while discussing this in the context of deep neural networks\. Related to this, we also discuss our prior, which is based on trace\-class neural networks\.

### A\.1Model Formulation

We now define the deep neural network model as follows\. LetD∈ℕD\\in\\mathbb\{N\},\(n0,…,nD\)∈ℕD\+1\(n\_\{0\},\.\.\.,n\_\{D\}\)\\in\\mathbb\{N\}^\{D\+1\}withn0=nn\_\{0\}=nfor the input layer andnD=mn\_\{D\}=mwith output layer\. Let weights beAd∈ℝnd×nd−1A\_\{d\}\\in\\mathbb\{R\}^\{n\_\{d\}\\times n\_\{d\-1\}\}and biases bebd∈ℝndb\_\{d\}\\in\\mathbb\{R\}^\{n\_\{d\}\}ford∈\{1,…,D\}d\\in\\\{1,\.\.\.,D\\\}, we denoteθ:=\(\(A1,b1\),…,\(AD,bD\)\)∈Θ=⊗d=1D\{ℝnd×nd−1×ℝnd\}\\theta:=\(\(A\_\{1\},b\_\{1\}\),\.\.\.,\(A\_\{D\},b\_\{D\}\)\)\\in\\Theta=\\otimes\_\{d=1\}^\{D\}\\\{\\mathbb\{R\}^\{n\_\{d\}\\times n\_\{d\-1\}\}\\times\\mathbb\{R\}^\{n\_\{d\}\}\\\}andμ\\muis the hyperparameter related toθ\\theta\. The layer is defined by

g1​\(x,θ,μ\)\\displaystyle g\_\{1\}\(x,\\theta,\\mu\):=A1​x\+b1,\\displaystyle:=A\_\{1\}x\+b\_\{1\},gd​\(x,θ,μ\)\\displaystyle g\_\{d\}\(x,\\theta,\\mu\):=Adσnd−1\(gd−1\(x\)\)\+bd,d∈\{2,…,D−1\},\\displaystyle:=A\_\{d\}\\sigma\_\{n\_\{d\-1\}\}\(g\_\{d\-1\}\(x\)\)\+b\_\{d\},\\ \\ d\\in\\\{2,\.\.\.,D\-1\\\},g⁡\(x,θ,μ\)\\displaystyle g\(x,\\theta,\\mu\):=AD​σnD−1​\(gD−1​\(x\)\)\+bD,\\displaystyle:=A\_\{D\}\\sigma\_\{n\_\{D\-1\}\}\(g\_\{D\-1\}\(x\)\)\+b\_\{D\},whereσk​\(z\):=\(ν⁡\(z1\),…,ν⁡\(zk\)\)T\\sigma\_\{k\}\(z\):=\(\\nu\(z\_\{1\}\),\.\.\.,\\nu\(z\_\{k\}\)\)^\{T\}with ReLU activationν⁡\(z\)=max⁡\{0,z\}\\nu\(z\)=\\max\\\{0,z\\\}\.

In order to infer the parameterθ\\thetaand hyperparameterμ\\mu, we establish the Bayesian model of the DNNs, which is called Bayesian DNNs\. Given a proper priorπ¯\\bar\{\\pi\}onΘ\\Thetaand a priorπ0\\pi\_\{0\}onμ\\mu, the joint posterior distribution of the Bayesian DNNs\[[Chada et al\., 2025](https://arxiv.org/html/2610.19074#bib.bib16),[Neal, 2012](https://arxiv.org/html/2610.19074#bib.bib14)\]can be computed up to a constant of proportionality

π⁡\(θ\)∝ℒ⁡\(θ,μ\)​π¯​\(θ\|μ\)​π0​\(μ\),\\pi\(\\theta\)\\propto\\mathcal\{L\}\(\\theta,\\mu\)\\bar\{\\pi\}\(\\theta\|\\mu\)\\pi\_\{0\}\(\\mu\),\(A\.1\)whereℒ⁡\(θ,μ\)\\mathcal\{L\}\(\\theta,\\mu\)is the joint likelihood onΘ×𝖲\\Theta\\times\\mathsf\{S\}\.

Define the quantity of interest through a mapΨ\\Psiapplied to the output\. That is,

φ⁡\(x,θ,μ\):=\(Ψ∘g\)​\(x,θ,μ\)=Ψ⁡\(g⁡\(x,θ,μ\)\)\.\\varphi\(x,\\theta,\\mu\):=\(\\Psi\\circ g\)\(x,\\theta,\\mu\)=\\Psi\(g\(x,\\theta,\\mu\)\)\.We then aim to compute the posterior expectation of the quantity of interest,

π​\(φ\):=𝔼π​\[φ\]\.\\pi\(\\varphi\):=\\mathbb\{E\}\_\{\\pi\}\[\\varphi\]\.\(A\.2\)

### A\.2Trace\-class neural network priors

The neural network priors onθ\\theta, the weights and bias, we consider are the trace class neural networks, first proposed by[Sell and Singh \[2020\]](https://arxiv.org/html/2610.19074#bib.bib15)\. These priors were introduced to mimic Gaussian priorsπ¯0∼𝒩⁡\(0,C\)\\overline\{\\pi\}\_\{0\}\\sim\\mathcal\{N\}\(0,C\), for function\-space inverse problems, where a common way to simulate Gaussian random fields is to use the Karhunen\-Loève expansion

θ=∑j∈ℤ\+λj​ξj​Φj,ξj∼𝒩⁡\(0,I\),\\theta=\\sum\_\{j\\in\\mathbb\{Z\}^\{\+\}\}\\sqrt\{\\lambda\_\{j\}\}\\xi\_\{j\}\\Phi\_\{j\},\\qquad\\xi\_\{j\}\\sim\\mathcal\{N\}\(0,I\),\(A\.3\)
where\(λj,Φj\)j∈ℤ\+\(\\lambda\_\{j\},\\Phi\_\{j\}\)\_\{j\\in\\mathbb\{Z\}^\{\+\}\}are the associated eigenbasis of the covariance functionCC, and\{ξj\}j∈ℤ\+\\\{\\xi\_\{j\}\\\}\_\{j\\in\\mathbb\{Z\}^\{\+\}\}is Gaussian white noise\. For a detail description of the derivation and its application to stochastic numerical problems, we refer the reader to[Lord et al\. \[2014\]](https://arxiv.org/html/2610.19074#bib.bib17)\. An issue that can arise with using \([A\.3](https://arxiv.org/html/2610.19074#A1.E3)\) is that using such priors do not scale well with high\-dimensional problems\. This acts as the initial motivation for trace\-class neural network \(TNN\) priors, where the work of[Sell and Singh \[2020\]](https://arxiv.org/html/2610.19074#bib.bib15)provided a justification for such priors\.

The hyperparametersσA,σb,s\\sigma\_\{A\},\\sigma\_\{b\},ssignificantly influence the inductive bias of the prior, and should be tuned for optimal performance\. TNN priors arenotbased on the Gaussian processes, which have been popular in a variety of Bayesian\-related tasks\. Recent literature has understood GP Neural Networks \(GPNNs\) in the context of width and depth limits of neural networks\[[Cho and Saul, 2009](https://arxiv.org/html/2610.19074#bib.bib4),[Agrawal et al\., 2020](https://arxiv.org/html/2610.19074#bib.bib1),[Lee et al\., 2018](https://arxiv.org/html/2610.19074#bib.bib20)\]\.

![Refer to caption](https://arxiv.org/html/2610.19074v1/image.png)Figure 4:Trace class nerual network architecture\.

## Appendix BState\-of\-the\-Art Estimation Methods for Bayesian Neural Networks

In the context of the Bayesian neural network with fixed hyperparametersμ\\muandss, the parameter estimation methods introduced here can be utilized for computing the marginal distribution ofθ\\thetaconditional on the hyperparameterμ\\mu\.

### B\.1Sequential Monte Carlo estimator

First, we look at the sequential Monte Carlo \(SMC\) estimator\[[Beskos et al\., 2018](https://arxiv.org/html/2610.19074#bib.bib12),[Del Moral et al\., 2006](https://arxiv.org/html/2610.19074#bib.bib7)\]\. Consider the conditioned posterior distribution,π⁡\(θ\|μ\)\\pi\(\\theta\|\\mu\), in \([A\.1](https://arxiv.org/html/2610.19074#A1.E1)\) with the fixed hyperparameters, we define the target distribution asπ⁡\(θ\|μ\)=f⁡\(θ\|μ\)/Z\\pi\(\\theta\|\\mu\)=f\(\\theta\|\\mu\)/Z, whereZ=∫Θf⁡\(θ\|μ\)​𝑑θZ=\\int\_\{\\Theta\}f\(\\theta\|\\mu\)d\\thetaandf⁡\(θ\|μ\):=ℒ⁡\(θ\|μ\)​π¯​\(θ\)f\(\\theta\|\\mu\):=\\mathcal\{L\}\(\\theta\|\\mu\)\\bar\{\\pi\}\(\\theta\)\. Given a quantity of interestφ:Θ→ℝ\\varphi:\\Theta\\rightarrow\\mathbb\{R\}, for simplicity, we define

f⁡\(φ\|μ\):=∫Θφ⁡\(θ\|μ\)​f​\(θ\|μ\)​𝑑θ=f⁡\(1\)​π​\(φ\|μ\),f\(\\varphi\|\\mu\):=\\int\_\{\\Theta\}\\varphi\(\\theta\|\\mu\)f\(\\theta\|\\mu\)d\\theta=f\(1\)\\pi\(\\varphi\|\\mu\),wheref⁡\(1\)=∫Θf⁡\(θ\|μ\)​𝑑θ=Zf\(1\)=\\int\_\{\\Theta\}f\(\\theta\|\\mu\)d\\theta=Z\. So that the conditional target distribution isπ⁡\(φ\)=f⁡\(φ\)f⁡\(1\)\\pi\(\\varphi\)=\\frac\{f\(\\varphi\)\}\{f\(1\)\}\.

Defineh1,…,hJ−1h\_\{1\},\.\.\.,h\_\{J\-1\}byhj=fj\+1/fjh\_\{j\}=f\_\{j\+1\}/f\_\{j\}, wheref1=π¯,fJ=f=ℒ​π¯f\_\{1\}=\\bar\{\\pi\},f\_\{J\}=f=\\mathcal\{L\}\\bar\{\\pi\}, and defineπj=fj/Zj\\pi\_\{j\}=f\_\{j\}/Z\_\{j\}whereZj=∫Θfj​𝑑θZ\_\{j\}=\\int\_\{\\Theta\}f\_\{j\}d\\theta\. Forj=2,…,Jj=2,\.\.\.,J, we letfjf\_\{j\}define in a annealing scheme:

fj=f11−λj​fJλj=ℒ​\(θ\)λj​π¯​\(θ\),f\_\{j\}=f\_\{1\}^\{1\-\\lambda\_\{j\}\}f\_\{J\}^\{\\lambda\_\{j\}\}=\\mathcal\{L\}\(\\theta\)^\{\\lambda\_\{j\}\}\\bar\{\\pi\}\(\\theta\),where0=λ1<λ2<⋯<λJ=10=\\lambda\_\{1\}<\\lambda\_\{2\}<\\cdots<\\lambda\_\{J\}=1andλj\\lambda\_\{j\}will be chosen adaptively according to the effective sample size \(ESS\), see Example[B\.1](https://arxiv.org/html/2610.19074#A2.Thmexample1)\.

###### Example B\.1\(Adaptive tempering\)\.

In order to keep a sufficient diversity of sample population, we let the effective sample size be at leastE​S​Smin=N/2ESS\_\{\\min\}=N/2at each temperingλj−1\\lambda\_\{j\-1\}and use it to compute the next temperingλj\\lambda\_\{j\}\. Forjjth tempering, we have weight samples\{wj−1k,θj−1k\}k=1N\\\{w^\{k\}\_\{j\-1\},\\theta^\{k\}\_\{j\-1\}\\\}\_\{k=1\}^\{N\}, then the ESS is computed by

E​S​S=1∑k=1N\(wj−1k\)2,ESS=\\frac\{1\}\{\\sum\_\{k=1\}^\{N\}\(w^\{k\}\_\{j\-1\}\)^\{2\}\},wherewj−1k=ℒ​\(θj−1k\)λj−λj−1/∑k=1Nℒ​\(θj−1k\)λj−λj−1w^\{k\}\_\{j\-1\}=\\mathcal\{L\}\(\\theta^\{k\}\_\{j\-1\}\)^\{\\lambda\_\{j\}\-\\lambda\_\{j\-1\}\}/\\sum\_\{k=1\}^\{N\}\\mathcal\{L\}\(\\theta\_\{j\-1\}^\{k\}\)^\{\\lambda\_\{j\}\-\\lambda\_\{j\-1\}\}\. LetΔ​λ=λj−λj−1\\Delta\\lambda=\\lambda\_\{j\}\-\\lambda\_\{j\-1\}, the effective sample size can be presented as a function ofΔ​λ\\Delta\\lambda, ESS\(Δ​λ\)\(\\Delta\\lambda\)\. Using suitable root finding method, one can findΔ​λ∗\\Delta\\lambda^\{\*\}such thatE​S​S​\(Δ​λ∗\)=E​S​SminESS\(\\Delta\\lambda^\{\*\}\)=ESS\_\{\\min\}, then set the next temperingλj=λj−1\+Δ​λ∗\\lambda\_\{j\}=\\lambda\_\{j\-1\}\+\\Delta\\lambda^\{\*\}\.

Algorithm 2SMC sampler1:Generate

θ1i∼π¯\\theta^\{i\}\_\{1\}\\sim\\bar\{\\pi\}for

i=1,…,Ni=1,\.\.\.,Nand

Z1N=1Z^\{N\}\_\{1\}=1\.

2:for

j=2j=2to

JJdo

3:Store

ZjN=Zj−1N​1N​∑k=1Nℒ​\(θj−1k\)λj−λj−1Z^\{N\}\_\{j\}=Z^\{N\}\_\{j\-1\}\\frac\{1\}\{N\}\\sum\_\{k=1\}^\{N\}\\mathcal\{L\}\(\\theta\_\{j\-1\}^\{k\}\)^\{\\lambda\_\{j\}\-\\lambda\_\{j\-1\}\}\.

4:for

i=1i=1to

NNdo

5:Define

wji=ℒ​\(θj−1i\)λj−λj−1/∑k=1Nℒ​\(θj−1k\)λj−λj−1w^\{i\}\_\{j\}=\\mathcal\{L\}\(\\theta\_\{j\-1\}^\{i\}\)^\{\\lambda\_\{j\}\-\\lambda\_\{j\-1\}\}/\\sum\_\{k=1\}^\{N\}\\mathcal\{L\}\(\\theta\_\{j\-1\}^\{k\}\)^\{\\lambda\_\{j\}\-\\lambda\_\{j\-1\}\}\.

6:Resample: Select

Iji∼\{wj1,…,wjN\}I\_\{j\}^\{i\}\\sim\\\{w\_\{j\}^\{1\},\.\.\.,w\_\{j\}^\{N\}\\\}and let

θ^ji=θj−1Iji\\hat\{\\theta\}\_\{j\}^\{i\}=\\theta\_\{j\-1\}^\{I\_\{j\}^\{i\}\}\.

7:Mutate: Draw

θji∼ℳj​\(θ^ji,⋅\)\\theta\_\{j\}^\{i\}\\sim\\mathcal\{M\}\_\{j\}\(\\hat\{\\theta\}\_\{j\}^\{i\},\\cdot\)\.

8:endfor

9:endfor

Now letℳj\\mathcal\{M\}\_\{j\}forj=2,…,Jj=2,\\dots,Jbe Markov transition kernels such that\(πj​ℳj\)​\(d​θ\)=πj​\(d​θ\)\(\\pi\_\{j\}\\mathcal\{M\}\_\{j\}\)\(d\\theta\)=\\pi\_\{j\}\(d\\theta\), analogous to any suitable MCMC kernel\[[Geyer, 1992](https://arxiv.org/html/2610.19074#bib.bib32),[Robert et al\., 1999](https://arxiv.org/html/2610.19074#bib.bib11),[Cotter et al\., 2013](https://arxiv.org/html/2610.19074#bib.bib33),[Law, 2014](https://arxiv.org/html/2610.19074#bib.bib34)\]\. This operation must sufficiently decorrelate the samples, and as such, we typically define the MCMC kernelsℳj\{\\mathcal\{M\}\_\{j\}\}by a suitable number of steps of some basic MCMC kernel, where it may or may not be 1\.

One can observe that

πj​\(hj\)=1Zj​∫fj\+1​\(θ\)fj​\(θ\)​fj​\(θ\)​𝑑θ=Zj\+1Zj\.\\pi\_\{j\}\(h\_\{j\}\)=\\frac\{1\}\{Z\_\{j\}\}\\int\\frac\{f\_\{j\+1\}\(\\theta\)\}\{f\_\{j\}\(\\theta\)\}f\_\{j\}\(\\theta\)d\\theta=\\frac\{Z\_\{j\+1\}\}\{Z\_\{j\}\}\.Since we letπJ=π\\pi\_\{J\}=\\pibe the target distribution andZ1=1Z\_\{1\}=1, the normalization constant ofπ\\piis given by

Z=Z1​∏j=1J−1πj​\(hj\)=∏j=1J−1πj​\(hj\)\.Z=Z\_\{1\}\\prod\_\{j=1\}^\{J\-1\}\\pi\_\{j\}\(h\_\{j\}\)=\\prod\_\{j=1\}^\{J\-1\}\\pi\_\{j\}\(h\_\{j\}\)\.
From Algorithm[2](https://arxiv.org/html/2610.19074#alg2), we can define the following estimators:

πjN​\(hj\)\\displaystyle\\pi\_\{j\}^\{N\}\(h\_\{j\}\):=1N∑i=1Nhj\(θji\),j=1,…,J−1,\\displaystyle:=\\frac\{1\}\{N\}\\sum\_\{i=1\}^\{N\}h\_\{j\}\(\\theta\_\{j\}^\{i\}\),\\ j=1,\.\.\.,J\-1,\(B\.1\)ZN\\displaystyle Z^\{N\}:=∏j=1J−1πjN​\(hj\),\\displaystyle:=\\prod\_\{j=1\}^\{J\-1\}\\pi\_\{j\}^\{N\}\(h\_\{j\}\),fN​\(φ\)\\displaystyle f^\{N\}\(\\varphi\):=∏j=1J−1πjN​\(hj\)​πN​\(φ\)=ZN​πN​\(φ\)\.\\displaystyle:=\\prod\_\{j=1\}^\{J\-1\}\\pi\_\{j\}^\{N\}\(h\_\{j\}\)\\pi^\{N\}\(\\varphi\)=Z^\{N\}\\pi^\{N\}\(\\varphi\)\.Hence, the sequential Monte Carlo estimator for estimating the expectation of the quantity of interest is given by

φ^SMC:=πN​\(φ\)=fN​\(φ\)fN​\(1\)\.\\hat\{\\varphi\}\_\{\\text\{SMC\}\}:=\\pi^\{N\}\(\\varphi\)=\\frac\{f^\{N\}\(\\varphi\)\}\{f^\{N\}\(1\)\}\.\(B\.2\)

### B\.2Novel Multilevel sequential Monte Carlo estimators

Integrating a multilevel scheme into the single SMC estimator further decreases the variance in the estimator, yielding particularly more gains in high\-dimensional settings\. In this section, we restate three multilevel Monte Carlo methods: MLSMC\[[Law et al\., 2022](https://arxiv.org/html/2610.19074#bib.bib5)\], randomized MLSMC\[[Cotter et al\., 2022](https://arxiv.org/html/2610.19074#bib.bib27)\]ratio estimators, and the original serial MLSMC\[[Chada et al\., 2025](https://arxiv.org/html/2610.19074#bib.bib16)\]\.

#### B\.2\.1Multilevel Structure

There are two fundamentally different ways used here to couple across levels: the “top‐down” construction and the “bottom‐up”\.

- •”top\-down”: one starts with a fine‐level sample at levellland then constructs its coarse analogue at levell−1l\-1, typically by truncating parameters in the data, inputs, or network width \(as in \(randomized\) MLSMC estimators\)\.
- •”bottom\-up”: one first samples a coarse\-level parameterθl−1∈Θl−1\\theta\_\{l\-1\}\\in\\Theta\_\{l\-1\}, and then extends or propagates it to a fine\-level parameterθl∈Θl\\theta\_\{l\}\\in\\Theta\_\{l\}\. This direction is used later for the serial MLSMC estimator\.

In what follows, unless otherwise stated, the approximate couplingΠ¯l\\bar\{\\Pi\}\_\{l\}refers to the top\-down truncation coupling\.

##### Paralleization\.

This distinction also affects the available parallelization\. In the ”top\-down” construction in \(randomized\) MLSMC, the estimators at different levels can be run independently, and within each level the SMC particles can be propagated and evaluated in parallel between resampling steps\. By contrast, the ”bottom\-up” construction in serial MLSMC has a stronger dependence across levels, since the fine\-level sample is obtained by propagating information from the previous coarse level\. Hence this construction is generally less parallel across levels, although particle\-wise operations and independent replications can still be parallelized\.

##### Approximate Coupling\.

At each levelll, we intend to simulate a pair of samples\(θli,θl−1i\)\(\\theta\_\{l\}^\{i\},\\theta\_\{l\-1\}^\{i\}\)follows\(θli,θl−1i\)∼Πl​\(θl,θl−1\|μ\)\(\\theta\_\{l\}^\{i\},\\theta\_\{l\-1\}^\{i\}\)\\sim\\Pi\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\), whereΠl​\(θl,θl−1\|μ\)\\Pi\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)is the exact coupling distribution such that

∫ΘlΠl​\(θl,θl−1\|μ\)​d​θl=πl−1​\(θl−1\|μ\),∫Θl−1Πl​\(θl,θl−1\|μ\)​d​θl−1=πl​\(θl\|μ\)\.\\int\_\{\\Theta\_\{l\}\}\\Pi\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)d\\theta\_\{l\}=\\pi\_\{l\-1\}\(\\theta\_\{l\-1\}\|\\mu\),\\quad\\int\_\{\\Theta\_\{l\-1\}\}\\Pi\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)d\\theta\_\{l\-1\}=\\pi\_\{l\}\(\\theta\_\{l\}\|\\mu\)\.However, it is not a trivial approach to compute the exact coupling; instead, we introduce the approximate coupling technique\. Following the same procedure in rewriting the posterior distribution at the single level in SMC, we have

There are two fundamentally different ways used here to couple across levels\. The first runs “top‐down”: one starts with a fine‐level sample at levellland then constructs its coarse analogue at levell−1l\-1, typically by truncating parameters in the data, inputs, or network width \(as in \(randomized\) MLSMC estimators\)\. The second runs “bottom‐up”: one begins with a coarse‐level sample at levell−1l\-1and serially propagates it to levelll\(as in the serial MLSMC estimator\)\.

Since the dataxxis the same and the hyperparameters are fixed, we omit thexxandμ\\muin the notations below\.

#### B\.2\.2Multilevel sequential Monte Carlo estimators

Given a starting levelLmin≥1L\_\{\\text\{min\}\}\\geq 1and a maximum levelLmax≥LminL\_\{\\text\{max\}\}\\geq L\_\{\\text\{min\}\}, the approximation of the expectation \([A\.2](https://arxiv.org/html/2610.19074#A1.E2)\) is given by

π⁡\(φ\)≈∑l=LminLmaxπl​\(Δ​φ\)\\displaystyle\\pi\(\\varphi\)\\approx\\sum\_\{l=L\_\{\\text\{min\}\}\}^\{L\_\{\\text\{max\}\}\}\\pi\_\{l\}\(\\Delta\\varphi\):=∑l=LminLmaxπl​\(φl​\(θl\)−φl−1​\(θl−1\)\)≈∑l=LminLmax1Nl​∑i=1Nl\{φl​\(θli\)−φl−1​\(θl−1i\)\},\\displaystyle:=\\sum\_\{l=L\_\{\\text\{min\}\}\}^\{L\_\{\\text\{max\}\}\}\\pi\_\{l\}\(\\varphi\_\{l\}\(\\theta\_\{l\}\)\-\\varphi\_\{l\-1\}\(\\theta\_\{l\-1\}\)\)\\approx\\sum\_\{l=L\_\{\\text\{min\}\}\}^\{L\_\{\\text\{max\}\}\}\\frac\{1\}\{N\_\{l\}\}\\sum\_\{i=1\}^\{N\_\{l\}\}\\Big\\\{\\varphi\_\{l\}\(\\theta\_\{l\}^\{i\}\)\-\\varphi\_\{l\-1\}\(\\theta\_\{l\-1\}^\{i\}\)\\Big\\\},whereθli\\theta^\{i\}\_\{l\}are i\.i\.d\. samples at each levelll,NlN\_\{l\}is the number of samples at levelllandφLmin−1:=0\\varphi\_\{L\_\{\\text\{min\}\}\-1\}:=0\.

Here, we use the sequential Monte Carlo sampler, stated in Algorithm[2](https://arxiv.org/html/2610.19074#alg2), to sample from the approximate coupling distribution at each level, shown in Algorithm[16](https://arxiv.org/html/2610.19074#alg3.l16)\. We now define intermediate distributions based on the coupling distribution,FlF\_\{l\},𝐋l\{\\bf L\}\_\{l\}andΠ¯l\\bar\{\\Pi\}\_\{l\}\. DefineHl,1,…,Hl,J−1H\_\{l,1\},\.\.\.,H\_\{l,J\-1\}byHl,j=Fl,j\+1/Fl,jH\_\{l,j\}=F\_\{l,j\+1\}/F\_\{l,j\}whereFl,j=𝐋lλj​Π¯lF\_\{l,j\}=\{\\bf L\}\_\{l\}^\{\\lambda\_\{j\}\}\\bar\{\\Pi\}\_\{l\}with0=λ1<λ2<⋯<λJ=10=\\lambda\_\{1\}<\\lambda\_\{2\}<\\cdots<\\lambda\_\{J\}=1333Use the same adaptive tempering technique as in the SMC estimator shown in Example[B\.1](https://arxiv.org/html/2610.19074#A2.Thmexample1)\. And letℳl,j\\mathcal\{M\}\_\{l,j\}forj=1,…,Jj=1,\.\.\.,Jbe such that\(Πl,j​ℳl,j\)​\(d​θ\)=Πl,j​\(d​θ\)\(\\Pi\_\{l,j\}\\mathcal\{M\}\_\{l,j\}\)\(d\\theta\)=\\Pi\_\{l,j\}\(d\\theta\)\.

Forj=1,…,Jj=1,\.\.\.,J,

Πl,jNl​\(d​θ\):=1Nl​∑i=1Nlδθji​\(d​θ\),\\Pi\_\{l,j\}^\{N\_\{l\}\}\(d\\theta\):=\\frac\{1\}\{N\_\{l\}\}\\sum\_\{i=1\}^\{N\_\{l\}\}\\delta\_\{\\theta\_\{j\}^\{i\}\}\(d\\theta\),and then, we define

ZlNl:=∏j=1J−1Πl,jNl​\(Hl,j\),FlNl​\(d​θ\):=ZlNl​Πl,JNl​\(d​θ\)\.Z\_\{l\}^\{N\_\{l\}\}:=\\prod\_\{j=1\}^\{J\-1\}\\Pi\_\{l,j\}^\{N\_\{l\}\}\(H\_\{l,j\}\),\\ \\ F\_\{l\}^\{N\_\{l\}\}\(d\\theta\):=Z\_\{l\}^\{N\_\{l\}\}\\Pi\_\{l,J\}^\{N\_\{l\}\}\(d\\theta\)\.
1:Simulate

θl,1i∼π¯l\\theta\_\{l,1\}^\{i\}\\sim\\bar\{\\pi\}\_\{l\}for

i=1,…,Nli=1,\\ldots,N\_\{l\}and set

θl−1,1i=Tl​\(θl,1i\)\\theta\_\{l\-1,1\}^\{i\}=T\_\{l\}\(\\theta\_\{l,1\}^\{i\}\)\.

2:Set

Zl,1Nl=1Z\_\{l,1\}^\{N\_\{l\}\}=1
3:Compute the likelihoods

ℒl,1​\(θl,1i\),ℒl−1,1​\(θl−1,1i\)\\mathcal\{L\}\_\{l,1\}\(\\theta\_\{l,1\}^\{i\}\),\\mathcal\{L\}\_\{l\-1,1\}\(\\theta\_\{l\-1,1\}^\{i\}\)
4:Store the approximate coupling

Fl,1​\(θl,1i,θl−1,1i\)=Π¯l​\(θl,1i,θl−1,1i\)F\_\{l,1\}\(\\theta\_\{l,1\}^\{i\},\\theta\_\{l\-1,1\}^\{i\}\)=\\bar\{\\Pi\}\_\{l\}\(\\theta\_\{l,1\}^\{i\},\\theta\_\{l\-1,1\}^\{i\}\)
5:for

j=2,…,Jj=2,\.\.\.,Jdo

6:Store

Zl,jNl=Zl,j−1Nl​1Nl​∑p=1NlHl,j−1​\(θl,j−1p\)Z\_\{l,j\}^\{N\_\{l\}\}=Z\_\{l,j\-1\}^\{N\_\{l\}\}\\frac\{1\}\{N\_\{l\}\}\\sum\_\{p=1\}^\{N\_\{l\}\}H\_\{l,j\-1\}\(\\theta\_\{l,j\-1\}^\{p\}\)
7:for

i=1,…,Nli=1,\.\.\.,N\_\{l\}do

8:Define

wl,ji=Hl,j−1​\(θl,j−1i\)∑p=1NlHl,j−1​\(θl,j−1p\)w\_\{l,j\}^\{i\}=\\frac\{H\_\{l,j\-1\}\(\\theta^\{i\}\_\{l,j\-1\}\)\}\{\\sum\_\{p=1\}^\{N\_\{l\}\}H\_\{l,j\-1\}\(\\theta^\{p\}\_\{l,j\-1\}\)\}
9:Resample: Select

Il,ji∼\{wl,j1,…,wl,jNl\}I\_\{l,j\}^\{i\}\\sim\\\{w\_\{l,j\}^\{1\},\\ldots,w\_\{l,j\}^\{N\_\{l\}\}\\\}, and set

θ^l,ji=θl,j−1Il,ji\\hat\{\\theta\}\_\{l,j\}^\{i\}=\\theta\_\{l,j\-1\}^\{I\_\{l,j\}^\{i\}\},

θ^l−1,ji=θl−1,j−1Il,ji\\hat\{\\theta\}\_\{l\-1,j\}^\{i\}=\\theta\_\{l\-1,j\-1\}^\{I\_\{l,j\}^\{i\}\}\.

10:Mutate: Draw

θl,ji∼ℳl,j​\(θ^l,ji,⋅\)\\theta\_\{l,j\}^\{i\}\\sim\\mathcal\{M\}\_\{l,j\}\(\\hat\{\\theta\}\_\{l,j\}^\{i\},\\cdot\)and set

θl−1,ji=Tl​\(θl,ji\)\\theta\_\{l\-1,j\}^\{i\}=T\_\{l\}\(\\theta\_\{l,j\}^\{i\}\)\.

11:Compute the likelihoods

ℒl,j​\(θl,ji\)\\mathcal\{L\}\_\{l,j\}\(\\theta\_\{l,j\}^\{i\}\)and

ℒl−1,j​\(θl−1,ji\)\\mathcal\{L\}\_\{l\-1,j\}\(\\theta\_\{l\-1,j\}^\{i\}\)\.

12:Store the approximate coupling

Fl,j​\(θl,ji,θl−1,ji\)=Π¯l​\(θl,ji,θl−1,ji\)​𝐋l​\(θl,ji,θl−1,ji\)λjF\_\{l,j\}\(\\theta\_\{l,j\}^\{i\},\\theta\_\{l\-1,j\}^\{i\}\)=\\bar\{\\Pi\}\_\{l\}\(\\theta\_\{l,j\}^\{i\},\\theta\_\{l\-1,j\}^\{i\}\)\{\\bf L\}\_\{l\}\(\\theta\_\{l,j\}^\{i\},\\theta\_\{l\-1,j\}^\{i\}\)^\{\\lambda\_\{j\}\}
13:where

𝐋l​\(θl,ji,θl−1,ji\)=max⁡\{ℒl,j​\(θl,ji\),ℒl−1,j​\(θl−1,ji\)\}\{\\bf L\}\_\{l\}\(\\theta\_\{l,j\}^\{i\},\\theta\_\{l\-1,j\}^\{i\}\)=\\max\\\{\\mathcal\{L\}\_\{l,j\}\(\\theta\_\{l,j\}^\{i\}\),\\mathcal\{L\}\_\{l\-1,j\}\(\\theta\_\{l\-1,j\}^\{i\}\)\\\}
14:endfor

15:endfor

16:Define

θli=θl,Ji\\theta\_\{l\}^\{i\}=\\theta\_\{l,J\}^\{i\},

θl−1i=θl−1,Ji\\theta\_\{l\-1\}^\{i\}=\\theta\_\{l\-1,J\}^\{i\}\. Store

fl​\(θli\)=π¯l​\(θli\)​ℒl,J​\(θli\)f\_\{l\}\(\\theta\_\{l\}^\{i\}\)=\\bar\{\\pi\}\_\{l\}\(\\theta\_\{l\}^\{i\}\)\\mathcal\{L\}\_\{l,J\}\(\\theta\_\{l\}^\{i\}\),

fl−1​\(θl−1i\)=π¯l−1​\(θl−1i\)​ℒl−1,J​\(θl−1i\)f\_\{l\-1\}\(\\theta\_\{l\-1\}^\{i\}\)=\\bar\{\\pi\}\_\{l\-1\}\(\\theta\_\{l\-1\}^\{i\}\)\\mathcal\{L\}\_\{l\-1,J\}\(\\theta\_\{l\-1\}^\{i\}\), and

Fl​\(θli,θl−1i\)=Fl,J​\(θl,Ji,θl−1,Ji\)F\_\{l\}\(\\theta\_\{l\}^\{i\},\\theta\_\{l\-1\}^\{i\}\)=F\_\{l,J\}\(\\theta\_\{l,J\}^\{i\},\\theta\_\{l\-1,J\}^\{i\}\)\.

Algorithm 3Sequential Monte Carlo\[[Giles, 2008](https://arxiv.org/html/2610.19074#bib.bib22),[Giles, 2015](https://arxiv.org/html/2610.19074#bib.bib6)\]sampler for coupled estimation at levelll##### Self\-normalized estimator

The trivial strategy for constructing the MLSMC estimator is to first normalized itself at each level and then compute the telescope sum\. Then, the MLSMC self\-normalized estimator is given as

φ^LMLSN\\displaystyle\\hat\{\\varphi\}^\{\\text\{MLSN\}\}\_\{L\}=∑l=LminLmax\(FlNl​\(φl\)FlNl​\(1\)−Fl−1Nl​\(φl−1\)Fl−1Nl​\(1\)\)\\displaystyle=\\sum\_\{l=L\_\{\\text\{min\}\}\}^\{L\_\{\\text\{max\}\}\}\\Bigg\(\\frac\{F^\{N\_\{l\}\}\_\{l\}\(\\varphi\_\{l\}\)\}\{F^\{N\_\{l\}\}\_\{l\}\(1\)\}\-\\frac\{F^\{N\_\{l\}\}\_\{l\-1\}\(\\varphi\_\{l\-1\}\)\}\{F^\{N\_\{l\}\}\_\{l\-1\}\(1\)\}\\Bigg\)\(B\.3\)=∑l=LminLmax\(∑i=1Nlφl​\(θli\)​fl​\(θli\)Fl​\(θli,θl−1i\)∑i=1Nlfl​\(θli\)Fl​\(θli,θl−1i\)−∑i=1Nlφl−1​\(θl−1i\)​fl−1​\(θl−1i\)Fl​\(θli,θl−1i\)∑i=1Nlfl−1​\(θl−1i\)Fl​\(θli,θl−1i\)\)\.\\displaystyle=\\sum\_\{l=L\_\{\\text\{min\}\}\}^\{L\_\{\\text\{max\}\}\}\\Bigg\(\\frac\{\\sum\_\{i=1\}^\{N\_\{l\}\}\\varphi\_\{l\}\(\\theta\_\{l\}^\{i\}\)\\frac\{f\_\{l\}\(\\theta\_\{l\}^\{i\}\)\}\{F\_\{l\}\(\\theta\_\{l\}^\{i\},\\theta\_\{l\-1\}^\{i\}\)\}\}\{\\sum\_\{i=1\}^\{N\_\{l\}\}\\frac\{f\_\{l\}\(\\theta\_\{l\}^\{i\}\)\}\{F\_\{l\}\(\\theta\_\{l\}^\{i\},\\theta\_\{l\-1\}^\{i\}\)\}\}\-\\frac\{\\sum\_\{i=1\}^\{N\_\{l\}\}\\varphi\_\{l\-1\}\(\\theta\_\{l\-1\}^\{i\}\)\\frac\{f\_\{l\-1\}\(\\theta\_\{l\-1\}^\{i\}\)\}\{F\_\{l\}\(\\theta\_\{l\}^\{i\},\\theta\_\{l\-1\}^\{i\}\)\}\}\{\\sum\_\{i=1\}^\{N\_\{l\}\}\\frac\{f\_\{l\-1\}\(\\theta\_\{l\-1\}^\{i\}\)\}\{F\_\{l\}\(\\theta\_\{l\}^\{i\},\\theta\_\{l\-1\}^\{i\}\)\}\}\\Bigg\)\.

##### Ratio estimator

Besides the self\-normalized estimator, a ratio estimator is given below\. Letζ=\{1,φ\}\\zeta=\\\{1,\\varphi\\\}, Define the difference functionΔ​ζ:Θ→ℝm\\Delta\\zeta:\\Theta\\rightarrow\\mathbb\{R\}^\{m\}as

Δ​ζ=ζl​\(θl\)​fl​\(θl\)Fl​\(θl,θl−1\)−ζl−1​\(θl−1\)​fl−1​\(θl−1\)Fl​\(θl,θl−1\),\\Delta\\zeta=\\zeta\_\{l\}\(\\theta\_\{l\}\)\\frac\{f\_\{l\}\(\\theta\_\{l\}\)\}\{F\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\)\}\-\\zeta\_\{l\-1\}\(\\theta\_\{l\-1\}\)\\frac\{f\_\{l\-1\}\(\\theta\_\{l\-1\}\)\}\{F\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\)\},and define the MLSMC ratio estimator as

φ^LMLRE\\displaystyle\\hat\{\\varphi\}^\{\\text\{MLRE\}\}\_\{L\}=∑l=LminLmaxFlNl​\(Δ​φ\)max⁡\{∑l=LminLmaxFlNl​\(Δ​1\),Zmin\}\\displaystyle=\\frac\{\\sum\_\{l=L\_\{\\text\{min\}\}\}^\{L\_\{\\text\{max\}\}\}F^\{N\_\{l\}\}\_\{l\}\(\\Delta\\varphi\)\}\{\\max\\\{\\sum\_\{l=L\_\{\\text\{min\}\}\}^\{L\_\{\\text\{max\}\}\}F^\{N\_\{l\}\}\_\{l\}\(\\Delta 1\),\\ Z\_\{\\min\}\\\}\}\(B\.4\)=∑l=LminLmaxZlNl​\(1Nl​∑i=1Nlφl​\(θli\)​fl​\(θli\)Fl​\(θli,θl−1i\)−φl−1​\(θl−1i\)​fl−1​\(θl−1i\)Fl​\(θli,θl−1i\)\)max⁡\{∑l=LminLmaxZlNl​\(1Nl​∑i=1Nlfl​\(θli\)Fl​\(θli,θl−1i\)−fl−1​\(θl−1i\)Fl​\(θli,θl−1i\)\),Zmin\}\.\\displaystyle=\\frac\{\\sum\_\{l=L\_\{\\text\{min\}\}\}^\{L\_\{\\text\{max\}\}\}Z\_\{l\}^\{N\_\{l\}\}\\Big\(\\frac\{1\}\{N\_\{l\}\}\\sum\_\{i=1\}^\{N\_\{l\}\}\\varphi\_\{l\}\(\\theta\_\{l\}^\{i\}\)\\frac\{f\_\{l\}\(\\theta\_\{l\}^\{i\}\)\}\{F\_\{l\}\(\\theta\_\{l\}^\{i\},\\theta\_\{l\-1\}^\{i\}\)\}\-\\varphi\_\{l\-1\}\(\\theta\_\{l\-1\}^\{i\}\)\\frac\{f\_\{l\-1\}\(\\theta\_\{l\-1\}^\{i\}\)\}\{F\_\{l\}\(\\theta\_\{l\}^\{i\},\\theta\_\{l\-1\}^\{i\}\)\}\\Big\)\}\{\\max\\\{\\sum\_\{l=L\_\{\\text\{min\}\}\}^\{L\_\{\\text\{max\}\}\}Z\_\{l\}^\{N\_\{l\}\}\\Big\(\\frac\{1\}\{N\_\{l\}\}\\sum\_\{i=1\}^\{N\_\{l\}\}\\frac\{f\_\{l\}\(\\theta\_\{l\}^\{i\}\)\}\{F\_\{l\}\(\\theta\_\{l\}^\{i\},\\theta\_\{l\-1\}^\{i\}\)\}\-\\frac\{f\_\{l\-1\}\(\\theta\_\{l\-1\}^\{i\}\)\}\{F\_\{l\}\(\\theta\_\{l\}^\{i\},\\theta\_\{l\-1\}^\{i\}\)\}\\Big\),Z\_\{\\min\}\\\}\}\.
It is worth noting that when computing the likelihoods in practice, it may occur computationally small number that the computer will recognize as 0\. However, the computer can store a much smaller number if one takes the logarithm of it\. So, here we store and run with the log\-likelihood instead of the likelihood and take the exponential of it when needed\. For the small number of IS weights, we work with the log\-likelihood for unnormalized weights and modify it with a constant\.

The practical settings for the MLSMC estimators are shown in Appendix[G\.1](https://arxiv.org/html/2610.19074#A7.SS1)\.

#### B\.2\.3Randomized multilevel sequential Monte Carlo ratio estimator

In randomized method, give a starting levelLminL\_\{\\text\{min\}\}as before, and instead of giving a maximum levelLmaxL\_\{\\text\{max\}\}in advance, we choose the sample of levels randomly from a distribution from numbers larger or equal thanLminL\_\{\\text\{min\}\}\. Let𝐩\{\\bf p\}be a probability measure onℙ=\{l∈ℕ:l≥Lmin\}\\mathbb\{P\}=\\\{l\\in\\mathbb\{N\}:l\\geq L\_\{\\text\{min\}\}\\\}, wherepl:=𝐩⁡\(l\)p\_\{l\}:=\{\\bf p\}\(l\)\. We simulateNNi\.i\.d samplesli∼𝐩l\_\{i\}\\sim\{\\bf p\}and define the number of samples at levelllbyNl=\#⁡\{i:li=l\}N\_\{l\}=\\\#\\\{i:l\_\{i\}=l\\\}for anyl∈ℙl\\in\\mathbb\{P\}\. Then, the posterior is

π⁡\(φ\)=∑l∈ℙπl​\(Δ​φ\)≈∑l∈ℙ1pl​N​∑i=1Nl\{φl​\(θli\)−φl−1​\(θl−1i\)\},\\displaystyle\\pi\(\\varphi\)=\\sum\_\{l\\in\\mathbb\{P\}\}\\pi\_\{l\}\(\\Delta\\varphi\)\\approx\\sum\_\{l\\in\\mathbb\{P\}\}\\frac\{1\}\{p\_\{l\}N\}\\sum\_\{i=1\}^\{N\_\{l\}\}\\Big\\\{\\varphi\_\{l\}\(\\theta\_\{l\}^\{i\}\)\-\\varphi\_\{l\-1\}\(\\theta\_\{l\-1\}^\{i\}\)\\Big\\\},whereθli\\theta^\{i\}\_\{l\}are i\.i\.d\. samples at each levelllandφLmin−1:=0\\varphi\_\{L\_\{\\text\{min\}\}\-1\}:=0\.

Following the same approximate technique and SMC sampler on coupling distributions in Algorithm[16](https://arxiv.org/html/2610.19074#alg3.l16), The randomized MLSMC ratio estimator is defined as

φ^rML\\displaystyle\\hat\{\\varphi\}^\{\\text\{rML\}\}=∑l∈ℙNlpl​N​FlNl​\(Δ​φ\)max⁡\{∑l∈ℙNlpl​N​FlNl​\(Δ​1\),Zmin\}\\displaystyle=\\frac\{\\sum\_\{l\\in\\mathbb\{P\}\}\\frac\{N\_\{l\}\}\{p\_\{l\}N\}F^\{N\_\{l\}\}\_\{l\}\(\\Delta\\varphi\)\}\{\\max\\\{\\sum\_\{l\\in\\mathbb\{P\}\}\\frac\{N\_\{l\}\}\{p\_\{l\}N\}F^\{N\_\{l\}\}\_\{l\}\(\\Delta 1\),Z\_\{\\min\}\\\}\}\(B\.5\)=∑l∈ℙZlNl​\(1pl​N​∑i=1Nlφl​\(θli\)​fl​\(θli\)Fl​\(θli,θl−1i\)−φl−1​\(θl−1i\)​fl−1​\(θl−1i\)Fl​\(θli,θl−1i\)\)max⁡\{∑l∈ℙZlNl​\(1pl​N​∑i=1Nlfl​\(θli\)Fl​\(θli,θl−1i\)−fl−1​\(θl−1i\)Fl​\(θli,θl−1i\)\),Zmin\},\\displaystyle=\\frac\{\\sum\_\{l\\in\\mathbb\{P\}\}Z\_\{l\}^\{N\_\{l\}\}\\Big\(\\frac\{1\}\{p\_\{l\}N\}\\sum\_\{i=1\}^\{N\_\{l\}\}\\varphi\_\{l\}\(\\theta\_\{l\}^\{i\}\)\\frac\{f\_\{l\}\(\\theta\_\{l\}^\{i\}\)\}\{F\_\{l\}\(\\theta\_\{l\}^\{i\},\\theta\_\{l\-1\}^\{i\}\)\}\-\\varphi\_\{l\-1\}\(\\theta\_\{l\-1\}^\{i\}\)\\frac\{f\_\{l\-1\}\(\\theta\_\{l\-1\}^\{i\}\)\}\{F\_\{l\}\(\\theta\_\{l\}^\{i\},\\theta\_\{l\-1\}^\{i\}\)\}\\Big\)\}\{\\max\\\{\\sum\_\{l\\in\\mathbb\{P\}\}Z\_\{l\}^\{N\_\{l\}\}\\Big\(\\frac\{1\}\{p\_\{l\}N\}\\sum\_\{i=1\}^\{N\_\{l\}\}\\frac\{f\_\{l\}\(\\theta\_\{l\}^\{i\}\)\}\{F\_\{l\}\(\\theta\_\{l\}^\{i\},\\theta\_\{l\-1\}^\{i\}\)\}\-\\frac\{f\_\{l\-1\}\(\\theta\_\{l\-1\}^\{i\}\)\}\{F\_\{l\}\(\\theta\_\{l\}^\{i\},\\theta\_\{l\-1\}^\{i\}\)\}\\Big\),Z\_\{\\min\}\\\}\},whereZminZ\_\{\\min\}is added here to prevent an unexpected computational small number\.

The practical settings for the rMLSMC estimator are shown in Appendix[G\.2](https://arxiv.org/html/2610.19074#A7.SS2)\.

#### B\.2\.4Serial multilevel sequential Monte Carlo estimator

1:Simulate

θ0i∼π¯0\\theta\_\{0\}^\{i\}\\sim\\bar\{\\pi\}\_\{0\}for

i=1,…,N0i=1,\.\.\.,N\_\{0\}
2:for

l=Lmin,…,Lmaxl=L\_\{\\text\{min\}\},\.\.\.,L\_\{\\text\{max\}\}do

3:Give

θl−1i\\theta\_\{l\-1\}^\{i\}for

i=1,…,Nl−1i=1,\.\.\.,N\_\{l\-1\}from previous level

l−1l\-1
4:Simulate only the new/missing dofs for the next level \(

ll\) from the prior sample, with

θli\|θl−1i∼π¯l\\theta\_\{l\}^\{i\}\|\\theta\_\{l\-1\}^\{i\}\\sim\\bar\{\\pi\}\_\{l\}for

i=1,…,Nli=1,\.\.\.,N\_\{l\}
5:Compute and store the likelihood

ℒl​\(θli\)\\mathcal\{L\}\_\{l\}\(\\theta\_\{l\}^\{i\}\)
6:for

i=1,…,Nli=1,\.\.\.,N\_\{l\}do

7:Define

wl,ji=ℒl​\(θli\)ℒl−1​\(θl−1i\)w\_\{l,j\}^\{i\}=\\frac\{\\mathcal\{L\}\_\{l\}\(\\theta\_\{l\}^\{i\}\)\}\{\\mathcal\{L\}\_\{l\-1\}\(\\theta^\{i\}\_\{l\-1\}\)\}
8:Resample: Select

Ili∼\{wl1,…,wlNl\}I\_\{l\}^\{i\}\\sim\\\{w\_\{l\}^\{1\},\.\.\.,w\_\{l\}^\{N\_\{l\}\}\\\}, and let

θ^li=θl−1Ili\\hat\{\\theta\}\_\{l\}^\{i\}=\\theta\_\{l\-1\}^\{I\_\{l\}^\{i\}\}
9:Mutate: Draw

θli∼ℳl​\(θ^li,⋅\)\\theta\_\{l\}^\{i\}\\sim\\mathcal\{M\}\_\{l\}\(\\hat\{\\theta\}\_\{l\}^\{i\},\\cdot\)
10:Compute and store the likelihood

ℒl​\(θli\)\\mathcal\{L\}\_\{l\}\(\\theta\_\{l\}^\{i\}\)
11:endfor

12:endfor

Algorithm 4Serial Sequential Monte Carlo samplerThe algorithm of the serial multilevel sequential Monte Carlo sampler is shown in Algorithm[12](https://arxiv.org/html/2610.19074#alg4.l12), which is taken from[Chada et al\. \[2025\]](https://arxiv.org/html/2610.19074#bib.bib16),[Beskos et al\. \[2017\]](https://arxiv.org/html/2610.19074#bib.bib8),[Beskos et al\. \[2018\]](https://arxiv.org/html/2610.19074#bib.bib12)\. This is a serial method in which the sampling in the current level is based on the previous level\. Here, give the level parameterL=\{Lmin,Lmax\}L=\\\{L\_\{\\text\{min\}\},L\_\{\\text\{max\}\}\\\}, we follow the same path as in MLSMC to generate the intermediate distributions and mutation kernel here\.

The serial MLSMC estimator is defined as

φ^sML\\displaystyle\\hat\{\\varphi\}^\{\\text\{sML\}\}=∑l=LminLmax\(∑i=1Nlφl​\(θli\)​ℒl​\(θli\)∑i=1Nlℒl​\(θli\)−1Nl​∑i=1Nlφl−1​\(θl−1i\)\),\\displaystyle=\\sum\_\{l=L\_\{\\text\{min\}\}\}^\{L\_\{\\text\{max\}\}\}\\Bigg\(\\frac\{\\sum\_\{i=1\}^\{N\_\{l\}\}\\varphi\_\{l\}\(\\theta\_\{l\}^\{i\}\)\\mathcal\{L\}\_\{l\}\(\\theta\_\{l\}^\{i\}\)\}\{\\sum\_\{i=1\}^\{N\_\{l\}\}\\mathcal\{L\}\_\{l\}\(\\theta\_\{l\}^\{i\}\)\}\-\\frac\{1\}\{N\_\{l\}\}\\sum\_\{i=1\}^\{N\_\{l\}\}\\varphi\_\{l\-1\}\(\\theta\_\{l\-1\}^\{i\}\)\\Bigg\),\(B\.6\)whereφLmin−1:=0\\varphi\_\{L\_\{\\text\{min\}\}\-1\}:=0\. The practical settings follow the same procedure as the MLSMC in Appendix[G\.1](https://arxiv.org/html/2610.19074#A7.SS1)\. However, the serial MLSMC is not parallelizable in level\-level since its procedure is serial\.

## Appendix CPMMH

1:Given

μ,θ\\mu,\\thetaand

pN​\(μ\)p^\{N\}\(\\mu\)
2:Propose

μ∗∼Q⁡\(⋅,μ\)\\mu^\{\*\}\\sim Q\(\\cdot,\\mu\)and sample

θ∗∼P\(⋅\|μ∗\)\\theta^\{\*\}\\sim P\(\\cdot\|\\mu^\{\*\}\)
3:Simulate

u∼U⁡\[\(0,1\)\]u\\sim U\[\(0,1\)\]and compute

α=min⁡\{π0​\(μ∗\)​Q​\(μ\|μ∗\)​pN​\(μ∗\)π0​\(μ\)​Q​\(μ∗\|μ\)​pN​\(μ\),1\}\.\\alpha=\\min\\Bigg\\\{\\frac\{\\pi\_\{0\}\(\\mu^\{\*\}\)Q\(\\mu\|\\mu^\{\*\}\)p^\{N\}\(\\mu^\{\*\}\)\}\{\\pi\_\{0\}\(\\mu\)Q\(\\mu^\{\*\}\|\\mu\)p^\{N\}\(\\mu\)\},\\ 1\\Bigg\\\}\.
4:If

u≤αu\\leq\\alpha, then

\(μ,θ,pN​\(μ\)\)=\(μ∗,θ∗,pN​\(μ∗\)\)\(\\mu,\\theta,p^\{N\}\(\\mu\)\)=\(\\mu^\{\*\},\\theta^\{\*\},p^\{N\}\(\\mu^\{\*\}\)\); otherwise,

\(μ,θ,pN​\(μ\)\)=\(μ,θ,pN​\(μ\)\)\(\\mu,\\theta,p^\{N\}\(\\mu\)\)=\(\\mu,\\theta,p^\{N\}\(\\mu\)\)\.

Algorithm 5PMMH at a given state\[[Andrieu et al\., 2010](https://arxiv.org/html/2610.19074#bib.bib29),[Andrieu and Roberts, 2009](https://arxiv.org/html/2610.19074#bib.bib30)\]We first describe the pseudo\-marginal Metropolis\-Hastings \(PMMH\) method for inferring the hyperparameterμ\\mu\. The target posterior distribution of this case is, given data𝒟\\mathcal\{D\},

π⁡\(μ\|𝒟\)∝π⁡\(𝒟\|μ\)​π0​\(μ\),\\pi\(\\mu\|\\mathcal\{D\}\)\\propto\\pi\(\\mathcal\{D\}\|\\mu\)\\pi\_\{0\}\(\\mu\),and one often found the posteriorπ⁡\(μ\|𝒟\)\\pi\(\\mu\|\\mathcal\{D\}\)might be analytically intractable or too complex to evaluate\. The pseudo\-marginal method then gives us an approach to computing it on an joint state spaceΘ×𝖲\\Theta\\times\\mathsf\{S\}, in which case we could computeπ⁡\(μ\|𝒟\)\\pi\(\\mu\|\\mathcal\{D\}\)by computing the marginal of the joint posterior distribution

π⁡\(θ,μ\|𝒟\)∝π⁡\(𝒟,d​θ\|μ\)​π0​\(μ\)∝f⁡\(θ\|μ\)​π0​\(μ\)\.\\pi\(\\theta,\\mu\|\\mathcal\{D\}\)\\propto\\pi\(\\mathcal\{D\},d\\theta\|\\mu\)\\pi\_\{0\}\(\\mu\)\\propto f\(\\theta\|\\mu\)\\pi\_\{0\}\(\\mu\)\.to get the target distribution as follows

π⁡\(μ\|𝒟\)=∫Θπ⁡\(𝑑θ,μ\|𝒟\)∝∫Θπ⁡\(𝒟,𝑑θ\|μ\)​π0​\(μ\)\.\\pi\(\\mu\|\\mathcal\{D\}\)=\\int\_\{\\Theta\}\\pi\(d\\theta,\\mu\|\\mathcal\{D\}\)\\propto\\int\_\{\\Theta\}\\pi\(\\mathcal\{D\},d\\theta\|\\mu\)\\pi\_\{0\}\(\\mu\)\.Then, the likelihoodp⁡\(μ\):=π⁡\(𝒟\|μ\)p\(\\mu\):=\\pi\(\\mathcal\{D\}\|\\mu\)can be computed by

p⁡\(μ\)=∫Θπ⁡\(𝒟,𝑑θ\|μ\)=∫Θf⁡\(𝑑θ\|μ\)\.p\(\\mu\)=\\int\_\{\\Theta\}\\pi\(\\mathcal\{D\},d\\theta\|\\mu\)=\\int\_\{\\Theta\}f\(d\\theta\|\\mu\)\.We often call this likelihoodp⁡\(μ\)p\(\\mu\)the marginal distribution since it computes the marginal of the joint likelihood\.

One can then run a SMC sampler with the selection functionshjh\_\{j\}in each tempering stepj=1,…​J−1j=1,\.\.\.J\-1to give an unbiased estimator of the marginal distribution

pN​\(μ\)=∏j=1J−11N​∑i=1Nhj​\(θji\|μ\),p^\{N\}\(\\mu\)=\\prod\_\{j=1\}^\{J\-1\}\\frac\{1\}\{N\}\\sum\_\{i=1\}^\{N\}h\_\{j\}\(\\theta\_\{j\}^\{i\}\|\\mu\),i\.e\.𝔼⁡\[pN​\(μ\)\]=p⁡\(μ\)\\mathbb\{E\}\[p^\{N\}\(\\mu\)\]=p\(\\mu\), by sampling from the following auxiliary joint, on the space\{1,…,N\}J×ΘN​J\\\{1,\.\.\.,N\\\}^\{J\}\\times\\Theta^\{NJ\},

P\(I1:J1:N,θ1:J1:N\|μ\)\\displaystyle P\(I\_\{1:J\}^\{1:N\},\\theta\_\{1:J\}^\{1:N\}\|\\mu\)=∏i=1Nπ¯​\(θ1I1i\|μ\)​∏j=1J−1ℳ⁡\(θj\+1Ij\+1i,θjIji\|μ\)​hj​\(θjIji\|μ\)∑k=1Nhj​\(θjk\|μ\)\.\\displaystyle=\\prod\_\{i=1\}^\{N\}\\bar\{\\pi\}\(\\theta\_\{1\}^\{I\_\{1\}^\{i\}\}\|\\mu\)\\prod\_\{j=1\}^\{J\-1\}\\mathcal\{M\}\(\\theta\_\{j\+1\}^\{I\_\{j\+1\}^\{i\}\},\\theta\_\{j\}^\{I\_\{j\}^\{i\}\}\|\\mu\)\\frac\{h\_\{j\}\(\\theta\_\{j\}^\{I\_\{j\}^\{i\}\}\|\\mu\)\}\{\\sum\_\{k=1\}^\{N\}h\_\{j\}\(\\theta\_\{j\}^\{k\}\|\\mu\)\}\.Moreover,pN​\(μ\)​π0​\(μ\)p^\{N\}\(\\mu\)\\pi\_\{0\}\(\\mu\)is an unbiased estimator of the joint distributionπ⁡\(𝒟,μ\)\\pi\(\\mathcal\{D\},\\mu\)which is proportional to our target distributionπ⁡\(μ\|𝒟\)\\pi\(\\mu\|\\mathcal\{D\}\)\. A general approach of the PMMH is given in Algorithm[4](https://arxiv.org/html/2610.19074#alg5.l4)\.

In the previous discussions, we fix the hyperparameter’s valueμ\\muand omit it in the computation\. We now expand the hyperparameterμ\\muterm in these expressions by writing joint variableφ⁡\(x,θ,μ\)\\varphi\(x,\\theta,\\mu\)\(or conditional variableφ⁡\(x,θ\|μ\)\\varphi\(x,\\theta\|\\mu\)\) and joint variableπ⁡\(θ,μ\)\\pi\(\\theta,\\mu\)\(or conditional variableπ⁡\(θ\|μ\)\\pi\(\\theta\|\\mu\)\)\. Also, sometimes we omit the termxxfor symbol simplicity since the data are the same\.

## Appendix DDetails of the S\-MLPMMH approximation

We provide further details of the empirical\-Bayes, Laplace, and MAP\-centered SBMC construction used in S\-MLPMMH\.

##### Hyperparameter transformation and empirical\-Bayes estimate\.

Consider the transformed hyperparameter variableη=log⁡μ\\eta=\\log\\mu, whereμ\>0\\mu\>0denotes the standard deviation hyperparameter\. We place a Gaussian hyper\-prior on the transformed variableη\\eta,

η∼𝒩⁡\(mη,ση2\),\\eta\\sim\\mathcal\{N\}\(m\_\{\\eta\},\\sigma\_\{\\eta\}^\{2\}\),which is equivalent to assigning a log\-normal prior toμ\\mu\. Given the mean and variance ofμ\\mu, the corresponding Gaussian parameters forη=log⁡μ\\eta=\\log\\muare

ση2=log⁡\(1\+Var⁡\(μ\)𝔼​\[μ\]2\),mη=log⁡\(𝔼⁡\[μ\]\)−12​ση2\.\\sigma\_\{\\eta\}^\{2\}=\\log\\left\(1\+\\frac\{\\operatorname\{Var\}\(\\mu\)\}\{\\mathbb\{E\}\[\\mu\]^\{2\}\}\\right\),\\qquad m\_\{\\eta\}=\\log\\left\(\\mathbb\{E\}\[\\mu\]\\right\)\-\\frac\{1\}\{2\}\\sigma\_\{\\eta\}^\{2\}\.
At each levelll, we first compute an empirical\-Bayes MAP estimate of the hyperparameter\. For each candidate value ofη\\eta, withμ=exp⁡\(η\)\\mu=\\exp\(\\eta\), we compute a MAP estimate of the model parametersθ\\thetaconditional onμ\\mu\. Around this conditional MAP estimate, we approximate the marginal likelihood using a diagonal Laplace approximation\. The level\-dependent empirical\-Bayes objective is then

Jl​\(η\)=log⁡π^l​\(𝒟∣μ=exp⁡\(η\)\)\+log⁡π0​\(μ\),J\_\{l\}\(\\eta\)=\\log\\hat\{\\pi\}\_\{l\}\\left\(\\mathcal\{D\}\\mid\\mu=\\exp\(\\eta\)\\right\)\+\\log\\pi\_\{0\}\(\\mu\),whereπ^l​\(𝒟∣μ=exp⁡\(η\)\)\\hat\{\\pi\}\_\{l\}\(\\mathcal\{D\}\\mid\\mu=\\exp\(\\eta\)\)is the diagonal Laplace approximation of the marginal likelihood at levelll, andπ0​\(μ\)\\pi\_\{0\}\(\\mu\)is the hyper\-prior\.

The resulting hyperparameter MAP center is

ηl,𝖬𝖠𝖯=arg⁡maxη​Jl​\(η\),μl,𝖬𝖠𝖯=exp⁡\(ηl,𝖬𝖠𝖯\)\.\\eta\_\{l,\\sf MAP\}=\\arg\\max\_\{\\eta\}J\_\{l\}\(\\eta\),\\qquad\\mu\_\{l,\\sf MAP\}=\\exp\(\\eta\_\{l,\\sf MAP\}\)\.\(D\.1\)

##### Outer Laplace approximation\.

After obtainingηl,𝖬𝖠𝖯\\eta\_\{l,\\sf MAP\}, we construct a local Laplace approximation to the hyperparameter posterior centered atηl,𝖬𝖠𝖯\\eta\_\{l,\\sf MAP\}\. Specifically, we estimate the local curvature of the empirical\-Bayes objective by

Hl,η=−∂2Jl​\(η\)∂η2\|η=ηl,𝖬𝖠𝖯,H\_\{l,\\eta\}=\-\\left\.\\frac\{\\partial^\{2\}J\_\{l\}\(\\eta\)\}\{\\partial\\eta^\{2\}\}\\right\|\_\{\\eta=\\eta\_\{l,\\sf MAP\}\},and use the corresponding Laplace varianceVl,η=Hl,η−1V\_\{l,\\eta\}=H\_\{l,\\eta\}^\{\-1\}\. The outer SBMC hyperparameter samples are then independently drawn from

ηk,l∼𝒩\(ηl,𝖬𝖠𝖯,Vl,η\),μk,l=exp\(ηk,l\),k=1,…,Ml\.\\eta\_\{k,l\}\\sim\\mathcal\{N\}\\left\(\\eta\_\{l,\\sf MAP\},V\_\{l,\\eta\}\\right\),\\qquad\\mu\_\{k,l\}=\\exp\(\\eta\_\{k,l\}\),\\qquad k=1,\\ldots,M\_\{l\}\.In practice,Hl,ηH\_\{l,\\eta\}is estimated using a finite\-difference approximation ofJlJ\_\{l\}aroundηl,𝖬𝖠𝖯\\eta\_\{l,\\sf MAP\}\.

##### Inner MAP\-centered SBMC approximation\.

For each sampled hyperparameterμk,l\\mu\_\{k,l\}, we then construct an inner MAP\-centered SBMC sampler\. The original inner prior is

π¯l=𝒩⁡\(0,μk,l2​𝖨\)\.\\bar\{\\pi\}\_\{l\}=\\mathcal\{N\}\\left\(0,\\mu\_\{k,l\}^\{2\}\\sf I\\right\)\.Conditional onμk,l\\mu\_\{k,l\}, we compute the MAP estimator

θl,𝖬𝖠𝖯=arg⁡maxθ​𝐋l​\(θl,θ∗l−1\)​π¯l​\(θl,θl−1\)\.\\theta\_\{l,\\sf MAP\}=\\arg\\max\_\{\\theta\}\{\\bf L\}\_\{l\}\\left\(\\theta\_\{l\},\\theta\*\{l\-1\}\\right\)\\bar\{\\pi\}\_\{l\}\\left\(\\theta\_\{l\},\\theta\_\{l\-1\}\\right\)\.The MAP estimate is recomputed for each hyperparameter sampleμk,l\\mu\_\{k,l\}\.

Following the SBMC framework, we then establish a new inner prior centered at the MAP estimate in \(\) and shrink its covariance by a tuning parametersθs\_\{\\theta\}\. The resulting prior is

π¯^l=𝒩⁡\(θl,𝖬𝖠𝖯,sθ​μk,l2​𝖨\)\.\\hat\{\\bar\{\\pi\}\}\_\{l\}=\\mathcal\{N\}\\left\(\\theta\_\{l,\\sf MAP\},s\_\{\\theta\}\\mu\_\{k,l\}^\{2\}\\sf I\\right\)\.In our numerical experiments, we usesθ=0\.1s\_\{\\theta\}=0\.1\. This MAP\-centered construction reduces the cost of the inner sampling stage while retaining the multilevel correction structure across successive levels\.

## Appendix ETheoretical results

In this section, we present our main results related to our proposed methodology, discussed in Section[3](https://arxiv.org/html/2610.19074#S3)\. Specifically, we first restate the theoretical results for MLSMC methods, provide a variance bound on the trace\-class neural network and a convergence theorem related to our joint posterior distribution of interest is followed by\.

We write𝖵=Θ×𝖲\\mathsf\{V\}=\\Theta\\times\\mathsf\{S\}with associatedσ−\\sigma\-algebra as𝒱\\mathcal\{V\}\. Let\(𝖶,𝒲\)\(\\mathsf\{W\},\\mathcal\{W\}\)be a measurable space such that𝖵∈𝖶\\mathsf\{V\}\\in\\mathsf\{W\}\. LetK:𝖶×𝒲→\[0,1\]K:\\mathsf\{W\}\\times\\mathcal\{W\}\\rightarrow\[0,1\]be any Markov kernel of invariant measureπ\\pisuch that, for anyB⊂𝒱B\\subset\\mathcal\{V\},

∫B×\(𝖶\\𝖵\)π⁡\(𝑑w\)=Π⁡\(B\)\.\\int\_\{B\\times\(\\mathsf\{W\}\\backslash\\mathsf\{V\}\)\}\\pi\(dw\)=\\Pi\(B\)\.
### E\.1Assumptions

##### MLMC bound assumptions\.

The assumption for boundedness of the functionfi,j​\(θ\|μ\)f\_\{i,j\}\(\\theta\|\\mu\)is given in Assumption[E\.1](https://arxiv.org/html/2610.19074#A5.Thmass1)and for rate decay in parameters required for the application of the MLMC is given in Assumption[E\.2](https://arxiv.org/html/2610.19074#A5.Thmass2)\.

###### Assumption E\.1\.

LetJ∈ℕJ\\in\\mathbb\{N\}be given\. For eachj∈\{0,…,J\}j\\in\\\{0,\\dots,J\\\}there exists someC\>0C\>0such that for all\(l,θ,μ\)∈ℕ×Θ×𝖲\(l,\\theta,\\mu\)\\in\\mathbb\{N\}\\times\\Theta\\times\\mathsf\{S\}

C−1<Zl,fl,j​\(θ\|μ\)≤C\.C^\{\-1\}<Z\_\{l\},\\ f\_\{l,j\}\(\\theta\|\\mu\)\\leq C\\,\.

###### Assumption E\.2\.

For anyφ∈ℬb​\(Θl×Θl−1×S\)∩Lip​\(Θl×Θl−1×𝖲\)\\varphi\\in\\mathcal\{B\}\_\{b\}\(\\Theta\_\{l\}\\times\\Theta\_\{l\-1\}\\times S\)\\cap\\text\{Lip\}\(\\Theta\_\{l\}\\times\\Theta\_\{l\-1\}\\times\\mathsf\{S\}\)there existc1,c2,c3,β,α,γ\>0c\_\{1\},c\_\{2\},c\_\{3\},\\beta,\\alpha,\\gamma\>0such that the following holds, under the joint posteriorΠl\\Pi\_\{l\},

- \(B\)𝔼⁡\[\|φl​\(θl,μ\)​fl​\(θl\|μ\)Fl​\(θl,θl−1\|μ\)−f⁡\(φ⁡\(θ,μ\)\)\|\]=:Bl≤c1​2−α​l\\mathbb\{E\}\\Big\[\\big\|\\varphi\_\{l\}\(\\theta\_\{l\},\\mu\)\\frac\{f\_\{l\}\(\\theta\_\{l\}\|\\mu\)\}\{F\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\}\-f\(\\varphi\(\\theta,\\mu\)\)\\big\|\\Big\]=:B\_\{l\}\\leq c\_\{1\}2^\{\-\\alpha l\};
- \(V\)𝔼⁡\[\|φl​\(θl,μ\)​fl​\(θl\|μ\)Fl​\(θl,θl−1\|μ\)−φl−1​\(θl−1,μ\)​fl−1​\(θl−1\|μ\)Fl​\(θl,θl−1\|μ\)\|2\]=:Vl≤c2​2−β​l\\mathbb\{E\}\\Big\[\\big\|\\varphi\_\{l\}\(\\theta\_\{l\},\\mu\)\\frac\{f\_\{l\}\(\\theta\_\{l\}\|\\mu\)\}\{F\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\}\-\\varphi\_\{l\-1\}\(\\theta\_\{l\-1\},\\mu\)\\frac\{f\_\{l\-1\}\(\\theta\_\{l\-1\}\|\\mu\)\}\{F\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\}\\big\|^\{2\}\\Big\]=:V\_\{l\}\\leq c\_\{2\}2^\{\-\\beta l\};
- \(C\)COST⁡\(Πl​\(φl\),Πl−1​\(φl−1\)\)=:Cl≤c3​2γ​l\{\\rm COST\}\(\\Pi\_\{l\}\(\\varphi\_\{l\}\),\\Pi\_\{l\-1\}\(\\varphi\_\{l\-1\}\)\)=:C\_\{l\}\\leq c\_\{3\}2^\{\\gamma l\}\.

Assumption[E\.2](https://arxiv.org/html/2610.19074#A5.Thmass2)is related to the standard assumptions required in MLMC, where one has assumptions on the “weak” error where \(B\) denotes the bias, the “strong” error where \(V\) denotes the variance and the costs hence \(C\), where the cost is related to the cost of functional evaluations, and not time taken to complete the experiment\. Fix the value of the hyperparameterμ\\mu, Assumption[E\.1](https://arxiv.org/html/2610.19074#A5.Thmass1)and[E\.2](https://arxiv.org/html/2610.19074#A5.Thmass2)falls into the original assumptions needed in MLMC\.

##### PMMH assumption\.

The assumption for transition in the PMMH is given in Assumption[E\.3](https://arxiv.org/html/2610.19074#A5.Thmass3)\.

###### Assumption E\.3\.

There exists aξ∈\(0,1\)\\xi\\in\(0,1\)andν∈𝒫⁡\(𝖶\)\\nu\\in\\mathcal\{P\}\(\\mathsf\{W\}\)such that, for eachw∈𝖶w\\in\\mathsf\{W\},φ∈ℬb​\(𝖶\)∩Lip​\(𝖶\)\\varphi\\in\\mathcal\{B\}\_\{b\}\(\\mathsf\{W\}\)\\cap\\text\{Lip\}\(\\mathsf\{W\}\),

∫𝖶φ⁡\(w′\)​K​\(w,d​w′\)≥ξ​∫Wφ⁡\(v\)​ν​\(𝑑v\)\.\\int\_\{\\mathsf\{W\}\}\\varphi\(w^\{\\prime\}\)K\(w,dw^\{\\prime\}\)\\geq\\xi\\int\_\{W\}\\varphi\(v\)\\nu\(dv\)\.

The verification of Assumption[E\.2](https://arxiv.org/html/2610.19074#A5.Thmass2)and[E\.3](https://arxiv.org/html/2610.19074#A5.Thmass3)above are given in Appendix[H\.5](https://arxiv.org/html/2610.19074#A8.SS5)\.

### E\.2MLSMC results

We restated theoretical results for MLSMC with self\-normalized estimator444It is worth noting that only the ratio estimator has a rigorous theoretical result in the MIMC case \(an extension from the MLMC\), where the self\-normalized does not\., MLSMC with ratio estimator, randomized MLSMC estimator and serial MLSMC estimator, which are shown in Theorem[E\.1](https://arxiv.org/html/2610.19074#A5.Thmtheorem1),[E\.2](https://arxiv.org/html/2610.19074#A5.Thmtheorem2),[E\.3](https://arxiv.org/html/2610.19074#A5.Thmtheorem3)and[E\.4](https://arxiv.org/html/2610.19074#A5.Thmtheorem4), respectively\.

###### Theorem E\.1\.

Assume Assumptions[E\.1](https://arxiv.org/html/2610.19074#A5.Thmass1)and[E\.2](https://arxiv.org/html/2610.19074#A5.Thmass2)with fixedμ\\muandβ\>γ\\beta\>\\gamma\. Then for anyε\>0\\varepsilon\>0and suitableφ:Θ→ℝ\\varphi:\\Theta\\rightarrow\\mathbb\{R\}, it is possible to choose a level paramemterL=\{Lmin,Lmax\}L=\\\{L\_\{\\text\{min\}\},L\_\{\\text\{max\}\}\\\}, whereLminL\_\{\\text\{min\}\}be the starting level andLmaxL\_\{\\text\{max\}\}be the maximum level, such that for someC\>0C\>0

𝔼⁡\[\(φ^LMLSN−π⁡\(φ\)\)2\]≤C​ε2,\\mathbb\{E\}\[\(\\widehat\{\\varphi\}^\{\\rm MLSN\}\_\{L\}\-\\pi\(\\varphi\)\)^\{2\}\]\\leq C\\varepsilon^\{2\}\\,,andCOST⁡\(φ^LMLSN\)≤C​ε−2\{\\rm COST\}\(\\widehat\{\\varphi\}^\{\\rm MLSN\}\_\{L\}\)\\leq C\\varepsilon^\{\-2\}, the canonical rate\. The estimatorφ^LMLSN\\widehat\{\\varphi\}^\{\\rm MLSN\}\_\{L\}is defined in equation \([B\.3](https://arxiv.org/html/2610.19074#A2.E3)\)\.

###### Proof\.

The proof is given in[Moral et al\. \[2017\]](https://arxiv.org/html/2610.19074#bib.bib13)\. ∎

###### Theorem E\.2\.

Assume Assumptions[E\.1](https://arxiv.org/html/2610.19074#A5.Thmass1)and[E\.2](https://arxiv.org/html/2610.19074#A5.Thmass2)with fixedμ\\muandβ\>γ\\beta\>\\gamma\. Then for anyε\>0\\varepsilon\>0and suitableφ:Θ→ℝ\\varphi:\\Theta\\rightarrow\\mathbb\{R\}, it is possible to choose a level paramemterL=\{Lmin,Lmax\}L=\\\{L\_\{\\text\{min\}\},L\_\{\\text\{max\}\}\\\}, whereLminL\_\{\\text\{min\}\}be the starting level andLmaxL\_\{\\text\{max\}\}be the maximum level, such that for someC\>0C\>0

𝔼⁡\[\(φ^LMLRE−π⁡\(φ\)\)2\]≤C​ε2,\\mathbb\{E\}\[\(\\widehat\{\\varphi\}^\{\\rm MLRE\}\_\{L\}\-\\pi\(\\varphi\)\)^\{2\}\]\\leq C\\varepsilon^\{2\}\\,,andCOST⁡\(φ^LMLRE\)≤C​ε−2\{\\rm COST\}\(\\widehat\{\\varphi\}^\{\\rm MLRE\}\_\{L\}\)\\leq C\\varepsilon^\{\-2\}, the canonical rate\. The estimatorφ^LMLRE\\widehat\{\\varphi\}^\{\\rm MLRE\}\_\{L\}is defined in equation \([B\.4](https://arxiv.org/html/2610.19074#A2.E4)\)\.

###### Proof\.

The proof is given in[Law et al\. \[2022\]](https://arxiv.org/html/2610.19074#bib.bib5)\. ∎

###### Theorem E\.3\.

Assume[E\.1](https://arxiv.org/html/2610.19074#A5.Thmass1)and[E\.2](https://arxiv.org/html/2610.19074#A5.Thmass2)with fixedμ\\muandβ\>γ\\beta\>\\gamma\. Then, for suitableφ:Θ→ℝ\\varphi:\\Theta\\rightarrow\\mathbb\{R\}, give a starting levelLminL\_\{\\text\{min\}\}, it is possible to choose a probability distribution𝐩\{\\bf p\}onℙ=\{l∈ℕ:l≥Lmin\}\\mathbb\{P\}=\\\{l\\in\\mathbb\{N\}:l\\geq L\_\{\\text\{min\}\}\\\}such that for someC\>0C\>0,

𝔼⁡\[\(φ^rML−π⁡\(φ\)\)2\]≤C​ε2,\\mathbb\{E\}\[\(\\widehat\{\\varphi\}^\{\\rm rML\}\-\\pi\(\\varphi\)\)^\{2\}\]\\leq C\\varepsilon^\{2\}\\,,and expectedCOST⁡\(φ^rMI\)≤C​ε−2\{\\rm COST\}\(\\widehat\{\\varphi\}^\{\\rm rMI\}\)\\leq C\\varepsilon^\{\-2\}, i\.e\. the canonical rate\. The estimatorφ^rMI\\widehat\{\\varphi\}^\{\\rm rMI\}is defined in equation \([B\.5](https://arxiv.org/html/2610.19074#A2.E5)\)\.

###### Proof\.

The proof is given in[Cotter et al\. \[2022\]](https://arxiv.org/html/2610.19074#bib.bib27)\. ∎

###### Theorem E\.4\.

Assume Assumptions[E\.1](https://arxiv.org/html/2610.19074#A5.Thmass1)and[E\.2](https://arxiv.org/html/2610.19074#A5.Thmass2)with fixedμ\\muandβ\>γ\\beta\>\\gamma\. Then for anyε\>0\\varepsilon\>0and suitableφ:Θ→ℝ\\varphi:\\Theta\\rightarrow\\mathbb\{R\}, it is possible to choose a level paramemterL=\{Lmin,Lmax\}L=\\\{L\_\{\\text\{min\}\},L\_\{\\text\{max\}\}\\\}, whereLminL\_\{\\text\{min\}\}be the starting level andLmaxL\_\{\\text\{max\}\}be the maximum level, such that for someC\>0C\>0

𝔼⁡\[\(φ^LsML−π⁡\(φ\)\)2\]≤C​ε2,\\mathbb\{E\}\[\(\\widehat\{\\varphi\}^\{\\rm sML\}\_\{L\}\-\\pi\(\\varphi\)\)^\{2\}\]\\leq C\\varepsilon^\{2\}\\,,andCOST⁡\(φ^LsML\)≤C​ε−2\{\\rm COST\}\(\\widehat\{\\varphi\}^\{\\rm sML\}\_\{L\}\)\\leq C\\varepsilon^\{\-2\}, the canonical rate\. The estimatorφ^LsML\\widehat\{\\varphi\}^\{\\rm sML\}\_\{L\}is defined in equation \([B\.6](https://arxiv.org/html/2610.19074#A2.E6)\)\.

###### Proof\.

The proof is given in[Beskos et al\. \[2017\]](https://arxiv.org/html/2610.19074#bib.bib8),[Beskos et al\. \[2018\]](https://arxiv.org/html/2610.19074#bib.bib12)\. ∎

### E\.3Trace\-class neural network priors

We begin this section by providing a result on the strong convergence, i\.e\. a bound on the variance related to the TNN prior\. In particular, as shown in[Chada et al\. \[2025\]](https://arxiv.org/html/2610.19074#bib.bib16), for a trace‐class neural network prior with smoothness parameters≥1/2s\\geq 1/2, the truncation error of the network output satisfies the decay rate in \([E\.1](https://arxiv.org/html/2610.19074#A5.E1)\)\. If the quantity of interest is obtained by applying a Lipschitz mapΨ\\Psito the network output, the same rate transfers directly toφl=Ψ∘gl\\varphi\_\{l\}=\\Psi\\circ g\_\{l\}, see the Proposition[E\.1](https://arxiv.org/html/2610.19074#A5.Thmprop1)\.

###### Proposition E\.1\.

Assume that for allz∈ℝz\\in\\mathbb\{R\},σ⁡\(z\)≤\|z\|\\sigma\(z\)\\leq\|z\|\. Then fors≥1/2s\\geq 1/2, the priors defined as in \([2\.1](https://arxiv.org/html/2610.19074#S2.E1)\) are trace\-class\. Furthermore, letx∈𝖷x\\in\\mathsf\{X\}andμ∈𝖲\\mu\\in\\mathsf\{S\}, and consider the TNNgl​\(x,θl,μ\)g\_\{l\}\(x,\\theta\_\{l\},\\mu\)truncated atnl=2ln\_\{l\}=2^\{l\}terms, forl∈ℕl\\in\\mathbb\{N\}, with the limit denotedg⁡\(x,θ,μ\)=liml→∞gl​\(x,θl,μ\)g\(x,\\theta,\\mu\)=\\lim\_\{l\\rightarrow\\infty\}g\_\{l\}\(x,\\theta\_\{l\},\\mu\), as described above\. Then there areC′​\(x\),C⁡\(x\)<\+∞C^\{\\prime\}\(x\),\\ C\(x\)<\+\\inftysuch that

∫Θ×𝖲\|gl​\(x,θl,μ\)−g⁡\(x,θ,μ\)\|2​π¯​\(𝑑θ\)​π0​\(𝑑μ\)≤C′​2−\(2​s−1\)​l\.\\int\_\{\\Theta\\times\\mathsf\{S\}\}\|g\_\{l\}\(x,\\theta\_\{l\},\\mu\)\-g\(x,\\theta,\\mu\)\|^\{2\}\\bar\{\\pi\}\(d\\theta\)\\pi\_\{0\}\(d\\mu\)\\leq C^\{\\prime\}2^\{\-\(2s\-1\)l\}\.\(E\.1\)ifΨ\\Psiis Lipschitz continuous, then naturally we have

∫Θ×𝖲\|φl​\(x,θl,μ\)−φ⁡\(x,θ,μ\)\|2​π¯​\(𝑑θ\)​π0​\(𝑑μ\)≤C​2−\(2​s−1\)​l\.\\int\_\{\\Theta\\times\\mathsf\{S\}\}\|\\varphi\_\{l\}\(x,\\theta\_\{l\},\\mu\)\-\\varphi\(x,\\theta,\\mu\)\|^\{2\}\\bar\{\\pi\}\(d\\theta\)\\pi\_\{0\}\(d\\mu\)\\leq C2^\{\-\(2s\-1\)l\}\.

###### Proof\.

The proof of the first part of Proposition[E\.1](https://arxiv.org/html/2610.19074#A5.Thmprop1)can be found in[Chada et al\. \[2025\]](https://arxiv.org/html/2610.19074#bib.bib16)\. The second part can be derived instantly since

\|φl−φ\|=\|Ψ⁡\(gl\)−Ψ⁡\(g\)\|≤L​\|gl−g\|,L\>0,\|\\varphi\_\{l\}\-\\varphi\|=\|\\Psi\(g\_\{l\}\)\-\\Psi\(g\)\|\\leq L\|g\_\{l\}\-g\|,\\quad L\>0,whereLLis the Lipschtiz constant\. ∎

### E\.4Supporting results

###### Theorem E\.5\.

Assume Assumptions[E\.1](https://arxiv.org/html/2610.19074#A5.Thmass1),[E\.2](https://arxiv.org/html/2610.19074#A5.Thmass2)\. For anyφ:Θ×𝖲→ℝd\\varphi:\\Theta\\times\\mathsf\{S\}\\rightarrow\\mathbb\{R\}^\{d\}, there exists aC<\+∞C<\+\\inftysuch that

𝔼⁡\[\(ΠlN,Ml​\(Δ​φl\)−Πl​\(Δ​φl\)\)2\]≤C​2−β​lMl\.\\mathbb\{E\}\\Bigg\[\\Big\(\\Pi\_\{l\}^\{N,M\_\{l\}\}\(\\Delta\\varphi\_\{l\}\)\-\\Pi\_\{l\}\(\\Delta\\varphi\_\{l\}\)\\Big\)^\{2\}\\Bigg\]\\leq\\frac\{C2^\{\-\\beta l\}\}\{M\_\{l\}\}\.\(E\.2\)

###### Proof\.

Proof sees Appendix[H\.2](https://arxiv.org/html/2610.19074#A8.SS2)\. ∎

###### Lemma E\.1\.

Leta,b,c,d,A,B,C,D∈ℝa,b,c,d,A,B,C,D\\in\\mathbb\{R\}withA,B,C,DA,B,C,Dnon\-zero then;

aA−\\displaystyle\\frac\{a\}\{A\}\-bB−\(cC−dD\)=\[a−b−\(c−d\)\]A−b⁡\[A−B−\(C−D\)\]A​B\+1A​C​\[C−A\]​\[c−d\]\\displaystyle\\frac\{b\}\{B\}\-\\Bigg\(\\frac\{c\}\{C\}\-\\frac\{d\}\{D\}\\Bigg\)=\\frac\{\[a\-b\-\(c\-d\)\]\}\{A\}\-\\frac\{b\[A\-B\-\(C\-D\)\]\}\{AB\}\+\\frac\{1\}\{AC\}\[C\-A\]\[c\-d\]\(E\.3\)−1A​B​\(b−d\)​\(C−D\)\+dC​B​D​\(B−D\)​\(C−D\)\+dA​C​B​\(A−C\)​\(C−D\)\.\\displaystyle\-\\frac\{1\}\{AB\}\(b\-d\)\(C\-D\)\+\\frac\{d\}\{CBD\}\(B\-D\)\(C\-D\)\+\\frac\{d\}\{ACB\}\(A\-C\)\(C\-D\)\.

###### Proof\.

The Lemma is from Lemma C\.3\. in[Jasra et al\. \[2017\]](https://arxiv.org/html/2610.19074#bib.bib31)and its proof is given in Appendix[H\.3](https://arxiv.org/html/2610.19074#A8.SS3)for completeness\. ∎

Definevk=\(θk,l,μk,l\)v^\{k\}=\(\\theta\_\{k,l\},\\mu\_\{k,l\}\)as the relevant variables ofwkw^\{k\}from the MCMC kernel, and defining

φ~l​\(vk\):=\(φl​\(θk,l,μk,l\)​fl​\(θk,l\|μk,l\)Fl​\(θk,l,θk,l−1\|μk,l\)−φl−1​\(θk,l−1,μk,l\)​fl−1​\(θk,l−1\|μk,l\)Fl​\(θk,l,θk,l−1\|μk,l\)\),\\tilde\{\\varphi\}\_\{l\}\(v^\{k\}\):=\\Big\(\\varphi\_\{l\}\(\\theta\_\{k,l\},\\mu\_\{k,l\}\)\\frac\{f\_\{l\}\(\\theta\_\{k,l\}\|\\mu\_\{k,l\}\)\}\{F\_\{l\}\(\\theta\_\{k,l\},\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)\}\-\\varphi\_\{l\-1\}\(\\theta\_\{k,l\-1\},\\mu\_\{k,l\}\)\\frac\{f\_\{l\-1\}\(\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)\}\{F\_\{l\}\(\\theta\_\{k,l\},\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)\}\\Big\),and the mapψ:wk→vk\\psi:w^\{k\}\\rightarrow v^\{k\}\.

###### Proposition E\.2\.

Assume Assumptions[E\.1](https://arxiv.org/html/2610.19074#A5.Thmass1),[E\.3](https://arxiv.org/html/2610.19074#A5.Thmass3),[E\.2](https://arxiv.org/html/2610.19074#A5.Thmass2)\. Suppose that\{Wi\}i\\\{W^\{i\}\\\}\_\{i\}is a Markov chain with the Markov kernelKK, andW1∼πW^\{1\}\\sim\\pi\. Then there exists aβ\>0\\beta\>0such that for anyφ∈ℬb​\(Θ×𝖲\)∩Lip​\(Θ×𝖲\)\\varphi\\in\\mathcal\{B\}\_\{b\}\(\\Theta\\times\\mathsf\{S\}\)\\cap\\text\{Lip\}\(\\Theta\\times\\mathsf\{S\}\)there exists aC\>0C\>0such that

𝔼⁡\[\(1Ml​∑k=1Mlφ~l​\(Vk\)−Πl​\(φ~l\)\)2\]≤C​2−β​lMl,\\mathbb\{E\}\\Bigg\[\\Big\(\\frac\{1\}\{M\_\{l\}\}\\sum\_\{k=1\}^\{M\_\{l\}\}\\tilde\{\\varphi\}\_\{l\}\(V^\{k\}\)\-\\Pi\_\{l\}\(\\tilde\{\\varphi\}\_\{l\}\)\\Big\)^\{2\}\\Bigg\]\\leq\\frac\{C2^\{\-\\beta l\}\}\{M\_\{l\}\},whereVk=\(θk,l,μ\)V^\{k\}=\(\\theta\_\{k,l\},\\mu\)is the relevant variables ofWkW^\{k\}\.

###### Proof\.

Proof sees the Appendix[H\.4](https://arxiv.org/html/2610.19074#A8.SS4)\. ∎

## Appendix FAdditional Numerical Experiments

We begin by presenting a comparative study of several MLSMC methods for approximating the target conditional posterior distribution\. Next, we conduct a numerical verification on a synthetic regression benchmark to confirm both the regularity and complexity rate of the MLPMMH method\. Furthermore, we justify the TNN architecture on the MNIST subset through controlled MAP training experiments\. The COST is measured as the time \(arithmetic\) complexity, i\.e\. the number of basic floating\-point operations needed in the coupled level\.

### F\.1Problem setting

In the case that data is defined in continuous space, such as the regression problem, we haveY=ℝmY=\\mathbb\{R\}^\{m\}, then we can assume the data follows

yi=g⁡\(xi,θ,μ\)\+εi,ε∼i\.i\.d\.𝒩m​\(0,Σi\),y\_\{i\}=g\(x\_\{i\},\\theta,\\mu\)\+\\varepsilon\_\{i\},\\ \\varepsilon\\mathrel\{\\overset\{\\mathclap\{\\text\{\\tiny i\.i\.d\.\}\}\}\{\\sim\}\}\\mathcal\{N\}\_\{m\}\(0,\\Sigma\_\{i\}\),\(F\.1\)wherei=\{1,…,Ntr\}i=\\\{1,\.\.\.,N\_\{\\text\{tr\}\}\\\}and𝒩m​\(0,Σi\)\\mathcal\{N\}\_\{m\}\(0,\\Sigma\_\{i\}\)is themm\-dimensional Gaussian distribution with mean zero and covariance matrixΣi\\Sigma\_\{i\}\.

Then the likelihood for the regression problem \([F\.1](https://arxiv.org/html/2610.19074#A6.E1)\) can be computed as

ℒ⁡\(θ,μ\)=∏i=1Ntrϕm​\(yi,g⁡\(xi,θ,μ\),Σi\),\\mathcal\{L\}\(\\theta,\\mu\)=\\prod\_\{i=1\}^\{N\_\{\\text\{tr\}\}\}\\phi\_\{m\}\(y\_\{i\};g\(x\_\{i\},\\theta,\\mu\),\\Sigma\_\{i\}\),whereψm​\(yi,g⁡\(xi,θ,μ\),Σi\)\\psi\_\{m\}\(y\_\{i\};g\(x\_\{i\},\\theta,\\mu\),\\Sigma\_\{i\}\)is themm\-dimensional Gaussian density of meang⁡\(xi,θ,μ\)g\(x\_\{i\},\\theta,\\mu\), covariance matrixΣi\\Sigma\_\{i\}evaluated atyiy\_\{i\}\.

In the discrete data set, like the classification problem, we have𝖸=\{1,…,m\}\\mathsf\{Y\}=\\\{1,\.\.\.,m\\\}, then we instead define the so\-calledsoftmaxfunction as

hk​\(x,θ,μ\)=exp⁡\(gk​\(x,θ,μ\)\)∑j=1mexp⁡\(gj​\(x,θ,μ\)\),k∈𝖸,h\_\{k\}\(x,\\theta,\\mu\)=\\frac\{\\exp\(g\_\{k\}\(x,\\theta,\\mu\)\)\}\{\\sum\_\{j=1\}^\{m\}\\exp\(g\_\{j\}\(x,\\theta,\\mu\)\)\},\\ k\\in\\mathsf\{Y\},\(F\.2\)and defineh⁡\(x,θ,μ\)=\(h1​\(x,θ,μ\),…,hm​\(x,θ,μ\)\)h\(x,\\theta,\\mu\)=\(h\_\{1\}\(x,\\theta,\\mu\),\.\.\.,h\_\{m\}\(x,\\theta,\\mu\)\)as a categorical distribution onmmoutcomes based on dataxx\. Then we assume thatyi∼i\.i\.d\.h⁡\(xi\)y\_\{i\}\\mathrel\{\\overset\{\\mathclap\{\\text\{\\tiny i\.i\.d\.\}\}\}\{\\sim\}\}h\(x\_\{i\}\)fori=\{1,…,Ntr\}i=\\\{1,\.\.\.,N\_\{\\text\{tr\}\}\\\}\.

The likelihood for the classification problem \([F\.2](https://arxiv.org/html/2610.19074#A6.E2)\) can be computed as

ℒ\(θ,μ\)=∏i=1Ntr∏k=1mhk\(xi,θ,μ\)𝟙\[yi=k\],\\mathcal\{L\}\(\\theta,\\mu\)=\\prod\_\{i=1\}^\{N\_\{\\text\{tr\}\}\}\\prod\_\{k=1\}^\{m\}h\_\{k\}\(x\_\{i\},\\theta,\\mu\)^\{\\mathds\{1\}\_\{\[y\_\{i\}=k\]\}\},where𝟙\[yi=k\]=1\\mathds\{1\}\_\{\[y\_\{i\}=k\]\}=1ifyi=ky\_\{i\}=kand00otherwise\.

#### F\.1\.1MNIST dataset

Consider the MNIST handwritten–digit dataset, which contains a total of60 00060\\,000training and10 00010\\,000test images of digits00–99\. Denote the firstNt​rN\_\{tr\}images from the training be the train dataset as\{\(xi,yi\)\}i=1Ntr⊂ℝ784×\{0,…,9\}\\\{\(x\_\{i\},\\,y\_\{i\}\)\\bigr\\\}\_\{i=1\}^\{N\_\{\\text\{tr\}\}\}\\subset\\mathbb\{R\}^\{784\}\\times\\\{0,\\dots,9\\\}and firstNtestN\_\{\\text\{test\}\}images from test be the test dataset as\{\(xjtest,yjtest\)\}j=1Ntest⊂ℝ784×\{0,…,9\}\\\{\(x\_\{j\}^\{\\text\{test\}\},\\,y\_\{j\}^\{\\text\{test\}\}\)\\bigr\\\}\_\{j=1\}^\{N\_\{\\text\{test\}\}\}\\subset\\mathbb\{R\}^\{784\}\\times\\\{0,\\dots,9\\\}, with every image vectorised and standardized to lie in\[0,1\]784\[0,1\]^\{784\}\.

### F\.2Comparison among MLSMC methods

We compare three MLSMC methods on the synthetic regression example using the Bayesian TNN model with fixed hyperparametersμ=\(σA,σb\)=\(1,1\)\\mu=\(\\sigma\_\{A\},\\sigma\_\{b\}\)=\(1,1\)\. The number of neural layers isD=3D=3, the input neural network width is then set ton=10n=10, and the output neural network width is set tom=1m=1\. The TNN prior uses a rate decay ofs=2s=2, and the pCN kernel is applied with1010mutation steps in three MLSMC methods\.

A regularity test is conducted to estimate the decay of bias and variance for normalized estimator differences between levels under the posterior distribution\. The results are shown in Figure[5](https://arxiv.org/html/2610.19074#A6.F5), supporting the decay rate assumed in Assumption[E\.2](https://arxiv.org/html/2610.19074#A5.Thmass2)with fixedμ\\mufor the convergence analysis\. Starting from levell=4l=4, the decay of biases and weights aligns with the expected rate\. In the MLSMC setting for TNN models, the bias decay is uniformly bounded due to Jensen’s inequality, ensuring at leastα=β/2\\alpha=\\beta/2, where the strong convergence rate isβ=2​s−1=3\\beta=2s\-1=3withs=2s=2\. It is important to note that the first level does not need to follow the rate, as it is treated as a single level\.

TakingLmin=3L\_\{\\text\{min\}\}=3as the starting level, the convergence test for the three MLSMC methods, \(r\)MLSMC and serial MLSMC, is presented in Figure[6](https://arxiv.org/html/2610.19074#A6.F6)\. Given a target accuracy levelε\\varepsilon, the maximum levelLmaxL\_\{\\max\}and the number of samples per levelNlN\_\{l\}in both MLSMC and serial MLSMC are followed the computation in Section[G\.1](https://arxiv.org/html/2610.19074#A7.SS1)\. For rMLSMC, level indiceslland correspondingNlN\_\{l\}are selected based on the index probability distribution, as described in Section[G\.2](https://arxiv.org/html/2610.19074#A7.SS2)\. The reference solution is computed at accuracyε=10−4\\varepsilon=10^\{\-4\}using MLSMC\.

![Refer to caption](https://arxiv.org/html/2610.19074v1/tnns_mean_norm.png)\(a\)
![Refer to caption](https://arxiv.org/html/2610.19074v1/tnns_var_norm.png)\(b\)

Figure 5:Regularity test for different MLSMC on Bayesian regression model with fixedμ\\mu:Empirical values for strong and weak convergence parameters in normalized estimator\. All plots use3030realizations with1000010000samples in each index, as the maximum levelLmax=7L\_\{\\max\}=7and starting level asLmin=3L\_\{\\text\{min\}\}=3at each experiment\. \(a\) Logarithm plot with base22for mean of difference and reference line with slope−3/2\-3/2; \(b\) Logarithm plot with base22for variance of difference and reference line with slope−3\-3\.![Refer to caption](https://arxiv.org/html/2610.19074v1/tnns_r_MLSMC.png)Figure 6:Complexity test for different MLSMC methods on Bayesian regression model:We use3030realisations to compute the MSE for each experiment and use the rate of convergence to compare the different methods\. Rate of regression: MLSMCSN\(red line\) is−0\.962\-0\.962, MLSMCRE \(yellow line\) is−0\.916\-0\.916, rMLSMC \(purple line\) is−0\.885\-0\.885and Serial MLSMC \(blue line\) is−0\.870\-0\.870\.
### F\.3Numerical justification for MLPMMH

![Refer to caption](https://arxiv.org/html/2610.19074v1/MLPMMH_reg_mean.png)\(a\)
![Refer to caption](https://arxiv.org/html/2610.19074v1/MLPMMH_reg_var.png)\(b\)

Figure 7:Regularity test for MLPMMH on Bayesian regression model:Empirical values for strong and weak convergence parameters in normalized estimator\. All plots use5050realizations withMl=100M\_\{l\}=100samples in each index, as the maximum levelLmax=8L\_\{\\max\}=8and starting level asLmin=1L\_\{\\text\{min\}\}=1at each experiment, with5050burn\-in samples\. For the MLSMC sampler,N=100N=100samples are used with an HMC kernel of one leapfrog step and1010mutation steps\. \(a\) Logarithm plot with base22for mean of difference and reference line with slope−3/2\-3/2; \(b\) Logarithm plot with base22for variance of difference and reference line with slope−3\-3\.The plots for autocorrelation functions ofσA\\sigma\_\{A\}andσb\\sigma\_\{b\}are given in Figure[8](https://arxiv.org/html/2610.19074#A6.F8), suggesting that a reasonable number of burn\-in samples in the MCMC in MLPMMH is5050\.

![Refer to caption](https://arxiv.org/html/2610.19074v1/mu_acc_l1.png)\(a\)l=1l=1:σA\\sigma\_\{A\}\(left\) andσb\\sigma\_\{b\}\(right\)
![Refer to caption](https://arxiv.org/html/2610.19074v1/mu_acc_l2.png)\(b\)l=2l=2:σA\\sigma\_\{A\}\(left\) andσb\\sigma\_\{b\}\(right\)
![Refer to caption](https://arxiv.org/html/2610.19074v1/mu_acc_l3.png)\(c\)l=3l=3:σA\\sigma\_\{A\}\(left\) andσb\\sigma\_\{b\}\(right\)
![Refer to caption](https://arxiv.org/html/2610.19074v1/mu_acc_l4.png)\(d\)l=4l=4:σA\\sigma\_\{A\}\(left\) andσb\\sigma\_\{b\}\(right\)
![Refer to caption](https://arxiv.org/html/2610.19074v1/mu_acc_l5.png)\(e\)l=5l=5:σA\\sigma\_\{A\}\(left\) andσb\\sigma\_\{b\}\(right\)
![Refer to caption](https://arxiv.org/html/2610.19074v1/mu_acc_l6.png)\(f\)l=6l=6:σA\\sigma\_\{A\}\(left\) andσb\\sigma\_\{b\}\(right\)
![Refer to caption](https://arxiv.org/html/2610.19074v1/mu_acc_l7.png)\(g\)l=7l=7:σA\\sigma\_\{A\}\(left\) andσb\\sigma\_\{b\}\(right\)
![Refer to caption](https://arxiv.org/html/2610.19074v1/mu_acc_l8.png)\(h\)l=8l=8:σA\\sigma\_\{A\}\(left\) andσb\\sigma\_\{b\}\(right\)

Figure 8:Autocorrelation functions forσA\\sigma\_\{A\}andσB\\sigma\_\{B\}on Bayesian regression model:All plots use5050realizations withMl=200M\_\{l\}=200samples in each indexllas the maximum levelLmax=8L\_\{\\max\}=8and starting level asLmin=1L\_\{\\text\{min\}\}=1at each experiment, with5050burn\-in samples\. For the MLSMC sampler,N=100N=100samples are used with an HMC kernel of one leapfrog step and1010mutation steps\.
### F\.4Empirical verification of TNN prior structure

We use a fully connected neural network \(D=3D=3\) with input widthn=784n=784, the hidden neural network widthn1=128n\_\{1\}=128\(equivalent tol=7l=7\), and the output neural network widthm=10m=10\. We compare three training regimes:

- •NoReg: standard training with no regularization, corresponding to pure data fitting;
- •Weight\-decay: MAP training with an isotropic Gaussian prior, corresponding to standardL2L\_\{2\}regularization;
- •Width\-decay: MAP training with the TNN\-style prior, where weights are regularized according to their expected variance under the TNN prior\. In particular, V​a​r​\[wi​j\]∝\{σ2​i−s,for the input\-to\-hidden layer,σ2​\(i​j\)−s,for hidden\-to\-hidden layer,σ2​j−s,for the hidden\-to\-output layer,Var\[w\_\{ij\}\]\\propto\\begin\{cases\}\\sigma^\{2\}i^\{\-s\},&\\text\{for the input\-to\-hidden layer\},\\\\ \\sigma^\{2\}\(ij\)^\{\-s\},&\\text\{for hidden\-to\-hidden layer\},\\\\ \\sigma^\{2\}j^\{\-s\},&\\text\{for the hidden\-to\-output layer\},\\end\{cases\}whereσ=σA\\sigma=\\sigma\_\{A\}for weights andσ=σb\\sigma=\\sigma\_\{b\}for biases\.

All models were initialized using the PyTorch default initialization and trained with the same Adam optimiser \(lr=10−310^\{\-3\}\) for2020epochs555This is the sufficient number of epochs to ensure convergence and avoid overfit for the model\.for NoReg, Weight\-decay, Width\-decay respectively\.

### F\.5Practical application of MLPMMH

We first verified that the PCA decay in the featured matrix of the MINST10 training set follows a similar decay as in the TNN\-style weight\-decay as in Figure[9](https://arxiv.org/html/2610.19074#A6.F9)\. The eigenvalues decay approximately linearly on the log\-log scale, with fitted slope−0\.615\-0\.615\. Interpreting this as a decay of the formj−s/2j^\{\-s/2\}gives the empirical values=1\.23s=1\.23\.

![Refer to caption](https://arxiv.org/html/2610.19074v1/MNIST_pca_decay.png)Figure 9:PCA featured matrix on the whole MNIST training dataset:We then numerically justify the PCA hierarchy in levelllwith Gaussian prior decays in bias and variance as required in Assumption[E\.2](https://arxiv.org/html/2610.19074#A5.Thmass2)for the convergence of MLPMMH, see Figure[10](https://arxiv.org/html/2610.19074#A6.F10)\. It shows the decay of the prior predictive increments in bias and variance\. For a NN\-based experiment, the cost rateγ=1\\gamma=1\. From the Figure[10](https://arxiv.org/html/2610.19074#A6.F10), we can observe the variance decay rateβ=1\.408\\beta=1\.408, which is close to the theoretical value2​s−12s\-1\. Also note thatβ\>1=γ\\beta\>1=\\gamma, hence satisfying all the conditions for the convergence of MLPMMH\.

![Refer to caption](https://arxiv.org/html/2610.19074v1/prior_bias_decay.png)\(a\)
![Refer to caption](https://arxiv.org/html/2610.19074v1/prior_var_decay.png)\(b\)

Figure 10:Regularity test for MLPMMH on MNIST dataset w\.r\.t prior:Empirical values for strong and weak convergence parameters under prior distributionπ0\\pi\_\{0\}\. All plots use100100realizations withMl=100M\_\{l\}=100samples in each index as the maximum levelLmax=8L\_\{\\max\}=8and starting level asLmin=0L\_\{\\text\{min\}\}=0at each experiment\. \(a\) Logarithm plot with base22for mean of difference and reference line with slope−0\.690\-0\.690; \(b\) Logarithm plot with base22for variance of difference and reference line with slope−1\.408\-1\.408\.

## Appendix GPractical settings in MLSMC methods

### G\.1Multilevel sequential Monte Carlo estimators

First, the empirical value ofLminL\_\{\\text\{min\}\}is determined by the regularity results, where theBBandVVneeds to decay followed the expected rate in Assumption[E\.2](https://arxiv.org/html/2610.19074#A5.Thmass2)with fixedμ\\mufor each levelLmin<l≤LmaxL\_\{\\text\{min\}\}<l\\leq L\_\{\\text\{max\}\}\. Suppose we know the parameter values and the constants, we can compute the empirical values forLmaxL\_\{\\text\{max\}\}andMlM\_\{l\}if we want the MSE to be smaller thanε2\\varepsilon^\{2\}, i\.e\. the variance and the bias terms are all smaller thanε22\\frac\{\\varepsilon^\{2\}\}\{2\}, using the complexity theorem of MLSMC \(Theorem[E\.2](https://arxiv.org/html/2610.19074#A5.Thmtheorem2)\):

- •The value ofLmaxL\_\{\\text\{max\}\}: If\|𝔼⁡\[φl−φ\]\|≤c1​hlα\|\\mathbb\{E\}\[\\varphi\_\{l\}\-\\varphi\]\|\\leq c\_\{1\}h\_\{l\}^\{\\alpha\}andhl=hlh\_\{l\}=h^\{l\}, we have the bias term as 𝔼⁡\[φL−φ\]=∑l=L\+1∞𝔼⁡\[φl−φl−1\]=𝔼⁡\[φL−φL−1\]h−α−1\.\\mathbb\{E\}\[\\varphi\_\{L\}\-\\varphi\]=\\sum\_\{l=L\+1\}^\{\\infty\}\\mathbb\{E\}\[\\varphi\_\{l\}\-\\varphi\_\{l\-1\}\]=\\frac\{\\mathbb\{E\}\[\\varphi\_\{L\}\-\\varphi\_\{L\-1\}\]\}\{h^\{\-\\alpha\}\-1\}\.\(G\.1\)In order to ensure𝔼​\[φL−φ\]2≤ε22\\mathbb\{E\}\[\\varphi\_\{L\}\-\\varphi\]^\{2\}\\leq\\frac\{\\varepsilon^\{2\}\}\{2\}, we need 𝔼​\[φL−φL−1\]2\(h−α−1\)2≤c12​hL2​α\(h−α−1\)2≤ε22\\displaystyle\\frac\{\\mathbb\{E\}\[\\varphi\_\{L\}\-\\varphi\_\{L\-1\}\]^\{2\}\}\{\(h^\{\-\\alpha\}\-1\)^\{2\}\}\\leq\\frac\{c\_\{1\}^\{2\}h^\{2\\alpha\}\_\{L\}\}\{\(h^\{\-\\alpha\}\-1\)^\{2\}\}\\leq\\frac\{\\varepsilon^\{2\}\}\{2\}\(G\.2\)⇒\\displaystyle\\RightarrowL≥1α​logh−1⁡\(2​c1ε⁡\(h−α−1\)\)\.\\displaystyle L\\geq\\frac\{1\}\{\\alpha\}\\log\_\{h^\{\-1\}\}\\Big\(\\frac\{\\sqrt\{2\}c\_\{1\}\}\{\\varepsilon\(h^\{\-\\alpha\}\-1\)\}\\Big\)\.We then takeLmax=⌈1α​logh−1⁡\(2​c1ε⁡\(h−α−1\)\)⌉L\_\{\\text\{max\}\}=\\Big\\lceil\\frac\{1\}\{\\alpha\}\\log\_\{h^\{\-1\}\}\\Big\(\\frac\{\\sqrt\{2\}c\_\{1\}\}\{\\varepsilon\(h^\{\-\\alpha\}\-1\)\}\\Big\)\\Big\\rceil\.
- •The value ofMlM\_\{l\}: From the observedLminL\_\{\\text\{min\}\}and estimatedLmaxL\_\{\\text\{max\}\}above, we can estimate the number of samples at each level, which is given by Ml=⌈2​ε−2​Vl/Cl​\(∑l=LminLmaxVl​Cl\)⌉,M\_\{l\}=\\Big\\lceil 2\\varepsilon^\{\-2\}\\sqrt\{V\_\{l\}/C\_\{l\}\}\\Big\(\\sum\_\{l=L\_\{\\text\{min\}\}\}^\{L\_\{\\text\{max\}\}\}\\sqrt\{V\_\{l\}C\_\{l\}\}\\Big\)\\Big\\rceil,\(G\.3\)to ensure the variance term is smaller thanε22\\frac\{\\varepsilon^\{2\}\}\{2\}, where⌈⋅⌉\\lceil\\cdot\\rceildenotes the smallest integer greater than the number inside\.

When the Assumption[E\.2](https://arxiv.org/html/2610.19074#A5.Thmass2)is not satisfied, it leads to a sub\-canonical case, stated in Remark[G\.1](https://arxiv.org/html/2610.19074#A7.Thmrem1)\. Moreover, since each level is independent, we can parallelize the algorithm for all the levels to accelerate the sampler\.

###### Remark G\.1\(Sub\-canonical regime in MLSMC\)\.

In the sub\-canonical case, where the variance decays more slowly than the cost grows, i\.e\.β≤γ\\beta\\leq\\gamma, the decay rate violates in the Assumption[E\.2](https://arxiv.org/html/2610.19074#A5.Thmass2)\. The MLSMC estimators incur an additional exponent in theε\\varepsilon–COST complexity\. Specifically, for levelll, if the bias decays like𝒪⁡\(2−α​l\)\\mathcal\{O\}\(2^\{\-\\alpha l\}\)and the cost grows likeO⁡\(2γ​l\)O\(2^\{\\gamma l\}\), then to achieveMSE≤ε2\\mathrm\{MSE\}\\leq\\varepsilon^\{2\}, one obtains

COST=𝒪⁡\(ε−2−\(γ−β\)/α\),\\mathrm\{COST\}=\\mathcal\{O\}\\bigl\(\\varepsilon^\{\-2\-\(\\gamma\-\\beta\)/\\alpha\}\\bigr\),rather than the canonicalO⁡\(ε−2\)O\(\\varepsilon^\{\-2\}\)whenβ\>γ\\beta\>\\gamma\.

### G\.2Randomized multilevel sequential Monte Carlo estimators

Again, the empirical value ofLminL\_\{\\text\{min\}\}is determined by the regularity results, where theBBandVVneeds to decay followed the expected rate in Assumption[E\.2](https://arxiv.org/html/2610.19074#A5.Thmass2)with fixedμ\\mufor each levelLmin<l≤LmaxL\_\{\\text\{min\}\}<l\\leq L\_\{\\text\{max\}\}\. Proposition 3 in[Law et al\. \[2022\]](https://arxiv.org/html/2610.19074#bib.bib5)states the requirement of

∑l∈ℙ1pl​2−l​β≤C\\sum\_\{l\\in\\mathbb\{P\}\}\\frac\{1\}\{p\_\{l\}\}2^\{\-l\\beta\}\\leq C\(G\.4\)to let𝔼⁡\[\(∑l∈ℙNlpl​N​FlNl​\(Δ​ζ\)−f⁡\(ζ\)\)2\]\\mathbb\{E\}\\big\[\\big\(\\sum\_\{l\\in\\mathbb\{P\}\}\\frac\{N\_\{l\}\}\{p\_\{l\}N\}F^\{N\_\{l\}\}\_\{l\}\(\\Delta\\zeta\)\-f\(\\zeta\)\\big\)^\{2\}\\big\]be of𝒪⁡\(N−1\)\\mathcal\{O\}\(N^\{\-1\}\), so that the convergence rate in Theorem[E\.3](https://arxiv.org/html/2610.19074#A5.Thmtheorem3)\. Additionally, we let the distribution follow from a simple constrained minimization of the expected cost for one sample

∑l∈ℙpl​2l​γ≤C\.\\sum\_\{l\\in\\mathbb\{P\}\}p\_\{l\}2^\{l\\gamma\}\\leq C\.\(G\.5\)Sincepl\>0p\_\{l\}\>0, one finds that

pl∝2−l\(γ\+β\)/2p\_\{l\}\\propto 2^\{\-l\(\\gamma\+\\beta\)/2\}\(G\.6\)is one suitable choice for probabilities to satisfy the two conditions stated in \([G\.4](https://arxiv.org/html/2610.19074#A7.E4)\) and \([G\.5](https://arxiv.org/html/2610.19074#A7.E5)\)\. Then define the probability distribution𝐩\{\\bf p\}onℙ\\mathbb\{P\}with𝐩⁡\(l\)=:pl\>0\{\\bf p\}\(l\)=:p\_\{l\}\>0\. Consider drawingN/NminN/N\_\{\\rm min\}i\.i\.d\. samplesli∼𝐩l\_\{i\}\\sim\{\\bf p\}, and define the number of samplesNl=Nmin\#i\{li=l\}N\_\{l\}=N\_\{\\rm min\}\\\#i\\\{l\_\{i\}=l\\\}\. Like in the MLSMC, since each level is independent, we can parallel the algorithm for all the levels to accelerate this sampler as well\. It is worth noting that when the Assumption[E\.2](https://arxiv.org/html/2610.19074#A5.Thmass2)is not satisfied, it leads to a sub\-canonical case like in the MLSMC, stated in Remark[G\.2](https://arxiv.org/html/2610.19074#A7.Thmrem2)\.

###### Remark G\.2\(Sub\-canonical regime in rMLSMC\[[Chada et al\., 2021](https://arxiv.org/html/2610.19074#bib.bib37)\]\)\.

Similarly, in the randomized MLSMC estimator, even with randomization that eliminates all bias, the fundamental rate‐limiting factor remains the mismatch between variance decay and cost growth\. Let the sampling probabilities be

pl∝2−β​l,p\_\{l\}\\propto 2^\{\-\\beta l\},which is the asymptotically optimal “borderline” choice, so that the estimator achieves a sub‐canonical rate with a logarithmic penalty:

COST=𝒪⁡\(ε−2​\|log⁡ε\|q\),\\mathrm\{COST\}=\\mathcal\{O\}\\bigl\(\\varepsilon^\{\-2\}\\,\\lvert\\log\\varepsilon\\rvert^\{q\}\\bigr\),for anyq\>2q\>2\.

## Appendix HProofs of Theoretical results

The proofs of the various results in the paper, along with restatements of the results and verification of the assumptions, are presented here\.

### H\.1Proof relating to Theorem[4\.1](https://arxiv.org/html/2610.19074#S4.Thmtheorem1)

Let us assume Assumptions[E\.1](https://arxiv.org/html/2610.19074#A5.Thmass1),[E\.3](https://arxiv.org/html/2610.19074#A5.Thmass3)and[E\.2](https://arxiv.org/html/2610.19074#A5.Thmass2)withβ\>γ\\beta\>\\gamma\. Then, for suitableφ:Θ×𝖲→ℝ\\varphi:\\Theta\\times\\mathsf\{S\}\\rightarrow\\mathbb\{R\}, there exists aC<\+∞C<\+\\inftysuch that for anyε\>0\\varepsilon\>0one can choose\(L=\{Lmin,Lmax\},\{Ml\}l=LminLmax\)\(L=\\\{L\_\{\\text\{min\}\},L\_\{\\text\{max\}\}\\\},\\\{M\_\{l\}\\\}\_\{l=L\_\{\\text\{min\}\}\}^\{L\_\{\\text\{max\}\}\}\)such that

𝔼⁡\[\(φ^LMLPMMH−π⁡\(φ\)\)2\]≤C​ε2,\\mathbb\{E\}\[\(\\hat\{\\varphi\}\_\{L\}^\{\\mathrm\{MLPMMH\}\}\-\\pi\(\\varphi\)\)^\{2\}\]\\leq C\\varepsilon^\{2\},and expected COST\(φ^LMLPMMH\)≤C​ε−2\(\\hat\{\\varphi\}\_\{L\}^\{\\mathrm\{MLPMMH\}\}\)\\leq C\\varepsilon^\{\-2\}, i\.e\. the canonical rate\. The estimator is defined in equation \([3\.1](https://arxiv.org/html/2610.19074#S3.E1)\)\.

###### Proof\.

We first define

Δ​Πl​\(φl\):=Πl​\(Δ​φl\)=Πl​\(φl−φl−1\),\\Delta\\Pi\_\{l\}\(\\varphi\_\{l\}\):=\\Pi\_\{l\}\(\\Delta\\varphi\_\{l\}\)=\\Pi\_\{l\}\(\\varphi\_\{l\}\-\\varphi\_\{l\-1\}\),whereΠl\\Pi\_\{l\}is the exact coupling\. Note that

Δ​Πl​\(φl\)=πl​\(φl\)−πl−1​\(φl−1\)\.\\Delta\\Pi\_\{l\}\(\\varphi\_\{l\}\)=\\pi\_\{l\}\(\\varphi\_\{l\}\)\-\\pi\_\{l\-1\}\(\\varphi\_\{l\-1\}\)\.Noting the independence between differences, we derive

𝔼⁡\[\(φ^LMLPMMH−π⁡\(φ\)\)2\]\\displaystyle\\mathbb\{E\}\[\(\\hat\{\\varphi\}\_\{L\}^\{\\mathrm\{MLPMMH\}\}\-\\pi\(\\varphi\)\)^\{2\}\]=\(πL​\(φL\)−π⁡\(φ\)\)2\+𝔼⁡\[\(φ^LMLPMMH−πL​\(φL\)\)2\]\\displaystyle=\\Big\(\\pi\_\{L\}\(\\varphi\_\{L\}\)\-\\pi\(\\varphi\)\\Big\)^\{2\}\+\\mathbb\{E\}\\Bigg\[\\Big\(\\hat\{\\varphi\}\_\{L\}^\{\\mathrm\{MLPMMH\}\}\-\\pi\_\{L\}\(\\varphi\_\{L\}\)\\Big\)^\{2\}\\Bigg\]\(H\.1\)=\(πL​\(φL\)−π⁡\(φ\)\)2\+𝔼⁡\[\(∑l=0LΔ​ΠlN,Ml​\(φl\)−Δ​Πl​\(φl\)\)2\]\\displaystyle=\\Big\(\\pi\_\{L\}\(\\varphi\_\{L\}\)\-\\pi\(\\varphi\)\\Big\)^\{2\}\+\\mathbb\{E\}\\Bigg\[\\Big\(\\sum\_\{l=0\}^\{L\}\\Delta\\Pi\_\{l\}^\{N,M\_\{l\}\}\(\\varphi\_\{l\}\)\-\\Delta\\Pi\_\{l\}\(\\varphi\_\{l\}\)\\Big\)^\{2\}\\Bigg\]=\(πL​\(φL\)−π⁡\(φ\)\)2\+∑l=0L𝔼⁡\[\(Δ​ΠlN,Ml​\(φl\)−Δ​Πl​\(φl\)\)2\]\\displaystyle=\\Big\(\\pi\_\{L\}\(\\varphi\_\{L\}\)\-\\pi\(\\varphi\)\\Big\)^\{2\}\+\\sum\_\{l=0\}^\{L\}\\mathbb\{E\}\\Bigg\[\\Big\(\\Delta\\Pi\_\{l\}^\{N,M\_\{l\}\}\(\\varphi\_\{l\}\)\-\\Delta\\Pi\_\{l\}\(\\varphi\_\{l\}\)\\Big\)^\{2\}\\Bigg\]\+∑q≠l=0L𝔼\[\(ΔΠqN,Mq\(φq\)−ΔΠq\(φq\)\)\(ΔΠlN,Ml\(φl\)−ΔΠl\(φl\)\)\]\.\\displaystyle\+\\sum\_\{q\\neq l=0\}^\{L\}\\mathbb\{E\}\\Bigg\[\\Big\(\\Delta\\Pi\_\{q\}^\{N,M\_\{q\}\}\(\\varphi\_\{q\}\)\-\\Delta\\Pi\_\{q\}\(\\varphi\_\{q\}\)\\Big\)\\Big\(\\Delta\\Pi\_\{l\}^\{N,M\_\{l\}\}\(\\varphi\_\{l\}\)\-\\Delta\\Pi\_\{l\}\(\\varphi\_\{l\}\)\\Big\)\\Bigg\]\.In order for us to process, we need to consider three terms, i\.e\. the bias term, the variance term and the covariance term\. We letβ=\(2​s−1\)\\beta=\(2s\-1\)\. Let us consider first the bias term, where a direct application of Assumption[E\.2](https://arxiv.org/html/2610.19074#A5.Thmass2)gives

\(πL​\(φL\)−π⁡\(φ\)\)2\\displaystyle\\Big\(\\pi\_\{L\}\(\\varphi\_\{L\}\)\-\\pi\(\\varphi\)\\Big\)^\{2\}=\(∫Θ×𝖲\(φL​\(θL,μ\)−φ⁡\(θ,μ\)\)​π​\(𝑑θ,𝑑μ\)\)2\\displaystyle=\\Bigg\(\\int\_\{\\Theta\\times\\mathsf\{S\}\}\(\\varphi\_\{L\}\(\\theta\_\{L\},\\mu\)\-\\varphi\(\\theta,\\mu\)\)\\pi\(d\\theta,d\\mu\)\\Bigg\)^\{2\}≤C​2−2​α​L≤C​2−β​L\.\\displaystyle\\leq C2^\{\-2\\alpha L\}\\leq C2^\{\-\\beta L\}\.\(H\.2\)
Next, we consider the Variance term, which by Theorem[E\.5](https://arxiv.org/html/2610.19074#A5.Thmtheorem5)we have

𝔼⁡\[\(Δ​ΠlN,Ml​\(φl\)−Δ​Πl​\(φl\)\)2\]≤C​2−β​lMl\.\\mathbb\{E\}\\Bigg\[\\Big\(\\Delta\\Pi\_\{l\}^\{N,M\_\{l\}\}\(\\varphi\_\{l\}\)\-\\Delta\\Pi\_\{l\}\(\\varphi\_\{l\}\)\\Big\)^\{2\}\\Bigg\]\\leq\\frac\{C2^\{\-\\beta l\}\}\{M\_\{l\}\}\.
By the Cauchy\-Schwarz inequality, we have for the covariance term

𝔼⁡\[\(Δ​ΠqN,Mq​\(φq\)−Δ​Πq​\(φq\)\)​\(Δ​ΠlN,Ml​\(φl\)−Δ​Πl​\(φl\)\)\]\\displaystyle\\mathbb\{E\}\\Bigg\[\\Big\(\\Delta\\Pi\_\{q\}^\{N,M\_\{q\}\}\(\\varphi\_\{q\}\)\-\\Delta\\Pi\_\{q\}\(\\varphi\_\{q\}\)\\Big\)\\Big\(\\Delta\\Pi\_\{l\}^\{N,M\_\{l\}\}\(\\varphi\_\{l\}\)\-\\Delta\\Pi\_\{l\}\(\\varphi\_\{l\}\)\\Big\)\\Bigg\]≤𝔼⁡\[\(Δ​ΠqN,Mq​\(φq\)−Δ​Πq​\(φq\)\)2\]​𝔼⁡\[\(Δ​ΠlN,Ml​\(φl\)−Δ​Πl​\(φl\)\)2\]\\displaystyle\\leq\\sqrt\{\\mathbb\{E\}\\Bigg\[\\Big\(\\Delta\\Pi\_\{q\}^\{N,M\_\{q\}\}\(\\varphi\_\{q\}\)\-\\Delta\\Pi\_\{q\}\(\\varphi\_\{q\}\)\\Big\)^\{2\}\\Bigg\]\}\\sqrt\{\\mathbb\{E\}\\Bigg\[\\Big\(\\Delta\\Pi\_\{l\}^\{N,M\_\{l\}\}\(\\varphi\_\{l\}\)\-\\Delta\\Pi\_\{l\}\(\\varphi\_\{l\}\)\\Big\)^\{2\}\\Bigg\]\}≤C′​2−β​qMq​C′′​2−β​lMl≤C​2−β\(q\+l\)/2Mq​Ml\.\\displaystyle\\leq\\sqrt\{\\frac\{C^\{\\prime\}2^\{\-\\beta q\}\}\{M\_\{q\}\}\}\\sqrt\{\\frac\{C^\{\\prime\\prime\}2^\{\-\\beta l\}\}\{M\_\{l\}\}\}\\leq C\\frac\{2^\{\-\\beta\(q\+l\)/2\}\}\{\\sqrt\{M\_\{q\}M\_\{l\}\}\}\.
Then, combining the results above, we have

\([H\.1](https://arxiv.org/html/2610.19074#A8.E1)\)≤C′​2−β​L\+C′′​∑l=0L2−β​lMl\+C′′′​∑q≠l=0L2−β\(q\+l\)/2Mq​Ml\.\\displaystyle\\eqref\{eqn:mse1\}\\leq C^\{\\prime\}2^\{\-\\beta L\}\+C^\{\\prime\\prime\}\\sum\_\{l=0\}^\{L\}\\frac\{2^\{\-\\beta l\}\}\{M\_\{l\}\}\+C^\{\\prime\\prime\\prime\}\\sum\_\{q\\neq l=0\}^\{L\}\\frac\{2^\{\-\\beta\(q\+l\)/2\}\}\{\\sqrt\{M\_\{q\}M\_\{l\}\}\}\.\(H\.3\)
In order to ensure the MSE is𝒪⁡\(ε2\)\\mathcal\{O\}\(\\varepsilon^\{2\}\), then all three terms in \([H\.3](https://arxiv.org/html/2610.19074#A8.E3)\) need to be𝒪⁡\(ε2\)\\mathcal\{O\}\(\\varepsilon^\{2\}\)\. Following from the first term in \([H\.3](https://arxiv.org/html/2610.19074#A8.E3)\), it will be satisfied to letL∝2​\|log⁡\(ε\)\|/βL\\propto 2\|\\log\(\\varepsilon\)\|/\\beta\. For the the second term in \([H\.3](https://arxiv.org/html/2610.19074#A8.E3)\), suppose we minimize total COST=∑l=0L2−γ​l​Ml=\\sum\_\{l=0\}^\{L\}2^\{\-\\gamma l\}M\_\{l\}subject to∑l=0L2−β​lMl=𝒪⁡\(ε2\)\\sum\_\{l=0\}^\{L\}\\frac\{2^\{\-\\beta l\}\}\{M\_\{l\}\}=\\mathcal\{O\}\(\\varepsilon^\{2\}\), we solve the number of samples at each levelM0,M1,…,MLM\_\{0\},M\_\{1\},\.\.\.,M\_\{L\}\. This gives

Ml∝ε−2KL2−\(β\+γ\)l/2,M\_\{l\}\\propto\\varepsilon^\{\-2\}K\_\{L\}2^\{\-\(\\beta\+\\gamma\)l/2\},\(H\.4\)whereKL=∑l=0L2−\(β−γ\)l/2K\_\{L\}=\\sum\_\{l=0\}^\{L\}2^\{\-\(\\beta\-\\gamma\)l/2\}\. This gives the total cost of𝒪⁡\(ε−2\)\\mathcal\{O\}\(\\varepsilon^\{\-2\}\)\. Finally, we need to check that the third term is of𝒪⁡\(ε2\)\\mathcal\{O\}\(\\varepsilon^\{2\}\):

∑q≠l=0L2−β\(q\+l\)/2Mq​Ml=∑q≠l=0L2−\(β−γ\)\(q\+l\)/4KLε2\.\\sum\_\{q\\neq l=0\}^\{L\}\\frac\{2^\{\-\\beta\(q\+l\)/2\}\}\{\\sqrt\{M\_\{q\}M\_\{l\}\}\}=\\sum\_\{q\\neq l=0\}^\{L\}2^\{\-\(\\beta\-\\gamma\)\(q\+l\)/4\}K\_\{L\}\\varepsilon^\{2\}\.\(H\.5\)Since we assumeβ\>γ\\beta\>\\gamma, the quantity∑q≠l=0L2−\(β−γ\)\(q\+l\)/4KL<∞\\sum\_\{q\\neq l=0\}^\{L\}2^\{\-\(\\beta\-\\gamma\)\(q\+l\)/4\}K\_\{L\}<\\inftyis finite, and finally

∑q≠l=0L2−β\(q\+l\)/2Mq​Ml=O⁡\(ε2\)\.\\sum\_\{q\\neq l=0\}^\{L\}\\frac\{2^\{\-\\beta\(q\+l\)/2\}\}\{\\sqrt\{M\_\{q\}M\_\{l\}\}\}=O\(\\varepsilon^\{2\}\)\.Therefore, this completes the proof\.

∎

### H\.2Proof relating to Theorem[E\.5](https://arxiv.org/html/2610.19074#A5.Thmtheorem5)

Assume Assumptions[E\.1](https://arxiv.org/html/2610.19074#A5.Thmass1),[E\.2](https://arxiv.org/html/2610.19074#A5.Thmass2)\. For anyφ:Θ×𝖲→ℝd\\varphi:\\Theta\\times\\mathsf\{S\}\\rightarrow\\mathbb\{R\}^\{d\}, there exists aC<\+∞C<\+\\inftysuch that

𝔼⁡\[\(ΠlN,Ml​\(Δ​φl\)−Πl​\(Δ​φl\)\)2\]≤C​2−β​lMl\.\\mathbb\{E\}\\Bigg\[\\Big\(\\Pi\_\{l\}^\{N,M\_\{l\}\}\(\\Delta\\varphi\_\{l\}\)\-\\Pi\_\{l\}\(\\Delta\\varphi\_\{l\}\)\\Big\)^\{2\}\\Bigg\]\\leq\\frac\{C2^\{\-\\beta l\}\}\{M\_\{l\}\}\.\(H\.6\)
###### Proof\.

Writing out \([H\.6](https://arxiv.org/html/2610.19074#A8.E6)\), we have

𝔼\[\\displaystyle\\mathbb\{E\}\\Bigg\[\(1Ml​∑k=1Mlφl​\(θk,l,μk,l\)​fl​\(θk,l\|μk,l\)Fl​\(θk,l,θk,l−1\|μk,l\)1Ml​∑k=1Mlfl​\(θk,l\|μk,l\)Fl​\(θk,l,θk,l−1\|μk,l\)−1Ml​∑k=1Mlφl−1​\(θk,l−1,μk,l\)​fl−1​\(θk,l−1\|μk,l\)Fl​\(θk,l,θk,l−1\|μk,l\)1Ml​∑k=1Mlfl−1​\(θk,l−1\|μk,l\)Fl​\(θk,l,θk,l−1\|μk,l\)CLOSE\\displaystyle\\Bigg\(\\frac\{\\frac\{1\}\{M\_\{l\}\}\\sum\_\{k=1\}^\{M\_\{l\}\}\\varphi\_\{l\}\(\\theta\_\{k,l\},\\mu\_\{k,l\}\)\\frac\{f\_\{l\}\(\\theta\_\{k,l\}\|\\mu\_\{k,l\}\)\}\{F\_\{l\}\(\\theta\_\{k,l\},\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)\}\}\{\\frac\{1\}\{M\_\{l\}\}\\sum\_\{k=1\}^\{M\_\{l\}\}\\frac\{f\_\{l\}\(\\theta\_\{k,l\}\|\\mu\_\{k,l\}\)\}\{F\_\{l\}\(\\theta\_\{k,l\},\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)\}\}\-\\frac\{\\frac\{1\}\{M\_\{l\}\}\\sum\_\{k=1\}^\{M\_\{l\}\}\\varphi\_\{l\-1\}\(\\theta\_\{k,l\-1\},\\mu\_\{k,l\}\)\\frac\{f\_\{l\-1\}\(\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)\}\{F\_\{l\}\(\\theta\_\{k,l\},\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)\}\}\{\\frac\{1\}\{M\_\{l\}\}\\sum\_\{k=1\}^\{M\_\{l\}\}\\frac\{f\_\{l\-1\}\(\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)\}\{F\_\{l\}\(\\theta\_\{k,l\},\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)\}\}−\(fl​\(φl​\(θl,μ\)\)Zl−fl−1​\(φl−1​\(θl−1,μ\)\)Zl−1\)\)2\]\.\\displaystyle\-\\bigg\(\\frac\{f\_\{l\}\(\\varphi\_\{l\}\(\\theta\_\{l\},\\mu\)\)\}\{Z\_\{l\}\}\-\\frac\{f\_\{l\-1\}\(\\varphi\_\{l\-1\}\(\\theta\_\{l\-1\},\\mu\)\)\}\{Z\_\{l\-1\}\}\\bigg\)\\Bigg\)^\{2\}\\Bigg\]\.Now we simplify the notation, let

a\\displaystyle a=1Ml​∑k=1Mlφl​\(θk,l,μk,l\)​fl​\(θk,l\|μk,l\)Fl​\(θk,l,θk,l−1\|μk,l\),\\displaystyle=\\frac\{1\}\{M\_\{l\}\}\\sum\_\{k=1\}^\{M\_\{l\}\}\\varphi\_\{l\}\(\\theta\_\{k,l\},\\mu\_\{k,l\}\)\\frac\{f\_\{l\}\(\\theta\_\{k,l\}\|\\mu\_\{k,l\}\)\}\{F\_\{l\}\(\\theta\_\{k,l\},\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)\},A\\displaystyle A=1Ml​∑k=1Mlfl​\(θk,l\|μk,l\)Fl​\(θk,l,θk,l−1\|μk,l\),\\displaystyle=\\frac\{1\}\{M\_\{l\}\}\\sum\_\{k=1\}^\{M\_\{l\}\}\\frac\{f\_\{l\}\(\\theta\_\{k,l\}\|\\mu\_\{k,l\}\)\}\{F\_\{l\}\(\\theta\_\{k,l\},\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)\},b\\displaystyle b=1Ml​∑k=1Mlφl−1​\(θk,l−1,μk,l\)​fl−1​\(θk,l−1\|μk,l\)Fl​\(θk,l,θk,l−1\|μk,l\),\\displaystyle=\\frac\{1\}\{M\_\{l\}\}\\sum\_\{k=1\}^\{M\_\{l\}\}\\varphi\_\{l\-1\}\(\\theta\_\{k,l\-1\},\\mu\_\{k,l\}\)\\frac\{f\_\{l\-1\}\(\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)\}\{F\_\{l\}\(\\theta\_\{k,l\},\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)\},B\\displaystyle B=1Ml​∑k=1Mlfl−1​\(θk,l−1\|μk,l\)Fl​\(θk,l,θk,l−1\|μk,l\),\\displaystyle=\\frac\{1\}\{M\_\{l\}\}\\sum\_\{k=1\}^\{M\_\{l\}\}\\frac\{f\_\{l\-1\}\(\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)\}\{F\_\{l\}\(\\theta\_\{k,l\},\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)\},c\\displaystyle c=fl​\(φl​\(θl,μ\)\),C=Zl,d=fl−1​\(φl−1​\(θl−1,μ\)\),D=Zl−1\.\\displaystyle=f\_\{l\}\(\\varphi\_\{l\}\(\\theta\_\{l\},\\mu\)\),\\quad C=Z\_\{l\},\\quad d=f\_\{l\-1\}\(\\varphi\_\{l\-1\}\(\\theta\_\{l\-1\},\\mu\)\),\\quad D=Z\_\{l\-1\}\.
By Lemma[E\.1](https://arxiv.org/html/2610.19074#A5.Thmlem1)and Jensen’s inequality, we have

\(aA−CLOSE\\displaystyle\\Bigg\(\\frac\{a\}\{A\}\-OPENbB−\(cC−dD\)\)2=\(\[a−b−\(c−d\)\]A−b⁡\[A−B−\(C−D\)\]A​B\+1A​C​\[C−A\]​\[c−d\]CLOSE\\displaystyle\\frac\{b\}\{B\}\-\\Bigg\(\\frac\{c\}\{C\}\-\\frac\{d\}\{D\}\\Bigg\)\\Bigg\)^\{2\}=\\Bigg\(\\frac\{\[a\-b\-\(c\-d\)\]\}\{A\}\-\\frac\{b\[A\-B\-\(C\-D\)\]\}\{AB\}\+\\frac\{1\}\{AC\}\[C\-A\]\[c\-d\]OPEN−1A​B​\(b−d\)​\(C−D\)\+dC​B​D​\(B−D\)​\(C−D\)\+dA​C​B​\(A−C\)​\(C−D\)\)2\\displaystyle\-\\frac\{1\}\{AB\}\(b\-d\)\(C\-D\)\+\\frac\{d\}\{CBD\}\(B\-D\)\(C\-D\)\+\\frac\{d\}\{ACB\}\(A\-C\)\(C\-D\)\\Bigg\)^\{2\}≤6\[\(\[a−b−\(c−d\)\]A\)2\+\(b⁡\[A−B−\(C−D\)\]A​B\)2\+\(1A​C\[C−A\]\[c−d\]\)2\\displaystyle\\leq 6\\Bigg\[\\Bigg\(\\frac\{\[a\-b\-\(c\-d\)\]\}\{A\}\\Bigg\)^\{2\}\+\\Bigg\(\\frac\{b\[A\-B\-\(C\-D\)\]\}\{AB\}\\Bigg\)^\{2\}\+\\Bigg\(\\frac\{1\}\{AC\}\[C\-A\]\[c\-d\]\\Bigg\)^\{2\}\+\(1A​B\(b−d\)\(C−D\)\)2\+\(dC​B​D\(B−D\)\(C−D\)\)2\+\(dA​C​B\(A−C\)\(C−D\)\)2\]\.\\displaystyle\+\\Bigg\(\\frac\{1\}\{AB\}\(b\-d\)\(C\-D\)\\Bigg\)^\{2\}\+\\Bigg\(\\frac\{d\}\{CBD\}\(B\-D\)\(C\-D\)\\Bigg\)^\{2\}\+\\Bigg\(\\frac\{d\}\{ACB\}\(A\-C\)\(C\-D\)\\Bigg\)^\{2\}\\Bigg\]\.\(H\.7\)
Now, we need to prove that all six terms in \([H\.7](https://arxiv.org/html/2610.19074#A8.E7)\) are𝒪⁡\(2−β​lMl\)\\mathcal\{O\}\(\\frac\{2^\{\-\\beta l\}\}\{M\_\{l\}\}\)\. By Proposition[E\.2](https://arxiv.org/html/2610.19074#A5.Thmprop2), the first term is verified directly and the second term is also verified by lettingφ==1\\varphi==1\.

By Assumption[E\.2](https://arxiv.org/html/2610.19074#A5.Thmass2)\(V\) and \(B\) withα≤β/2\\alpha\\leq\\beta/2, we have

𝔼⁡\[\(b−d\)2\]≤C​2−β​lMl,𝔼⁡\[\(A−C\)2\]≤C​2−β​lMl,\\mathbb\{E\}\[\(b\-d\)^\{2\}\]\\leq C\\frac\{2^\{\-\\beta l\}\}\{M\_\{l\}\},\\quad\\mathbb\{E\}\[\(A\-C\)^\{2\}\]\\leq C\\frac\{2^\{\-\\beta l\}\}\{M\_\{l\}\},withφ==1\\varphi==1\(the same for𝔼⁡\[\(B−D\)2\]≤C​2−β​lMl\\mathbb\{E\}\[\(B\-D\)^\{2\}\]\\leq C\\frac\{2^\{\-\\beta l\}\}\{M\_\{l\}\}\)\. Since both

fl​\(θk,l\|μk,l\)Fl​\(θk,l,θk,l−1\|μk,l\)≤1,fl−1​\(θk,l−1\|μk,l\)Fl​\(θk,l,θk,l−1\|μk,l\)≤1,\\frac\{f\_\{l\}\(\\theta\_\{k,l\}\|\\mu\_\{k,l\}\)\}\{F\_\{l\}\(\\theta\_\{k,l\},\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)\}\\leq 1,\\quad\\frac\{f\_\{l\-1\}\(\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)\}\{F\_\{l\}\(\\theta\_\{k,l\},\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)\}\\leq 1,by the construction of coupled distributionFl​\(θk,l,θk,l−1\|μk,l\)F\_\{l\}\(\\theta\_\{k,l\},\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\), both

A=1Ml​∑k=1Mlfl​\(θk,l\|μk,l\)Fl​\(θk,l,θk,l−1\|μk,l\),B=1Ml​∑k=1Mlfl−1​\(θk,l−1\|μk,l\)Fl​\(θk,l,θk,l−1\|μk,l\),A=\\frac\{1\}\{M\_\{l\}\}\\sum\_\{k=1\}^\{M\_\{l\}\}\\frac\{f\_\{l\}\(\\theta\_\{k,l\}\|\\mu\_\{k,l\}\)\}\{F\_\{l\}\(\\theta\_\{k,l\},\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)\},\\quad B=\\frac\{1\}\{M\_\{l\}\}\\sum\_\{k=1\}^\{M\_\{l\}\}\\frac\{f\_\{l\-1\}\(\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)\}\{F\_\{l\}\(\\theta\_\{k,l\},\\theta\_\{k,l\-1\}\|\\mu\_\{k,l\}\)\},are bounded\. Noting thatc,C,d,Dc,C,d,Dare not random, and they are all bounded following Assumption[E\.1](https://arxiv.org/html/2610.19074#A5.Thmass1)\. Hence, the last four terms are bounded by𝒪⁡\(2−β​lMl\)\\mathcal\{O\}\(\\frac\{2^\{\-\\beta l\}\}\{M\_\{l\}\}\)\.

∎

### H\.3Proof relating to Lemma[E\.1](https://arxiv.org/html/2610.19074#A5.Thmlem1)

Leta,b,c,d,A,B,C,D∈ℝa,b,c,d,A,B,C,D\\in\\mathbb\{R\}withA,B,C,DA,B,C,Dnon\-zero then;

aA−\\displaystyle\\frac\{a\}\{A\}\-bB−\(cC−dD\)=\[a−b−\(c−d\)\]A−b⁡\[A−B−\(C−D\)\]A​B\+1A​C​\[C−A\]​\[c−d\]\\displaystyle\\frac\{b\}\{B\}\-\\Bigg\(\\frac\{c\}\{C\}\-\\frac\{d\}\{D\}\\Bigg\)=\\frac\{\[a\-b\-\(c\-d\)\]\}\{A\}\-\\frac\{b\[A\-B\-\(C\-D\)\]\}\{AB\}\+\\frac\{1\}\{AC\}\[C\-A\]\[c\-d\]\(H\.8\)−1A​B​\(b−d\)​\(C−D\)\+dC​B​D​\(B−D\)​\(C−D\)\+dA​C​B​\(A−C\)​\(C−D\)\.\\displaystyle\-\\frac\{1\}\{AB\}\(b\-d\)\(C\-D\)\+\\frac\{d\}\{CBD\}\(B\-D\)\(C\-D\)\+\\frac\{d\}\{ACB\}\(A\-C\)\(C\-D\)\.
###### Proof\.

Expand the terms in the right\-hand side of \([H\.8](https://arxiv.org/html/2610.19074#A8.E8)\),

\(aA−bA−cA\+dA\)−\(b​AA​B−b​BA​B−b​CA​B\+b​DA​B\)\+\(C​cA​C−C​dA​C−A​cA​C\+A​dA​C\)\\displaystyle\\Bigg\(\\frac\{a\}\{A\}\-\\frac\{b\}\{A\}\-\\frac\{c\}\{A\}\+\\frac\{d\}\{A\}\\Bigg\)\-\\Bigg\(\\frac\{bA\}\{AB\}\-\\frac\{bB\}\{AB\}\-\\frac\{bC\}\{AB\}\+\\frac\{bD\}\{AB\}\\Bigg\)\+\\Bigg\(\\frac\{Cc\}\{AC\}\-\\frac\{Cd\}\{AC\}\-\\frac\{Ac\}\{AC\}\+\\frac\{Ad\}\{AC\}\\Bigg\)−\(b​CA​B−b​DA​B−d​CA​B\+d​DA​B\+\)\+\(d​B​CC​B​D−d​B​DC​B​D−d​D​CC​B​D\+d​D2C​B​D\)\\displaystyle\-\\Bigg\(\\frac\{bC\}\{AB\}\-\\frac\{bD\}\{AB\}\-\\frac\{dC\}\{AB\}\+\\frac\{dD\}\{AB\}\+\\Bigg\)\+\\Bigg\(\\frac\{dBC\}\{CBD\}\-\\frac\{dBD\}\{CBD\}\-\\frac\{dDC\}\{CBD\}\+\\frac\{dD^\{2\}\}\{CBD\}\\Bigg\)\+\(d​A​CA​C​B−d​A​DA​C​B−d​C2A​C​B\+d​C​DA​C​B\)\\displaystyle\+\\Bigg\(\\frac\{dAC\}\{ACB\}\-\\frac\{dAD\}\{ACB\}\-\\frac\{dC^\{2\}\}\{ACB\}\+\\frac\{dCD\}\{ACB\}\\Bigg\)=\(aA−bA−cA\+dA\)−\(bB−bA−b​CA​B\+b​DA​B\)\+\(cA−dA−cC\+dC\)\\displaystyle=\\Bigg\(\\frac\{a\}\{A\}\-\\frac\{b\}\{A\}\-\\frac\{c\}\{A\}\+\\frac\{d\}\{A\}\\Bigg\)\-\\Bigg\(\\frac\{b\}\{B\}\-\\frac\{b\}\{A\}\-\\frac\{bC\}\{AB\}\+\\frac\{bD\}\{AB\}\\Bigg\)\+\\Bigg\(\\frac\{c\}\{A\}\-\\frac\{d\}\{A\}\-\\frac\{c\}\{C\}\+\\frac\{d\}\{C\}\\Bigg\)−\(b​CA​B−b​DA​B−d​CA​B\+d​DA​B\)\+\(dD−dC−dB\+d​DC​B\)\+\(dB−d​DC​B−d​CA​B\+d​DA​B\)\\displaystyle\-\\Bigg\(\\frac\{bC\}\{AB\}\-\\frac\{bD\}\{AB\}\-\\frac\{dC\}\{AB\}\+\\frac\{dD\}\{AB\}\\Bigg\)\+\\Bigg\(\\frac\{d\}\{D\}\-\\frac\{d\}\{C\}\-\\frac\{d\}\{B\}\+\\frac\{dD\}\{CB\}\\Bigg\)\+\\Bigg\(\\frac\{d\}\{B\}\-\\frac\{dD\}\{CB\}\-\\frac\{dC\}\{AB\}\+\\frac\{dD\}\{AB\}\\Bigg\)=aA−bB−cC\+dD=aA−bB−\(cC−dD\)\.\\displaystyle=\\frac\{a\}\{A\}\-\\frac\{b\}\{B\}\-\\frac\{c\}\{C\}\+\\frac\{d\}\{D\}=\\frac\{a\}\{A\}\-\\frac\{b\}\{B\}\-\\Bigg\(\\frac\{c\}\{C\}\-\\frac\{d\}\{D\}\\Bigg\)\.∎

### H\.4Proof relating to Proposition[E\.2](https://arxiv.org/html/2610.19074#A5.Thmprop2)

Assume Assumptions[E\.1](https://arxiv.org/html/2610.19074#A5.Thmass1),[E\.3](https://arxiv.org/html/2610.19074#A5.Thmass3),[E\.2](https://arxiv.org/html/2610.19074#A5.Thmass2)\. Suppose that\{Wi\}i\\\{W^\{i\}\\\}\_\{i\}is a Markov chain with the Markov kernelKK, andW1∼πW^\{1\}\\sim\\pi\. Then there exists aβ\>0\\beta\>0such that for anyφ∈ℬb​\(Θ×𝖲\)∩Lip​\(Θ×𝖲\)\\varphi\\in\\mathcal\{B\}\_\{b\}\(\\Theta\\times\\mathsf\{S\}\)\\cap\\text\{Lip\}\(\\Theta\\times\\mathsf\{S\}\)there exists aC\>0C\>0such that

𝔼⁡\[\(1Ml​∑k=1Mlφ~l​\(Vk\)−Πl​\(φ~l\)\)2\]≤C​2−β​lMl,\\mathbb\{E\}\\Bigg\[\\Big\(\\frac\{1\}\{M\_\{l\}\}\\sum\_\{k=1\}^\{M\_\{l\}\}\\tilde\{\\varphi\}\_\{l\}\(V^\{k\}\)\-\\Pi\_\{l\}\(\\tilde\{\\varphi\}\_\{l\}\)\\Big\)^\{2\}\\Bigg\]\\leq\\frac\{C2^\{\-\\beta l\}\}\{M\_\{l\}\},whereVk=\(θk,l,μ\)V^\{k\}=\(\\theta\_\{k,l\},\\mu\)is the relevant variables ofWkW^\{k\}\.

###### Proof\.

This proof is similar as Proposition A\.1 in[Jasra et al\. \[2018\]](https://arxiv.org/html/2610.19074#bib.bib9), where we repeat a similar procedure\.

𝔼\[\(1Ml∑k=1Mlφ~l\(Vk\)−Πl\(φ~l\]\)\)2\]=𝔼\[\(1Ml∑k=1Mlf\(Wk\)\)2\],\\displaystyle\\mathbb\{E\}\\bigg\[\\bigg\(\\frac\{1\}\{M\_\{l\}\}\\sum\_\{k=1\}^\{M\_\{l\}\}\\tilde\{\\varphi\}\_\{l\}\(V^\{k\}\)\-\\Pi\_\{l\}\(\\tilde\{\\varphi\}\_\{l\}\]\)\\bigg\)^\{2\}\\bigg\]=\\mathbb\{E\}\\bigg\[\\bigg\(\\frac\{1\}\{M\_\{l\}\}\\sum\_\{k=1\}^\{M\_\{l\}\}f\(W^\{k\}\)\\bigg\)^\{2\}\\bigg\],wheref⁡\(w\)=φ~l∘ψ⁡\(w\)−πl​\(φ~l∘ψ\)=φ~l​\(v\)−Πl​\(φ~l\)f\(w\)=\\tilde\{\\varphi\}\_\{l\}\\circ\\psi\(w\)\-\\pi\_\{l\}\(\\tilde\{\\varphi\}\_\{l\}\\circ\\psi\)=\\tilde\{\\varphi\}\_\{l\}\(v\)\-\\Pi\_\{l\}\(\\tilde\{\\varphi\}\_\{l\}\)\. Then, we have

𝔼⁡\[\(1Ml​∑k=1Mlf⁡\(Wk\)\)2\]\\displaystyle\\mathbb\{E\}\\bigg\[\\bigg\(\\frac\{1\}\{M\_\{l\}\}\\sum\_\{k=1\}^\{M\_\{l\}\}f\(W^\{k\}\)\\bigg\)^\{2\}\\bigg\]=1Ml2​∑p,q=1Ml⟨f,K\|p−q\|​f⟩\\displaystyle=\\frac\{1\}\{M\_\{l\}^\{2\}\}\\sum\_\{p,q=1\}^\{M\_\{l\}\}\\langle f,K^\{\|p\-q\|\}f\\rangle=1Ml​‖f‖22\+2Ml​∑p=1Ml−1\(Ml−p\)​⟨f,Kp​f⟩\\displaystyle=\\frac\{1\}\{M\_\{l\}\}\\left\\lVert f\\right\\rVert\_\{2\}^\{2\}\+\\frac\{2\}\{M\_\{l\}\}\\sum\_\{p=1\}^\{M\_\{l\}\-1\}\(M\_\{l\}\-p\)\\langle f,K^\{p\}f\\rangle≤1Ml​‖f‖22\+2Ml2​∑p=1Ml−1\(Ml−p\)​‖Kp‖2​‖f‖22\\displaystyle\\leq\\frac\{1\}\{M\_\{l\}\}\\left\\lVert f\\right\\rVert\_\{2\}^\{2\}\+\\frac\{2\}\{M\_\{l\}^\{2\}\}\\sum\_\{p=1\}^\{M\_\{l\}\-1\}\(M\_\{l\}\-p\)\\left\\lVert K^\{p\}\\right\\rVert\_\{2\}\\left\\lVert f\\right\\rVert\_\{2\}^\{2\}≤1Ml​‖f‖22\+2Ml​∑p=1∞‖Kp‖2​‖f‖22\.\\displaystyle\\leq\\frac\{1\}\{M\_\{l\}\}\\left\\lVert f\\right\\rVert\_\{2\}^\{2\}\+\\frac\{2\}\{M\_\{l\}\}\\sum\_\{p=1\}^\{\\infty\}\\left\\lVert K^\{p\}\\right\\rVert\_\{2\}\\left\\lVert f\\right\\rVert\_\{2\}^\{2\}\.\(H\.9\)By Assumption[E\.3](https://arxiv.org/html/2610.19074#A5.Thmass3), we have

\([H\.9](https://arxiv.org/html/2610.19074#A8.E9)\)≤1Ml​‖f‖22​\(1\+1−ξξ\)≤2​‖f‖22Ml​ξ\.\\displaystyle\\eqref\{eqn:f\_sum\}\\leq\\frac\{1\}\{M\_\{l\}\}\\left\\lVert f\\right\\rVert\_\{2\}^\{2\}\\big\(1\+\\frac\{1\-\\xi\}\{\\xi\}\\big\)\\leq\\frac\{2\\left\\lVert f\\right\\rVert\_\{2\}^\{2\}\}\{M\_\{l\}\\xi\}\.\(H\.10\)Since‖f‖22=𝔼⁡\[\(φ~l​\(Vk\)−Π⁡\(φ~l\)\)2\]\\left\\lVert f\\right\\rVert\_\{2\}^\{2\}=\\mathbb\{E\}\[\(\\tilde\{\\varphi\}\_\{l\}\(V^\{k\}\)\-\\Pi\(\\tilde\{\\varphi\}\_\{l\}\)\)^\{2\}\]and Assumption[E\.2](https://arxiv.org/html/2610.19074#A5.Thmass2)\(V\), we have

‖f‖22=𝔼⁡\[\(φ~l​\(Vk\)−Π⁡\(φ~l\)\)2\]≤Vl≤C​2−β​l\.\\displaystyle\\left\\lVert f\\right\\rVert\_\{2\}^\{2\}=\\mathbb\{E\}\[\(\\tilde\{\\varphi\}\_\{l\}\(V^\{k\}\)\-\\Pi\(\\tilde\{\\varphi\}\_\{l\}\)\)^\{2\}\]\\leq V\_\{l\}\\leq C2^\{\-\\beta l\}\.∎

### H\.5Verification of Assumptions

We now want to verify Assumption[E\.3](https://arxiv.org/html/2610.19074#A5.Thmass3)and Assumption[E\.2](https://arxiv.org/html/2610.19074#A5.Thmass2)\. Let us begin with Assumption[E\.3](https://arxiv.org/html/2610.19074#A5.Thmass3)\.

There exists aξ∈\(0,1\)\\xi\\in\(0,1\)andν∈𝒫⁡\(𝖶\)\\nu\\in\\mathcal\{P\}\(\\mathsf\{W\}\)such that, for eachw∈𝖶w\\in\\mathsf\{W\},φ∈ℬb​\(𝖶\)∩Lip​\(𝖶\)\\varphi\\in\\mathcal\{B\}\_\{b\}\(\\mathsf\{W\}\)\\cap\\text\{Lip\}\(\\mathsf\{W\}\),

∫𝖶φ⁡\(w′\)​K​\(w,d​w′\)≥ξ​∫Wφ⁡\(v\)​ν​\(𝑑v\)\.\\int\_\{\\mathsf\{W\}\}\\varphi\(w^\{\\prime\}\)K\(w,dw^\{\\prime\}\)\\geq\\xi\\int\_\{W\}\\varphi\(v\)\\nu\(dv\)\.
###### Proof\.

The proposal kernelψ\\psiof PMMH is

ψ\(w,dw∗\)=P\(θ1:M,l\|μq,l∗\)Q\(dμq,l∗\|μk,l\)Fq,l​\(θq,l,θq,l−1\|μq,l∗\)∑i=1MFi,l​\(θi,l,θi,l−1\|μq,l∗\),\\psi\(w,dw^\{\*\}\)=P\(\\theta\_\{1:M,l\}\|\\mu^\{\*\}\_\{q,l\}\)Q\(d\\mu^\{\*\}\_\{q,l\}\|\\mu\_\{k,l\}\)\\frac\{F\_\{q,l\}\(\\theta\_\{q,l\},\\theta\_\{q,l\-1\}\|\\mu^\{\*\}\_\{q,l\}\)\}\{\\sum\_\{i=1\}^\{M\}F\_\{i,l\}\(\\theta\_\{i,l\},\\theta\_\{i,l\-1\}\|\\mu^\{\*\}\_\{q,l\}\)\},\(H\.11\)wherew=\(θk,l,μk,l\)w=\(\\theta\_\{k,l\},\\mu\_\{k,l\}\),w∗=\(θq,l∗,μq,l∗\)w^\{\*\}=\(\\theta^\{\*\}\_\{q,l\},\\mu^\{\*\}\_\{q,l\}\)and

P\(θ1:M,l\|μ\)=∏k=1Mπ¯\(θk,l\)Fk,l​\(θk,l,θk,l−1\|μ\)∑i=1MFi,l​\(θi,l,θi,l−1\|μ\)\.P\(\\theta\_\{1:M,l\}\|\\mu\)=\\prod\_\{k=1\}^\{M\}\\bar\{\\pi\}\(\\theta\_\{k,l\}\)\\frac\{F\_\{k,l\}\(\\theta\_\{k,l\},\\theta\_\{k,l\-1\}\|\\mu\)\}\{\\sum\_\{i=1\}^\{M\}F\_\{i,l\}\(\\theta\_\{i,l\},\\theta\_\{i,l\-1\}\|\\mu\)\}\.The transition kernelKKis

K⁡\(w,d​w∗\)=ψ⁡\(w,d​w∗\)​α​\(w,w∗\)\+δw​\(d​w∗\)​Q​\(w\)\\displaystyle K\(w,dw^\{\*\}\)=\\psi\(w,dw^\{\*\}\)\\alpha\(w,w^\{\*\}\)\+\\delta\_\{w\}\(dw^\{\*\}\)Q\(w\)where the acceptance probabilityα⁡\(w,w∗\)\\alpha\(w,w^\{\*\}\)is

α⁡\(w,w∗\)=min⁡\{1,π0​\(μq,l∗\)​Q​\(μk,l\|μq,l∗\)​pk,lN​\(μk,l∗\)π0​\(μq,l∗\)​Q​\(μq,l∗\|μk,l\)​pk,lN​\(μk,l\)\},\\displaystyle\\alpha\(w,w^\{\*\}\)=\\min\\Bigg\\\{1,\\frac\{\\pi\_\{0\}\(\\mu^\{\*\}\_\{q,l\}\)Q\(\\mu\_\{k,l\}\|\\mu^\{\*\}\_\{q,l\}\)p\_\{k,l\}^\{N\}\(\\mu\_\{k,l\}^\{\*\}\)\}\{\\pi\_\{0\}\(\\mu^\{\*\}\_\{q,l\}\)Q\(\\mu^\{\*\}\_\{q,l\}\|\\mu\_\{k,l\}\)p\_\{k,l\}^\{N\}\(\\mu\_\{k,l\}\)\}\\Bigg\\\},and the rejection probabilityQ⁡\(w\)Q\(w\)is

Q⁡\(w\)=1−∫w∗ψ⁡\(w,d​w∗\)​α​\(w,w∗\)\.\\displaystyle Q\(w\)=1\-\\int\_\{w^\{\*\}\}\\psi\(w,dw^\{\*\}\)\\alpha\(w,w^\{\*\}\)\.
Assumeπ:𝖲→ℝ\+\\pi:\\mathsf\{S\}\\rightarrow\\mathbb\{R\}\_\{\+\}andQ:𝖲2→ℝ\+Q:\\mathsf\{S\}^\{2\}\\rightarrow\\mathbb\{R\}\_\{\+\}are continuous and strictly positive and Assumption[E\.1](https://arxiv.org/html/2610.19074#A5.Thmass1)are satisfied, we have

C1:=infw,w∗α⁡\(w,w∗\)\>0,\\displaystyle C\_\{1\}:=\\inf\_\{w,w^\{\*\}\}\\alpha\(w,w^\{\*\}\)\>0,infwψ⁡\(w,d​w∗\)\\displaystyle\\inf\_\{w\}\\psi\(w,dw^\{\*\}\)≥\{minμk,l,μq,l∗Q\(μq,l∗\|μk,l\)\}1𝖲\(μq,l∗\)dμq,l∗P\(θ1:M,l\|μq,l∗\)Fq,l​\(θq,l,θq,l−1\|μq,l∗\)∑i=1MFi,l​\(θi,l,θi,l−1\|μq,l∗\)\\displaystyle\\geq\\big\\\{\\min\_\{\\mu\_\{k,l\},\\mu^\{\*\}\_\{q,l\}\}Q\(\\mu^\{\*\}\_\{q,l\}\|\\mu\_\{k,l\}\)\\big\\\}1\_\{\\mathsf\{S\}\}\(\\mu^\{\*\}\_\{q,l\}\)d\\mu^\{\*\}\_\{q,l\}P\(\\theta\_\{1:M,l\}\|\\mu^\{\*\}\_\{q,l\}\)\\frac\{F\_\{q,l\}\(\\theta\_\{q,l\},\\theta\_\{q,l\-1\}\|\\mu^\{\*\}\_\{q,l\}\)\}\{\\sum\_\{i=1\}^\{M\}F\_\{i,l\}\(\\theta\_\{i,l\},\\theta\_\{i,l\-1\}\|\\mu^\{\*\}\_\{q,l\}\)\}=:C2​ψ​\(d​w∗\),\\displaystyle=:C\_\{2\}\\psi\(dw^\{\*\}\),whereC2=Leb\(𝖲\)minμk,l,μq,l∗Q\(μq,l∗\|μk,l\]\)\>0C\_\{2\}=\\text\{Leb\}\(\\mathsf\{S\}\)\\min\_\{\\mu\_\{k,l\},\\mu^\{\*\}\_\{q,l\}\}Q\(\\mu^\{\*\}\_\{q,l\}\|\\mu\_\{k,l\}\]\)\>0with a probability measureψ⁡\(d​w∗\)\\psi\(dw^\{\*\}\)\. Thus, we have

K⁡\(w,d​w∗\)≥C1​C2​ψ​\(d​w∗\)\.\\displaystyle K\(w,dw^\{\*\}\)\\geq C\_\{1\}C\_\{2\}\\psi\(dw^\{\*\}\)\.Hence the Assumption[E\.3](https://arxiv.org/html/2610.19074#A5.Thmass3)is satisfied\. ∎

Next, we verify the bounds for the \(B\) and \(V\) in Assumption[E\.2](https://arxiv.org/html/2610.19074#A5.Thmass2)\.

For anyφ∈ℬb​\(Θl×Θl−1×S\)∩Lip​\(Θl×Θl−1×𝖲\)\\varphi\\in\\mathcal\{B\}\_\{b\}\(\\Theta\_\{l\}\\times\\Theta\_\{l\-1\}\\times S\)\\cap\\text\{Lip\}\(\\Theta\_\{l\}\\times\\Theta\_\{l\-1\}\\times\\mathsf\{S\}\)there existc1,c2,c3,β,α,γ\>0c\_\{1\},c\_\{2\},c\_\{3\},\\beta,\\alpha,\\gamma\>0such that the following holds, under the joint posteriorΠl\\Pi\_\{l\},

- \(B\)𝔼⁡\[\|φl​\(θl,μ\)​fl​\(θl\|μ\)Fl​\(θl,θl−1\|μ\)−f⁡\(φ⁡\(θ,μ\)\)\|\]=:Bl≤c1​2−α​l\\mathbb\{E\}\\Big\[\\big\|\\varphi\_\{l\}\(\\theta\_\{l\},\\mu\)\\frac\{f\_\{l\}\(\\theta\_\{l\}\|\\mu\)\}\{F\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\}\-f\(\\varphi\(\\theta,\\mu\)\)\\big\|\\Big\]=:B\_\{l\}\\leq c\_\{1\}2^\{\-\\alpha l\};
- \(V\)𝔼⁡\[\|φl​\(θl,μ\)​fl​\(θl\|μ\)Fl​\(θl,θl−1\|μ\)−φl−1​\(θl−1,μ\)​fl−1​\(θl−1\|μ\)Fl​\(θl,θl−1\|μ\)\|2\]=:Vl≤c2​2−β​l\\mathbb\{E\}\\Big\[\\big\|\\varphi\_\{l\}\(\\theta\_\{l\},\\mu\)\\frac\{f\_\{l\}\(\\theta\_\{l\}\|\\mu\)\}\{F\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\}\-\\varphi\_\{l\-1\}\(\\theta\_\{l\-1\},\\mu\)\\frac\{f\_\{l\-1\}\(\\theta\_\{l\-1\}\|\\mu\)\}\{F\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\}\\big\|^\{2\}\\Big\]=:V\_\{l\}\\leq c\_\{2\}2^\{\-\\beta l\};
- \(C\)COST⁡\(Πl​\(φl\),Πl−1​\(φl−1\)\)=:Cl≤c3​2γ​l\{\\rm COST\}\(\\Pi\_\{l\}\(\\varphi\_\{l\}\),\\Pi\_\{l\-1\}\(\\varphi\_\{l\-1\}\)\)=:C\_\{l\}\\leq c\_\{3\}2^\{\\gamma l\}\.

###### Proof\.

Bl\\displaystyle B\_\{l\}=∫Θ×𝖲\|fl​\(θl\|μ\)Fl​\(θl,θl−1\|μ\)​φl​\(θl,μ\)−f⁡\(θ\|μ\)​φ​\(θ,μ\)\|​π​\(𝑑θ,𝑑μ\)\\displaystyle=\\int\_\{\\Theta\\times\\mathsf\{S\}\}\\Bigg\|\\frac\{f\_\{l\}\(\\theta\_\{l\}\|\\mu\)\}\{F\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\}\\varphi\_\{l\}\(\\theta\_\{l\},\\mu\)\-f\(\\theta\|\\mu\)\\varphi\(\\theta,\\mu\)\\Bigg\|\\pi\(d\\theta,d\\mu\)=∫Θ×𝖲\|fl​\(θl\|μ\)Fl​\(θl,θl−1\|μ\)​φl​\(θl,μ\)−f⁡\(θ\|μ\)​φ​\(θ,μ\)\|​ℒ​\(θ\|μ\)​π¯​\(𝑑θ\)​π0​\(𝑑μ\)×∫Θ×𝖲ℒ⁡\(θ\|μ\)​π¯​\(𝑑θ\)​π0​\(𝑑μ\)\\displaystyle=\\int\_\{\\Theta\\times\\mathsf\{S\}\}\\Bigg\|\\frac\{f\_\{l\}\(\\theta\_\{l\}\|\\mu\)\}\{F\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\}\\varphi\_\{l\}\(\\theta\_\{l\},\\mu\)\-f\(\\theta\|\\mu\)\\varphi\(\\theta,\\mu\)\\Bigg\|\\mathcal\{L\}\(\\theta\|\\mu\)\\bar\{\\pi\}\(d\\theta\)\\pi\_\{0\}\(d\\mu\)\\times\\int\_\{\\Theta\\times\\mathsf\{S\}\}\\mathcal\{L\}\(\\theta\|\\mu\)\\bar\{\\pi\}\(d\\theta\)\\pi\_\{0\}\(d\\mu\)=∫Θ×𝖲\|fl​\(θl\|μ\)Fl​\(θl,θl−1\|μ\)φl\(θl,μ\)−fl​\(θl\|μ\)Fl​\(θl,θl−1\|μ\)φ\(θ,μ\)\+fl​\(θl\|μ\)Fl​\(θl,θl−1\|μ\)φ\(θ,μ\)\\displaystyle=\\int\_\{\\Theta\\times\\mathsf\{S\}\}\\Bigg\|\\frac\{f\_\{l\}\(\\theta\_\{l\}\|\\mu\)\}\{F\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\}\\varphi\_\{l\}\(\\theta\_\{l\},\\mu\)\-\\frac\{f\_\{l\}\(\\theta\_\{l\}\|\\mu\)\}\{F\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\}\\varphi\(\\theta,\\mu\)\+\\frac\{f\_\{l\}\(\\theta\_\{l\}\|\\mu\)\}\{F\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\}\\varphi\(\\theta,\\mu\)−f⁡\(θ\|μ\)​φ​\(θ,μ\)\|ℒ⁡\(θ\|μ\)​π¯​\(𝑑θ\)​π0​\(𝑑μ\)×∫Θ×𝖲ℒ⁡\(θ\|μ\)​π¯​\(𝑑θ\)​π0​\(𝑑μ\)\.\\displaystyle\-f\(\\theta\|\\mu\)\\varphi\(\\theta,\\mu\)\\Bigg\|\\mathcal\{L\}\(\\theta\|\\mu\)\\bar\{\\pi\}\(d\\theta\)\\pi\_\{0\}\(d\\mu\)\\times\\int\_\{\\Theta\\times\\mathsf\{S\}\}\\mathcal\{L\}\(\\theta\|\\mu\)\\bar\{\\pi\}\(d\\theta\)\\pi\_\{0\}\(d\\mu\)\.we defineC′=maxθ∈Θ⁡ℒ⁡\(θ\|μ\)C^\{\\prime\}=\\max\_\{\\theta\\in\\Theta\}\\mathcal\{L\}\(\\theta\|\\mu\)andC′=∫Θ×𝖲ℒ⁡\(θ\|μ\)​π¯​\(𝑑θ\)​π0​\(𝑑μ\)C^\{\\prime\}=\\int\_\{\\Theta\\times\\mathsf\{S\}\}\\mathcal\{L\}\(\\theta\|\\mu\)\\bar\{\\pi\}\(d\\theta\)\\pi\_\{0\}\(d\\mu\)which results

Bl\\displaystyle B\_\{l\}≤C′​C′′​∫Θ×𝖲\|fl​\(θl\|μ\)Fl​\(θl,θl−1\|μ\)​φl​\(θl,μ\)−fl​\(θl\|μ\)Fl​\(θl,θl−1\|μ\)​φ​\(θ,μ\)\|​π¯​\(𝑑θ\)​π0​\(𝑑μ\)\\displaystyle\\leq C^\{\\prime\}C^\{\\prime\\prime\}\\int\_\{\\Theta\\times\\mathsf\{S\}\}\\Bigg\|\\frac\{f\_\{l\}\(\\theta\_\{l\}\|\\mu\)\}\{F\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\}\\varphi\_\{l\}\(\\theta\_\{l\},\\mu\)\-\\frac\{f\_\{l\}\(\\theta\_\{l\}\|\\mu\)\}\{F\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\}\\varphi\(\\theta,\\mu\)\\Bigg\|\\bar\{\\pi\}\(d\\theta\)\\pi\_\{0\}\(d\\mu\)\+C′∫Θ×𝖲\|fl​\(θl\|μ\)Fl​\(θl,θl−1\|μ\)φ\(θ,μ\)−f\(θ\|μ\)φ\(θ,μ\)\|ℒ\(θ\|μ\)π¯\(dθ\)π0\(dμ\)\\displaystyle\+C^\{\\prime\}\\int\_\{\\Theta\\times\\mathsf\{S\}\}\\Bigg\|\\frac\{f\_\{l\}\(\\theta\_\{l\}\|\\mu\)\}\{F\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\}\\varphi\(\\theta,\\mu\)\-f\(\\theta\|\\mu\)\\varphi\(\\theta,\\mu\)\\Bigg\|\\mathcal\{L\}\(\\theta\|\\mu\)\\bar\{\\pi\}\(d\\theta\)\\pi\_\{0\}\(d\\mu\)≤C′​C′′​∫Θ×𝖲\|fl​\(θl\|μ\)Fl​\(θl,θl−1\|μ\)\|​\|φl​\(θl,μ\)−φ⁡\(θ,μ\)\|​π¯​\(𝑑θ\)​π0​\(𝑑μ\)\\displaystyle\\leq C^\{\\prime\}C^\{\\prime\\prime\}\\int\_\{\\Theta\\times\\mathsf\{S\}\}\\Bigg\|\\frac\{f\_\{l\}\(\\theta\_\{l\}\|\\mu\)\}\{F\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\}\\Bigg\|\|\\varphi\_\{l\}\(\\theta\_\{l\},\\mu\)\-\\varphi\(\\theta,\\mu\)\|\\bar\{\\pi\}\(d\\theta\)\\pi\_\{0\}\(d\\mu\)\+C′∫Θ×𝖲\|fl​\(θl\|μ\)Fl​\(θl,θl−1\|μ\)−f\(θ\|μ\)\|\|φ\(θ,μ\)\|ℒ\(θ\|μ\)π¯\(dθ\)π0\(dμ\)\.\\displaystyle\+C^\{\\prime\}\\int\_\{\\Theta\\times\\mathsf\{S\}\}\\Bigg\|\\frac\{f\_\{l\}\(\\theta\_\{l\}\|\\mu\)\}\{F\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\}\-f\(\\theta\|\\mu\)\\Bigg\|\|\\varphi\(\\theta,\\mu\)\|\\mathcal\{L\}\(\\theta\|\\mu\)\\bar\{\\pi\}\(d\\theta\)\\pi\_\{0\}\(d\\mu\)\.Let us defineC′′′=maxθ∈Θ,μ∈𝖲⁡\|φ⁡\(θ,μ\)\|C^\{\\prime\\prime\\prime\}=\\max\_\{\\theta\\in\\Theta,\\mu\\in\\mathsf\{S\}\}\|\\varphi\(\\theta,\\mu\)\|andC~\\tilde\{C\}to be a collection of various constants\. Furthermore since

\|fl​\(θl\|μ\)Fl​\(θl,θl−1\|μ\)\|≤1,π⁡\(fl​\(θl\|μ\)Fl​\(θl,θl−1\|μ\)\)=π⁡\(fl​\(θ\|μ\)\),\\Bigg\|\\frac\{f\_\{l\}\(\\theta\_\{l\}\|\\mu\)\}\{F\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\}\\Bigg\|\\leq 1,\\quad\\pi\\bigg\(\\frac\{f\_\{l\}\(\\theta\_\{l\}\|\\mu\)\}\{F\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\}\\bigg\)=\\pi\(f\_\{l\}\(\\theta\|\\mu\)\),we have

Bl\\displaystyle B\_\{l\}≤C′​C′′​∫Θ×𝖲\|φl​\(θl,μ\)−φ⁡\(θ,μ\)\|​π¯​\(𝑑θ\)​π0​\(𝑑μ\)\\displaystyle\\leq C^\{\\prime\}C^\{\\prime\\prime\}\\int\_\{\\Theta\\times\\mathsf\{S\}\}\|\\varphi\_\{l\}\(\\theta\_\{l\},\\mu\)\-\\varphi\(\\theta,\\mu\)\|\\bar\{\\pi\}\(d\\theta\)\\pi\_\{0\}\(d\\mu\)\+C′C′′′∫Θ×𝖲\|fl​\(θl\|μ\)Fl​\(θl,θl−1\|μ\)−f\(θ\|μ\)\|ℒ\(θ\|μ\)π¯\(dθ\)π0\(dμ\)\\displaystyle\+C^\{\\prime\}C^\{\\prime\\prime\\prime\}\\int\_\{\\Theta\\times\\mathsf\{S\}\}\\Bigg\|\\frac\{f\_\{l\}\(\\theta\_\{l\}\|\\mu\)\}\{F\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\}\-f\(\\theta\|\\mu\)\\Bigg\|\\mathcal\{L\}\(\\theta\|\\mu\)\\bar\{\\pi\}\(d\\theta\)\\pi\_\{0\}\(d\\mu\)≤C~​\(∫Θ×𝖲\|φl​\(θl,μ\)−φ⁡\(θ,μ\)\|​π¯​\(𝑑θ\)​π0​\(𝑑μ\)\+\|fl​\(θl\|μ\)−f⁡\(θ\|μ\)\|​ℒ​\(θ\|μ\)​π¯​\(𝑑θ\)​π0​\(𝑑μ\)\)\\displaystyle\\leq\\tilde\{C\}\\Bigg\(\\int\_\{\\Theta\\times\\mathsf\{S\}\}\|\\varphi\_\{l\}\(\\theta\_\{l\},\\mu\)\-\\varphi\(\\theta,\\mu\)\|\\bar\{\\pi\}\(d\\theta\)\\pi\_\{0\}\(d\\mu\)\+\|f\_\{l\}\(\\theta\_\{l\}\|\\mu\)\-f\(\\theta\|\\mu\)\|\\mathcal\{L\}\(\\theta\|\\mu\)\\bar\{\\pi\}\(d\\theta\)\\pi\_\{0\}\(d\\mu\)\\Bigg\)≤C~​\(∫Θ×𝖲\|φl​\(θl,μ\)−φ⁡\(θ,μ\)\|​π¯​\(𝑑θ\)​π0​\(𝑑μ\)\+\|φl​\(θl\|μ\)−φ⁡\(θ\|μ\)\|​π¯​\(𝑑θ\)​π0​\(𝑑μ\)\)\\displaystyle\\leq\\tilde\{C\}\\Bigg\(\\int\_\{\\Theta\\times\\mathsf\{S\}\}\|\\varphi\_\{l\}\(\\theta\_\{l\},\\mu\)\-\\varphi\(\\theta,\\mu\)\|\\bar\{\\pi\}\(d\\theta\)\\pi\_\{0\}\(d\\mu\)\+\|\\varphi\_\{l\}\(\\theta\_\{l\}\|\\mu\)\-\\varphi\(\\theta\|\\mu\)\|\\bar\{\\pi\}\(d\\theta\)\\pi\_\{0\}\(d\\mu\)\\Bigg\)Applying Proposition[E\.1](https://arxiv.org/html/2610.19074#A5.Thmprop1)on it and lettingα≤\(2​s−1\)/2\\alpha\\leq\(2s\-1\)/2≤C​2−α​l,\\displaystyle\\leq C2^\{\-\\alpha l\},where we have assumed thatfl​\(θ\|μ\)f\_\{l\}\(\\theta\|\\mu\)is Lipschitz w\.r\.tφl​\(θ,μ\)\\varphi\_\{l\}\(\\theta,\\mu\)with Lipschitz constantGG,f∞=ff\_\{\\infty\}=fandC′=maxθ∈Θ⁡ℒ⁡\(θ\|μ\)C^\{\\prime\}=\\max\_\{\\theta\\in\\Theta\}\\mathcal\{L\}\(\\theta\|\\mu\)\. Now, for the variance term, which we denote asVlV\_\{l\}, we proceed as follows:

Vl\\displaystyle V\_\{l\}=∫Θ2×𝖲\(fl​\(θl\|μ\)Fl​\(θl,θl−1\|μ\)φl\(θl,μ\)−fl−1​\(θl−1\|μ\)Fl​\(θl,θl−1\|μ\)φl−1\(θl−1,μ\)\)2×\\displaystyle=\\int\_\{\\Theta^\{2\}\\times\\mathsf\{S\}\}\\Bigg\(\\frac\{f\_\{l\}\(\\theta\_\{l\}\|\\mu\)\}\{F\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\}\\varphi\_\{l\}\(\\theta\_\{l\},\\mu\)\-\\frac\{f\_\{l\-1\}\(\\theta\_\{l\-1\}\|\\mu\)\}\{F\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\}\\varphi\_\{l\-1\}\(\\theta\_\{l\-1\},\\mu\)\\Bigg\)^\{2\}\\times𝐋l​\(θl,θl−1\|μ\)​Π¯l​\(d⁡\(θl,θl−1\)\)​π0​\(𝑑μ\)​\(∫Θ2×𝖲𝐋l​\(θl,θl−1\|μ\)​Π¯l​\(d⁡\(θl,θl−1\)\)​π0​\(𝑑μ\)\)−1\\displaystyle\{\\bf L\}\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\\bar\{\\Pi\}\_\{l\}\(d\(\\theta\_\{l\},\\theta\_\{l\-1\}\)\)\\pi\_\{0\}\(d\\mu\)\\Bigg\(\\int\_\{\\Theta^\{2\}\\times\\mathsf\{S\}\}\{\\bf L\}\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\\bar\{\\Pi\}\_\{l\}\(d\(\\theta\_\{l\},\\theta\_\{l\-1\}\)\)\\pi\_\{0\}\(d\\mu\)\\bigg\)^\{\-1\}≤C′​∫Θ2\(fl​\(θl\|μ\)​φl​\(θl,μ\)−fl−1​\(θl−1\|μ\)​φl−1​\(θl−1,μ\)\)2​Π¯l​\(d⁡\(θl,θl−1\)\)​π0​\(𝑑μ\)\\displaystyle\\leq C^\{\\prime\}\\int\_\{\\Theta^\{2\}\}\\Bigg\(f\_\{l\}\(\\theta\_\{l\}\|\\mu\)\\varphi\_\{l\}\(\\theta\_\{l\},\\mu\)\-f\_\{l\-1\}\(\\theta\_\{l\-1\}\|\\mu\)\\varphi\_\{l\-1\}\(\\theta\_\{l\-1\},\\mu\)\\Bigg\)^\{2\}\\bar\{\\Pi\}\_\{l\}\(d\(\\theta\_\{l\},\\theta\_\{l\-1\}\)\)\\pi\_\{0\}\(d\\mu\)where we define

C′=infθ∈Θ2𝐋l​\(θl,θl−1\|μ\)Fl​\(θl,θl−1\|μ\)2×\(∫Θ2×𝖲𝐋l​\(θl,θl−1\|μ\)​Π¯l​\(d⁡\(θl,θl−1\)\)​π0​\(𝑑μ\)\)−1,C^\{\\prime\}=\\inf\_\{\\theta\\in\\Theta^\{2\}\}\\frac\{\{\\bf L\}\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\}\{F\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)^\{2\}\}\\times\\Bigg\(\\int\_\{\\Theta^\{2\}\\times\\mathsf\{S\}\}\{\\bf L\}\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\\bar\{\\Pi\}\_\{l\}\(d\(\\theta\_\{l\},\\theta\_\{l\-1\}\)\)\\pi\_\{0\}\(d\\mu\)\\bigg\)^\{\-1\},Vl\\displaystyle V\_\{l\}=C′​∫Θ×𝖲\(fl​\(θ\|μ\)​φl​\(θ,μ\)−fl−1​\(θ\|μ\)​φl−1​\(θ,μ\)\)2​π¯l​\(𝑑θ\)​π0​\(𝑑μ\)\\displaystyle=C^\{\\prime\}\\int\_\{\\Theta\\times\\mathsf\{S\}\}\(f\_\{l\}\(\\theta\|\\mu\)\\varphi\_\{l\}\(\\theta,\\mu\)\-f\_\{l\-1\}\(\\theta\|\\mu\)\\varphi\_\{l\-1\}\(\\theta,\\mu\)\)^\{2\}\\bar\{\\pi\}\_\{l\}\(d\\theta\)\\pi\_\{0\}\(d\\mu\)=C′​∫Θ×𝖲\(fl​\(θ\|μ\)​φl​\(θ,μ\)CLOSE\\displaystyle=C^\{\\prime\}\\int\_\{\\Theta\\times\\mathsf\{S\}\}\\Big\(f\_\{l\}\(\\theta\|\\mu\)\\varphi\_\{l\}\(\\theta,\\mu\)OPEN−fl​\(θ\|μ\)​φl−1​\(θ,μ\)\+fl​\(θ\|μ\)​φl−1​\(θ,μ\)−fl−1​\(θ\|μ\)​φl−1​\(θ,μ\)\)2​π¯l​\(d​θ\)​π0​\(d​μ\)\\displaystyle\-f\_\{l\}\(\\theta\|\\mu\)\\varphi\_\{l\-1\}\(\\theta,\\mu\)\+f\_\{l\}\(\\theta\|\\mu\)\\varphi\_\{l\-1\}\(\\theta,\\mu\)\-f\_\{l\-1\}\(\\theta\|\\mu\)\\varphi\_\{l\-1\}\(\\theta,\\mu\)\\Big\)^\{2\}\\bar\{\\pi\}\_\{l\}\(d\\theta\)\\pi\_\{0\}\(d\\mu\)≤C′​∫Θ×𝖲\(fl​\(θ\|μ\)​φl​\(θ,μ\)−fl​\(θ\|μ\)​φl−1​\(θ,μ\)\)2​π¯l​\(𝑑θ\)​π0​\(𝑑μ\)\\displaystyle\\leq C^\{\\prime\}\\int\_\{\\Theta\\times\\mathsf\{S\}\}\(f\_\{l\}\(\\theta\|\\mu\)\\varphi\_\{l\}\(\\theta,\\mu\)\-f\_\{l\}\(\\theta\|\\mu\)\\varphi\_\{l\-1\}\(\\theta,\\mu\)\)^\{2\}\\bar\{\\pi\}\_\{l\}\(d\\theta\)\\pi\_\{0\}\(d\\mu\)\+C′∫Θ×𝖲\(fl\(θ\|μ\)φl−1\(θ,μ\)−fl−1\(θ\|μ\)φl−1\(θ,μ\)\)2π¯l\(dθ\)π0\(dμ\)\\displaystyle\+C^\{\\prime\}\\int\_\{\\Theta\\times\\mathsf\{S\}\}\(f\_\{l\}\(\\theta\|\\mu\)\\varphi\_\{l\-1\}\(\\theta,\\mu\)\-f\_\{l\-1\}\(\\theta\|\\mu\)\\varphi\_\{l\-1\}\(\\theta,\\mu\)\)^\{2\}\\bar\{\\pi\}\_\{l\}\(d\\theta\)\\pi\_\{0\}\(d\\mu\)≤C′​∫Θ×𝖲\(fl​\(θ\|μ\)\)2​\(φl​\(θ,μ\)−φl−1​\(θ,μ\)\)2​π¯l​\(𝑑θ\)​π0​\(𝑑μ\)\\displaystyle\\leq C^\{\\prime\}\\int\_\{\\Theta\\times\\mathsf\{S\}\}\(f\_\{l\}\(\\theta\|\\mu\)\)^\{2\}\(\\varphi\_\{l\}\(\\theta,\\mu\)\-\\varphi\_\{l\-1\}\(\\theta,\\mu\)\)^\{2\}\\bar\{\\pi\}\_\{l\}\(d\\theta\)\\pi\_\{0\}\(d\\mu\)\+C′∫Θ×𝖲\(φl−1\(θ,μ\)\)2\(fl\(θ\|μ\)−fl−1\(θ\|μ\)\)2π¯l\(dθ\)π0\(dμ\)\\displaystyle\+C^\{\\prime\}\\int\_\{\\Theta\\times\\mathsf\{S\}\}\(\\varphi\_\{l\-1\}\(\\theta,\\mu\)\)^\{2\}\(f\_\{l\}\(\\theta\|\\mu\)\-f\_\{l\-1\}\(\\theta\|\\mu\)\)^\{2\}\\bar\{\\pi\}\_\{l\}\(d\\theta\)\\pi\_\{0\}\(d\\mu\)Sincefl​\(θ\|μ\)≤Fl​\(θl,θl−1\|μ\)≤C′f\_\{l\}\(\\theta\|\\mu\)\\leq F\_\{l\}\(\\theta\_\{l\},\\theta\_\{l\-1\}\|\\mu\)\\leq C^\{\\prime\}, and we defineC′′=maxθ∈Θ,μ∈𝖲⁡φl−1​\(θ,μ\)C^\{\\prime\\prime\}=\\max\_\{\\theta\\in\\Theta,\\mu\\in\\mathsf\{S\}\}\\varphi\_\{l\-1\}\(\\theta,\\mu\)

Vl\\displaystyle V\_\{l\}≤1C′​∫Θ×𝖲\(φl​\(θ,μ\)−φl−1​\(θ,μ\)\)2​π¯l​\(𝑑θ\)​π0​\(𝑑μ\)\\displaystyle\\leq\\frac\{1\}\{C^\{\\prime\}\}\\int\_\{\\Theta\\times\\mathsf\{S\}\}\(\\varphi\_\{l\}\(\\theta,\\mu\)\-\\varphi\_\{l\-1\}\(\\theta,\\mu\)\)^\{2\}\\bar\{\\pi\}\_\{l\}\(d\\theta\)\\pi\_\{0\}\(d\\mu\)\+C′\(C′′\)2∫Θ×𝖲\(fl\(θ\|μ\)−fl−1\(θ\|μ\)\)2π¯l\(dθ\)π0\(dμ\)\\displaystyle\+C^\{\\prime\}\(C^\{\\prime\\prime\}\)^\{2\}\\int\_\{\\Theta\\times\\mathsf\{S\}\}\(f\_\{l\}\(\\theta\|\\mu\)\-f\_\{l\-1\}\(\\theta\|\\mu\)\)^\{2\}\\bar\{\\pi\}\_\{l\}\(d\\theta\)\\pi\_\{0\}\(d\\mu\)Assume thatfl​\(θ\|μ\)f\_\{l\}\(\\theta\|\\mu\)is Lipschitz w\.r\.tφl​\(θ,μ\)\\varphi\_\{l\}\(\\theta,\\mu\)with Lipschitz constantGG≤C~​\(∫Θ×𝖲\(φl​\(θ,μ\)−φl−1​\(θ\)\)2​π¯l​\(𝑑θ\)​π0​\(𝑑μ\)​\(φl​\(θ,μ\)−φl−1​\(θ,μ\)\)2​π¯l​\(𝑑θ\)​π0​\(𝑑μ\)\)\\displaystyle\\leq\\tilde\{C\}\\Bigg\(\\int\_\{\\Theta\\times\\mathsf\{S\}\}\(\\varphi\_\{l\}\(\\theta,\\mu\)\-\\varphi\_\{l\-1\}\(\\theta\)\)^\{2\}\\bar\{\\pi\}\_\{l\}\(d\\theta\)\\pi\_\{0\}\(d\\mu\)\(\\varphi\_\{l\}\(\\theta,\\mu\)\-\\varphi\_\{l\-1\}\(\\theta,\\mu\)\)^\{2\}\\bar\{\\pi\}\_\{l\}\(d\\theta\)\\pi\_\{0\}\(d\\mu\)\\Bigg\)Applying the Proposition[E\.1](https://arxiv.org/html/2610.19074#A5.Thmprop1)on it and lettingβ≤2​s−1\\beta\\leq 2s\-1≤C​2−β​l\.\\displaystyle\\leq C2^\{\-\\beta l\}\.∎

## Appendix IFuture research directions

For future work, we list several potential projects to consider below\.

- •The examples given in this work were to firstly verify the theory, and also the provide a justification\. As we have shown this, it would be of interest to consider more challenging datasets\. intensive and high\-dimensional examples\. Such examples could include bigger datasets such as IMDB or the CIFAR datasets\. If dimensionality is a problem, one could exploit PCA\.
- •A useful direction would be to consider the complexity in terms of the depth of a neural network\. We have only done this, in terms of the width, motivated by the original paper[Sell and Singh \[2020\]](https://arxiv.org/html/2610.19074#bib.bib15)\. Doing so would provide a fuller understanding of its impact algorithmically\.
- •Related to the above point, a natural question is to ask how it would compare to generic state\-of\-the\-art MCMC methods like HMC\. This has been answered, in comparison with parallel\-SMC, but not in the context of our multilevel trace\-class neural network setup\. To do so, one could run SMC samplers withHamiltorchmutations, and thus provide a comparison\. This goes beyond this work, and will be done in a future paper\.
- •A final direction would be consider the use of such multilevel\-trace class priors in the context of non\-fully Bayesian methods\. An example of this would be unadjusted Langevin algorithms, which have proven to be successful, which include recent multilevel algorithms[Chada et al\. \[\]](https://arxiv.org/html/2610.19074#bib.bib25),[Paulin et al\. \[2025\]](https://arxiv.org/html/2610.19074#bib.bib24)\. This would involve a comparison between SMC/MCMC methodologies, as tested in this work, and ULA\-algorithms\.