Multilevel neural simulation-based inference
Authors: Yuga Hikida
Ayush Bharti
Niall Jeffrey
Francois-Xavier Briol
Venue: NeurIPS 2025 poster ()
Keywords: Simulation-based Inference
Multilevel Monte Carlo
Neural posterior estimation
Neural likelihood estimation
Primary area: probabilistic_methods
TL;DR: We propose a multi-fidelity simulation-based inference method based on multilevel Monte Carlo for computationally costly simulators.

Abstract: Neural simulation-based inference (SBI) is a popular set of methods for Bayesian inference when models are only available in the form of a simulator. These methods are widely used in the sciences and engineering, where writing down a likelihood can be significantly more challenging than constructing a simulator. However, the performance of neural SBI can suffer when simulators are computationally expensive, thereby limiting the number of simulations that can be performed. In this paper, we propose a novel approach to neural SBI which leverages multilevel Monte Carlo techniques for settings where several simulators of varying cost and fidelity are available. We demonstrate through both theoretical analysis and extensive experiments that our method can significantly enhance the accuracy of SBI methods given a fixed computational budget.

Multilevel neural simulation-based inference
Yuga Hikida Aalto University yuga. hikida@aalto..ff
Ayush Bharti Aalto University ayush. bharti charti@aalto. @aalto..a.......................
Jeffree
François-Xavier Briol University College London f. briol0ucl... acc.acc...
University College London n. j.fffeeyuuuc..c ac.
Abstract
Neural simulation-based inference (SBI) is a popular set of methods for Bayesian inference when models only available in the form of a simulator. These methods are widely used in the sciences and engineering, where writing down a likelihood can be significantly more challenging than constructing simulator. However, the performance of neural SBI can suffer when simulators are computationally expensive, thereby limiting the number of simulations that be performed. In this paper, we propose novel approach to neural SBI which leverages multilevel Monte Carlo techniques for settings where several simulators of varying cost and fidelity are available. We demonstrate through both theoretical analysis and extensive experiments that method significantly enhance the accuracy of SBI methods given a fixed computational budget.
Introduction
Simulation-basee inference (SBI) [Cranmer et al., 2020] is set of methods used to estimate parame- ters of complex models for which the likelihood is intractable but simulating data is possible. is particularly useful in fields such cosmology [Jeffrey al., 2021], epidemiology [Kypraios et al., 2017], ecology [Beaumont, 2010], synthetic biology [Lintusaari et al., 2017], and teleconmunications engineering [Bharti et al., 2022a], where models describe intricate physical or biological processes such as galaxy formation, spread of diseases, the interaction of cells, or propagation of radio signals.
For long time, SBI was dominated by methods such as approximate Bayesian computation (ABC) [Beaumont et al., 2002, Beaumont, 2019], which compared summary statistics of simulations and of the observed data. However, SBI methods using neural networks to approximate likelihoods [Papamakarios al., 2019, Lueckmann et al., 2019, Boelts et al., 2022], likelihood ratios [TThmmas al., 2022, Durkan et al., 2020, Hermans et al., 2020], or posterior distributions [Papamakarios and Murray, 2016, Lueckmann et al., 2017, Greenberg al., 2019, Radev al., 2022] are now quickly becoming the preferred approach. These neural SBI methods are often favoured because they allow for amortisation [Zammit- Manniin et al., 2024], meaning that they require a large offline cost to train the neural network, but once the network is trained, the method can rapidly infer parameters for new observations or different priors without requiring additional costly simulations. This particularly useful when the simulator is computationally expensive, it reduces the need to repeatedly run simulations for each new inference task, making the overall process significantly less costly.
Nevertheless, the initial training phase of neural SBI methods typically still requires a large number of simulations, preventing their application on computationally expensive (and often more realistic) models which can take hours of compute time for simulating a single data- point. Examples include most tsunami [Behrens and Dias, 2015], wind farm [Kirby et al., 2023], nuclear fusion [Hoppe et al., 2021] and cosmology [Jeffrey et al., 2025] simulators.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

One avenue to mitigate this issue is multi-fidelity methods [Peherstorfer et al., 2018]: we often have access a sequence of simulators with increasing computational cost and aceu- racy which we may be able to use to refine existing methods based on single simulator. This setting is quite common. One example simulators requiring the numerical solution of ordinary, partial, or stochastic differential equa- tions, where the choice of mesh size or step- size affects both the accuracy of the solution and the computational cost. Another example arises when modelling complex physical, chemical, or biological processes, where low-fidelity simula- tors can be obtained by neglecting certain aspects of the system. This is exemplified in Figure 1, where a high-fidelity cosmological simulation in- cludes baryonic astrophysics and thus appears smoother than the low-fidelity simulation.
10R
10a
Urendessry
139
10-1
Low-fidelity
High-fidelig
Figure I: Low- and high-fidelite cosmological simulations from the CAMELS data [Villaescusa- Navarro et al., 2023] studied in Section 5.4.
The key idea behind multi- fultii-fieelity methods is to leverage cheaper, less accurate simulations to supplement the more expensive, high-fidelity simulations, ultimately improving efficiency without sacrificing accuracy. This has been popular in the emulation literature since the seminal work of Kennedy and O'Hagan [2000], but applications to SBI are much more recent and include Jasra et al. [2019], Warne et al. [2018], Prescott and Baker [2020, 2021], Warne et al. [2022], Prescott et al. (2024), who proposed multi-fidelity versions of ABC. More recently, [Krouglova et al., 2025] also proposed a multi-ffdeliin method to enhance neural SBI approximations of the posterior based on transfer learning. However, their method is not supported by theoretical guarantees.
In this paper, we propose a novel fulti-fidelity method which is broadly applicable to neural SBI methods. Taking neural likelihood and neural posterior estimation the main case studies, show that our approach is able to significantly reduce the computational cost of the initial training through the use of multilevel Monte Carlo [Giles, 2015, Jasra et al., 2020] estimates of the training objective. The approach has strong theoretical guarantees; our main result (Theorem 1) directly links the reduction in computational cost to the accuracy of the low-fidelity simulators, and we demonstrate (in Theorem 2) how to best balance the number of simulations at each fidelity level in the process, Our extensive experiments on models from finance, synthetic biology, and cosmology also demonstrate the significant computational advantages provided by our method.
2 Background
We first recall the basics of SBI methods and the related works reducing computational cost in Section 2.1, then provide a brief introduction to multilevel Monte Carlo Section 2.2.
2.1 Simulation-based inference
Let (Falaee be parametric family of distributions on some space X [Rdx parameterised by o θ Rde We assume that this model is available in the form of computer code, i.e. simulator, where simulating from Po is straightforward, but the likelihood function p(- [0] associated with Po is intractable. Simulators can be characterised by a pair (U, Go), where U is distribution (typically simple, such uniform randomness, and G. 24 + X' a (deterministic) parametric map called the generator. Simulating Pg can be achieved by first simulating u U, and then applying the generator x (u(u). In this paper, we consider Bayesian inference for the parameters
Gaussian) space U C dd which captures all of the
this simulator-basee model given independent and identically distributed (iid) data 2 {x; {- 2 collected from some data-generating process. Specifically, we are interested in approximating the posterior with density (0 ff.m) x I p(>? (0), where (0) is the prior density. As introduced below, this be achieved via neural SBI method which approximates the likelihood posterior.
Neural likelihood estimation (NLE). NLEs [Papamakarios et al., 2019, Lueckmann et al., 2019, Boelts et al., 2022, Radev et al., 2023a] are extensions of the synthetic likelihood approach [Wood, 2010, Price et al., 2018] that use flexible conditional density estimators, typically normalising flows

---

[Rezende and Mohamed, Papamakaris et al., 2021], as surrogate models for the likelihood
function associated with Pg. The surrogate conditional density 9s NE X' Q + [0, oc) where T6 (- 0) is a density function for each 0 E ② and & E Φ S I denotes its learnable parameters, is trained by minimising the negative log-likelihoog with respect to & on simulated samples. More precisely, let {00,, be the training data such that 8 ~ π are realisations from the prior and 2; ~ Ie, are realisations from the simulator. Then éMc == arg min ENLE (M)),, where (NLE is an 'CC
empirical (Monte Carlo) estimate of negative expected log-density:
(Ø) ;= -Bgwa [EewEo
NLE (x go
[0]] - E Di log To
NLE (2; 0.) =; "MC
(M)
Once the surrogate likelihood is trained, Markov chain Monte Carlo (MCMC) or variational infer- ence methods are used to sample from the (approximate) posterior distribution TNLE (0) == xx (x 0) (0). NLEs can therefore be regarded as being partially amortised- the surrogate likelihood need not be trained for a new observed dataset, however, MCMC needs to be carried out again to obtain the new posterior.
m 9 Ai
For a computationally costly simulator, bottleneck, which affects the accuracy of estimating the expected loss. Thus, estimating the loss accurately with fewer samples is key to handling costly simulators.
note that obtaining training samples can become
Neural posterior estimation (NPE). Instead of learning a surrogate likelihood, NPEs learn a mapping p(0 x) from the data to the posterior using conditional density estimators. These are often based on mixture density networks [Bishop. 1994, Papamakarios and Murray, 2016] or normalising flows [Dinh et al., 2014, Papamakarios et al., 2017, Radev al., 2022]. Similar to NLE, the conditional density NPE ② x X + [0, 00) is trained by minimising the negative log likelihood with respect to ó using data {0,, x i.m generated by first sampling from the prior 0, ~ π and then the simulator Ti:m,i (I1,, ................................ Im.i) Po::
NPE (0) = -EEas Emrm- To
log
NPE (0 (mmm) d0
2 - 1 E.
log 9 1
NPE 11:m.a) =:
NPE MC (Ø)
Once @ is estimated, the NPE posterior is obtained as TNPE (0) =f..) = NPE (0) x..)) Although training 4% NPE incurs an upfront cost, this is a one-time cost as approximate posteriors for new observed datasets are obtained by simple forward pass of æf.m. through the trained networks, making NPEs fully amortised (in contrast with the partial amortisation of NLE). Similarly to NLE, the computationally costly step in NPE is the generation of training samples from running expensive simulators. Note that both of and d NPE usually include a summary function (often architecturally implicit in NPE). This is helpful when X' is high-dimensional or the number of observations m is large [Alsing et al., 2018, Radev et al., 2022], and for NPE it allows conditioning sizes. Recently, alternative training objectives for NPE based on flow matching [Wildberger et al., 2023], diffusion [Geffner et al., 2023, Sharrock et al., 2024, Gloeckler al., 2024], and consistency models [Schmitt et al., 2024b] have been proposed.
datasets of different
Related work. We briefly note that beyond the aforementioned multi-fidelity methods, other works also aim to reduce the computational cost of SBI. In the context of ABC, adaptive sampling of the posterior using either sequential Monte Carlo techniques [Sisson et al., 2007, Beaumont et al., 2009, D. Moral et al., 201 or Gaussian process surrogates [Gutmann and Corander, 2016, Meeds and Welling, 2014] has been a popular approach. A similar approach has been applied to neural SBI [Papamakarios and Murray, 2016, Lueckmann et al., 2017, Greenberg et al., 2019, Papamakarios et al., 2019, Hermans et al., 2020, Durkan et al., 2020], where sequential training schemes are employed to reduce the number of calls to the simulator. Other works have tackled this problem using cost-aware sampling [Bharti al., 2025], side-stepping high- high-dimensional estimation [Jeffrey and Wandelt, 2020], early stopping of simulations [Prangle, 2016], dependent simulations [Niu et al., 2023, Bharti et al., 2023), expert-in-the-logg methods [Bharti et al., 2022b), self-consistency properties [Schmitt et al., 2024a], parallelisation of computations [Kulkarni and Moritz, 2023], and the Markovian structure of certain simulators [Gloeckler et al., 2025]. Our proposed method, introduced in Section 3, can be combined with all of these compute-efficient SBI methods and is hence complementary to them.
3

---

2.2 Multilevel Monte Carlo method
Consider some square-integrabbe function f 2 -+ R and distribution μe on domain 2 C Rd2. We consider the task of estimating [=)]] A first approach is standard Monte Carlo (MC) [Robert and Casella, 2000, Owen, 2013], which yields the following estimator: 1/n Li=1 (zz), where 21:-..,. An ~ J.. The roo--mean-squared error (RMSE) of this estimator converges at a rate ((n-//) where the rate constant is controlled by Var [f)] [Owen, 2013, Ch. 2]. In cases where f is expensive to evaluate f is expensive sample from, the RMSE can therefore be relatively large.
This issue can be mitigated by multilevel Monte Carlo (MLMC), which was first proposed by Heinrich [2001], Giles [2008], and more recently reviewed in Giles [2015], Jasra et al. [2020]. Suppose we have a sequence of qquare-integrable functions fi : 2 - R for I E {0, 1......,, L} which are approximations of and which are ordered such that fu = f. and both the cost of evaluation Cj and the accuracy (or fidelity) of fi increase with 1. In that case, MLMC consists of expressing [x)]
through a telescoping sum and approximating each term through MC based on samples 2 ~ p for i= l...... n4 and I E {0, 1,....... L)::
E Eg [f]]] = Ezzz [f(())] + Et. Eewn [f((z) fi-1
L 10g
Σ fo(2)) + +
1. E
fi-1 (i))
Note that this can also be thought of as using the low-fidelity functions as approximate control variates. By carefully balancing the number samples no. ................................ Rr, according to the costs Co.... Cr. and variances Var[]o(=)], V=r - fo(=)]..... (r.(Sr.(=) fr. [=)] at each level, one can show that MLMC can significantly improve on the accuracy of MC given fixed computational budget. For this reason, MLMC has found numerous applications in statistics and machine learning. including for optimisation [Asi et al., 2021, Hu et al., 2023, Yang et al., 2024], sampling [Dodwell et al., 2019, Jasra et al., 2020], variational inference [Fujisawa and Sato, 2021. Shi and Cornish, 2021], probabilistic numerics [Li et al., 2023, Chen et al., 2025] and the design of experiments [Goda et al., 2020, 2022].
3 Methodology
We now present NLE and NPE versions of our approach, termed multilevel-NLE and multilevel-NPE.
MLMC for NLE and NPE, Recall that we are performing inference for a simulator (Ge, U) with prior π the parameter 0 O. We reparameterise the NLE and the NPE objective in terms of u instead of æ (akin to the reparametrisation trick in variational inference), and express them more broadly using an arbitrary loss é + - IR (to represent either @NLE or ENPE) and an arbitrary function fø : tm x θ + R (to represent either f or SNPE) as:
{(0)= == Eg um [fo(Urm, 0)]:
where for NLE we have for (NLE(U, 0) log ALE (Go(u) 00) = log NLE (x 0) with m 1, and for NPE we have FNPB (u1mm 0) == - log o (0) Go(u))... G.(um)) Hereafter, we present our methods using f and in order to avoid duplication.
log & (P 11:m).
The MC estimator of the loss {()) is given by (M)()) == 1 E'= fa 0.). Note that the variance of this MC estimator depends on the number of iid samples n. Hence, a small n owing to a computationally expensive simulator will lead to a poor estimator for the loss. Now suppose that we have access to a sequence of generators (u), G² ((u) with varying fidelity levels. Then, for
I = O. L - 1, we can define a corresponding sequence of functions f%. 4m x O + R: NLE. (ui=m 0) = - log co &
NLE (t) (), and 5 pNPE. (11::m 0) == - log g (P (G(t))...,GGC
such that C%
Go(u) and f (u1:m; 0) == (u 9). This sequence of functions gives evaluations of the log conditional density at evaluations of the (approximate) simulator. Recall that the larger the value of I, the more accurate (and computationally expensive) such simulations will tend to be. At this point,
re-express the objective using a telescoping sum
(()) = Eg urom ~ [fo (U1m, 0)1 = Eg ~ [f (u::m: })]
E [ [8 0) + 1==1 Bo w [fa(uim, 0) (uim(w) (1)

---

Equation (I) follows by adding then subtracting some terms, and using linearity of expectations. Now suppose that we can simulate from each of these approximate simulators to obtain:
{ 0% u/ um Go/ (u),a) Gog ( (
where 0g ~ u u U,
for 1,..., ny and I ....., L. We can then these simulations to approximate each term in the telescoping sum through Monte Carlo estimator as follows:
....,
no
12
(()) MMLMC (0)
(u 0.)
(f(wimmss (u 03) (2)
T1
VIL
(d)
hi(d)
This corresponds to an (unbiased) MLMC estimator of our original objective in (I). The first term, ho(o), approximates f()), but is biased since it uses fg (i.e. the lowest fidelity simulator) rather than f (the highest fidelity simulator). The additional terms, h1 (ő)..... .u.))......... (o), correct this bias by estimating the expected difference between the objectives at consecutive fidelity levels.
Within each term (o), the functions f and f/=1 are evaluated the same samples u' and 6', i.c., they are seed-matched. This ensures that f (u1..m, 02) and fll-(ullt (u11.m,: i 0,) are coupled and highly correlated, which leads to a reduction in variance [Owen, 2013, Ch. 8] since
Var[h
(3)
- (Varl(' (u (m 0,] + (ul 0D 2Cov (1' 0,), (u o
and the covariance will be large. This covariance will be particularly large the more similar the two functions l and f are, and therefore expect [ør] be smallest in those settings. We also immediately see that without seed-matching, the covariance will be small and the variance will be large, highlighting why seed-matching essential for MLMC.
Computational cost, Let Ct be the computational cost of sampling one from the 1th level generator Go and evaluating fg and recall that Co C1 <...< CC, Then, the cost of iLMC
=0(n (noCo Et=r mi (C) + CL))
Cost(fMLMc(0)) no. TL) O
while that of the MC estimator is Cost(Emc(o):n)) = ((CL)). Using solely the high-fidelity generator (as is customary in SBI) would require large in order reasonably estimate A thus increasing the total cost. However, with multiple lower-fidelity generators available, we can have a different number of simulated samples per level (i.e. we can take ni0 * ni * ... * nt.), and can select ro...... .,, n. such that ni < f1-1- This allows us to take a much larger number of samples from the cheaper (or low Ci) approximations of the simulator, and much smaller number of samples from the expensive (or high Ci) approximations of the simulator, making MLMC particularly attractive for reducing the total computational cost of simulation in neural SBI.
Extensions. There are several straightforward extensions of our approach which are not covered above so as to not overload notation. Firstly, each simulator could have its own base measure U', which could be defined on spaces ( of different dimensions. This is not a problem since could simply consider U to be the tensor product measure and 24 to be the corresponding tensor product space, in which case all equations above remain valid. Similarly, the parameter space may differ across simulators. However, to ensure the best possible performance, it will still be essential to seed-match random numbers where there is overlap; see Owen [2013, Ch. 8] for more details on the use of common random numbers, and Section 5 for a study of this issue for multilevel neural SBI.
Secondly, although our discussion has pertained to multilevel-NLL and multilevel-NPE so far, it straightforward to extend the MLMC approach to other neural SBI methods such as neural ratio estimation [Hermans et al., 2020, Durkan et al., 2020, Miller et al., 2022], score-based NPE [Geffner al., 2023], flow-matching NPE [Wildberger et al., 2023), and GAN-based NPE [Ramesh al., 2022] since these are all based on objectives which can be expressed as MC estimators.
5

---

Optimisation. Gradient-based optimisation of the MLMC objective (o) can be challenging due to the 'conflicting' gradients which appear in consecutive terms Vohi(() and V Voh4+1(0). More precisely, the term V (o) always contains
Algorithm MLMC gradient adjustment
Input V Voho(o), Is''t 5'5
3 -
0
I! Rescaling the gradients:
for 1or1110 1==10
11c2-=12 112 1-1,-
c d- end for
110 65l11 Sp
ch++ + = R2 E of ( O)),
V Toh(o)) + Σ1= (Stt + +50 +1-11-)
which approximates V VoElAA whilst the term V Vohi+1 always contains
11 Projecting the gradients (only when confticting)
V oho(o)) V Vohe(o) < then
E
V hho(l) + Voh((d)
(chold)-oohhll
( H
V (pl)(()
ML+1
V voho(0) + V Voho(o) Vahalo) vahe))
V Voho(o)
which approximates -V In
else V To)o(d) t (o)o (ohc(0) t (ohc(Ø)
the safiee limit, Vollif and [[)) cancel out, but this is not the case for sa and Ch since we are typically work- ing with only small number of expensive
end if Output V VoLMLMC()) + voho(o) + Voho(d)
simulations. When naively applying standard gradient-based optimisation methods on this loss, we observe that the training dynamics is typically dominated by only one of these two quantities until approaching stationarity, at which point the conflicting gradients lead to unstable updates and, ultimately, often cause divergence.
To mitigate this issue, we use a combination of gradient adjustments summarised in Algorithm 1. Firstly, we rescale ct and c/ ensure that they have comparable norms and that their difference remains small and stable. Secondly,
apply the gradient projection technique of Liu al. [2020], projecting the gradients of (o) and he(6) := Σi= by (o) onto each other's normal planes to reduce the impact of conflicting gradients. We observe empirically that combining these two techniques significantly improves the stability of the optimisation throughout the training and leads to better
performance; see Appendix B.6 for detailed comparison.
Theory
We now present our main theoretical results. Theorem 1 expresses the variance of each of the terms in the telescoping sum approximation (see (2)) as a function of the number of simulations per level, and the magnitude of the difference in generators between consecutive levels.
We say that f is log-concave if it has density of the form exp(->(2)) for some convex function 2 2 + IR. We recall that for E [1, oo), d, d' E N and a non-empty, open, connected set 2 C IRdd the space of vector-valued r-integrable functions with respect to a probability distribution is given by ('(x)= == {g: Z + Id ( == (ʃz [=(=))])(ddz))" < co). For T E N. the corresponding Sobolev space of vector-valued functions of smoothness r is given by W (n) := 2 1 D is the weak derivative operator corresponding to Finally, recall that function g is locally -mopthh if its gradient is locally Lipschitz continuous with Lipschitz constant Krip > O; i.e. for all 2, z' in some open set of Z, we have that [|Vg(=)= )) Vg(2')ll2 S Kripll= z'112. For simplicity, we will write t6 (#1m 8) for the conditional density model used for either NLE (in which case NLE
{g
Rd [[9llw->(()) = jalsr D (1) < oo), where for a multi-index a E N
Lo (1 0) = L (21 (0)) and NPE (in which case q (Z1mm 0) = 1s
NPE (0) m1:m)).
Theorem 1. Let & E and suppose the following assumptions hold:
(Al) The prior π and the base measure IU are log-concave distributions.
(A2) The generators satisfy [["|llwnnraxop < St for E 10,,1...,,1)
(A3) log qe is continuously differentiable, locally (0) smooohh and satisfies the growth condition V/ log To (Zim, 0)ll S (o)( K grow (0) > O.
||o:||2 + |101|2 + 1) for some

---

Then,for E fl...................... L D and Kolo)....,KK (δ) > 0 independent of no..... ..... R.L. we have that:
Var [ar]) - Ko(Ø) (IG)) (I. +
K((Ø) | - G'-11 GI-1I
and Var[hi(0)) -
n2
See Appendix A.I for the proof, and we emphasise that the variance is over both parameters 0 and noise, but conditional on Φ. (A1) requires log-concavity, which is strong condition, but this is only required of the prior and the base measure. However, in SBI these tend to be simple distributions such as Gaussians or uniforms, which satisfy this assumption; see Saumard and Wellner [2014) for how to verify log-concavity in practice. (A2) relatively mild: it asks that Go, one derivative in 0 and u, and for these generators and their derivatives to have fourth moment. This will hold when the simulators are used to define sufficiently well-behaved data-generating processes, but will be violated for sufficiently heavy-tailed distributions (e.g. the Cauchy). Finally, (A3) is mild; holds when the gradient of the log conditional density estimator is Lipschitz continuous, which is for example the case when the conditional density is twice continuously differentiable and the Hessian is bounded (see c.g. Lemma 2.3 of Wright and Recht [2022]). also holds for models such mixtures of Gaussians, which are used in mixture density networks and for normalising flows with sufficiently regular transformations; see Table in Liang et al. [2022] for some known Lipschitz continuous transformations. There two key implications of Theorem 1. The first is bound on the variance of the MC objective:
G have least
Var (+)) S K(ő) (IG1Ii +1) n2
which is obtained by noticing that the bound on Var[ho[@]] is simply a bound on a Monte Carlo objective. From this, we immediately notice that the bound will be large whenever the high-fidelity simulator is expensive and n is small, or whenever the simulator is complex as measured in this Sobolev norm. The second implication is a bound on the variance of the MLMC objective:
Var [6ML
Ko(d) ni
(1IG [[G"I (XXU( +1) + Σ
Ki(0) ||'' d'-'Ilwaa(xxx)
(5)
In order make this bound small, we need make each term small. Typically, be large, but this will be counterbalanced by taking to be large. For the higher-fidelit levels, the number of samples ny will typically smaller, but this will be mounte-bbllance, by the fact that ||G' Gi II we can directly use (5) to get an indication of how to select ,
will
(xx)( is small whenever Gd11 is a good approximation of G¹. As we now show,
Theorem 2. Suppose the assumptions of Theorem I hold. Then, the values of no..... un which minimise the upper bound on VarlEMLMC(() (o)] in (5) given a fixed computational budget; i.e. for Cost((MMc((o): To. S Cuudger for some Ohudger > O, are given by
Chadges
VCo VIGGiI (III(EX + ni x 1
VC++CL-1 ||G' -G' -I''' (ll'IIlwa((xx))
Obudget
no x
See Appendix A.2 for the proof. This result provides useful intuition on how to select the number of samples per level. For instance, consider the case L 1, i.c., two levels. If G is known to be a good approximation of G it makes sense to allocate a large budget to generating low-fidelity simulations while only a small number of high-fidelin simulations would be sufficient. On the other hand, if Go and G differ substantially, allocating a larger budget to high-fidelity simulations makes sense, despite the higher cost, it can help accurately capture this difference. Although Theorem 2 provides intuition, it may be hard to obtain the optimal A...... nr exactly since it requires computing quantities which are often unknown. For example, computing Sobolev norms can be challenging, and the results implicitly depend on e through the constants Ko()) (o))). (...... ..... (x) of Theorem (which are unlikely to be tight). This is a common limitation of MLMC theory: see [Giles, 2015, Sec. 2-3).
Before concluding this section, we note that our theory focussed on the variance of (MLMC (+)), but another important quantity for gradient-based optimisation will be the variance of the gradient V (oCMLMM(()) of this objective. turns out that similar results are straightforward to prove for this quantity under very minor modifications of the assumptions, see Appendix A.3.
7

---

Empirical Coverage
NPE (kipk))
KLD
10P
NLPD
ML NLE (200)
NIE (kos))
Demsity
0.2
NLE (high)
NLE (low)
1071
0.1
-10
10
300
100
100
0.0
==
Confidence Level
(a) ML-NLE
(b) ML-NPE
(c) NLE comparison
(d) NPE comparison
Figure 2: Performance of our ML-NLE and ML-NPE method on the g-and-k example. KL-- divergence (L) between the estimated and the (almost) exact density for ML-NLE under different high-fidelity samples 721. We compare it with NLE (low) trained on only low-fidelity data (n 10*) and NLE (high), trained on only high-fidelity data (n 300). (b) Negative log-posterios density (NLPD L) for ML-NPE, NPE (low) with n 10 and NPE (high) with 100. (c) One instance of learned densities using NLE. (d) Empirical coverage plot for ML-NPE, NPE (low), and NPE (high).
5 Numerical Experiments
We compare the performance of our multilevel version of NLE and NPE, termed ML-NLE and ML-NPE, respectively, against their standard counterpart with the MC loss. We use the sbi library (Teeerr-CCanter Cantero et al., 2020] implementation for NLE and NPE, see Appendix for the details. The code reproduce experiments is available http://gihuuboommmuugakkkkllitit
5.1 The g-and-k distribution: an illustrative example
We first consider the g-and-k distribution [Prangle, 2020] an illustrative example. This is very flexible univariate distribution that is defined via its quantile function and has four parameters, controlling the mean, variance, skewness, and kurtosis respectively, making it challenging for SBI methods. It does not typically have a low-fidelity simulator, so we construct one through a Taylor approximation of the quantile function, see Appendix B.1 for details. This makes it slightly contrived example, but the fact that the g-and-k allows for an efficient approximation of the likelihood will make it particularly convenient to study the performance of NLE-based methods.
We fix number of low-fidelity samples (ri- = 10 for ML-NLE and ni = 100 for ML-NPE) and vary the number of high-fidelity samples n to asses the improvement in performance of our methods as nı increases. For MI-NLE, we compute the Kullback-Leeblle divergence (KID) between estimated conditional density and numerical approximation of the likelihood. For ML-NPE, use the negative log peg-pstterio density (NLPD) of the true 0 under the estimated posterior density the metric. The results in Figure 2a and 2b show that the multilevel versions of NLE and NPE perform better than their standard counterparts (with MC loss) using just a fraction of the high-fidelity samples. Unsurprisingly. the performance of MLE-NLE and ML-NPE improves n1 increases.
In Figure 2c, show an example of the NLE densities learned, with additional examples in Appendix B.I. Our ML-NLE method is able to approximate the almost exact g-and-k density the best. The coverage plot in Figure 2d shows that ML-NPE yields slightly conservative posteriors as opposed to the overconfident posteriors obtained from NPE trained on either all low- or all high-fidelig data.
5.2
Ornstein-Uhlenben process: a popular financial model
Our next experiment involves the Omsteen-Uhllnbeekk (OU) process-a stochastic differential equa- tion model commonly used in financial analysis [Minenna, 2003], which in our case outputs 100-dimensional Markovian time-series and has three parameters. The process is known to converge to stationary Gaussian distribution, which we use as the low-fidelity simulator, see Appendix B.3. The example is reproduced from Krouglova et al. [2025], who used it for benchmarking their transfer learning approach to NPE (TL-NPE). TL-NPE first trains an NPE network on low-fidelity simula- tions, and then uses a small set of high-fidelity data to refine the network parameters until a stopping criterion is met. We implement TL-NPE with n4 = 1100 and n; (10, 100), and compare against {10, 100). These values are picked to keep the simulation
ML-NPE with no 1000 and n
8

---

NLPD
XLD
(a) 10
(b) 100
Figure 3: Performance of ML-NPE and TL-NPE measured by NLPD (L) and KL divergence (L) with different number of high-fidelity samples n1 the OU process. (a) When ni 10, ML-NPE outperforms on both metric. (b) When nı 100, TL-NPE outperforms on both metric. Note that for visualisation purpose, removed outliers, see Appendix B.3 for the result with all the data.
budget the same for both methods. TL-NPE uses the default early stopping criterion from the sbi library Tiejero-Cantero et al., 2020], which terminates training if the validation loss increases for 20 epochs. Once the criterion is met, the network parameters achieving the lowest validation loss are selected. This criterion is used for both the low- and the nigh-fidelity training stages. We use 20% of data as validation set for early stopping.
We run both methods 20 times and compute the NLPD of the ground-truth parameter under the report the KLD between the posterior
estimated posterior across 500 test points. Additionally, obtained using TL-NPE or ML-NPE and reference NPE posterior trained with n = 10,000 high- fidelity simulations. Results aggregated across the 20 runs and all test points, and reported in Figure We observe that both the methods achieve comparable performance: for nı = 10, ML-NPE performs slightly better both metrics, whereas for n 100, TL-NPE slightly outperforms ML-NPE on both.
We additionally investigate a challenging scenario where the high- and low-fidelity simulators have differing dimensionalities of 8, in which our method exhibits some limitations, possibly due to instability during optimisation; see Appendix B.3.1 for details.
5.3 Toggle-switch model: a Systems Biology example
2.0
We now consider the toggle-switch model [Bonassi al., 2011, Bonassi and West, 2015]. This model describes the interaction between two genes over time, and has seven parameters and scalar observation at the end of a time interval. Simulators typically use time-discretisation, and the total number of time-steps T acts as fidelity parameter: running the model with large 7' incurs larger computational cost but leads to accurate simulations, while smaller T' leads to cheap but inaccurate samples. Thus, this model illustrates a setting with more than two fidelity levels (by taking T to be more than values). We take Ty 50, Ti 80, and Te = 300 to be the number of steps for the three fidelity levels with no 104, nı 500, and n2 100. Guided by the intuition from Theorem 2, we allocate a large budget to estimating the difference between the low- and the cndiium-fidelity simulators, leveraging prior knowledge that this difference substantial. See Appendix B.4.1 for results under alternative budget allocations. Note that in this case each fidelity level has different base measure of dimension 27' - 1; however, there still sufficient seed-matching (see Appendix B.4 for details), and therefore variance reduction. thanks to MLMC. We compare our ML-NLE method with the standard NLE trained on samples from either the low- (n 12, 060), the medium- (n 7537), or the high-fieeity simulator (n 2010). The number of training data n in each case is selected so as to match the total computational cost of simulation between ML-NLE and NLE, see Appendix B.4. Figure reports the maximum mean discrepancy (MMD) [Gretton et al., 2012] between 500 samples from the learned conditional densities and 500 samples from the high-fidelity simulator across for 5000 different parameter values. We observe that ML-NLE performs better than all the NLE baselines for the same computational cost.
1.5
MMD
1.0
0.5
0.0
ML-NLE
NLE (high)
NLE (med.)
NLE (low)
Figure 4: MMD (4) across 5000 parameter values for NLEE and MI-NLE.

---

5.4 Cosmological Simulations
We now consider a cosmological simulator using the CAMELS suite (Villaescusa-Navees et al., 2021, 2023]-one of the most computationally intensive cosmological simulations to date-which comprises both low- and high-fidelity data. These are state-of-the-an simulations being used with real- world data. Developing surrogate, multilevel techniques for cosmology has become key research focus [Chartier et al., 2021, Chartier and Wandelt, 2022]. Our task is to infer a standard cosmological target parameter using 39-dimensional power spectra of cosmology data, see Appendix B.5 for an example. The low-fidelity simulations are gravity-only N-body simulations, whose physical behaviour is controlled only by the parameter and some Gaussian fluctuations of the initial conditions. The high-fidell fideliiy hydrodynamic simulations have additional physics, controlled by an additional five parameters. We include these additional parameters as part of the U/ space, constituting case of partial common random numbers between the low- and high-fidelity simulators, similar to Section 5.3.
Here, the high-fidelit simulations can be orders of magnitude (> x x100) slower to generate than the low-fidelity ones, making this representa- tive problem for our method. Assuming we only have access to n n1 20 figh-fideiity sim- ulations, we wish to ascertain the improvement in inference accuracy by including 1000 low- fidelity simulations (i.e. no 980) using our ML-NPE method. To that end, we measure the NLPD and empirical coverage of the estimated posteriors for 980 test data. The result in Fig- ure 5 shows that ML-NPE performs better than standard NPE for both the metries. Standard NPE tends to produce overconfident posteriors, while ML-NPE yields calibrated or underconfident posteriors for most confidence levels. Thus, including low-fidelity samples using our method leads to better inference outcomes. Before concluding, we note that some recent papers demonstrating the potential of multi-fidelity SBI methods in cosmology appeared around the same time as our paper [Saoulis al., 2025, Thiele et al., 2025]. These more in-depth studies clearly highlight the potential for impact of advanced multi-fidelity SBI methods.
Empirical Coverage
ML-NPE
NPE
NLPD
0.5
Overconfidest
0.5 1
ML-NPE
NPE
Confidence Level
Figure 5: NLPD & and empirical coverage of NPE and ML-NPE for the cosmological inference task.
6
Conclusion
This paper demonstrated how reduce the cost of SBI using MLMC, but could more broadly be seen as a way to perform multilevel training of conditional density estimators in scenarios where data from different sources with different accuracy levels needs to be combined. Our method can be readily applied to scenarios with more than two fidelity levels. It is also particularly appealing since it is complementary to other compute-efficient SBI methods. For example, Tatsuoka al. [2025] recently proposed to train an NPE network on low-fidelity data, and to then use the resulting posterior approximation to guide sampling from the high-fidelity simulator. This approach could easily be combined with our method.
In terms of limitations, our method involves gradient adjustments during optimisation, and we did observe a minor increase of roughly 15%-20% in training time compared to that of standard SBI methods (see Appendix B.2). However, this is not significant issue training time is usually negligible compared to costly high- fidelity simulations. Another limitation of our approach is that is not applicable in cases where seed seed-matching of the low- and high- high-fidelity simulators is not possible. This was not an issue in any of the examples we encountered, but limit its applicability in some cases.
Acknowledgments
The authors are grateful to Sam Power, Tim Sullivan and David Warne for helpful discussions, and to the authors of Krouglova et al. [2025] for sharing their code and identifying a bug in our implementa- tion of their method in a preprint version of this paper. YH and AB were supported by the Research Council of Finland grant no. 362534. FXB was supported by the EPSRC grant [EP/Y022300/1 NJ was supported by the FRC-selected UKRI Frontier Research Grant EP/Y03015X/1 and by the Simons Collaboration on Learning the Universe.
10

---

References
Alsing, B. Wandelt, and S. Feency. Massive optimal data compression and density estimation for scalable, likelihood-fri free inference in cosmology. Monthly Notices of the Royal Astronomical Society, 477(3))2874-288) 2885, 2018. 3
H. Asi, Y. Carmon, A. Jambulapati, Y. Jin, and A. Sidford. Stochastic bias-reduced gradient methods. In Advances in Neural Information Processing Systems, volume 34, pages 10810-10822, 2021.
M. A. Beaumont. Approximate Bayesian computation in evolution and ecology. Annual Review of Ecology, Evolution, and Systematics, 41(1):379-406 2010.
M. A. Beaumont. Approximate Bayesian computation. Annual Review of Statistics and Its Application, 6(1):379-403, 2019. 1
M. Beaumont, W. Zhang, and D. J. Balding. Approximate Bayesian computation in population genetics. Genetics, (((4))2222-2035 2002. 1
M. A. Beaumont, J-M. Cornuet, J-M. Marin, and C. P. Robert. Adaptive approximate Bayesian computation. Biometrika, 96(4):983-99 2009.
Behrens and F. Dias. New computational methods in tsunami science. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 373(20553):01403322 2015. 1
A. Bharti, F-X. Briol, and T. Pedersen. A general method for calibrating stochastic radio channel models with kernels. IEEE Transactions on Antennas and Propagation, 7066):3986-4001, 2022a.
A. Bharti, L. Filstroff, and S. Kaski. Approximate Bayesian computation with domain expert in the loop. In Proceedings of the 39th International Conference Machine Learning, volume 162, pages 1893-1905, 2022b. 3
A. Bharti, M. Naslidnyk. O. Key, S. Kaski, and F-X. Briol. Optimally-weighted estimators of the maximum mean discrepancy for likelihood-free inference. In Proceedings of the 40th International Conference on Machine Learning, volume 202, pages 2289-2312, 2023. 3
A. Bharti, D. Huang, S. Kaski, and F-X. Briol. Cost-aware simulation-basee inference. In Proceedings of 28th International Conference on Artificial Intelligence and Statistics, volume 258, pages 28 28-36,7 2025. 3
C. M. Bishop. Mixture density nelworks. Aston University, 1994. 3, 29
S. G. Bobkov. Isoperimetric and analytic inequalities for log-concave probability measures. The Annals of Probability, 1903-199-1 1999. 17
J. Boelts, J-M. Lueckmann. R. Gao, and H. Macke. Flexible and efficient simulation inference for models of decision-making Elife, 11:e77220, 2022. 1,2
F. V. Bonassi and M. West. Sequential Monte Carlo with adaptive weights for approximate Bayesian computation. Bayesian Analysis, 10(1), 2015. 9
F. V. Bonassi, L. You, and M. West. Bayesian learning from marginal data in bionetwork models.
Statistical applications in genetics and molecular biology, 10(1). 2011. 9
N. Chartier, B. Wandelt, Akrami, and F. Villaescusa-Navarre CARPool: Fast, accurate computation of large-scale structure statistics by pairing costly and cheap cosmological simulations. Monthly Notices of the Royal Astronomical Society, 503(2):189--1114 2021. ISSN 13652966. 10, 30
Nicolas Chartier and Benjamin Wandelt. CARPool covariance: Fast, unbiased covariance estimation for large-scale structure observables. Monthly Notices of the Royal Astronomical Society, 509(2):2220-222) 2022. ISSN 13652966. doi: 10. 00.1093//mra///tab3007.. 10, 30
Z. Chen, M. Naslidnyk, and F-X. Briol. 1828, 2025, 4
Nested expectations with kernel quadrature.

---

K. Cranmer, J. Brehmer, and G. Louppe. The frontier of simulation basod inference. Proceedings of the National Academy of Sciences, 117(48)30055530000 2020.
P. D. Moral, A. Doucet, and A. Jasra. An adaptive sequential Monte Carlo method for approximate Bayesian computation. Statistics and Computing, 2205::009-1120, 2011.
Laurent Dinh, David Krueger, and Yoshua Bengio. Nice: Non-linear independent components estimation. arXiv preprint arXiv: 2014. 3
T. J. Dodwell, C. Ketelsen, R. Scheichl, and A. L. Teckentrup. Multilevel Markov chain Monte Carlo. SIAM Review, 61(3):509-544 2019.
C. Durkan, A. Bekasov, I. Murray, and G. Papamakarios. Neural spline flows. Advances in Neural Information Processing Systems, volume 32, pages 7511-7522, 2019. 26
C. Durkan, I. Murray, and Papamakarios. On contrastive learning for likelihood-free inference. In International Conference Machine Learning, volume 119, pages 2771-2781, 2020. 3,
M. Fujisawa and I. Sato. Multilevel Monte Carlo variational inference. Journal of Machine Learning Research, 22(278): 44, 2021. 4
M. Gatti al. Dark Energy Survey Year results: Simulation-bass cosmological inference with wavelet harmonics, scattering transforms, and moments of weak lensing mass maps. II. cosmological results. Physical Review D, 1666004633 2025. 30
T. Geffner, G. Papamakarios, and A. Mnih. Compositional score modeling for simulation-bail based inference. In Proceedings of the 40th International Conference on Machine Learning, volume 202, pages 1 -11981 2023. 3,5
M. B. Giles. Multilevel Monte Carlo path simulation. Operations Research, 56(3):607-617, 2008.
M. B. Giles. Multilevel Monte Carlo methods. Acta Numerica, 24:259-328, 2015. 2, 7,7
M. Gloeckler, M. Deistler, C. D. Weilbach, F. Wood, and J. H. Macke. All-in-one simulation based inference. In Proceedings of the 41st International Conference on Machine Learning, volume 235, pages 15735-15766, 2024. 3
M. Gloeckler, Toyota, Fukumizu, and H. Macke. Compositional simulation-based inference for time series. In The Thirteenth International Conference Learning Representations, 2025. 3
Goda, Hironaka, and T. Iwamoto. Multilevel Monte Carlo estimation of expected information gains. Stochastic Analysis and Applications, 38(4))581-600 2020. 4
Goda, T. Hironaka, W. Kitade, and A. Foster. Unbiased MLMC stochastic gradient-based
optimization of Bayesian experimental designs. SIAM Journal on Scientific Computing, 44(1): A286-A311, 2022.
D. Greenberg, M. Nonnenmacher, and J. H. Macke. Automatic posterior transformation for likelihood- free inference. In Proceedings of the 36th International Conference on Machine Learning, pages 2404-2414, 2019, 1,
A Gretton, K Borgwardt, MJ Rasch, and B Scholkopf. kernel two-sample test. Journal Machine Learning Research, 13:723-773, 2012. 9
M. U. Gutmann and J. Corander. Bayesian optimization for likelihood-free inference of simulator- based statistical models. Journal of Machine Learning Research, 17(125):1-47, 2016.
S. Heinrich. Multilevel Monte Carlo methods. In Large-Scale Scientific Computing, volume 24, pages 68-67. Springer, 2001. 4
J. Hermans, V. Begy, and G. Louppe. Likelihood-free MCMC with amortized approximate ratio estimators. In Proceedings of the 37th International Conference on Machine Learning, pages 4239 4239-4248, 2020. 3,
12

---

M. Hoppe, O. Embreus, and T. Fülóp. Dream: A fluid-kinetic framework for tokamak disruption runaway electron simulations. Computer Physics Communications, 268:108098, 2021.
Y. Hu, J. Wang, Y. Xie, A. Krause, and D. Kuhn. Contextual stochastic bilevel optimization. In seirth- Conference on Neural Information Processing Systems, 2023. 4
A. Jasra, S. Jo, D. Nott, C. Shoemaker, and R. Tempone. Multilevel Monte Carlo in approximate Bayesian computation. Stochastic Analysis and Applications, 37(3):346-36/4 2019.
A. Jasra, K. Law, and C. Suciu, Advanced Multilevel Monte Carlo Methods. Internationaa Statistical Review, 88(3):548-55 579599 2020. 2, 4
N Jeffrey and B. D. Wandelt. Solving high- -igensiinnni parameter inference: marginal posterior densities & Moment Networks. Third Workshop on Machine Learning and the Physical Sciences, NeurlPS 2020, art. arXiv:201 1.22001100000 1.05991, 2020. 3
N. Jeffrey, J. Alsing, and F. Lanusse. Likelihood free inference with neural compression of DES SV weak lensing map statistics. Monthly Notices of the Royal Astronomical Society, 501(1):954 969, 2021. 1
N. Jeffrey et al. Dark energy survey year 3 results: likelihood-free, simulation-based wCDM inference with neural compression of weak-lensing map statistics. Monthly Notices of the Royal Astronomical Society, 536(2):1003-1322, 2025. 1, 30
M. C. Kennedy and A. O' O'Hagan. Predicting the output from a complex computer code when fast approximations are available. Biometrika, 87:1-13, 2000. 2
O. Key, A. Gretton, F.-X. Briol, and Fernandez. Composite goodness boodness-of-fit tests with kernels. Journal of Machine Learning Research, 26(51):1-60, 2025. 28
D. P. Kingma and J. Ba. Adam: A method for stochastic optimization. In International Conference for Learning Representations, 2015. 25
A. Kirby, F-X. Briol, T. D. Dunstan, and T. Nishino. Data-driven modelling of turbine wake interactions and flow resistance in large wind farms. Wind Energy, 36(9):875-1011, 2023.
A. N. Krouglova, H. R. Johnson, B. Confavreux, M. Deistler, and P. J. Gonçalves. Multifidelity simulation-based inference for computationally expensive simulators. arXiv:2502.08416, 2025. 8, 10
S. Kulkarni and C. A. Moritz. Improving effectiveness of simulation- based inference in the massively parallel regime. IEEE Transactions on Parallel and Distributed Systems, 34(4): 100::1000-110 110011144 2023.
T. Kypraios, P. Neal, and D. Prangle. A tutorial introduction to Bayesian inference for stochastic epidemic models using approximate Bayesian computation. Mathematical Biosciences, 287:42-53, 2017. 1
K. Li, D. Giles, T. Karvonen, S. Guillas, and F-X. Briol. Multilevel Bayesian quadrature. In International Conference on Artificial Intelligence and Statistics, pages 1845-1868, 2023.
Liang, L. Hodgkinson, and M. W. Mahoney. Fat-tailed variational inference with anisotropic tail adaptive flows. In Proceedings of the 39th International Conference on Machine Learning, volume 162, pages 13257-13270, 2022.
Lintusaari, M. U. Gutmann, Dutta, S. Kaski, and J. Corander. Fundamentals and recent developments in approximate Bayesian computation. Systematic Biology, 66:66-82, 2017.
A. J. Z. Liu, J-S. Denain, K. Gimpel, S. Sidor, S. Levine, and P. Abbeel. Gradient surgery for multi-task learning. In Advances in Neural Information Processing Systems, volume 33, pages 5824-5836, 2020. 6
J-M. Lueckmann, P. J. Gonçalves, G. Bassetto, K. Öcal, M. Nonnenmacher, and J. H. Macke. Flexible statistical inference for mechanistie models of neural dynamics. In Advances in Neural Information Processing Systems, page 1289-1299, 2017. 3,3
13

---

J-M. Lueckmann, G. Bassetto, T. Karaletsos, and J. H. Macke. Likelihood-free inference with emulator networks. In Symposium on Advances Approximate Bayesian Inference, pages 32-53, 2019. 1,
Meeds and M. Welling. GPS-ABC: Gaussian process surrogate approximate Bayesian computation. In Proceedings of the Thirlieth Conference on Uncertainty in Artificial Intelligence, page 593-602, 2014. 3
B. K. Miller, C. Weniger, and P. Forré. Contrastive neural ratio estimation. In Advances in Neural Information Processing Systems, 2022.
M. Minenna. The detection of market abuse financial markets: A quantitative approach. Quaderni di finanza, 54, 2003. 8
Z. Niu, J. Meier, and F-X Briol. Discrepancy- inference for intractable generative models using quasi-Monte Carlo. Electronic Journal of Statistics, 17(1):411-1456, 2023.
A. B. Owen. Monte Carlo theory, methods and examples. https://artower su. /ntprtowww...uddmmmiii///////// 2013. 5,
G. Papamakarios and I. Murray. Fast e -free inference of simulation models with Bayesian conditional density estimation. In Advances in Neural Information Processing Systems, pages 1036-1044 2016. 1,
G. Papamakarios, D. Sterratt, and I. Murray. Sequential neural likelihood: Fast likelihood-free inference with autoregressive flows. In Proceedings of the Twenty-Second International Conference Artificial Intelligence and Statistics, pages 837-848, 2019, 2,
G. Papamakarios, E. Nalisnick, S. Rezende, D. J.and Mohamed, and B. Lakshminarayanan. Normal- izing flows for probabilistic modeling and inference. Journal of Machine Learning Research, 22 (57):1-64, 2021. 3
George Papamakarios, Theo Pavlakou, and lain Murray. Masked autoregressive flow for density estimation. In Advances in neural information processing systems, volume 30, 2017.
S. Paszke, A.and Gross, S. Chintala, E. Chanan, G.and Yang, Z. DeVito, Z. Lin, A. Desmaison,
L. Antiga, and A. Lerer. Automatic differentiation in pytorch. In NIPS W, 2017. 25
B. Peherstorfer, K. Willcox, and M. Gunzburger. Survey of multifidelity methods in uncertainty propagation, inference, and optimization. SIAM Review, 60(3):550 5910 2018.
D. Prangle. Lazy ABC. Statistics and Computing, 26: 2016. 26
D. Prangle. gk: An R Package for the g-and-k and Generalised g-and-h Distributions. The R Journal, 12(1):7, 2020.
T. P. Prescott and R. E. Baker. Multifidelity approximate Bayesian computation. SIAM-ASA Journal on Uncertainty Quantification, ((1):114-13, 2020.
T. P. Prescott and R. Baker. Multifidelity approximate Bayesian computation with sequential Monte Carlo parameter sampling. SIAM-ASA Journal on Uncertainty Quantification, 9(2):7888817, 2021.
T. P. Prescott, D. J. Warne, and R. E. Baker. Efficient multifidelity bkklihoodd-ree Bayesian inference with adaptive computational resource allocation. Journal of Computational Physics, 496:112577, 2024. 2
W. Press. Numerical recipes 3rd edition: The art of scientific computing. Cambridge university press, 2007. 26
L. Price, C. C. Drovandi, Lee, and Nott. Bayesian synthetic likelihood. Journal of Computational and Graphical Statistics, 27(1):1-11, 2018. 2
14

---

S. T. Radev, U. K. Mertens, A. Voss, L. Ardizzone, and U. Köthe. Bayesflow: Learning complex stochastic models with invertible neural networks. IEEE Transactions on Neural Networks and Learning Systems, 33(4): 1452-1466, 2022. 1,3
S. T. Radev, M. Schmitt, V. Pratz, U. Picchini, U. Köthe, and P-C. Bürkner. Jana: Jointly amortized neural approximation of complex bayesian models. In Uncertainty Artificial Intelligence, pages 1695-1706, 2023a. 2
S. T. Radev, M. Schmitt, L. Schumacher, L. Elsemüller, V. Pratz, Y. Schälte, U. Köthe, and P-C. Bürkner. Bayesflow: Amortized bayesian workflows with neural networks. Journal of Open Source Software, 8(89):5702. 2023b. 25
P. Ramesh, J-M. Lueckmann, J. Boelts, Á. Ceeere,,ntero, D. S. Greenberg, P. J. Goncalves, and J. H. Macke. GATSBI: Generative adversarial training for dimulation-basee inference. In International Conference on Learning Representations, 2022. 5
G. D. Rayner and H. L. MacGillivray. Numerical maximum likelihood estimation for the g-and-k and generalized g-and-h distributions. Statistics and Computing, 12(1))57-75, 2002. 25
D. Rezende and S. Mohamed. Variational inference with normalizing flows. International conference on machine learning, pages 1530-1538, 2015. 3
C. P. Robert and G. Casella. Monte Carlo Statistical Methods. Springer, 2000.
A. A. Saoulis, D. Piras, N. Jeffrey, A. Spurio-Mancini, A. M. G. Ferreira, and B. Joachimi. Transfer learning for multifidelity simulation-basee inference in cosmology. wrXiv:2505.2/215. 2025. 10
A. Saumard and J. A. Wellner. Log-concavity and strong log-concavit; A review. Statistics Surveys, 8:45-114, 2014. 7, 17, 19
M. Schmitt, D. R. Ivanova, D. Habermann, U. Köthe, P-C. Bürkner, and S. T. Radev. Leveraging self-consistency for data-efficient amortized Bayesian inference. In Proceedings of the 41st International Conference on Machine Learning, 2024a. 3
M. Schmitt, V. Pratz, U. Köthe, P-C. Bürkner, and S. Radev. Consistency models for scalable and fast simulation-based inference. Advances in Neural Information Processing Systems, 37: 126908-126945, 2024b. 3
I.. Sharmok,, J. Simons, S. Liu, and M. Beaumont. Sequential neural score estimation: Likelihood-free inference with conditional score based diffusion models. In Proceedings of the 41st International
Conference on Machine Learning, volume 235, pages 44565-44602, 2024. 3
Y. Shi and R. Cornish. On multilevel monte carlo unbiased gradient estimation for deep latent variable models. In Proceedings of The 24th International Conference on Artificial Intelligence and Statistics, volume 130, pages
2021. 4
S. A. Sisson, Y. Fan, and Mark M. Tanaka. Sequential Monte Carlo without likelihoods. Proceedings of the National Academy of Sciences,
2007. 3
C. Tatsuoka, M. Yang. D. Xiu, and G. Zhang. Multi-fidelity parameter estimation using conditional diffusion models. arXiv:2504.0184 2025. 10
A. Tejero-Cantero, J. Boelts, M. Deistler, J-M. Lueckmann, C. Durkan, P. J. Gonçalves, D. S. Greenberg, and H. Macke. sbi: toolkit for simulation-based inference. Journal of Open Source Software, 5(52):2505, 2020. 8, 9, 25
Leander Thiele, Adrian E. Bayer, and Naoya Takeishi. Simulation-efffcient cosmological inference with multi-fidelity SBI. 2025. 10
O. Thomas, R. Dutta, J. Corander, S. Kaski, and M. U. Gutmann. Likelihood-fiee inference by ratio estimation. Bayesian Analysis, 17(1):1331, 2022.
Villaescusa-Navarr al. The CAMELS Project: Cosmology and Astrophysics with Machine- learning Simulations. The Astrophysical Journal, 915(1):71, 2021. 10, 30
15

---

F. Villaescusa-Navee Navarro et al. The CAMELS Project: Public Data Release. The Astrophysical Journal Supplement Series, 265(2):54, 2023. 2, 10, 30
Virtanen et al. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods, 17:261-272, 2020, 25
D. J. Warn, R. E. Baker, and M. Simpson. Multilevel rejection sampling for approximate bayesian computation. Computational Statistics Data Analysis, 124:71-86, 2018. 2
D. J. Warne, T. P. Prescott, R. E. Baker, and M. J. Simpson. Multifidelity multilevel Monte Carlo to accelerate approximate Bayesian parameter inference for partially observed stochastic processes. Journal of Computational Physics, 469: 1143,, 2022. 2
J. B. Wildberger, M. Dax, S. Buchholz, S. R. Green, J. H. Macke, and B. Schölkopf. Flow matching for scalable simulation-basee -based inference. In Thirty-seventh Conference on Neural Information Processing Systems, 2023. 3,
S. N. Wood. Statistical inference for noisy nonlinear ecological dynamic systems. Nature, 466(7310): 1102-1104, 2010.
J. Wright and B. Recht. Optimization for Data Analysis. Cambridge University Press, 2022.
S. Yang, V. Zankin, M. Balandat, S. Scherer, K. Carlberg, N. Walton, and K. J. H. Law. Accelerating look-ahead in Bayesian optimization: Multilevel Monte Carlo is all you need. In International Conference on Machine Learning, pages 56722-56748, 2024. 4
A. Fammit-Mangion, M. Sainshury-Dake and R. Huser. Neural methods for amortized inference. Annual Review of Statistics and Its Application, 2024.
16

---

Supplemental Material
The appendix is arranged as follows: Appendix A contains the proofs of the theoretical results presented in Section Appendix B consists of the experimental details and additional results.
A Proof of theoretical results
A.1 Proof of Theorem
Proof. Note that here, and throughout the rest of this proof, we will simplify the notation. We will drop the subscript for variances and expectations, and these should all be understood as being u, U for E {l...., m and ~ x. We will also use the variable & to denote a vector containing V1am and 0, so that z E mddd Finally, we will write I and W114 to denote the spaces L'(TT x U) and W1 (T x U).
The proof will be structured as follows. We will express the variance of each term using the variance of individual samples. We will then use a Poincaré-type inequality to bound the variance of the objective in terms of the expected squared norm of its gradient, then will upper bound the norm using only terms which depend on constants and terms expressing how well G approximates G'.
Since each term in the MLMC expansion is an MC estimator (and therefore based on independent samples), can express the variances of each term follows:
Var [ho[ = Var
f (zi)
no i=1
- Var [f8 (z)
Var
(z.)
(6)
where we the fact that the variance of sum of independent random variables is the sum of variances. Similarly for the other levels,
Var [+(())]
= Var [1 1 (=) for E (1...., L.).
1 TH
For the first step. we use version of Poincaré-type inequality for log-concave measures due to [Bobkov, 1999] (note that simpler statement and additional discussion is provided in Proposition 10.1 (b) of Saumard and Wellner [2014] for the case of strongly log log-coecave measures). This result shows that for any log- log-concave measure there exists K > 0 such that for any sufficiently regular integrand f Z + where C Rddz, Var[.(=2]] < Ex-i [ [ Applying this Poincaré inequality to each term of the learning objective (i.e. to the terms in (6) and (7)), we get
Var [)r[)](=)] S Kroin
(8)
Var (f(()) f/- (x) S K Poin
- -
(x)
(9)
for 1 E f1,.... D}. Here, we emphasise again that the expectation is over 2, which encompasses both 0 and U1mm- and the vector Vaf((2) E Rduu+da for any E {1, {1...., L}. This result requires the joint distribution of the prior æ and times the base measure UJ to be log-concave, which holds since the product of log-concave densities also log-coneave (since the sum of convex functions is convex) and we have assumed that the prior and base measures are independent and separately log-concave through Assumption (AI).
To simplify notation, we now introduce the vector-valued function g'(2) = g'(0, Upm) = (u,) ,..., G% (Wm))T TT so that d (x1mmm 0) qo (0,, = (z))) We will now derive first bound, which looks at the first term in the MLMC expansion. To do we simplify
17

---

(8) as follows
f|1 = ||V: log ás (q
E
|| log (0 (z) [
E | [2 (2) log do (9 (z)) [
(10)
S E
| 16 (=) (=)
log T0 (g (=)
(11)
< E [ [
E | do (o ()))))
(12)
Here, we have that (10) follows due to the chain rule, (11) follows due to the definition of the matrix 2-norm, (12) follows due to the Cauchy- Schwwaz inequality of expectations.
We now bound each term in (12) separately. For the first term, we get
V
(2)
V 4:
(2)
Vog"(2)
(13)
i=1
m
m
V
12,))11 +
V (u,) 3 IV
(14)
Σ
(1g)) +
v (u, +
(15)
m
du
m
de
Mgh G2
(2y))
Vo, G)(u,)1 de
(16)
==
where (13) and (14) follow from the fact that the two norm squared of matrix is less than the sum of the two norm squared of sub-matrices constructed through rows and columns, (15) follows by noticing that ||V de and V/V,, G2t whenever i # j. and (16) follows similarly to (13) and (14). Taking squares and an expectation, we get:
E [[| 20 [2)
27
m de
<
E
I n.J Ga +
Vo ((g))11 de
(17)
k=
< (mdd + mde 1)F
uea G3(
en Go (u,)) + de
(18)
< (mdu + mde + 11 (m
Ii + d2)
(19)
< K grad (110 +1)
(20)
where (18) follows from (Enn a;) <n En_ norm and the fact that 81 .i....; Uyn have the same distribution and hence the same expectation, and (20) follows by grouping constants together.
= az. (19) follows from the definition of the Sobolev
18

---

We now move on to bounding the second term in (12). Since we assumed in Assumption (A3) that v log 9s satisfies a linear growth condition, we must have that
( Kgow (Ø)
(KSOn
m Σ
E V (2) SE
(ug(2))
110112 +
(21)
m
(m + 2) K grow (+))
Σ
[G)( + 11011
(22)
(m + 2) K grow (0)* (melG + [101 1)
(23)
kooon (0) [Iwaa]+
(24)
Here, (21) uses the growth condition, (22) holds by applying (2) a;)² < n2 Dm of (23) follows
from the definition of Sobolev norm, and (24) uses the fact that E [1||]]] is upper bounded by constant since the expectation against 9, which log-concave distribution and hence all its moments are finite (see Section 5.1. of Saumard and Wellner [2014]).
Combining (6), (8), (12), (20), and (24) therefore gives:
Var [o)] - - Kpoin grad (IG*))ima + 1) Kocorc (119)) + +1)²
1
no
Ko(Ø) (IG" +1)
(25)
where Ko(Ø) is used combine all of the constants. This now concludes the first part of our results, which bounds the variance of the first term in the MLMC telescoping sum.
We can now derive a similar bound for the other terms (i.e. to simplify (9)). We do this by first only considering the norm inside the expectation:
v.f((=) -V.f
(x)
(z)
V log do (g'(2)) (V, log (g'
V , log ās (g) V. log (g (2))l1
V V (2) lo ( (z)) V (z) log ds (g (z) 11, (26)
V (=)'(=) log To (gf (2)) v log (z))) 1| (z))
(Vig((>) (2) v V v
log To
(27)
V (= v 1og ( (z))
v is {2))
V V (z)
v log ( (z))
(28)
Here, (26) follows from the chain rule, (27) follows by adding and subtracting the term (2) log q (1' and using the triangle inequality, and (28) follows from the Cauchy- Schwarz inequality of the 2-norm. Squaring both sides of this inequality and taking expectations,
19

---

then obtain:
-=fd())-Voftt -
v=)
(2g((2) V (2) v Vlog (o' (z) 1/)
E 2 V
(=)) V Vlog (z)) |
(o (x)) 1
( - V
(29)
V (=) [ T [
2B
v log is (z))
2E -
veg' (x)
v log do (z))
(30)
where (29) follows from that fact that for a, E R, we have (a + b)2 < 2a² + 25², and (30) follows from the Schwarz inequality for expectations. To conclude this proof, notice that the derivation from (13) to (16) can be modified by replacing Go by G gives:
E Vig'
< Kgra (1o'"liu++ +1) s Kgnu (Sf + +1.
(31)
where the last inequality holds thanks to Assumption (A2). Similarly, replacing G by G - Gl-1 and following the derivations from (13) to (20) gives
Gl-1
E (ig
V V=g
(=)]2
S K grad ||G' G II
(32)
We notice that we lose the additive term since the last columns of the matrix Vag((z) form a zero matrix. This is because
Vzg'-1(=)
V_g' (2) (x)
(33)
V.Gu(un) V.G. 9 (1)) ... V.Gi
V.G. (um)
V.O V.0
(u1) V.G.
(w))
... VoGg(um))
VoG (um)
V V 0
V.G.(u)) V.G. (mn) V.G. V.Ga (um)
0
(un) V.G. (1n) VoG (um)
VoG (um)
The non-zero lower-right block, i.e., Va0 in (14), the cause of the additive term in (20) and the final result (the upper-right block is always zero since v.0 0), which now cancelled out. Additionally, we could replace G by G!-1 in the the bound from (22) to (24) in order to get:
E
V log ao (g' (x))
< Kscone (x) (I ll +1)
(34)
< Kscore (0) (St + 1)
(35)
where once again used Assumption (A2).
20

---

Finally, we split the following expression to make use of our local Lipschitz property:
E V log ds (g'(z)) V log do (9
(z)) (x))
(36)
= E || log ds (z)) - Vlogū de (g' (=)) |1 o
2 $]
x IP
+ E
V log is (z)) v log de ( ( |
<o]
- 9 (z)112 < 6)
x IP
(37)
log do (9 (z)))) |1 + V ds (d (3)1 [1ld(() - -
+ (rip(o) 'E |1o'(
- o (=) x1
(38)
K wore (0) + Stt + 2) + Kip(o)*) E (2) -==
(39)
K soore aor (0.)) - G²-1 - I
(40)
Here, (37) follows due to the law of total expectation. For (38). the bound on the first term follows due to ||a - bl2 < 8([|a]12 + []b] and Markov's inequality, and the bound on the second term follows due to the local-Lipschitz condition (i.e. Assumption (A3)) and the fact that probability is always upper bounded by 1. Then, (39) follows using (35). Finally, (40) follows by grouping all constants together and noting that
E (1))(()) -
=E
(ui) a² (w:
(41)
S m
(u.) - G (u:) S ml|GG G Wl
(42)
Combining all of the above (i.e. (7), (9), (28), (31), (32), (35), and (39)), we end up with
Var [hr[o)] - Koin
2
K grnd (S) + 1) K score- (0, (0. G² -i- Wr
n4
+ K 3/2
(1|G* - -0-'--11 1
(0) (st + )
grad
Ki(0)
G G1-1112
W1.4
74
where Ki(Ø) combines all constants. This proves our second result and therefore concludes our proof.
A.2 Proof of Theorem
Proof. We first recall both the cost and the variance of our estimator, modifying our notation slightly emphasise the number of samples no,..., The total cost of this method is given by
Cost(Emr (0): no. ..., nn.) Co +
Σ mc(Ci + Cr-1). 1=1
where C is the cost of evaluating fl at level 1. In addition, we also have the following upper bound on the variance using Theorem 1:
Var 1MMMM (#); no....,... nz] < Ko(#)
(| +
Ki(
||'
21

---

where once again we write || I|wia to denote the norm solve the following optimisation problem:
|| (lx)) Overall, we would like to
:.n..... == arg min (11........) such that Cost(EM)Mc(0)) no. nr. S Chadgets
Var (0); no nr.
where Chudges is the computational budget. We relax the problem slightly by minimising the upper bound on the variance instead, thus the problem can be expressed in the following Lagrangian form:
_no.....,nr n,, v)
= Ko(Ø)
IIG +1) +
Σ 1=1
L
Ki(Ø)
G G
W1.4
no
+ V (Cost(EMLMc(()) ,o,,,, Cbungea)
Ka(Φ)
(16 (1liiaa1))))) +1) +
K((Φ)
||G - G/-1112 W/1,
Yo
n
nn Σ
mi(C + C--1) Chuudee
This can be solved by setting all partial derivatives with respect to (no...., n1) and v equal to zero and solving the associated system of equations. Firstly, we get:
aC(no..... V..V) = (o( (IG"lI + +1) no + vCa =0
(43)
Ono
Ko(Ø)
(I69))ma+ +1)
n/ =
(44)
vCa
and for 1 (1...., L}, we have:
...(oo,...,,,..... -K(ld)
2 Il ni 2 +v((+++++))
Oni
Ki(O)
ni
v(C C1-1) |G' |G' G/-
(45)
Finally, taking the partial derivative with respect to confirms that our constraint is active (i.e. we on the boundary of the feasible region):
aC(no.... OL,,o...., no +
L.
Σ
nu mC Ci-1)
Chudget
dv
Coudget noCo + Σ +
(1(C + CC-)).
(46)
1=1
22

---

Plugging the results of (44) and (45) into (46), we get:
Chudget Do +
(C + Ct-1)
1=1
Ko(Ø)
(IG + +1) x Co
VCD
Ki(Ø)
C(-+ 16' - GG' (CI + C--1)
Ko(o)Co (10%)/oza +1)
(G -G
which gives
( Kkoto)Co (I6") + +1) + Ele +
2
G'-11|
(47)
budget
We can
the expression for that we obtained in (47) obtain simplified expression for Th (using (44) and (45)):
no
Chudget
IIGIIiia +1,
Chudget
|| 1
Iww I for orlle E (1,..., L..
ng
Ti.
VC,
VC++
This completes our proof.
A.3 Extension of Theorem 1 to the gradient
We provide an upper bound on the variance for each element of the gradient V VoLMLMC(O). that is, on each partial derivative Ve, (Ø) for E f dø}. The partial derivatives are given by
V By (d)M)(0) = Vo, ho()) +
V o, (1())
1
no
n1
(o
Ve 01)
n4
Theorem 3. Let Φ E C Rd and suppose the following assumption hold in addition to A1-A2 in theorem 1;
(A3') Ve, log qo is continuously differentiable, locally Kipp (m)-shhoth and satisfies the growth condition |[V V log (m 0)l12 I K (d) m (L))()), KJ groe (b) for all j 1,...,d..
||0/|22 + ||0122 + 1) for some
Then, for E (1....,.. L}. (o) Kg(() > -.. (..e,........................... ......,d)) independent of nu, n. we have that:
Var Voghho(o) v Ki(0)
+ +1),
Ki (e)
and Var [()) <
||G' - Gl -1
23

---

Proof. The variance of each term can be expressed as
Var V b [o[oo] Var V (o
(48)
Var V VV
(49)
Assume that Vo, So is sufficiently regular, applying Poincaré inequality to (48) and (49) (as in (8)
and (9)) gives:
Var Var[oo, Vo, & Kunin E
V V fo
(50)
Var ( (y V L5 (z) K Poin
V VoVoddd(>) - V. Vo
(51)
where the expectation over z. We can simplify (50) as
V 18
E
V log (o")
(52)
=E
V 1
(2)
log to (
(53)
9 22)
VV log (o'
(54)
Here (53) is due to the chain rule and (54) follows (11) - (12). The bound for the first expectation is given by (20). For the second expectation, have:
E VV log ( < Ki (Ø) (I +
(55)
by Assumption (A3') and following (21) (24). Combining (48), (50), (54), (20), and (55) gives:
Kg(0)
(IG") +1)
Var 6, ho())] I
(56)
no
where (o) a constant depending on and j, independent of no- Now, we derive an upper bound on Var [V (o)) Following (26) (30), have:
V Vo V.V Vo, (z) E
(57)
E || (+)11)) & (B VV log (o'())) 2
VV ay log Te (=))
E [=) v v=00'( e || Øg log io $
-
(58)
Bound for the first expectation given by (31) and the third expectation given by (32). For the forth expectation, from (55) (replacing 0 with 1) and Assumption (A2), we have:
E VV Φ3 log go (g' (z)) < Kimme (1) (S + +1)
(59)
24

---

For the second expectation, following (36)-(40), and using our local Lipschitz property, we get:
|| o, log ig (z)) VV % log de (o (=))]1)
(60)
K acooo-din 6) || - G w/1
(61)
where K2 score (6) combines all the constant, independent of ny and the difference in the simulators. Finally, applying the bounds for the each expectation; (31), (61), (32), and (59), to (58), and combining with (49) gives:
Poin [K] Agaa (St (St + 1) (K3 scor-diff (0)) || G' li
Var oy hu(o))
n4
+ K grad
(IG'- - -C' -l- (Kamm (o)) ($)) 1)
(62)
KP (0)
GIt
(63)
n4
where Kr combines all the constant. Combining (56) and (63) gives bound on the variance of the partial derivative:
(d)
k((0) ||'
Var dy (c)) S
+1)
(ir (64)
n4
1=1
B Experimental details & additional results
For all the experiments, we use Mae M4 CPU with 16 GB memory for the training of neural networks. Running all the experiments roughly takes half a day. For optimisation, we batch gradient descent with the Adam optimiser [Kingma and Ba, 2015] using Pytorch [Paszke et al., 2017]. For the construction of the conditional density estimator, et al., 2020]. We also use the BayesFlow package [Radev et al., 2023b] for some of the figures. For each experiment, we fix the number of epochs unless stated otherwise. We use the same setting and stopping criterion for the NPE and the NLE baseline as in the default implementation of the SBI package. We set the learning rate to 10 4 for all the experiments.
the SBI package [Tejero-Cantero
B.1 The g-and-k distribution
Simulator setup. The g g-and-k distribution can be written in simulator form as follows:
Go(u) = 01 + 0g 1 0.8 (1 --epp-003s) (xp(-03z(w)) ) (1 2(1)2) log((.) 2(u),
1 exp((Ab2(u))
2(u) = (u) = /2erf (2u - 1),
u ~~ Unif([0, 1]),
where 4-1 is the quantile function of the standard normal distribution and erf the inverse of enror function. We set the prior distribution to be tensor product of marginal distributions on each parameter: 01, f22 8g- ~ Unif([0. [0, 3) and 0₂ ~ Unif([O, ex.(). Note that we always resort to an approximation method for the evaluation of erf eof (.). For the high-fidelity simulator, we use an accurate approximation implemented in Scipy [Virtanen et al., 2020). For the low-fidelity simulator,
we use Taylor expansion of erf up to the third order:
Zlow (u) = /2ert (2u - 1), erf =
π
u
2
12
An advantage of the g-and-k distribution is that its density can be approximated numerically almost exactly. Following Rayner and MacGillivray [2002], we first numerically obtain Fe(x) G (x)
25

---

by solving numerically for 2; Ga(()) - 0 using root-solving algorithm [Press, 2007]', then we obtain the density function by taking p0x 0) Fz = 0Ff(z)/Oz. For this second step. typically use finite difference approximation. Seed-matching is simply done by using same the and 0.
Neural network details
NLE: We the neural spline flow (NSF) [Durkan et al., 2019]. We pick 10 bins, span of [-7, 7] and 1 coupling layer since we only have one dimensional input. The conditioner for the NSF is multilayer perceptron neural network (MLP) with hidden layers of 50 units, and 10% dropout, trained for 10, 000 epochs.
NPE:
We use NSF with bins, span of [-3,3] and 3 coupling layers. The conditioner for the NSF is a MLP with 2 hidden layers of 50 units, and 10% dropout, trained for 800 epochs. For each true parameter values, we produce m 1000 iid samples and obtain quantile based 4-dimensional summary statistics following Prangle [2016].
Evaluations details
NLE: We calculate the forward KL-divergence between the almost exact density and the approxi- mated density over 2000 equidistant points in [-30, 30).
NPE: We calculate NLPD over 500 simulations. Empirical coverage is estimated using average value 500 simulated datasets, where we draw 2000 posterior samples from each simulations. We then calculate - 3-highest posterior density credible interval, where 3 is 101 equidistant points between and 1.
B.2 Training time comparison
We report training time per epoch for multilevel versions of NLE and NPE and their standard counterparts for the g-and-k experiment; see Appendix B.2. For both methods, we picked n = 2000 for the standard NLE/NPE with MC loss and no 1000, nı 500 for ML-NLE/ML-NPPE with MLMC loss such that the total number samples are the same. Each network is trained independently 10 times assess uncertainty, with training budget of 500 epochs for NLE and 100 epochs for NPE.
Method
Training Time per Epoch (s x x10
ML-NLE NLE
1.82 (+0.05)
1.58 (+0.06)
ML-NPE NPE
8.04 (#0.17)
6.56 (40.14)
Table 1: Average training time per epoch (standard deviation in gray).
B.3
Ornstein-Uhlenbeed process
Simulator setup. Given T 10 (total time), At 0.1 (time step), 2o = 2.0 (initial value), 0 = 17, a] N = [ 100 (number of steps), the high- and the low-fidelity OUP simulators defined as follows:
High-fidelity:
For / O....N-
Je+l = Iz + AT,
Ar, = 7/1 10--x) x.) + ot, u, ~(0,1
Note that since Go(u) (u) is quantile function, Fo(x) G7 (x) is cumulative distribution function.
26

---

MI-NPE (100) PC T
Estimand
BE
BC SE
(Almest) ense
High dideity
0.2
Lne fideleyy
1SE
MI-NPE (50)
10P
R 0,50
BC SILE
0.0
ML-NPE (1O)
Estimatdd
PC Rile
F-BR22
2% BAL
(a) Data
(b) ISE comparison
0,50
(higho
0:22
Emirard
PLM
0.4
0.50
NPE (low)
Emmmadd
LL0
0.00
Trae
Trx
Truo
Tim R4.
(c) Density comparison
(d) NPE performance
Figure 6: Additional results for the g-and-k experiment. (a) Histogram of 1000 seed-matched samples from low- and high- fideliicyiity simulator. (b) ISE (integrated squared error) for NLE (L): sum of the squared distance between (almost) exact density and approximated density over 2000 points in the interval [-30, 30). The performance of ML-NLE improves nı increases. (c) Approximated density and (almost) exact density for NL.E across four simulations. The first row of (c) shows that we get better approximation of the density when combining low and fidelity samples than using them separately. The second row of (c) shows the improvement of the performance of our ML-NLE as nı increases. These results demonstrate the effectiveness of our method. It measures how well the ground truth value is captured by the median of approximated posterior. The x-axis shows the ground truth parameter value and the y-axis shows the median of the posterior distribution with the median absolute deviation, i.e. +ediinn(1 points. Here r and R² denotes the correlation coefficient (f) and the coefficient of determination (f) between the ground truth values and the estimated median, respectively. The dashed diagonal line indicates perfect recovery of the true parameter. Using MI.MC leads significant improvements in the parameter recovery, which become more pronounced as increases.
(d) Recovery plot for NPE:
median(0)). shown around the
Low-fidelity:
F1.N = u +
u1:N ~ 1))
We set the prior distribution ~ Unif( (0i,( x [0.1, 3.0j x [0.1,0.6) The resulting is 100-dimensional time series. Seed-matching simply done by using the same u and 6.
Neural network details. We use NSF with 2 bins, span of [_ and 2 coupling layers. The conditioner for the NSF is MLP with hidden layer of 32 units, and 10% dropout. For ML-NPE, we trained the network for 500 epochs. For TL-NPE, we set high maximum number of epochs and used on early stopping. The stopping criterion is based on the validation loss computed on 20% of the data, with patience parameter controlling how many epochs to wait before stopping. For each true parameter values, produce statistics. The subsamples are taken at equal interval in log space such that we take 0, 3, 10, 31, 99th data points.
1 sample and use 5 representative data points a summary
Evaluations details. We train each MI-NPP and TL-NPE 20 times. For each trained network, we compute the NLPD and KL divergence against reference posterior across 500 simulated datasets.
27

---

NPE rePfrroor
10000)
JC
M
Es Baiimuod
Indes
01
10
0.1
06
(a)
(b)
Istimatdd
Entinsed
(d)
ML- (nu
1000
TL-PE (oy F 0,07
1003
Estimateed
0227
Exiinted
J 0.53
0.4
(e)
Figure 7: Additional results for the OUP experiment. (a) Example of one seed-matched sample from high and low- fow-fidelity simulators. (b) Recovery plot for the reference NPE posterior with n = 10 (f)-(f) Recovery plots for ML-NPE and TL-NPE for nı 10 and n1 100.
101
NL2D
-500
KLDD
200
-75
ML-NPE TL-PPE
-00
(a) vi
(b) ml
Figure 8: Figure 8 with all the data, without removing outliers. Outliers identified using IQR method with factor 4, which removes only extreme outliers.
B.3.1 Experiment with different dimension of
We now conduct an experiment explore a failure mode of our method. We include an additional parameter in the high-fidelity simulator, termed "initial value" 04 zo uniform((0, 4)). instead of fixing it at zip 2.0 as we did previously. For this experiment, we picked no 1, 000 and M1 100. We additionally trained the conditional density estimator using only 100 data from the figh-iidlityy simulator. All the other settings remain the same.
Table 2: Mean and standard deviation of NLPD and KLD.
ML-NPE NPE (high only)
NLPD
-0.18
-1.21
(0.04)
(0.53)
We observe that it notably underperforms in NLPD compared to using only high-fidelity data. This may be due to our training objective being more unstable and therefore more sensitive large differences in simulator outputs introduced by the additional parameter
B.4 Toggle switch model
Simulator setup. We follow the implementation by Key al. [2025]. Given parameters 0
[on, (22, 31 L2, p1 o, 7] T we can sample from the simulator by Npp + ur, no paa//////// /u??), where N denotes the truncated normal distribution R,, Here, ur is given by the discretised-tin time equation
28

---

0006
MILNEE NER-Digh)
0109
0.002
(BOK)
0001
KURO
0.04
0.002
0100
200
1000 3:30
730
5000 1.400
20
623
(a)
(b)
Figure 9: Additional results for the Toggle-switch experiment. (a) Histograms of 1000 cetdmmachhed samples from the high-, the medium-, and the low-fidelity simulators. Note that the difference between the low- and the medium-fidelll simulator larger than the difference between the medium- and the high-fidelity simulator. This suggests that should take n1 n22. (b) Approximated density and samples from high-fidelity simulator across four simulations. ML-NLE demonstrates superior performance compared to the NLE baselines trained exclusively high-, medium-, and low-fidelity data with the same total simulation budget. In particular, the density estimated by ML-NLE aligns more closely with high-fidelity simulations, highlighting the improved accuracy of method.
as follows:
Fort = -O.... 1:
u++1 (H 0.5),
Pu,t = 2t +
a1
(1 +1+ 0.03u:)
(1 (1+
U++1 NHHo,, 0.5), Mv,t = +
0/2
-11 + 0.03c;)
(1 + u4
We set the initial state uo [0.01, 5] x [0.01, 5] x [250, 450) x [0.01, 0.5] [0.01, 0.4)). The cost of the simulator increases linearly with I", and T + co leads to an exact simulation. Different choice of leads to the simulator with different fidelity levels, having different cost and precision. The expression for the cost of simulating data for the MLMC and MC is given by
= 10 and prior distribution 0 ~ Unif([0.01, 50] X [0.01, 50) x
Σ mi(i +T--i) 1=1
CMI-NLE
noTo
ONLE cTn.
This yields a total data generation cost for ML-NLE with (no, n1, n2) = (104, 500, 100), and (To, Ti, T2) (50, 80, 300) CMLMC 603000, assuming unit cost 1. We then allocate the same computational budget to NLE at different fidelities. This results in the following sample sizes: 2010, n1 CML.NLE/Ti 7537, and no CMI-NLE/
n2 = DCCL-LLL//T
12060.
Note that / varies with the fidelity level can complicate seed-matching across fidelities. To address this, we define unified domain for the noise 44 = - 2/ shared across all fidelity levels, and assume that for lower fidelities (smaller I), certain components of the input are simply disregarded. For matchin, we share and Uy O u1-1 1 where uy U C U, Vl.
simulator, as is given by 2 This discrepancy
Neural network details. use Gaussian mixture density network [Bishop, 1994] with two mixture components. We use MLP with hidden layers and 20 units to estimate the parameters: means, variances, and mixture weights, trained for 10, 000 epochs.
Evaluations details. densities, and calculate MMD over 5000 simulations. The length scale estimated by the median heuristic.
We sample 500 from both high-fidelity simulator and each approximated
29

---

B.4.1 Experiment with different allocation of samples n4
We now repeat the same experiment, but with different allocations of Ro. R1, n12, keeping the total computational cost roughly the same. Apart from the choice of ni 10, 000, n 500, n2 100 that we use in Section 5.3 (option A) which is guided by our theory, we include two other alternatives: (B) Ro 9260, 7i = 200, n2 = 300 and (C) Ri: R1, 122 = 1077. Option B places more budget on learning the difference between the medium- and the high-fidelity simulator (opposite of option A), and option C allocates equal number of samples for all the terms. The results in Table 3 indicate that option A is the best performing, indicating performance can be improved by careful assignment of the computational budget led by insights from our theory.
Table 3:
Mean and standard deviation of MMD. The results for option A, only high, medium, and low remains the same in Section 5.3.
Option A
Option B
Option C
only high
only medium
only low
MMD (
0.16
0.33
0.29 (0.28)
0.43
0.37 0.59 (0.65)
(0.23)
(0.26)
(0.29)
(0.45)
B.5
Cosmological simulations
Simulator setup. We use the CAMELS simulation suite [Villaescusa-Nant Navarro et al., 2021, 2023], benchmark dataset for machine learning in astrophysics, to study multi-fidelity simulation- bbssd inference in cosmology. The simulation is one of the most expensive cosmological suites ever run; a small fraction of the next generation are being run the UK DIRAC HPC Facility with 15M CPUh.
The dataset includes both low-fidelity (gravity-only N-body) and high-fidelity (hydrodynamic) simu- lations of 25 Mpc/h / cosmological volumes. The original inference task involves two cosmological parameters, 0- (Iem, rs), where Ω. the matter density and og the amplitude of fluctuations, with mock power spectrum measurements P(k) used to infer the parameters. However, due to the limited availability of both high- and low-fidelity simulations, we focus on inferring only dg and treat Ωum as part of the nuisance parameter. We found that attempting to jointly infer both parameters led to poor performance (expected due to physical ag- degeneracy) even when using 90% of the high- aig-liielitty data.
In this setup, low- owwfiidlitty simulations governed solely by 0 and the initial condition seed while high-fidelitt simulations additionally incorporate complex astrophysical processes (c.g., feedback), modelled via four extra parameters. These high-fidelity simulations are significantly more compu- tationally expensive- often more than 100x slower to generate than low-fidelity ones. Given these limitations, cosmological analyses often rely on conservative data cuts exclude small-scale modes (e.g., 2 0.1, h/ Mpc), where simulation inaccuracies are most pronounced [e.g. Jeffrey et al., 2025, Gatti et al., 2025]. Multi-fidelity approaches aim mitigate such constraints by leveraging both inexpensive and high-fidelity simulations to improve the accuracy of cosmological inference.
Note that the idea of combining low- and high- fidelity simulations has previously been proposed in cosmology; see for example the work of Chartier et al. [2021], Chartier and Wandelt [2022] who use low-fidelity simulations to construct approximations of quantities of interest which can be used as control variates. This work differs from our proposed approach in that it targets quantities such as means and covariances, rather than f the training objective of neural SBI.
Neural network details. We NSF with 3 bins, span of [-3,3] and 3 coupling layers. The conditioner for the NSF is an MI.P with 2 hidden layer of 30 units, and 10% dropout, trained for 400 epochs.
Evaluations details. We use the same setting as the g-and-k experiment for NPE except that we use 980 test simulations for evaluation.
B.6 Ablation study of the gradient adjustment technique
We conduct experiments to evaluate the effect of the different gradient adjustment techniques we employ on the performance of our method. We train both MI-NPPE and MI-NLEE using (i) our gradient adjustment approach which involves both rescaling and projection, (ii) only rescaling, (iii)
30

---

[P03) (.66 15tps]7
80
ML-NPE
NPE
1.0
Jf
= 0.20
1.0
R
0,.00
Fsstiated
0..37
Estimated
High lidelity
60
fidelity
0.8
0.8
lidelity
Low libelity noise
0.6
0.6
0.6
0.8
1.0
0.6
0.8
1.0
8 (h Mpc
03g
(b)
(c)
Figure 10: Additional results for the cosmological simulator experiment. (a) seed-matched sample from high and low-fidelity simulators. The one with noise used train the neural networks. (b)-(c) Recovery plot of ML-NPE and NPE. Adding low-fidelity samples lead better recovery of the parameter.
only projection, and (iv) no gradient adjustment (standard training). We then compare the performance of NPE and NLE on the g-and-k and NLE on the toggle switch experiment. Other than the choice of the gradient adjustment, the training method, hyperparameters, and evaluation methods remain the same. The results are shown in Table 4.
Table Mean and standard deviation of the metrics. For g-and-k NPE, no 1000, ni 100,
1000 and for g-and-k NLE, no 10 7, 300, m as in Section 5.1. 1. For toggle switch NLE, no 10000, 1 500, 722 100 Section 5.3.
(i) both
(ii) only rescaling
(iii) only projection
(iv) standard
g-and-k NPE (NLPD L)
-0.30
-0.21
-0.11
-0.13
(0.31)
(0.25)
(0.11)
(0.27)
g-and-k NLE (KLD L)
0.22
0.21
0.24
0.24
(0.47)
(0.45)
(0.46)
(0.28)
Toggle-switch NLE (MMD 4)
0.26
0.50
0.35
0.59
(0.25)
(0.37)
(0.27)
(0.31)
We observe that when the MLMC loss diverges under standard training, as is the case for the toggle switch experiment and g-and-k with NPE, our gradient adjustment approach that combines both gradient rescaling and projection yields the best results. In the case of NLE on the g-and-k simulator, the loss does not diverge and standard training is sufficient. In such cases, the gradient adjustment yields similar results as standard training, albeit with an increase in the variance of the metric. Therefore, we suggest optimising using our gradient adjustment approach when using ML-NLE or ML-NPE avoids need first detect whether the loss diverging or not.
We further compare the training losses with and without gradient adjustment. With the gradient adjustment, the loss diverges, whereas with it, the loss curve exhibits convergence; Figure 11. To clarify the reason for this behaviour, we also provide an illustration of our gradient scaling procedure in Figure 12.
We that the optimisation requires some form of regularisation, and the one used one possible choice. In our experiments, the gradient adjustment approach performed the best among the options tried (e.g., regularising the contribution the difference term when it dominates the first term, or penalising high variance in the difference term).
31

---

hu(o)
1.0
0.5
0.0
1000
1000
1000 Epoch
1000
Epoch
Epoch
Epoch
(a) Loss without gradient adjustment
(
(u(o)
+.+
S
0,-
1.5
1.0
0.5
0.0
1000
Epoch
Epoch
Epuch
Epoch
Loss with gradient adjustment
Figure 11: Comparison of training losses with and without gradient adjustment. (a) Without gradient adjustment, the contribution of 20 begins to dominate the loss around epoch 1000. The resulting conflict between gradient components leads to unstable optimisation, as evidenced by strong fluctua- tions in the loss and eventual divergence. (b) With gradient adjustment, all components contribute more stably, and the overall loss decreases steadily eventually, indicating convergence. Loss functions shown for the g-and-k experiment with NLE.
hu
ho
0.
So
6
(a)
(b)
(c)
Figure 12: Illustration of the gradient scaling with two levels. The coloured arrows indicate gradient of the correction term, total gradient, and scaled gradient respectively. (a) Ideal case: ||V (o)) remains small and works as a correction to Vaho(o). (b) In the later stage of training, Vahu(o) and
<) c.+ diminishes and 0. starts to dominate the optimisation. (c) Gradient scale adjustment: We scale
such that ||5, ||c 11c--1 and ||V (o) remains small throughout the training as intended.
32

---

NeurIPS Paper Checklist
1. Claims
Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Question:
Answer: [Yes]
Justification: In Section 4, we provide theoretical analysis of our method, and in Section 5, we empirically show sample efficiency of our method compared to standard NPE/NLE and the related work.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: Limitations of our method are discussed in Section 6.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Answer: [Yes]
Justification: Assumptions and theoretical results are provided in Section 4 and complete proof is provided in Appendix A.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yss]
Justification:
The loss function is provided in (2), training method is provided in Section 3.
All the experimental details are provided Appendix
5. Open access to data and code
Question: Does the paper provide open access the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: Detailed experimental settings, including all hyperparameters, are provided in Appendix B. Reproducible code, along with the simulator and evaluation data, is included in the submission.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: All the hyperparameters for the training and the evaluation methods are provided in Appendix B.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: We report performance as box and violin plots over multiple runs and/or multiple independently trained networks in Section 5, All the details are in Appendix B.
8. Experiments compute resources
33

---

For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification: is provided the first paragraph of Appendix B.
9. Code of ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics https http:p:///nnuriie http//eerrippp.eeeeeee https///ccppppcce/ccccccccceeeea
Answer: [Yes]
Justification: All data used in this paper is synthetically generated through simulation, posing no ethical or privacy concerns.
Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: We discuss the potential positive societal impacts of our method in Section 1, and not aware of any significant negative impacts.
11. Safeguards
Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (c.g., pretrained language models, image generators, scraped datasets)?
Answer: [NA]
Justification: We do not use data or models with high risk for misuse.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g.. code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: All software used in this work is open-source and has been properly cited.
13. New assets
Question: Are new assets introduced in the paper well documented and the documentation provided alongside the assets?
Answer: [Yes]
Justification: Reproducible code is provided, along with instructions to replicate all experi- ments presented in the paper.
14. Crowdsoureing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification: The paper does not involve any such experiments.
15. Institutional review board (IRB) approvals equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
34

---

Justification: The paper does not involve any such experiments.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, originality of the research, declaration is not required.
Answer: [NA]
Justification: The core method development in this research does not involve LLMs as any important, original, non-standard component.
35


Reviews:
[{"review_id": "WEyBjzjGKS", "reviewer": "Reviewer_kZeL", "summary": "This paper extends existing amortized inference methods with a multilevel Monte Carlo (MLMC) estimator. It scopes its contributions through the lens of simulator based inference (SBI), and integrates MLMC both in neural likelihood estimators (NLEs) and neural posterior estimators (NPEs) to improve on sampling budgets for expensive high fidelity simulators. It achieves this by levering cheaper lower fidelity simulators to draw more training samples for the neural artifacts given the same computational budget. Besides introducing the two more efficient gradient estimators it introduces heuristics for allocation simulation efforts across fidelities, and covers a critical set of gradient surgery techniques to stabilize training. It also provides a theoretical convergence analysis of the introduces estimators. The method is empirically validated on a range of SBI benchmarks: g-and-k, Ornstein-Ornstein–Uhlenbeck, a toggle-switch model, and the CAMELs cosmology dataset.", "questions": "1. Why was the paper not targeted on an expensive real world simulator experiment? The papers establishing the Toggle-Switch model (sec. 5.3) and Cosmological Simulation (sec 5.4) are from the SBI machine learning literature and samples can be simulated sufficiently quickly even with high fidelity in them. While they can be made more expensive it is less clear whether the more subtle assumptions about alignment between fidelities in the paper translate into a real world setting, and it would have helped to highlight those. The biggest gains would be demonstratable with larger scale simulations from respective sciences, e.g. a full cell simulator or a more expensive astrophysics simulator, in which single samples can take hours or days to be computed. The SBI literature often suffers from picking very simplified simulator experiments that are of little practical interest, yet this paper aims to address one of the biggest practical limitations of SBI based on amortized inference. It would have been sufficient to cover such an experiment and two smaller ones such as g-and-k and maybe the OU-process.\n\n2. Similarly, why is the evaluation in 5.4 for the Cosmological experiment not more convincing? Figure 5 left shows very high variance of ML-NPE and the empirical coverage (right) is very hard to interpret without understanding the specific simulator well. I would expect a more thorough evaluation and clearer gains here.\n\n3. Krouglova et al. 2025 is taken as the main baseline, but as is pointed out in the discussion their approach is somewhat orthogonal. Why was no other multi-level MCMC method or baseline compared?", "limitations": "yes", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "# Strengths\n\n1. The paper is a straightforward transfer of multi-level estimators to an amortized inference setting as it is common in SBI.\n\n2. The problem is clear and one of the most important limitations of deep leanring based methods, as they are fairly sample inefficient and are not applicable to expensive simulators, fundamentally limiting the applications of state-of-the-art SBI methods.\n\n3. The framework is established with some generality for both NPE and NLE.\n\n4. The experiments demonstrate the benefits for different forms of low-fidelity approximations.\n\n\n# Weaknesses\n\nOverall the paper would benefit from being more like a recipe for the application of multi-level methods within its ML-NPE and ML-NLE settings.\n\n\n## 1. Implicit Structural Fidelity Assumptions\n\nMLMC’s success hinges on structural similarity between simulators at adjacent fidelity levels:\n\n* Outputs must be comparable and aligned for the same parameters and randomness.\n* Simulators must support shared seeds or correlated randomness.\n\nThese assumptions are not clearly articulated in the main text, and their failure modes are not explored. It is not fully clear whether the results only apply to these experiments because they were chosen in this way, or whether they generalize to a large classes of existing simulators.\n\n## 2. Gradient Surgery Requirements\n\nWhile the MLMC loss is principled, the paper heavily relies on gradient surgery heuristics—especially balancing and orthogonalization of gradients between fidelity levels—to make training stable. These are:\n\n* Presented with no theoretical justification.\n* Disconnected from the variance bounds derived in the paper.\n* Not evaluated via ablation studies, despite being empirically essential.\n\nThis creates a major gap between the theory and practice of the paper. \n\n## 3. Vacuous Variance Bounds\n\nThe paper presents formal bounds on loss variance in Section 4, but:\n\n* This does not translate directly to gradient variance, which governs its optimization.\n* The optimization dynamics are dominated by gradient imbalance and instability, which are not treated theoretically.\n\nThis leaves the theoretical results largely disconnected from the practical concerns addressed through heuristic gradient adjustments. While the contributions of bounds is appreciated their relationship to the rest of the paper should be made clearer.", "was_revised": true, "final_justification": "The paper is an important contribution to simulator based inference. My concerns have been addressed in the rebuttal. The remaining limitations are the need for gradient surgery, and it is unclear how generally it is applicable.", "rebuttal": "Thank you very much for your effort and time in writing a detailed review.\n\n>Implicit Structural Fidelity Assumptions: Outputs must be comparable and aligned for the same parameters and randomness. These assumptions are not clearly articulated in the main text, and their failure modes are not explored.\n\nWe agree that these assumptions should have been more clearly stated, and we will make this discussion more prominent in the method section of the camera-ready version.\n\nMLMC requires coupling between the data from different fidelity levels in terms of having common random numbers and parameters. In cases where this is not possible, MLMC is not applicable, which is a limitation. However, our method is still applicable when there are partial common random numbers between the different fidelities, as is the case with the toggle switch example (dimension of $u$ is $2T + 1$, where $T$ denotes the fidelity parameter). Similarly, our method can also be applied when the low and high-fidelity simulators only partially share the parameters. However, estimating posterior of parameter that is only present in the high-fidelity simulator can be challenging since we only use high-fidelity data in the correction terms of the MLMC loss. We explore this failure mode in the following additional experiment.\n\nWe include an additional parameter in the high fidelity simulator, termed \"initial value\" $\\theta_4 = x_0 \\sim \\text{uniform}([0, 4])$ in Ornstein–Uhlenbeck process (experiment 5.2) with $n_1 = 100$, instead of fixing it at $x_0 = 2.0$. We additionally trained the conditional density estimator using only $n = 100$ data from the high fidelity simulator. All the other settings remain the same. \n\n\n|             |     ML-NPE          |      TF-NPE         |   NPE (high only)    |\n|:-----------:|:-------------------:|:-------------------:|:--------------------:|\n| **NLPD ↓**  |  $\\mathbf{-0.18}$  $(0.04)$  |  $-0.89$  $(0.52)$  |  $\\mathbf{-1.21}$  $(0.53)$  |\n| **KLD ↓**   |  $\\mathbf{0.98}$  $(0.08)$   |  $1.24$  $(0.40)$   |  $1.05$  $(0.14)$            |\nTable 1: Mean and standard deviation of the metrics.\n\nWe observe that although our method shows a slight improvement in KLD, it notably underperforms in NLPD compared to using only high-fidelity data. This may be due to our training objective being more unstable and therefore more sensitive to large differences in simulator outputs introduced by the additional parameter $x_0$. We will include this result and discuss the different assumptions in depth in the main text. \n\n---\n>Gradient Surgery Requirements\n\nWe agree that the approach of the optimisation is heuristic without theory, but theory for gradient surgery would be a major project in itself and is therefore beyond the scope of this paper. However, as suggested, we conducted ablation studies on the effect of gradient adjustment and the performance which now provide much more detail on this approach. \n\nWe trained both ML-NPE and ML-NLE using (i) our gradient adjustment approach which involves both rescaling and projection, (ii) only rescaling, (iii) only projection, and (iv) no gradient adjustment (standard training). We then compared the performance of NPE and NLE on the g-and-k and NLE on the toggle switch experiment. Other than the choice of the gradient adjustment, all the training method / hyperparameters, and evaluation methods remain the same. The results are shown in the table below.\n\n\n|                                |        (i) both         |    (ii) only rescaling    |    (iii) only projection    |       (iv) standard        |\n|:------------------------------:|:-----------------------:|:-------------------------:|:---------------------------:|:--------------------------:|\n| **g-and-k NPE (NLPD ↓)**       |  $\\mathbf{-0.30}$  $(0.31)$  |    $-0.21$  $(0.25)$       |     $-0.11$  $(0.11)$        |     $-0.13$  $(0.27)$       |\n| **g-and-k NLE (KLD ↓)**        |     $0.22$  $(0.47)$       |  $\\mathbf{0.21}$  $(0.45)$  |      $0.24$  $(0.46)$        |      $0.24$  $(0.28)$       |\n| **Toggle-switch NLE (MMD ↓)**  |  $\\mathbf{0.26}$  $(0.25)$  |     $0.50$  $(0.37)$        |      $0.35$  $(0.27)$        |      $0.59$  $(0.31)$       |\n\nTable 2: Mean and standard deviation of the metrics. For g-and-k NPE, $n_0 = 1000$, $n_1 = 100$, $m = 1000$, for g-and-k NLE, $n_0 = 10^4$, $n_1 = 300$ as in Section 5.1, for toggle-switch NLE, $n_0 = 10000$, $n_1 = 500$, $n_2 = 100$ as in Section 5.3.\n\nWe observe that when the MLMC loss diverges under standard training, as is the case for the toggle switch experiment and g-and-k with NPE, our gradient adjustment approach that combines both gradient rescaling and projection  yields the best results.  In the case of NLE on the g-and-k simulator, the loss does not diverge and standard training is sufficient. In such cases, the gradient adjustment yields similar results as standard training, albeit with an increase in the variance of the metric. Therefore, we suggest optimising using our gradient adjustment approach when using ML-NLE or ML-NPE as it avoids the need to first detect whether the loss is diverging or not. We will naturally include these results and the associated discussion in the main text.\n\nWe note that the optimisation requires some form of regularization, and the one we used is one possible choice. According to our experiments, gradient adjustment performed the best among the options we tried (e.g., regularizing the contribution of the difference term when it dominates the first term, or penalizing high variance in the difference term). We will clarify this in the main text.\n\n \n---\n\n>The paper presents formal bounds on loss variance in Section 4, but This does not translate directly to gradient variance, which governs its optimization.\n\nThank you for the insightful point; we completely agree with you that a bound on the gradient would be interesting, and thankfully this is not very difficult to derive from the existing proof technique. This only requires replacing (A3) so that it holds for the gradient of the log-conditional density ($\\nabla_\\phi \\log q_\\phi^\\text{NLE}(x | \\theta)$ for NLE and $\\nabla_\\phi \\log q_\\phi^\\text{NPE}(\\theta | x_1,\\dots x_m)$) rather than the log conditional density itself. The rest of the proof is almost identical. We will now include this additional result in the paper.\n\n---\n\n>Why was the paper not targeted on an expensive real world simulator experiment?\n\nThank you for the suggestion. Actually, the CAMELS simulations suite we used in Section 5.4 is one of the most expensive cosmological suites ever run (a small fraction of the next generation are being run on the UK DIRAC HPC Facility with 15M CPUh). These are real state-of-the-art physics simulations being used for SBI analyses with real-world (observational) data; see [4]. This was clearly poorly explained, and we will remedy this for the camera-ready version. \n\n---\n\n>Similarly, why is the evaluation in 5.4 for the Cosmological experiment not more convincing? Figure 5 left shows very high variance of ML-NPE and the empirical coverage (right) is very hard to interpret without understanding the specific simulator well. I would expect a more thorough evaluation and clearer gains here.\n\n We agree that our analysis of the results from the cosmological simulations is limited, primarily because of the low amount of data, but we did use a state-of-the-art simulator. We meant our results as illustrations that combining multifidelity simulators with our method leads to better performance, and we refer the readers to recent works [3-4] for a more in-depth illustration of the benefits of multifidelity in this context. \n\nRegarding the variance of ML-NPE, actually it is lower than that of NPE (ML-NPE: $0.78$, NPE: $0.98$). It probably seems like the variance is higher due to the long tail of NLPD for ML-NPE, which is due to a few outliers ($<1$\\% of the test dataset). Removing those outliers reduces variance of ML-NPE to $0.60$, however, we reported the results with the outliers for full transparency. \n\n---\n\n>Krouglova et al. (2025) is taken as the main baseline, but as is pointed out in the discussion their approach is somewhat orthogonal. Why was no other multi-level MCMC method or baseline compared?\n\n   For amortized SBI methods, Krouglova et al. 2025 is the only existing multi-fidelity method, although [5] came out just as we were submitting this paper to NeurIPS (this is also a two-stage approach like Krouglova et al. 2025). Other alternatives include the multi-fidelity ABC methods referenced in the paper. However, such a comparison would not lead to conclusions about the different multi-fidelity approaches as ABC (or any other sampling based SBI method) is known to yield very different posteriors than NPE/NLE for the same amount of data. Moreover, it is not amortised, so computing performance metrics over multiple observed datasets is computationally expensive. If the reviewer has any particular method in mind then we would be happy to add some comparison.\n\n---\n**References:**\n\n[1] Krouglova, Anastasia N., et al. \"Multifidelity simulation-based inference for computationally expensive simulators.\" arXiv preprint arXiv:2502.08416 (2025).\n\n[2] Yu, Tianhe, et al. \"Gradient surgery for multi-task learning.\" Advances in neural information processing systems 33 (2020): 5824-5836.\n\n[3] Thiele, Leander, Adrian E. Bayer, and Naoya Takeishi. \"Simulation-Efficient Cosmological Inference with Multi-Fidelity SBI.\" arXiv preprint arXiv:2507.00514 (2025).\n\n[4] Lovell, Christopher C., et al. \"Learning the Universe: Cosmological and Astrophysical Parameter Inference with Galaxy Luminosity Functions and Colours.\" arXiv preprint arXiv:2411.13960 (2024).\n\n[5] Tatsuoka, Caroline, et al. \"Multi-fidelity parameter estimation using conditional diffusion models.\" arXiv preprint arXiv:2504.01894 (2025).\n\n---\n\nThank you again for your careful assessment. Please let us know if there is anything else we can clarify further."}, {"review_id": "yZF49YBh2Z", "reviewer": "Reviewer_HeBn", "summary": "This paper proposes an approach to more simulation-efficient SBi in the case that multiple simulators with varying cost/fidelity are available, similar to Krouglova et al. The approach is grounded in multilevel Monte Carlo and is applicable to NLE and NPE. In various comparisons, the method is competitive with or outperforms baselines.", "questions": "1. Can you decouple the contributions of the two gradient modifications used? How does optimization go if you use only one of them?\n2. Given some rough estimates of cost and fidelity, can you provide any heuristics on how to choose the number of simulations for each simulator? Can you show how performance varies as you vary the allocation of a fixed computational budget to e.g., two simulators, one high- and one low-fidelity?", "limitations": "Yes", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strengths:** The proposed method is theoretically grounded and quite straightforward. It performs well in all of the evaluations shown without requiring any additional hyperparameters. The method can be applied to cases where there are more than 2 varying-fidelity simulators and it can also be combined with other existing approaches, so it is quite broadly applicable. The text itself is very clearly written and easy to understand.\n\n**Weaknesses:**\n- The question of how to allocate a fixed computational budget to simulations remains unaddressed (though it is arguably beyond the scope of this work). The authors note that computing the optimal samples for each level is not tractable, but some empirical results/guidance on this question would still be useful in my opinion.\n- In a few of the experiments (Fig 2a, Fig 5), the proposed method achieves better mean/median scores but also has longer tails with poor performance, suggesting the multilevel objective can introduce higher performance variability across samples than the baselines.\n- Though simple, the proposed loss poses problems with optimization that require some gradient hacks.", "was_revised": true, "final_justification": "The proposed method is relatively simple, theoretically grounded, addresses an important problem, and it works. My concerns/questions about the gradient modifications and general heuristics for effectively allocating simulation budget were well addressed by the authors' response. I don't have a good sense of how limiting the need for seed matching and structure alignment across fidelities is, as mentioned by reviewer kZeL, but the work is otherwise thorough and clearly presented.", "rebuttal": "Thank you for your constructive feedback and the time and effort applied in reviewing our manuscript. We address your comments and questions below.\n\n>The question of how to allocate a fixed computational budget to simulations remains unaddressed (though it is arguably beyond the scope of this work). The authors note that computing the optimal samples for each level is not tractable, but some empirical results/guidance on this question would still be useful in my opinion.\n   \n Agreed. Though our theory does not provide explicit results on how to select the $n_0, \\dots, n_L$ due to some intractable constants in Theorem 2, it can still provide some useful intuition on how to choose the number of samples at each level. \n\nFor instance, consider a two-level simulator. According to our theory, if the low-fidelity simulator is known to be a good approximation of the high-fidelity one, it makes sense to allocate large budget to generating low-fidelity data, while only a small amount of high-fidelity data would be sufficient. On the other hand, if the outputs differ significantly, allocating more budget for high-fidelity data makes more sense, despite the higher cost, to accurately capture the differences. This principle is formalised by our theory and can be extended to simulators with arbitrary fidelity levels. In the toggle switch simulator, we allocated a large budget to approximating the difference between the low- and medium-fidelity simulators, based on prior knowledge that this difference is substantial and worth spending large budget. This lead to improved performance, as we demonstrate with an additional experiment below.\n\n---\n\n>Can you show how performance varies as you vary the allocation of a fixed computational budget to e.g., two simulators, one high- and one low-fidelity?\n\nWe conducted an experiment on the toggle switch simulator with different allocation of $n_0, n_1, n_2$, keeping the total computational cost roughly the same. Apart from the choice of $n_0 = 10,000, n_1 = 500, n_2 = 100$ that we used in the paper (option A) and is guided by our theory, we also include two other alternatives: (B) $n_0 = 9260, n_1 = 200, n_2 = 300$ and (C) $n_0, n_1, n_2 = 1077$. Option B places more budget on learning the difference between the medium and high-fidelity simulator (opposite of option A), and option C allocates equal number of samples for all the terms. The results below indicate that option A is the best performing, indicating performance can be improved by careful assignment of computational budget led by insights from our theory.\n\n\n|                |    Option A           |    Option B         |    Option C         |    only high       |    only medium     |    only low        |\n|:--------------:|:--------------------:|:-------------------:|:-------------------:|:------------------:|:------------------:|:------------------:|\n| **MMD ($\\downarrow$)** | $\\mathbf{0.16}$  $(0.23)$ |  $0.33$  $(0.26)$   |  $0.29$  $(0.28)$   |  $0.43$  $(0.29)$  |  $0.37$  $(0.45)$  |  $0.59$  $(0.65)$  |\n\nTable 1: Mean and standard deviation of MMD. The results for Option A, only high, medium, and low remain the same as in the paper.\n\n---\n\n>Can you decouple the contributions of the two gradient modifications used? How does optimization go if you use only one of them?\n\nGood point. As discussed with reviewer kZeL, we have conducted an ablation study on the gradient adjustment. We trained both ML-NPE and ML-NLE using (i) our gradient adjustment approach which involves both rescaling and projection, (ii) only rescaling, (iii) only projection, and (iv) no gradient adjustment (standard training). We then compared the performance of NPE and NLE on the g-and-k and NLE on the toggle switch experiment. Other than the choice of the gradient adjustment, all the training method / hyperparameters, and evaluation methods remain the same. The results are shown in the table below.\n\n|                                |        (i) both         |    (ii) only rescaling    |    (iii) only projection    |       (iv) standard        |\n|:------------------------------:|:-----------------------:|:-------------------------:|:---------------------------:|:--------------------------:|\n| **g-and-k NPE (NLPD ↓)**       |  $\\mathbf{-0.30}$  $(0.31)$  |    $-0.21$  $(0.25)$       |     $-0.11$  $(0.11)$        |     $-0.13$  $(0.27)$       |\n| **g-and-k NLE (KLD ↓)**        |     $0.22$  $(0.47)$       |  $\\mathbf{0.21}$  $(0.45)$  |      $0.24$  $(0.46)$        |      $0.24$  $(0.28)$       |\n| **Toggle-switch NLE (MMD ↓)**  |  $\\mathbf{0.26}$  $(0.25)$  |     $0.50$  $(0.37)$        |      $0.35$  $(0.27)$        |      $0.59$  $(0.31)$       |\n\nTable 2: Mean and standard deviation of the metrics. For g-and-k NPE, $n_0 = 1000$, $n_1 = 100$, $m = 1000$, for g-and-k NLE, $n_0 = 10^4$, $n_1 = 300$, $m = 1$ as in Section 5.1, for toggle-switch NLE, $n_0 = 10000$, $n_1 = 500$, $n_2 = 100$ as in Section 5.3.\n\nWe observe that when the MLMC loss diverges under standard training, as is the case for the toggle switch experiment and g-and-k with NPE, our gradient adjustment approach that combines both gradient rescaling and projection yields the best results. In the case of NLE on the g-and-k simulator, the loss does not diverge and standard training is sufficient. In such cases, the gradient adjustment yields similar results as standard training, albeit with an increase in the variance of the metric. We will include this in the limitations. Therefore, we suggest optimising using our gradient adjustment approach when using ML-NLE or ML-NPE as it avoids the need to first detect whether the loss is diverging or not. \n\n---\n>In a few of the experiments (Fig 2a, Fig 5), the proposed method achieves better mean/median scores but also has longer tails with poor performance, suggesting the multilevel objective can introduce higher performance variability across samples than the baselines.\n\n Yes, this is true from what we observe. For Figure 2a, higher variance is not necessarily associated with poor performance. Upper tail of the box plot roughly remains the same for all $n, n_1$. This means that as we increase the simulation budget, performance improves consistently—not just on average, but also without any drop in the worst-case performance. \n\nIn the case of Figure 5, the longer tail is due to a few outliers. In fact, our method has narrower $50$\\%, $90$\\%, $95$\\%, empirical interval without removing outliers. \nRegarding the variance of ML-NPE, actually it is lower than that of NPE (ML-NPE: $0.78$, NPE: $0.98$). It probably seems like the variance is higher due to the long tail of NLPD for ML-NPE, which is due to a few outliers ($<1$\\% of the test dataset). Removing those outliers reduces variance of ML-NPE to $0.60$, however, we reported the results with the outliers for full transparency.\n\n|         |  2.5%       |  5%          |  25%         |  50%            |  75%            |  95%            |  97.5%          |\n|:-------:|:-----------:|:------------:|:------------:|:---------------:|:---------------:|:---------------:|:---------------:|\n| ML-NPE  |  $-2.07$    |  $-1.87$     |  $-1.41$     |  $\\mathbf{-0.97}$ |  $\\mathbf{-0.24}$ |  $\\mathbf{0.58}$  |  $1.08$         |\n| NPE     |  $\\mathbf{-2.38}$ |  $-\\mathbf{2.20}$ |  $\\mathbf{-1.47}$ |  $-0.30$         |  $0.32$          |  $0.73$          |  $\\mathbf{0.85}$ |\n\nTable 3: Quantile of NLPD ($\\downarrow$)\n\n\n|         |   50%         |   90%         |   95%         |\n|:-------:|:-------------:|:-------------:|:-------------:|\n| ML-NPE  | $\\mathbf{3.15}$ | $\\mathbf{2.45}$ | $\\mathbf{1.17}$ |\n| NPE     | $3.23$          | $2.92$          | $1.80$          |\n\nTable 4: Width of the empirical interval of NLPD ($\\downarrow$)\n\n---\n\nThank you again for your thorough evaluation of our paper. Please don’t hesitate to let us know if there’s anything else we can clarify."}, {"review_id": "CJJ0SEPgTd", "reviewer": "Reviewer_PW7W", "summary": "This paper addresses the challenge of reducing simulation costs in simulation-based inference (SBI). The authors propose a method that leverages the Multilevel Monte Carlo (MLMC) technique for SBI, enabling significant cost reduction by utilizing a sequence of simulations with varying levels of accuracy. Theoretical results support the correctness of the proposed method. The effectiveness of the approach is demonstrated using several simulators, including a real-world cosmology simulator.", "questions": "• P115: Regarding the phrase “making NPEs fully amortized” — it would be helpful to clarify what “amortized” means in this context, particularly in comparison to NLE. My understanding is that “amortized” refers to the fact that NPEs do not rely on classical sampling schemes such as MCMC (as used in NLE), thereby avoiding their associated drawbacks.\n\n• I believe Theorem 1 and 2 is a very important result as it shows how the upper bound of the variance depends on the properties of the simulator. My understanding is that the motivation for applying MLMC to SBI lies in its ability to reduce the variance of the loss compared to standard MC. In addition to these theorems, if it can be theoretically shown that the upper bound of MLMC is smaller than that of MC, it would be very interesting.\n\nRegarding the paragraph on the optimisation of the MLMC objective: I understand that this optimisation is a crucial aspect of implementing MLMC. In the current structure, the discussion is placed in Appendix A, but since even the core ideas of the method are important enough, I would hope to see a more detailed discussion included in the main text.", "limitations": "Yes.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 4, "originality": 2, "strengths_and_weaknesses": "Strength \n\n•  The paper is clearly written and well-structured. \n\n•  It is well-motivated by real-world challenges, addressing a major limitation in SBI: the high computational cost of simulations.\n\n•  In addition to extensive simulation results, the paper includes a detailed theoretical analysis. Some of the theoretical results also provide useful implementation insights, which further support the proposed method.\n\n•  As the authors mention, applying MLMC to modern problems such as SBI or probabilistic numerics is a timely and relevant direction.\n\nWeekness \n\n・The technical and methodological novelty of applying MLMC to SBI may be relatively straightforward.\n\nTaking the weekness into the consideration, the paper’s theoretical contributions—particularly the non-trivial conditions that offer valuable insight into implementation—and solutions to practical challenges, especially in optimization, are substantial enough to outweigh this weekness. Overall, I like this paper, which make me think that this paper deserve accepted in NeurIPS. I decide to assign a score of 5.  Although I am giving an Accept rating prior to the rebuttal phase, I hope that further revisions during the discussion phase will enhance the quality of the paper!", "was_revised": true, "final_justification": "This paper is well-written, and has a good contribution. I vote for acceptance.", "rebuttal": "Thank you for your positive comments and strong support for our work. We address your comments and questions below.\n\n>The technical and methodological novelty of applying MLMC to SBI may be relatively straightforward.\n\n  We agree that the application of multi-level Monte Carlo to NPE and NLE is straightforward. We see this as an advantage of our approach, since it makes it more likely that the method will be widely adopted. Also, as you pointed out, the extensive theory, experiments, and practical solutions hopefully provide strong contributions towards solving the computational issue of dealing with costly simulators.\n\n---\n\n>P115: Regarding the phrase “making NPEs fully amortized”—it would be helpful to clarify what “amortized” means in this context, particularly in comparison to NLE.\n\n We use the phrase \"fully amortized\" in the context of NPE to mean that obtaining the posterior distribution for a new observed dataset is instantaneous through a simple forward pass of the neural network, once the network is trained. This is to differentiate NPE from NLE, which can be considered to be \"partially amortized\": although the conditional density $q_\\phi(x \\mid \\theta)$ can be evaluated instantaneously as well for different $\\theta$ and $x$ values once trained, obtaining the posterior distribution still requires MCMC schemes. We will make this distinction explicit in the main text.\n\n---\n\n>I believe Theorem 1 and 2...if it can be theoretically shown that the upper bound of MLMC is smaller than that of MC, it would be very interesting.\n\n This is a very good point, thank you very much. In fact, our theory already allows for this since the bound for $\\text{Var}[h_0(\\phi)]$ in Theorem 1 is a bound for an MC-type objective for the high-fidelity simulator if we replace $n_0$ with $n_L$ and $G^0$ with $G^L$ (assuming the MC objective is evaluated only on high-fidelity simulations). By comparing such a bound with the MLMC bound in Equation 3, we can come up with conditions on $n_0, \\dots, n_L$ and $G^0, \\dots, G^L$ that lead to a smaller MLMC upper bound than that of MC. However, note that this is only a comparison of upper bounds, as opposed to a direct comparison of variances, which is typically out of reach of MLMC theory. We will clarify this interesting point in the camera-ready version.\n\n---\n\n>Regarding the paragraph on the optimisation...I would hope to see a more detailed discussion included in the main text.\n\n  Agreed. Due to space constraints, we initially had to move some of the discussion to the appendix, but we are planning on moving this back to the main text for the camera-ready version of our paper.\n\n---\n\nThank you again for your careful review of our paper. Please let us know if there is anything else we can clarify."}]

Author rebuttal:

Meta-review: The average rating is 5, with all reviewers unanimously recommending acceptance (5, 5, 5). We encourage the authors to implement the modifications discussed with the reviewers in the final version of the paper.

Recommendation: acceptance.
Decision: Accept (poster)