Multilevel Control Functional
Authors: Kaiyu Li
Yiming Yang
Xiaoyuan Cheng
Yi He
Zhuo Sun
Venue: ICLR 2026 Poster (2026)
Keywords: Variance Reduction
Monte Carlo
Primary area: probabilistic methods (Bayesian methods, variational inference, sampling, UQ, etc.)

Abstract: Control variates are variance reduction techniques for Monte Carlo estimators. They play a critical role in improving Monte Carlo estimators in scientific and machine learning applications that involve computationally expensive integrals. We introduce \emph{multilevel control functionals} (MLCFs), a novel and widely applicable extension of control variates that combines non-parametric Stein-based control variates with multi-fidelity methods. We show that when the integrand and the density are smooth, and when the dimensionality is not very high, MLCFs enjoy a faster convergence rate. We provide both theoretical analysis and empirical assessments on differential equation examples, including Bayesian inference for ecological models, to demonstrate the effectiveness of our proposed approach. Furthermore, we extend MLCFs for variational inference, and demonstrate improved performance empirically through Bayesian neural network examples.

Published as a conference paper at ICLR 2026
MULTILEVEL CONTROL FUNCTIONAL
Kaiyu Li
Yiming Yang? Xiaoyuan Cheng?
Yi He
Zhuo Sun 4.5+
1 COMAC Shanghai Aircraft Design and Research Institute
2 *Department of Statistical Science, University College London
Dynamicc Systems Lab, University College London
"School of Statistics and Data Science, Shanghai University of Finance and Economics
$ "Intiituee Big Data Research, Shanghai University of Finance and Economics
ABSTRACT
Control variates are variance reduction techniques for Monte Carlo estimators. They play critical role in improving Monte Carlo estimators in scientific and machine learning applications that involve computationally expensive integrals. We introduce multilevel control functionals (MLCFs), novel and widely applicable extension of control variates that combines non-parametric Stein-based control variates with multi-fidelity methods. We show that when the integrand and the density are smooth, and when the dimensionality is not very high, MLCFs enjoy faster convergence rate. We provide both theoretical analysis and empirical assessments on differential equation examples, including Bayesian inference for ecological models, to demonstrate the effectiveness of our proposed approach. Fur- thermore, we extend MLCFs for variational inference, and demonstrate improved performance empirically through Bayesian neural network examples.
INTRODUCTION
The paper focuses on the estimation of intractable integrals, where the integrands lack closed-form solutions expressions and are computationally expensive evaluate. The integrals are of the form
[[ = [ f(x)
where II is a distribution with a Lebesgue density on X C Rd and : - R is the integrand of interest. Assume that f is square-integrable, that is. II II oc. This is a common challenge in diverse applied fields such as finance (Glasserman 2004 Chen et al. 2024). aerospace engineering (Morio & Balesdent 2015 Geraci et al.] 2017), hazard ana asilysis (Geist & Parsons 2006 Dalbey et al.) 2008), medical physics (Rogers) 2006). among many others. This also frequently arises in statistics and machine learning, such as computing normalizing constants in probabilistic models (Atay-Kaya
& Massam, 2005 Ohsaka & Matsuoka, 2020 Chehab et al., 2023), performing Bayesian inference
(Golinski el al.] 2019), training energy-basee based models (Song & Kingma 202T) and variational inference
(Buchholz et al. 2018: Vahdat & Kautz 2020
Fujisawa & Sato 2021).
Motivation Monte Carlo (MC) Liu 2001 Rubinstein & Kroese 2016) is the most widely used approach for estimating the integrals defined in Equation (Ip. However, MC estimators tend to have high variance and slow convergence rates (Assaraf & Caffarel 1999 Oates et al., 2019). Meanwhile, when dealing with complex scientific models, sampling or evaluating the integrand can be computationally expensive, e.g., large-scale computer simulators or costly experiments (Sánchez) Linares al. 2016). To achieve the desired accuracy, the overall sampling and evaluation cost with standard MC can be prohibitive. One way to improve the efficiency of MC estimators is to reduce the variance of the integrand by control variates. Control variates (CVs) (Robert et 1999). including both parametrized CVs (Assaraf & Caffarel 1999; South al. 2022; Sun et al. 2023b Sun 2023) and non non-paramen con-parameterized CVs (Oates et al.] 20117 2019), are variance reduction techniques for MC estimators. This is achieved by designing and learning function g (known control variate)
well correlated to the integrand f. Another approach for computationally expensive functions is
"Corresponding Author. Correspondence to Zhuo Sun: sunzhuo@ma shufe.edu.co

---

Published as a conference paper at ICLR 2026
(a)
(b)
(c)
1.0
(MUM
0.04-
0.8
n(G)
0.03-
0.9
0.6
(MLCF)
0.02-
0.01
0.4
A 4-6 (ML
6.03-
0.2
5.4
0'a
1.0
0.0
0/
0.0
1/b
0.6-
(d)
0008
(e)
0.5
4-4 (HLM))
0.530
1.002
6-4 (ML
0.8 0000
0.4
Figure 1: Illustration Example. (a): fo. fi and f2 are coarse, medium and fine approximations to (b)cc)d):: Compared with MIMC, after applying MLCF, the green curves become much flatter and closer to [[ - ft-1] (red dotted lines) than the original Ji fi-1 (blue curves) used at each MLMC level. (e): Compared with CF, although applying CF to fz (purple curve) already reduces variance, MLCF leverages the multilevel structure. Thus, fine levels, c.g. f2 - fi (in (d)), themselves show smaller variance than f2 (in (e)) and MLCF further decreases this variance (in (d)). This demonstrates that MLCF reduces the variance significantly.
the use of multi-fidelity methods (Peherstorfer et al. 2018 Fernández-GGoding 2023). which have shown to be useful scheme by compensating high-cost function approximations with low-cost ones, thereby reducing the overall computational cost. This particular scheme is highly related to the framework multilevel Monte Carlo (MLMC) (Giles) 2008 2015), which designed for expensive integrands when cheaper approximations
be constructed at several levels, e.g. tsunami modeling (Sánchez-Linarg et al., 2016 Li 2024) Li et al. (025). MLMC uses a sequence of mulli--fddelin models to construct telescoping sum for the original integral in Equation ( The telescoping sum consists of the expectations of increments between successive multi-fidelity models, namely Si and fi-1, with indexing the hierarchy of accuracy levels. Given fixed computational budget, MLMC yields more accurate estimates than standard MC estimators (Giles) 2015). This gain in performance has led the widespread use of MLMC, not only in scientific applications but also as powerful tool for enhancing many areas in machine learning, such as variational inference (Fujisawa Sato, 2021). From the perspective of variance reduction for MC, MLMC treats the low-fidelity model J control variate for the high-fidelity model ft. This implies, however, that the variance of their differences not minimized. Meanwhile, the existing work of CVs does not fully explore and exploit the inherent property of these multi-fidelity models of which the 'precision' can be decided and altered. See Figure for an illustrative example, with the corresponding results in Appendix D.I
available
Motivated by the above insights and unsolved problems, we propose a broadly applicable extension of control variates that integrates the strengths of control variates and multi-fidelity methods, in a manner analogous to multilevel Monte Carlo (Giles) 2008 2015). In particular, we focus on non-parametric Stein-based control variates, control functionals (CFs), which naturally lead a novel class of estimators named Multilevel Control Functionals (MICFFs. By leveraging the multi-fidelity structure alongside derived variance bound allows us to assign the optimal sample sizes across all fidelity levels. Mcanwhile, as we shown in Section 3 and Section 4 MLCFs are widely applicable. It can be applied to inference problems under un-oormalized densities which are common Bayesian inference, and can also be extended to variational inference as shown in this work.
Limitations of Previous Techniques Although several related methods exist, their practical use often constrained by inherent limitations. For example, (Li et al., (023) combined Bayesian quadrature with MLMC, but the method was restricted to specific kernel-distributin pairs and was not suitable for complex distributions such as unnormalized distributions in Bayesian inference. Similarly, existing works multilevel control variates are tailored for specific cases. Nobile & Teseil (2013) used the solution to an auxiliary diffusion problem with smoothed coefficients as the control variate, which was only applicable specific partial differential equations, |Firbanks et al. (2017) used a low-rank representation for high-fidelity models to construct a control variate. /eracii et al. (2017) used a simplified physical model, which required additional expert knowledge. In contrast, the
2

---

Published as a conference paper at ICLR 2026
proposed MLCFs are broadly applicable and can be implemented without reliance on Nomain-speccife expertise.
In summary, the main contributions of this work are: (i) We propose novel variance reduction method, MLCFs, which leverages both the multi-fidelity structure of integrands and CFs, and is broadly applicable; (ii) We provide theoretical variance bounds as well as theoretical optimal sample sizes for MLCFs in Section 3 (iii) We extend MLCFs for variational inference in Section 3 and Section [ (iv) We demonstrate that MLCFFs are widely applicable through a series of carefully designed experiments in Section 4
Our Contributions
2
BACKGROUND
In this section, we briefly review existing constructions for Stein-based control variates, multi-fideii models and multilevel Monte Carlo methods.
2.1 CONTROL VARIATES
Monte Carlo and Control Variates We assume f is in C (II)
:= -
This assumption is often required as the variance of f. V[f] [[ - is bounded OOtes
(TIS)),
ct al. 2017 South et al. 2022). Given evaluations of the integrand ff at n independent and identically distributed (1.i.)) realisations {r,) from I1, the Monte Carlo estimator of Equation ( is
LM = E f(x)).
The above estimator follows central limiting theorem: [i[Icclf--iij - ((OVV( Thus, the convergence rate of the estimator is often determined by the sample size n and the variance V[f]. often requires a large number of function evaluations to achieve the desired accuracy. Similar results hold for Markov Chain Monte Carlo (MCMC) (Dellaportas & Kontoyiannis 2012 Alexopoulos ict al., (023) and quasi-Monte Carlo methods (Hickernel et al.; In this work, we will focus on the Monte Carlo case. One way to improve the performance of Monte Carlo (MC) estimators is to identify a function g E (1I) with known mean 1l[g] such that Vlf g] is small. Such is also known control variates (CVs). Finding g with known mean [[g] is then the first challenge. When II is relatively simple, ad-hoc methods such as Taylor expansions of the integrand f can be used (Paisley et al. 2012; Wang et al., 2013). While for more general and complex distributions, the ones often encountered in Bayesian inference, we can employ Stein's 's method (Anastasiou et al. 2023 to construct such functions which are also known as Stein-based control variates.
General Recipe of based Control Variates Stein-based control variates (Oates et al. 2017 al. 2022; South et al.) 2022] [Snn et al. b) are variance reduction tools for Monte Carlo integration. They also widely used in the cases when the density is unnormalized and when the samples are MCMC samples, which often appears in Bayesian inference. The first step is to construct candidate set G such that [[y] 0 for Vy G. This can be achieved by using Stein's operators Sn (e.g. the Langevin Stein operator); see (Anastasiou et al. 2023) for a detailed review. By using the zero-mean property, we have II[f- - 9] control variate E G with reduced variance, i.c., V[f g] = II[[f- - - II[f g]) < [[]]] (Oates) al. 2017 Zhu et al. 2019; South et al. 2022). Such effective control variate g is often learnt by minimizing the empirical (penalized) variance of Vif g] conditioning on m samples {2, ) and their function evaluations from all samples { available. Then, through estimating II gl
II for Vg G. The second step is to select an effective
with the remaining m function evaluations, we can get an estimate of I[[S] with reduced variance and improved accuracy, given by
flvlff= = 1 El (f(z,) g((x))
Control Functionals We consider non-parametric family of control variates, control functionals (CFs) (Oates et al., 2017 2019), which is designed for single Monte Carlo integration problems. is class of non-parann parammeriic Saeen-based control variates based reproducing kernel Hilbert spaces (RKHS). applies the Langevin Stein operator Sn onto vector-valued functions u E C'l (X') x ... x Cd (2)) which takes the form Sn (n]))(r)
V. u(xx) u(x)- V2 log n (x) where V. denotes the
3

---

Published as a conference paper at ICLR 2026
divergence operator and V denotes the gradient operator. Each component function u, : + R is constrained to belong to a Hilbert space H. Let Hy be the RKHS induced by a reproducing kemee k. The image of u := H₂ X X H₂ under Sn is RKHS g with kernel Kg (also known as a Stein kernel); see Equation (TO) in Appendix [A22 Oates et al. 2017 2019) used functional approximations
s(x) 3 Sn [u] (xx) where / and are selected by sol volving constraint least-square optimisation problem in G conditioning on m samples { and (f(z) (z+)) y The control functional takes the form of: 9m (x) s(x) [[[s]. The standard control functional estimator is given by
[F "2 [f] n-m 1 (f(x)) ko(X', (X., (xxx=)(((x)) - (2 -2 ko Ro
.x0 X²) 1(xx))
}1}}
where X = (#) Fmm T X1 = (Im+1, ,intt,,.,,,np Fn) Fn) (I(X))); = (xx), (Ko (bo = ko(z), 2.), for all i,j {l.... J............................... and (f(X*)): f(zm++), (Ko (oox') ko(=m+i, T.y. for all i E {1.....
m) and for all j {1.....,) ). A major drawback of control functional is the ((m²) computational cost, which can be alleviated by using stochastic optimization et al. 2019 Si et al., [2022] Sun et al.) 2023a b). Meanwhile, this would not be an severe issue in the setting considered in this work as such cost much smaller than the cost of the evaluation of integrand.
in (Zhu
2.2 MULTI-FIDELITY MODELS AND MULTILEVEL MONTE CARLO
Multi-fidelity models have been used to accelerate a wide range of algorithms and related applications, including uncertainty propagation, inference, and optimization. The main intuition behind multi- fidelity models is to employ cheaper and less accurate models with low computational cost (a.k.a. low fidelity models) to generate additional supplementary data for the expensive and more accurate models ka. high fidelity models). See (Peherstorfer et al. 2018) for detailed review.
Multilevel Monte Carlo (Gilss] 2008 2015) uses a hierarchy of approximations fo- fi...., f to fr f with increasing levels of accuracy and cost to estimate the integral of interest. The method can achieve a higher accuracy with a lower computational cost compared to MC using only the fo f. Given the sequence of approximations, MLMC sums up the estimates of the corrections with respect to the consecutive lower level and obtain the telescoping sum
[[f] [[fz]= Ei llf: fi-1),
where f= 0 to simplify the equations. MLMC estimates each of these integrals in the telescoping sum independently. At each level, MLMC uses a MC estimator to estimate II -ufi-Ja fi-] by drawing ii.d samples {I(( J from II and evaluating fi(().)) and fi-1(xq.)). Therefore, the unbiased MLMC estimator takes the form
AMLMc[ = Et AMclfi f(-1) -t Xy E -
From the view of variance reduction, fi-1 can be regarded as a control variate for fi for all levels. has shown that by carefully analysis of the number of samples assigned to each level, MLMC can largely improve efficiency over standard MC when the desired accuracy is fixed (Giles) 2015). Such MLMC methods have also shown successes in various fields including probabilistic numeric (Li et al. 2023), simulation-based inference for multiple simulators of various fidelity (Hikida el al. 2025 and variational inference (Fujisawa & Sato 2021
Shi & Cornish 2021
METHODOLOGY
We now present the proposed method multilevel control functionals (MLCFs). In general, a MLCFs estimator of TT[) in Equation (1) takes the form of,
[f E'e a- E. (f(t))) f fi-1(2u.))
(s((2(.)) - I[s;)),,
(3)
where 81 [[s]] is the control functional at each level I, a non-parametric Stein-based control variate as discussed in Section
Proposition 3.1. Given Xt, the associated score evaluations {V log (F M and the function evaluations (f((ma,)) - (1-1( 3h for I E (O...., L.), split it into two parts: X²
4

---

Published as a conference paper at ICLR 2026
(.. (1.m))) D and X] ()) ) together with their score evaluations and function evaluations. The MICCF estimator of [f] is unbiased and has the following form:
i == Elzo " [/i fc-i]
(4)
- 1 ((f((x)) ft- (x-)( -b where 0₁ 1 kx -Xx()))(())())))))))) fi- (X 11)
mi),
The proof is provided in Appendix B.1 It is also common in practice to use simplified estimator for each level; see Appendix [ for details. Although the simplified estimator is biased, it usually has superior mean squared error (Oates et al. 2019). We also provide practical methods for kernel hyperparameter selection in Appendix C.2
The proposed method, MI.CFs, is simple yet effective, widely applicable, and offers several advan- tages. (i) The MLCF estimator in Proposition i3 unbiased and achieves a fast convergence rate under mild assumptions. Moreover, users have the flexibility to modify the estimator, such as using control functionals only on selected low levels. (ii) The restriction on II can be relaxed. We only assume that π is smooth and x(x) > 0, so that the gradient of log π can be evaluated pointwise. In Bayesian statistics, we often only know π up to an unknown normalization constant due to the intractable marginal likelihood. (iii) The simplified MLCF estimator defined in Equation (TT) (in Appendix (B66) has no restrictions on how to generate samples. It can employ any experimental design to further improve efficiency. (iv) Implementing MI.CFs is simple and straightforward and does not require domain-specifi expertise from users.
Next, we provide theoretical amalysis of the variance of MLCF estimators, which is based on the proof of Theorem 1 of (Oates et al.] 2019). We will see that the convergence rate of MLCF is related to the smoothness of π and ft. We use (xt) ) to denote the set of measurable functions for which continuous partial derivatives exist X' up to order 9 E No. For & (X)), paak/Or, ...Ori, Ox', ... -'',,
a continuous function for all i,, ig. J1, :3 E {I....,d].
Assumptions Let 0.'' denote the boundary of X. We make following assumptions: (A1) X satisfies the interior cone condition; (A2) π E Catl( for a E No: (A3) π > 0 on X: (A4) V., log π E (2 II') for i L, .,. d for all distributions II' on X; (A5) (z))k((z.)
0 for
E 0x;; (A6) for each I E {O...., L). kn E Ch++ (x) for by No: (A7) fi, fi-1 H for every (I...., L}, where H', is a RKHS with the kernel k4, (z,x²) positive constant Cj and Stein kernel ky obtained by plugging k; into Equation 100 (A8) for each
C kg (z. x') with
I E {O..... L}. the fill-distance of the samples XP, h == supzex min; m....mm | T(,))) 2|
satisfies h4 S qm,,
1ldd for a constant q > 0.
Assumption A1 ensures that the domain is sufficiently regular so that the scattering samples can adequately cover the domain. A simple example of such domain is a hyperrectangle. Assumption A5 satisfied by a constructive approach, as demonstrated in Oates et al. (2019) and in the syn- thetic example of Section 4 Assumption A8 requires that X at each level used to construct the control functional quasi uniform. This condition can be satisfied by space-filling designs such as quasi-Monte Carlo sequences, Latin hypercube sampling, or other cow-discrepancy point sets. The rest of the assumptions ensure that the problem is well-defined.
Theorem 3.2. Suppose that the assumptions A1-8 hold and X are i.i.d at each level, when X. are sufficiently dense, the upper bound of the variance of MLCF estimator is given by
(rim,
T/d 11/ alley )²
V Vxg... x [[ S Et
(5)
AL-mi
where n. = min(a, (a, bi} and ri is constant independent of ft. fi-1 and data points.
The proof is provided in Appendix B.2
Here, 'sufficiently dense' means that the fill distance is less
than a certain threshold. This condition is common scattered data approximation theory (Wendland] 2004]. The mean squared error of MLCF is MSE (f]) = _x; xy. [[[MLCF[]]] 1[
V Xg xx..... [f + (Exg....x [[mc Ilff))2 Since MLCF is an unbiased estimator, MSB(TMM.rr()) = V(nnarr [f]]. If we assume that the proportion my mu/m is the same all levels,
5

---

Published as a conference paper at ICLR 2026
then each level, the convergence rate is O(n(-n/0)-1/2) Compare to the convergence rate of MLMC each level, which is O(n the convergence rate of MLCF is faster. The theoretical results show that MLCF converges fast when the dimensionality is small or moderate and both the integrand and the density smooth.
Theorem 3.3. Suppose that assumptions A1-A8 hold, mi/n = p and == 17 min (az, bi} do not depend on I. Then MLCF = -no ..... 11 fi-111ze_ /2/ mi) subject to Σt Cini = T for T > 0. The solution is
MLCF obtained by minimizing Etmormmm n// [|
MLCF = -Rrillli-- - fi- S1-111es, )= 2-422 W (O....,.)) L).
(6)
where R= T El
Cr
Cp (rellfe fV- न|त
2ndd
See Appendix B.3 for proof. According to the theorem, the higher the evaluation cost C₁ at level I, the fewer samples are assigned to that level. Since we expect function norm |l -f--1lla! to decrease with level, and Cj typically increases with level, Theorem 3.3 implies that the sample size decreases with level. Moreover, for larger т and lower dimensionality d, the allocation is less sensitive to [1 - f1-1| and C₁, meaning that the penalty higher levels is reduced.
The result is then used to compute the optimal sample size at each level for the synthetic example in Section 4 Since the constant rı is independent of fi. fi and data points, and the domain, d and are the same across all levels, the effect of rı could be normalized away. Alternatively, one may assume uniform bound 2 max r and use this in the proof of the theorem to get result without rj. The RKHS norm can be computed using reproducing property of reproducing kernels for the synthetic example. Formally, Kanagawa et al. (2018) Theorem 2.4) provides a general equation for the RKHS norm. practical applications, though functions may not exist the form of linear combinations of kernel functions, users can use based methods to estimate the RKHS norm, For example, Scharnhorst et al. (2022 Appendix A) provides an example of estimating the RKHS norm using randomly sampled data. There are also other data-driven approaches illustrated in Karvonen (2022): Tokmak et al. (025a
Plugging the optimal sample size from Theorem 3.3 into Theorem 3.2 yields,
2r
Vxg.....X X!..... [MLLr [f]l s ATT 2ngd
Σt
C
2r fi 11
Ilfi
where - 1-p 2r 2d,2 The detailed derivation of Equation is provided in Appendix B.4 For CF, there is only one level, which is the finest level I. with evaluation cost C. Hence,
Wx [ucelf] s AT 3nd 2 ll
We denote the two upper bounds by BMLCE and BeF, respectively. The behavior of BMLCE depends 2rtd
on how the term 77 Ifi
fi नं in BMLCF varies with I. If this term increases rapidly H4
with level I, level L dominates and BMLCE Fs A-T - Cr 2nrd l15e S1-1l1 Comparing with
BeF gives BMLCF/ BeF C llfz fL-1 11 /////// /C 21 Ilfel (1 Although Ct > C, when the evaluation cost of fL tends to be much higher than fi-1, the difference between Ct and C' is moderate. Mcanwhile, llfz fL-1llu4 much smaller than [[flllu4. Hence, BMLCE < Ber in this case. If the term decreases rapidly with level 1, the coarsest level 0 dominates and BALER RS A Ind
2rrd
2rid Co 11 leading to BAICF/ BeF = C. 2n
1 (C 2يي ||
Since Co < C and fllool < Il.zz lzzz, we again conclude BALCY < Be
Extensions: MLCFss for Variational Inference MLCFs be further extended for variational
inference (VD). To be precise, we consider the scenarios when the objective is to minimize the Kull- back-Leibler (KI.) divergence between the variational distribution {(=[A) and the posterior p(=

---

Published as a conference paper at ICLR 2026
(D denotes the observations) with respect to the parameters X. This minimization is equivalent to maximizing the evidence lower bound (ELBO): C()) = Ey(zjx)|ogg p(z, D) - log q A)]. We will focus the re- paraaetterized gradient estimator to optimize the objective function, where 2 T(r, A) is expressed as a deterministic transformation 7 of a noise variable x with distribution (3r). Fujisawa & Sato (2021) proposed the multilevel re-parameterized gradient (MLRG), which reduces the variance by recycling the previous parameters and gradients. At iteration L, the multilevel
A)
re-parameterized gradient (MLRG) has the following form: VALRG [()x.) = Ei=o [[ Sx, 1 where fx, (2) = Vx, log p Txr, 2r) Vx, log ((T Ax) with fx_, = 0 for notational convenience. Fujisawa & Sato (202)) used Monte Carlo to estimate MLRG and resulted mul- tilevel Monte Carlo re-parameterized gradient (MLMCRG) estimators, which took the form of V & (),,) [MMMM[/AA.] Σt [/M fx, .l. The proposed ML.FF then naturally re- sults in a novel series of estimators, named as multilevel control functional re-parameterized gradient
(MLCFRG) estimators:
MLLCER ((AL) fMLL [ax.cr[fx.]:= Eio - fx-ill
(8)
whose variance is controlled as shown in Theorem
3.2
Furthermore, we also provide an simplified
form for MICFRG estimators in Proposition 3.4
Proposition 3.4. The update of MI.CERG estimator under stochastic gradient descent update can be rewritten in a simplier form, given by,
X1+11 = Ar. + aL-1 (Ar. - - -sxe ilb
where on is the learning rate at the 1-th iteration. Under the same assumptions as Theorem 3.2 its variance is bounded by V[MMLCC [(XL)] S (FLmp 1 - fx / mm1).
See Appendix B.5 for proof. Note that this largely reduces the computation cost from O Eto ln to OdnE)) and the memory cost from Old Eip In,) (dd ) at iteration L where d is the number of parameters of the neural networks. The variance bound is conditional on Å Ar.-1 while the variance of full optimization trajectory is bounded by Theorem (3.2] When nt, is small, the extra time and space cost is controlled. Mcanwhile, it is also possible to further accelerate the method by implementing modern model parallelism (Shoeybi al. 2019), wrapping it as a new optimizer and using Py PyTorch foreach operators (Paszke et 2019) to avoid loops through tensors individually.
4 EXPERIMENTAL RESULTS
We now assess the performance of MLCFs through: (i) A synthetic experiment validates the effec- tiveness of the optimal sample sizes. (ii) A boundary-value ordinary differential equation (ODE) example shows that MLCFs can be further enhanced by incorporating experimental design tech- niques, and suggests that developing adaptive strategies would be a potential future direction. (iii) The Voltera system example shows that MLCFF can be used generally for Bayesian inference (e.g. when the target density is unnormalized). For both (ii) and (iii), the implementation of the other methods reviewed in Section [ is either very difficult or not feasible. neural networks illustrate the extension of MLCF to variational inference. Additional results on the illustration example are presented in in Appendix D. demonstrate that the performance of MLCFs aligns with the motivation and the intuition discussed in Section All of these examples illustrate the broad applicability and high efficiency of our method.
(iv) Examples of Bayesian
Synthetic Example The variation of the test-bed example in Oates et al. (2019) (Section 1) is used as an illustrative example to verify the effectiveness of optimal sample size derived in Theorem 3.3 II
the uniform distribution in X jo, 1]a. The kernel takes the form (a, x') = 6(z)6(z)) k) (a 2') such that the assumption A5 holds through Matefn 2.5 kernel Re x²) times a smooth function
5(x) = I. x1(1 - I.). In this example, f(xx)) are defined as
fi(x) Ei onkl, 21),
where the specific setups are provided in Appendix D.2 The optimal sample size RML.CF for MLCF and PMLMC for MLMC (see Appendix B.7 the expression) then computed.
7

---

Published as a conference paper at ICLR 2026
Under the same budget limit, compare the performance of MLCF with MMLCF, MLCF with RMLMC. MLMC with PMLMC and CF. The sample sizes are listed in Table 2 absolute error from 50 replications and visualize the result in Figure MLCF with both AMLCF and PMLMC significantly outperforms MLMC and CF. While MI.CF with TMmlcc slightly outperforms MLCF with RMLF when the sample size is small, this could be because UMLCE minimizes the upper bound of the variance of MLCF, rather than minimizing the variance directly. When the sample size high, MLCF with TMPCPP slightly outperforms MLCF with TMMMC. This is promising because, although computing PM may be difficult in various practice, implementing MLCF with AMLMC still yields good performance.
For each setting, we compute the
Boundary-value ODE The boundary-value ODE example can also be viewed as a one- dimensional elliptic partial differential equation, with random coefficient and random forcing:
_ (c{ x3 for E (0,
10"
with u(0) = (1) = where c(z) = + 21 31 (O, 0.2) and x2 ~ (0, 1). This exam- ple is a variation of the test case for MLMC in Section 7.1 Giles 2015). The integral of inter- est is [[[]] = -xx (x)dd, where x = (2),
10
X R23 f(x) = So u(z)dz is approximated with h - u(zz), where h is the step size and each u(2;) is obtained by solving the ODE with the finite difference method; see Appendix DD. for more details.
Butget
medum
high
Figure 2: Synthetic Example: Absolute integration error under a budget constraint (Y-axis log-scale).
To compare MCLF using quasi-Monte Carlo points (QMC), Latin hypercube sampling (LHS), and i..d points (ID) with MLMC, multilevel Bayesian quadrature (MLBQ) (Li et al.) 2023) and CF using i.i.d points, we repeat the experiment 100 times. Multilevel methods in this example utilize the optimal sample size for MLMC. Details about the sample size, evaluation cost at each level, and other relevant information can be found in Appendix D.3] As shown in Figure 3 under the same evaluation cost constraint, MLCF outperforms MLMC, MLBQ and CF. Figure designs can improve the performance of MLCF
also shows that experimental
1.2
Estimattor MLGF GMC
Estimetor
2.0
MLSF ID
MLDF L45
Absolute Enarr
B.
NLVO ID ID
1.5
MLCF ID
0.8
Absolute
MIVIC
1.0
MLBC
0.2
0.
0:25
0.51
0.77
0.30
0.91
1.52
Budgel
Budget ((8)
Figure 4: Bayesian Inference for Lotka- Volterra: Absolute integration error under a budget constraint.
Figure 3: Boundary-value ODE: Absolute in- tegration error under budget constraint.
Bayesian Inference forr Volterra We now consider to perform Bayesian inference for the Lotka-Voltere system (Lotka 1925 1927 Volterra 1927). which is also known the predator-prey model. The model usually uses a system of differential equations:
du (t)
z zndl (t) x2u1(()
di
duy(t) 2z (t) xt)2(t),
de
to describe the interaction between predator and its prey in an ecosystem. un (t) and (z) are the prey population and the predator population at time t E [0, s), for some s E R,. The initial conditions of the system un (0) 375 and u2 (0) zg. The observations un (1:) and u2(4,) obtained exhibit
8

---

Published as a conference paper at ICLR 2026
oog-nommal noise with independent standard deviation 27 and zş respectively, for all i {1....,,, m). We can re-parameterize I as in (Sun et al. 2023a b) such that the re- re-parameterized model has parameters 2 R With the Gaussian distribution priors assign on 2 and the observations, can construct the posterior distribution of ă. The quantity of interest [[]]] is the posterior expectation
of the average prey population over the time period between 0 and 8, i.e. [T Sus (())+(()))) where is the posterior probability distribution of i and f()) is the average prey population between O and with the model parameter ž. f(ā) = s S u1 (t)dt is approximated with (s) Σ ii (tz). where h the step size and each u1 (t,) is obtained by solving the differential equations numerically. The real-world dataset (Hewitt) 1921) consisting of the population of snowshoe hares (prey) and Canadian lynxes (predators) used as observations for our study. With the real-world observations, conduct Bayesian inference and use MCMC sampler tno-U-turn sampler) in Stan (Carpenter iet al. 2017) to obtain samples.
We compare (i) MLCF with MCMC points, (ii) MLMC framework with MCMC points (MLMCMC), (iii) CF with MCMC points and (iv) MCMC, We repeat the experiment 50 times. The sample size, sampling and evaluation cost at each level, and other details can be found in Appendix D.4 As shown in Figure under the same budget constraint, MLCF outperforms all other methods.
Variational Inference for Bayesian Neural Networks We further extend the proposed MLCF estimators for variational inference, i.e. using MLCFRG estimators for the gradient of the ELBO. In particular, we applied a Bayesian neural network regression model to the UCI wine-quality-red dataset as (Fujisawa & Sato 2021). The Bayesian neural network consists of 15-unit or 20-unit hidden layer with ReLU activations. For each weight Wi. we set its prior as Gaussian wo; ~ N N(Pur, Gw., and the response variable y N(o(1, {wo 1d o)) The model then have posterior of dl 392 (when the number of hidden units is 15) or d = 522 (when the number of hidden units is 20) and was applied to dataset sub-sampled from the wine-quality-reed dataset. The posterior approximated by deploying diagonal Gaussian distribution as the variational distribution. We compare the MLCFRG with MLMCRG, MLMC (with the number of levels being 2) and MC estimators. For MLCFRG, we use preconditioned
kernels (Oates et al. 2017); see details in Appendix C.1 The results, the average training ELBO/test log-likelihooo standard deviations, are illustrated in Figure 5 and Figure [
(averaged over runs) and associated For each run, the networks are initialized under identical settings. See Appendix D.5 for extra details and experiments. Empirical results show that the MLCFRG estimators tend outperform the other estimators they have shown faster convergence rate towards the optimal training ELBO and testing log-likelihood over number of iterations equivalently the number of training data.
-2000
-2000
-29900
-2000
-2000
15930-
-0000
-4000
-0000
-20000
MC
-25900
130 200
tendion
200
320
2000
302
300
2007
Figure 5: MLCF for Variation Inference of Figure 6: MLCF for Variation Inference of Bayesian Neural Networks: Hidden Dimension Size is 15 (num. of param. is 392).
Bayesian Neural Networks: Hidden Dimension Size is 20 (num. of param. is 522).
5 CONCLUSION
We introduced a generally applicable, flexible, and efficient method for estimating intractable integrals, MLCFs. MLCFs leverage the advantages of the multifidelity structure of integrands and CFs simultaneously, and thus: (i) provide reduced variance and faster convergence rates under mild conditions; (ii) are broadly applicable, e.g., suitable for complex and un un-ormaalized distributions, enabling their application in various areas; (iii) offer the flexibility be combined with experimental

---

Published as a conference paper at ICLR 2026
design techniques; (iv) have theoretically guaranteed variance bounds and optimal sample sizes across levels. MLCFs are evaluated across diverse scenarios, including Bayesian inference for the Lotka- Volterra system and variational inference for Bayesian neural networks, demonstrating the versatility of our method.
MLCFs also have several limitations. For instance, the computational cost scales cubically due to the inversion the kernel Gram matrix; the smoothness of the integrands and the density impacts the convergence rate; the method relies on the construction of hierarchy of fidelity levels and the coupling of coarse-fine level evaluations. There are a number of possible extensions. Firstly, on theoretical level, one could investigate the conditions on which MLCFs will lead to sufficient improvement in accuracy. Secondly, on an implementation level, future work will be needed to make MI.CERG more efficient. For instance, one can improve the empirical running speed MLCFRG by wrapping it an optimizer with hardware accelerated operators. Thirdly, on an algorithmic level, one could consider joint optimization of the way of sampling and the location of samples together with the cost minimization of MLCFs.
ACKNOWLEDGMENTS
Zhuo Sun is supported by Fundamental Research Funds for the Central Universities 20251 10590 of Shanghai University of Finance and Economics.
ETHICS STATEMENT AND REPRODUCIBILITY STATEMENT
This work raises no specific ethical concerns beyond standard practices in machine learning research. All methods, datasets, and hyperparameters described in detail, and core code is released in supplementary materials.
REFERENCES
Angelos Alexopoulos, Petros Dellaportas, and Michalis Titsias. Variance reduction for Metropolis- Hastings samplers. Statistics and Computing, 33(1):6, 2023.
Andreas Anastasiou, Alessandro Barp, François- François-Xavier Xavier Briol, Bruno Ebner, Robert Gaunt, Fatemeh Ghaderinezhad, Jackson Gorham, Arthur Gretton, Christophe Ley, Qiang Liu, et al. Stein's method meets computational statistics: review of some recent developments. Statistical Science, 38(1): 120-139, 2023.
Roland Assaraf and Michel Caffarel. Zero- Zero-variance principle for Monte Carlo algorithms. Physical Review Letters, 83(23):4682, 1999.
Aliye Atay-Kayis and Hélene Massam. A Monte Carlo method for computing the marginal likelihood in nondecomposable Gaussian graphical models. Biometrika, 92(2):317-335 2005.
Alexander Buchholz, Florian Wenzel, and Stephan Mandt. Quasi-Monte Carlo variational inference. In International Conference on Machine Learning, pp. 668 677. PMLR, 10-15 Jul 2018.
Bob Carpenter, Andrew Gelman, Matthew D Hoffman, Daniel Lee, Ben Goodrich, Michael Be- tancourt, Marcus Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell. Stan: A probabilistic programming language. Journal of Statistical Sofiware, 76(1), 2017.
Omar Chehab, Aapo Hyvarinen, and Andrej Risteski. Provable benefits of annealing for estimating normalizing constants: Importance Sampling, Noise-Contrastive Estimation, and beyond. Advances in Neural Information Processing Systems, 36:45945-45970, 2023.
Zonghao Chen, Masha Naslidnyk, Arthur Gretton, and Francois- Praver Briol. Conditional bayesian quadrature. In Uncertainty in Artificial Intelligence, pp. 648-684. PMLR, 2024.
K Dalbey, AK Patra, EB Pitman, MI Bursik, and MF Sheridan. Input uncertainty propagation methods and hazard mapping of geophysical mass flows. Journal of Geophysical Research: Solid Earth, 113(5)) 2008.
10

---

Published as a conference paper at ICLR 2026
Petros Dellaportas and Ioannis Kontoyiannis. Control variates for estimation based on reversible Markov chain Monte Carlo samplers. Joumal of the Royal Statistical Society Series B: Statistical Methodology, 74(1):133-611, 2012.
Hillary R Fairbanks, Alireza Doostan, Christian Ketelsen, and Gianluca laccarino. low-rank control variate for multilevel Monte Carlo simulation of high nigh-dimensional uncertain systems. Journal of Computational Physics, 341:121-139, 2017.
M Giselle Germádez-Godiin Review of multi-fidelltt models. Advances in Computational Science and Engineering. I(LLNL-JRNL-84799)) 2023.
Masahiro Fujisawa and Issei Sato. Multilevel monte carlo variational inference. Journal of Machine Learning Research, 22(278):1-44 2021.
Eric L Geist and Tom Parsons. Probabilistic analysis of tsunami hazards. Natural Hazards, 37(3): 277-314, 2006.
Gianluca Geraci, Michael S. Eldred, and Gianluca laccarino. A multifidelity multilevel monte carlo method for uncertainty propagation in aerospace applications. In 19th AIAA Non-Deterministic Approuches Conference, pp. 1951. American Institute of Aeronautics and Astronautics, 2017.
Michael B Giles. Multilevel Monte Carlo path simulation. Operations research, 56(3):607-617, 2008.
Michael B Giles. Multilevel Monte Carlo methods. Acta numerica, 24:259-328, 2015.
Paul Glasserman. Monte Carlo methods in financciaal engineering, volume 53. Springer. 2004.
Adam Golinski, Frank Wood, and Tom Rainforth. Amortized Monte Carlo integration. In Interna- tional Conference on Machine Learning, pp. 2309-2318, PMLR, 2019.
Charles Gordon Hewitt. The conservation of the wild life of Canada. New York: C. Scribner, 1921.
Fred J Hickernell, Christiane Lemieux, Art Owen, and Pierre L'Ecuyer. Control Variates for Quasi-Monte Carlo. Statistical Science, pp. 1-31, 2005.
Yuga Hikida, Ayush Bharti, Niall Jeffrey, and Trançois-Xavier Briol. Multilevel neural simulation- based inference. Advances in Neural Information Processing Systems, 2025.
Motonobu Kanagawa, Philipp Hennig, Dino Sejdinovic, and Bharath K Sriperumbudur. Gaus- sian processes and kernel methods: A review on connections and equivalences. arXiv preprint 1rX1v:1807.02582, 2018
Toni Karvonen. Error bounds and the asymptotic setting in kernel-based approximation. Dolomites Research Notes on Approximation, 15(RNNA Volume 15.3) 2022.
Kaiyu Li. Multilevel Methods for Monte Carlo Integration, with Applications to Tsunami Modelling. PhD thesis, UCL (University College London), 2024.
Kaiyu Li, Daniel Giles, Toni Karvonen, Serge Guillas, and François- Kaviorr Briol. Multilevel Bayesian quadrature. In International Conference on Artificial Intelligence and Statistics, pp. 1845-1868 PMLR, 2023.
Kaiyu Li, Dimitra Salmanidou. Devaraj Gopinathan, Mohammad Heidarzadch, and Serge Guillas. Uncertainty in manning's roughness coefficient in multilevel simulations of future tsunamis in sumatra. Proceedings of the Royal Society A, 48122366)20240637, 2025.
Jun S Liu. Monte Carlo strategies in scientific computing, volume 10. Springer, 2001.
Alfred J Lotka. Fluctuations in the abundance of species considered mathematically. Nature, 119
(2983):12-12, 1927.
Alfred James Lotka. Elements of physical biology. Williams & Wilkins, 1925.
11

---

Published as a conference paper at ICLR 2026
Jérôme Morio and Mathieu Balesdent. Estimation of rare event probabilities in complex aerospace and other systems: a practical approach. Woodhead publishing, 2015.
Fabio Nobile and Francesco Tesei. A multilevel Monte Carlo method with control variate for elliptic PDEs with log-normal coefficients. Stochastic Partial Differential Equations: Analysis and Computations, 3:398-444, 2015.
Chris J Oates, Mark Girolami, and Nicolas Chopin. Control functionals for Monte Carlo integration. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 79(3):695-71 2017.
Chris J Oates, Jon Cockayne, François-Xave Xavier Briol, and Mark Girolami. Convergence rates for class of estimators based on Stein's method. Bernoulli, 25(2):114--1155 2019.
Naoto Ohsaka and Tatsuya Matsuoka. On the (In)tractabilitt of computing normalizing constants for the product determinantal point processes. In International Conference on Machine Learning, pp. 7414-7423. PMLR, 2020.
John Paisley, David M Blei, and Michael I Jordan. Variational bayesian inference with stochastic search. In International Coference on International Conference on Machine Learning, pp. 1363- 1370, 2012.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al. Pytorch: An imperative style, high-performance deep learning library. Advances in Neural Information Processing Systems, 32, 2019.
Benjamin Peherstorfer, Karen Willcox, and Max Gunzburger. Survey of multifidelity methods in uncertainty propagation, inference, and optimization. SIAM Review, 60(3):550-591, 2018.
Christian P Robert, George Casella, and George Casella. Monte Carlo statistical methods, volume 2.
Springer, 1999.
DWO Rogers. Fifty years of monte carlo simulations for medical physics. Physics in Medicine & Biology, 51(3))R227, 2006.
Reuven Y Rubinstein and Dirk P Kroese. Simulation and the Monte Curlo method. John Wiley & Sons, 2016.
Carlos Sánchez- Linares, Mare de la Asunción, Manuel Castro, José M González- Vida, Jorge Macías, and Siddhartha Mishra. Uncertainty quantification in tsunami modeling using multi-level monte carlo finite volume method. Journal of Mathematics in Industry, 6(1):5, 2016.
Paul Scharnhorst, Emilio T Maddalena, Yuning Jiang, and Colin N Jones. Robust uncertainty bounds in reproducing kernel hilbert spaces: convex optimization approach. IEEE Transactions on Automatic Control, 68(5))2848-288) -2661, 2022.
Yuyang Shi and Rob Cornish. On multilevel Monte Carlo unbiased gradient estimation for deep latent variable models. In International Conference on Artificial Intelligence and Statistics, pp. 3925-3933. PMLR, 2021.
Mohammad Shoeybi, Mostofa Patwary, Raul Puri, Patrick LeGresley, Jared Casper, and Bryan Catan- zaro. Megatron-Ime Training multi multi-billion parameter language models using model parallelism. arXiv preprint 1rXiv:199.008053 2019.
Shijing Si, Chris Oates, Andrew B Duncan, Lawrence Carin, François-Xoe Xavier Briol, et al. Scalable control variates for Monte Carlo methods via stochastic optimization. In International Conference on Monte Carlo and Quasi-Mont Monte Carlo Methods Scientific Computing, pp. 205-221. Springer, 2022.
Yang Song and Diederik P Kingma. How to train your energy-based models. arXiv preprint arXiv:2101.03288, 2021.
Leah South, Toni Karvonen, Chris Nemeth, Mark Girolami, and Chris Oates. Semi-exact control functionals from Sard's method. Biometrika, 109(2)3311-67,, 2022.
12

---

Published as a conference paper at ICLR 2026
Zhuo Sun. Transfer learning in Monte Carlo Methods and Beyond. PhD thesis, UCL (University College London), 2023.
Zhuo Sun, Alessandro Barp, and François-Xavier Briol. Vector-valued control variates. In Interna- tional Conference on Machine Learning, pp. 32819-32846. PMLR, 2023a.
Zhuo Sun, Chris Oates, and François-Xavier Briol. Meta-learning control variates: Variance reduction with limited data. In Uncertainty in Artificial Intelligence, pp. 2047-2057. PMLR, 2023b.
Abdullah Tokmak, Christian Fiedler, Melanie Zeilinger, Sebastian Trimpe, and Johannes Köhler. Automatic nonlinear mpe approximation with closed-loop guarantees. IEEE Transactions Automatic Control, 2025a.
Abdullah Tokmak, Kiran G Krishnan, Thomas B Schön, and Dominik Baumann. Safe exploration in reproducing kernel hilbert spaces. In International Conference on Artificial Intelligncce and Statistics, pp. 784-792. PMLR, 2025b.
Arash Vahdat and Jan Kautz. NVAE: A deep hierarchical variational autoencoder. Advances in Neural Information Processing Systems, 33:1966:111 1996699--96799 2020.
Vito Volterra. Variazioni e flutuazioni del numero d'individui in specie animali conviventi, volume 2.
Societá anonima tipografica" Leonardo da Vinci", 1927.
Chong Wang, Xi Chen, Alexander J Smola, and Eric Xing. Variance reduction for stochastic gradient optimization. Advances in Neural Information Processing Systems, 26, 2013.
Holger Wendland. Scattered data approximation, volume 17. Cambridge university press, 2004.
Zhanxing Zhu, Ruosi Wan, and Mingjun Zhong. Neural control variates for variance reduction. ECML PKDD, pp. 533-547, 2019.
13

---

Published as a conference paper at ICLR 2026
Appendix
In Appendix A we provide relevant preliminaries such the expression of the Stein kernels. In
Appendix B we provide proofs of the theoretical results stated in the main text. In Appendix C we
provide more details on the implementation of the proposed MLCFs. In Appendix [ provide additional details for the experiments and extra experiments.
THE USE OF LARGE LANGUAGE MODELS (LLMs)
During the preparation of this manuscript, the authors used ChatGPT to polish the writing (e.g.. improving grammar, readability, and clarity). The content, technical contributions, and conclusions of the paper were developed entirely by the authors, who take full responsibility for all ideas and results presented.
A PRELIMINARIES
A.1 AN ALTERNATIVE LANGEVIN STEIN OPERATOR
The Langevin Stein operator can also be adapted to apply the derivative of twice differentiable scalar-valued functions u X' + R, in which case it is called the second-order Langevin Stein operator:
(1 Anu(x) Vou(z) V (())
where 4, -x V..
A.2 STEIN KERNELS
The Stein kernels from the first-order Langevin Stein operator have the following form,
(0( V. Vark(z, V log (x) Vor((E,z)) + V log Vzh(o,
vv (or) V Tog
(10)
where V. = (0000,,..,,,000.))
B PROOFS OF THEORETICAL RESULTS
B. PROOF OF PROPOSITION 3.1
Proof. The unbiasedness can be obtained by taking the expectation with respect to the distribution IT
of the ny random variables that constitute X7 for E {O..... L.J. Firstly, we have that
((()E)))))))))))))))(((()))))))) -
f(-1(x?)) - am = 0,
14

---

Published as a conference paper at ICLR 2026
due to the property of the Stein kernel kg that the Stein kernel ki satisfies fx ko (') for all x E X' OOtes et al. 2017: 2019). Then, we have
MLCF
E
m " [ft fi-1ll
1=0
E8
[(f(())) - fi-1(x/))
m
1=0
- fi-1 -
4
1 ({{(((((})) fx-1(x
11 m
1=0
- (P -1-
- fi-1(x9))
a
11lf fi-1l
1=0
I[[f].
B.2 PROOF OF THEOREM 3.2
Proof.
Following the proof of Theorem of (Oates et al. 2019) or Theorem 11.13 of (Wendland
2004 under assumptions A A1-7, there exists rf > 0 and h 0, for hy < hj,
(x(( fi-1 sx)2 S ri h2 llft fi-1lla,
for all x E X. Since X i.i at each level, combing the bound above, have
V x [ [[]]
Xg Ng X1
Σ
(f((sa)) fi-1(1(1))) - (s((x(,)) [[ss:)))
1=0 111 mi
1 VI fi-1 - 81]
mi
[[ - -i-- - 31)
II - fi-1 st
11 m1
[1 fi-1 st)
n1 m
IM-IM-
[1 - fi-1 - 11
(rthT 1lfi f1-111111
WA
T4 11-m
Under the assumption A8, and let = 4Th ri, we can then write
rimm, n/d IIf - fi- 11 )2
x xt [M
n11 my!
(rim) -n/d fr fi l1ee. )²
711 m
15

---

Published as a conference paper at ICLR 2026
B.3 PROOF OF THEOREM 3.3
Proof of Theurem 3.3 Suppose that < mi/nt p < 1, then m mi/p The sample sizes
Ng
MLLC
MLCF ,..., RL.
that minimize the variance of MLCF in Theorem 3.2 with the overall cost constraint
are
nu
MLCF
MLCF
argmin
A,my
/(n
mi)
s.t.
Σ 1=0
Cm==
=T==
"g.
1=0
where A₁ = (rilfi fi-ille, ja. For convenience, we solve mo
ML.CF
ML.CP and ML.CF can be
m.
easily obtained by taking n nM mi MLM.F// /p. The optimisation problem above can be solved by using Lagrange multipliers. For some À > 0, we define
FMaLCr(mo..... mL,
Σ
A,m
# -1-1
1(
1) A(T Σ Cume
1=0
Differentiating (mo.. A) with respect to mo...., m., and setting the equations equal to gives
2T (- 1)Am, 2/ -/- -1) - -1)-11) 0
AC1
- 1)
-7-22
mi (p(2r/d + 1)A1
for lE L and Σ Cpmp/p-
By plugging the first equation into the second, we get
XCr
-1-11) 1212
Σ Cr
10 T
p(2r/d+
1)Ar
1=0
T
Eor
(1
1 1 11) -1)
Cr
(-
/p
1)Ar
U=
Plugging this last expression for À into expression for my, get
MLCF -T
C1
-
2r422
Σ Cr
C₁
[ -11
27422
10
mi!
p(27/d +
1)A1
p(27 d
1)Ai
1=-0
(C)
-FrIza
Cr.
-2-12a
-1
=pT
Σ
Cr
A
Av
1==0
Ir
=pT(rillft - fi-1l1
C 3गकc8
Σ
C
2r+22
(rellfrr fv- )
P=0
forl E (0, .0. L.).
16

---

Published as a conference paper at ICLR 2026
B.4 DERIVATION OF EQUATION
7
According the Theorem 3.2 and Theorem 3.3
(nim, 1/d I fi-1 7899 )2
V Xg Xg [f] 1
14 mm
1=0
rimi
27/0-1 -1.f
51-111
1/p-1
r((p(rrlllfi - -1-111111
A Cj
2542d
2zg
11
fi 1111
1/p-1
rpp # R
2n ri
- 2d
2-
[|st
fi 1
1
1=0
With R =T Ek- Cj
22+42
2rad
(ril.ir - Sv ----- )F2
-1
and 2 maxt=
V X1 x [Mn.cr/]]
2r
r2p 25 7
Σ =0
Cj
2r 22
||
C
2F
[|
ft
H
74
= T
2rid
C
2rdd 1r
llft
f1-11 11-11117
740
2r
7!!
where A 1
27 27/d
1-p
P
B.5 PROOF OF PROPOSITION 3.4
Proof. Note that, the multilevel re-parametrized gradient (MLRG) at iteration T. can be written as,
MLRG C(Ar.) BiSan (x)) Σ (B)A, -s, (x)) E[JA, Sxn-1(x)
AD
MIRG L(Ar--)) (x) -Ax (x)
The multilevel re-parametrized control functional gradient (MLCFRG) iteration can be written as,
MLCCRGG ((AL) = fTmLCr(ff [facl
Σ
[ar - fx-al
1=0
(iellfx fn-.l) + Morlfaa fxumil
G (()x-1) + iIRRsfre fxemal
17

---

Published as a conference paper at ICLR 2026
Under stochastic gradient descent, the MLCFRG can then be simplified into the following update rule of X,
AL+1 dL -aL ( (()L-1) McrlSA, -fxc-il)
aL
MLCFRG C(AL-1) al (xterr(ss - fAi- (x)
aL-1
aL-1
A +
ar.
(Az X1-1) - on fNc-Il-
a1-1
Then, under the assumptions of Theorem 3.2 its variance is upper bounded by,
V VIV
[(\x.) a 2 Maz vC(Ai)
V aL
(Az XL-1) an fAc-1l
aL-1
=ViLCelfA - -ll (remy I SAL allme)2
nr m.
B.6 SIMPLIFIED MLCFs
In this section, we provide simplified version of MLCF estimators. The simplified MLCF estimator takes the form of,
[ = Σ
Litlfi fi-1l Σ (((X) - (1 Kg (Xi, Xi) -1
(1)
1=0
1=0
(11)
where X1 (G(.))).........)))))
B.7 OPTIMAL SAMPLE SIZE FOR MLMC
For completeness, we recall the optimal allocation of sample sizes for MLMC
MLMC
MLMC
70
n
Since MLMC is an unbiased estimator, the mean squared error (MSE) of MLMC is equal to its variance:
MSE((IMMM) = [[[ [I[])) = VILMMMM! [fll (ELLÍMLMC [f]] 2
Σ Vin
1=0
where Vi
MLMC are obtained by minimizing the variance of MLMC estimates under total computational cost budget T:
Thus, the optimal sample sizes ni
MLMC
m
MLMC 120
MLMC
argmin
Σ Ving
Σ P=0
Crny =
1=0
To solve this, introduce the Lagrange multipliers with A >> 0:
FMIMCCHo.-.. nL, A) Σ Vin X(T Σ Crne).
Differentiating (no RII, A) with respect to no...., nI, À and setting the derivatives equal yields expressions for and ny- Substituting the expression for À into that of ny gives the closed form expression:
V. C
MLMC = T my
Vv.Cr
for IE (O...., ....L),) L}.
18

---

Published as a conference paper at ICLR 2026
C IMPLEMENTATION DETAILS
In this section, provide implementation details which is helpful to use the proposed method in the main text. We provide kernels and their derivatives which are important for coding the Stein kernels. We then provide the approaches select the associated kernel hyperparameters.
C.1 KERNELS AND THEIR DERIVATIVES
Mateńn 2.5 Kernel
The Matefn 2.5 Kernel
((x,z')= = o2² v5v(0 (z) x²) 5(x x²)T (x x²))
v5y (x (x) x¹)
exp(
3X2
with amplitude o² and length-scale has derivatives given by
5
V5v (a - (2)
(a x²)
exp( v5 (x (3-'''((---))
Vok(x,)))
(x
x
x')(
+
3
12
A3
5
1
v5
x')(( x²)
v5 (x x2) x¹)
Vzr((z,z)) =
(x x')
+
eep(
3
12
A3
d (x 2²) x1))
502
1
v5y (x
x¹) (x x¹)
5
Vx Vzr(z,z'))
12
+
1
da
v5 "5V2==2''(zz x¹)
exp(
Cquared-Exponentian Kernel
The squared-exponential kernel (also known as Gaussian kernel)
|| x''13
k(xxx') exp(
2A
with lengthscale À > 0 has derivatives given by
(x x²)
Vak(z,z)
k(xxx')
(x x¹)
V Vxrk(z,2))
k( x'),
d
82
0 (x, xb))
x'7) k(z,
V. Verk(z,))
k(2, x¹) =
Σ
Oz' Ox Oxx
Ox's
1=1
j=1
(x; - 2²)
[ (x 2') (x) x²)
k(xx,')= =
k(2,z').
X2
Preconditioned Squared-Exponenteee Kernel A preconditioned squared-exponential kernel (Otess et al. 2017) is:
|| x x11
kxxxx))= (1+ a|(x113)(1 a|331113)
exp
2X²
19

---

Published as a conference paper at ICLR 2026
with lengthscale À > 0 and preconditioner parameter > O. This kernel has derivatives given by:
20xx
(x
x¹)
Vek(a,x') =
k(z,z',
1 + a: ||111|
X2
2nxt''
(x
(x-x') x') k(xxx'),
1 + a11111112
X2
d
02
d
[1 al|x11
-2azj
- (og 22)) z/)
V. Vork(z,
Σ
k(
Σ
or,
Ox'
12
j=1
-200
(2) - x
12 (xx,z') +
dzy k(xxx²)
+ ae a<121113
X2
3=1
-200,
(x3 =11]]
|
-2att
+ (B)-21)
Σ 3=1
12 x(),z') +
k(2
1 + a||2011
X2
1 + a 21128111
12
402xx Txd
20(x x²) Txx
2a(x xx)Tx
,(2,,0))
(1+
(1)++
21 + a: | | 11
X² + a a112:13))
(x - x') [x - x') J.
+ 12
1a
C.2 HYPER-PARAMER PARAMETER SELECTION
In this section, we discuss and present the way of selecting have hyperparameters which we will have to select. For example, the squared-exponente kernel will often have a length-scale or amplitude parameter, and these will have significant impact on the performance.
of kernels. Most kernels
We propose to select kernel hyperparameters through marginal likelihood objective by noticing the equivalence between the optimal Stein kernel-based control variates on the objectives in (Oates) et al. 2017 Sun al. 2023a) and the posterior mean of zero-mean Gaussian process model with covariance matrix ko (z. y) See (Oates et al. (017)) for a discussion the scalar-valued control variates case; see (Sun et al. 2023a) for discussion in the vector-valued control variates case. Denotes all kernel-hyperparametee by V v= (MI... }r), and we maximize the sum of the marginal likelihood,
Σ (Vi(x)) Si fi-1(xi)) (X)) (H,(X,, + Am (m(X))(((())))) (x)) 1=0
3¹ arg max-
-
+ log det(|tt (X,, XRiv) + Amml).
where (X,, X2; is a matrix with entries (b));;; (a zu.)): 2) with ky being a Stein reproducing kernel of the form in Equation (10) specialized to the distribution II with hyperparameter v. Equivalently, we can maximize the marginal likelihood for each fi fi-1 as follows,
11 = arg max (1(())) -t (K6, X(;K) + ) - (x))
+ log det (X,, X;; + Alm.l
for forlE E {1, L). We found that it performed well in our experiments. Note that the regularization parameter can also be selected through the marginal likelihood or cross validation. However, in practice are in an interpolation setting. Therefore, we choose À small as possible whilst still being large enough to guarantee numerically stable computation of the inverse of the matrix above.
Additionally, for kernels with hyperparameters (such as the gquared-exponent kernels) that can be regarded as 'lengthsccale' median heuristics are also be adopted as effective estimators. That is, using the median of the pairwise distances between all data points estimator for the length-scales.
20

---

Published as a conference paper at ICLR 2026
10
10"
Absolute Emor
10
10
Estimator
MLCF
MLMC
medium Budget T (s)
high
Figure 7: Illustration Example: Absolute integration error under a budget constraint (Y-axis log-
scale).
D ADDITIONAL DETAILS AND EXTRA EXPERIMENTS
In Section 4 the performance of MLCF being evaluated through empirical assessments. We used different probability distributions in these experiments including uniform, Gaussian and intractable posterior distributions. Although some the assumptions are not fulfilled in these experiments, we still use these examples to study the versatility of our method across variety of settings.
D.1 EXPERIMENTAL DETAILS FOR THE ILLUSTRATION EXAMPLE
Construction of Levels this example, we use the forward Fuler method to solve first-order ODE
-(xx==((a)) f(0)
[0.
The integral of interest is [[xx] = So f(z)dz. fo. f1 and f2 are constructed by using different step size h1 to approximate f(x) at grid points, and applying linear interpolation in between. MLCF, CF and MLMC are used estimate II[f] and the result from 50 replications is represented in Figure When using the same evaluation budget to estimate [ MLCF significantly outperforms MLMC. and CF that only evaluations the finest level in this example. The result is consistent with those shown in Figure !
At level we use the forward Euler method to approximate the solution fi(x) with step size hy by updating the value:
f(( f((x)) h=f((x;),
at {a, = ihg 1/b After solving f(;)) at grid point {z, = ihs 1/h we construct a continuous function f(()) by using linear interpolation
1//h
zi (f((#41) f((x)))
f((z) f((x.)
h4
for x E (x1,
Additional Details For the illustration example, we use step size ho 0.25, h 0.05 and h22 = 0.005. Since the computational time is very low for this example, we approximate the cost at level / by assuming C, is inversely proportional to the step size. The sample sizes used for MLCF, MLMC and CF under different budget constraints are presented in Table [] The example is implemented using a Matefn 2.5 kernel, with hyper-parametee tuned at each level. as discussed in Appendix C
21

---

Published as a conference paper at ICLR 2026
Table 1: Illustration example: Number of samples at level given budget constraint T.
T
2
CF
low
11
3
4
medium
14
high
23
17
5
8
Table 2: Synthetic Example: Number of samples level ( given budget constraint T.
setting
MMLC
MMMMCC
CF
T
1=
0
[ =
=
0
==1
1
1 2
I = 2
low
97
24
153
13
3
67
medium
145
36
10
229
4
101
high
193
48
14
305
26
5
134
D.2 EXPERIMENTAL DETAILS FOR THE SYNTHETIC EXAMPLE
Construction of Levels
In this synthetic example, fi(x) are defined
fo()) ko (Ir, =0);
f1 (x) aokk (x, 2n) + on k4 (I, 21); fz(x) = ank4 (z, 20) + ank (z, 21) + ark (r, z2)
where ki0 kt. and k2 are obtained by applying Stein operators to kin his and k2, respectively, and then adding positive constants co. C1, and cg.
For implementation, we set ao = 10, 0 3, 02 1, and choose 2o (0.1, 0.5), 21 0..7, 22 (0.1, 0.3). The amplitudes of k k4. k2 are 6, 4 and respectively. Their length-scales
vO.I, v0.2, and v0.4. The constants co. C1 and C2 are set to 1, 0.5 and 0.15.
Additional Details
To determine the optimal sample size, we compute norm of fi- fi-1 as follows
|| 1 = 2 (zo.20) & [|ʃ - sollig oi (11, & |l/2 f1 1x a² (z2, #2)-
Following the result of Theorem 3.3 we compute the optimal sample size MMLCF for MI.F when evaluation budget is limited. With the expression in Appendix B.7 we also compute the optimal
sample size PMLMC for MLMC. The sample size is listed in Table 2
D.3 EXPERIMENTAL DETAILS FOR THE ODE EXAMPLE
Construction of Levels We follow the settings described in Giles 2015 Li et al. 2023). For completeness, we provide details of the solver (finite difference approximation) The boundary-value ordinary differential equation (ODE) with random coefficient and random forcing is given by:
d (c(z)
du 50
for (0,1)
dz
dz (0)) x(1) =
where c(z) = =1 + 212. Expanding the equation gives
du
du
dz + (1 + 212)
50x² for E ().
d=
22

---

+ 2att))))))((((((((((( Published as a conference paper at ICLR 2026
Let u(2;) = u(th) for i fi,....(1- h)/h} with boundary conditions u(0) = =(1) 0. Using a finite difference approximation with step size 0, the left-hand side of the equation above is approximated as:
((2)) u(z; h) + (1 + x122)
a(2: + h) 2u(22) u12s h) 50x2
21
h
122
u(zi) u(22-1) (1 + ih) u u(2241) 2u(2;) + u(2;-1) 5022
21
h2
which, after rearrangement, can expressed as
21
h
haa
5022.
Incorporating the random coefficient and the random forcing, the level-/ approximation given by
1/b
fi(x) Σ
hru(ihr, x),
where u = au), (2h,, r).... ................................ u(1 h1, z)) solves the linear system
(a (i/Qi/ + L, (1/h))u: = 50x31.
Here 1 E R(1-he//h is vector of ones. The stiffness matrices Qi RRI1-h)/bxx(L-/))// and L E (B(-tt)//hsx((-aa)) hi)// are tridiagonal:
- -0)===(())--=-----1
and
=-2, (L)) (1))-1;;
Additional Details For the ODE example, we have evaluation cost at each level C (Co, C1 (1.22,3.57 11.89) (all measured in 10-3 seconds). Under the same evaluation cost constraint, compared (1) MLCF with Quasi-Monte Carlo points (QMC), (2) MLCF with Latin hypercube sampling (LHS), (3) MLCF with i.i.d points, (4) CF with i.i.d points, (5) MLMC with i..d points, (6) MLBQ with i.i.d points. The sample size is the optimal sample size for MLMC, which listed Table 3 In this example, we used squared-exper exponential kernels. The prper-parameters of the squared-exponential kernels are tuned independently each level as illustrated in Appendix C.2 The kernel hyper-parameters at each level are tuned by maximizing the marginal likelihood associated with each fi Лi-1-
Table 3: ODE example: Number of samples at level given budget constraint T',
T
=
0
[
1
[ =
2
CF
0.30
70
10
2
0.91s
209
31
5
45
1.52s
349
52
6
75
Extra Experiments: Effect of Kernels and Kernel Hyper-parameters We also study the effect of kernels for the ODE example. We use i.i.d samples for all settings in Figure & Under the same evaluation cost constraint, we compared (1) MLCF with the square exponential kernels (length-scale tuned by maximizing marginal likelihoods) (2) MLCF with the preconditioned square exponential kernels (length-scall tuned by maximizing marginal likelihoods) (3) MLCF with the square expo- nential kernels (length-scale selected by median heuristic) (4) MLMC (5) CF (length-scale tuned by maximizing marginal likelihoods). We still use the sample size listed in Table We find that setting the 'sength-scaae' of the square exponential kemnlss by median heuristic works well in practice.
23

---

Published as a conference paper at ICLR 2026
Estimator
MLGF SE
MLCF PSE
1.5
MH
Absolute
MLMC
1.0
CF
0.0
0.30
Budget T
Figure 8: Effect of Kernels and Kernel constraint (Y-axis log-scale). PSE stands for preconditioned square exponential kernels. SE stands for square exponential kernels. MH stands for median heuristic.
Absolute integration under budget
D.4 EXPERIMENTAL DETAILS FOR THE LOTKA- EXAMPLE
Construction of Levels The Lotka- Volterra system (Lotka 1925
1927
Voltema
1927) of ordinary
differential equations is given by:
du, (t)
21 21 (t) 1ru1 (t)
dwa(t) (2gu(())u2 (1) x4))(1),
de
for [0, s), with initial condition un (0) Iy and uy (0) Zg. To solve this system numerically, we use Stan (Carpenter et al., (017). Stan solves the ODE system for the times provided using the Dormand- Cormand-Prince algorithm, a 4th/5th order Runge-Kutte method. For the multilevel construction, we consider uniform time discretization with step size hy at level I. The grid points are {t; = ihy 1s/h,
3 At higher levels smaller values of hı are chosen. fi is then defined as fi (T) = (s) hu Σ Di (t.). which corresponds to an approximation of the time average of the prey population (t) over the interval (0, s).
Additional Details For the Lotka- Lotka-Volterra example, have sampling and evaluation cost each level C' = (Co, C1, C2) = (6.88, 34,41, 165.18) (all measured in 10- seconds). We use different step sizes for different levels. Under the same budget constraint, we compare (1) MLCF with MCMC points. (2) MLMC framework with MCMC points (MLMCMC), (3) CF with MCMC points, (4) MCMC. The sample size is listed in Table We used squared-exponential kernels this example, whose hyperparameters are tuned by maximizing marginal likelihood for fi fi-1 at each level.
Table 4: Lotka-Vou Volterra: Number of samples at level / given budget constraint T.
T
1=
0
[ -=
1 = 2
CF
MCMC
0.26 s
207
23
20
20
0.51s
413
47
40
40
0.77s
620
70
60
60
For completeness, recall the reparameterization of model parameters this section. The log-exp transform on model parameters as follows:
ce-parametrization and Priors
xy cxp((E)
i log(z;),
24

---

Published as a conference paper at ICLR 2026
for j in {1,.. 8), with Gaussian distribution priors: T4 ~ (0, 0.5²)
22 ~ (..3
T5. TB (log 10, 10)
37, By 19)
Then, the Lotka- Volterra system is
duy (t)
exp exp( exp($2) (1)
dt
duz(1) (t)
exp(@1)
dt
The observations are
(0) (ognormal cxxx(( exp(
(2( Lognormal exp (xp())))
(t) Lognormal(log(uu un (t)), exp(i7))
(2(t) Lognormal(log uy exppx))
D.5 EXPERIMENTAL DETAILS FOR VARIATIONAL INFERENCE OF BAYESIAN NEURAL NETWORKS
Construction of Levels In this case, the term level of MLCFs naturally corresponds to the term iteration of optimization of the Bayesian neural networks.
Additional Details To demonstrate the effectiveness of the proposed method for variational in- ference, we utilize two-layer Bayesian neural networks, with the middle hidden layer having some hidden units, and use ReLU the activation function. The prior of weights are set to be standard Gaussians N0, 1). We use Adam as the stochastic optimizer to training the variational parameters À with initial learning rate 10-4 and step-based decay function for the learning rate iteration I: 17/ gfoor(\/r) with drop parameter estimators, we use Monte Carlo samples to estimate the gradients during the optimization process while for MLCFRG estimators, only use 1 sample. The details of optimized training ELBO and test log- liglliiel are presented in Table and Table 6
250 and / 0.95. For MC, MLMC and MLMCRG
Table 5: Training ELBO and Test Log-likelihook (Hidden Dimension Size is 15)
Method MC 69779.5 (+1593.7)
Training ELBO
Test Log-llkeeihooo
-1683.71 (#552.1)
MLMC
-5214.09 (#949.3)
-1353.94 (+395.8)
MLMCRG
-2187.77 (=83.8)
-388. 42 (+0)5)
MLCERG
-2101.95 (+36.3)
-38.34 (+0.9)
Table 6: Training ELBO and Test Log-likelihood (Hidden Dimension Size is 20)
Method
Training ELBO
Test Lig-likel
MC
-3581.99 (57.
-771.04 (1553
MLMC
-3534.61
(+33.
-762.67 (+27.5)
MLMCRG
-2010.09
(=18.1)
-386.27 (=011)
MLCERG
-2014.54
(+20.2)
386.557 (+0.3)
Extra Experiments: Effect of Kernels We also study the effects kernels used in variational inference for BNNs. In Figure and Figure I0 we square-exponential kernels for the cases when the number of hidden states is 15 and 20, respectively. The hyperparameters kernels chosen by median heuristies. For MC, MLMC and MLMCRG estimators, estimate the gradients during the optimization process while for MLCFRG estimators, we only use sample.
5 Monte Carlo samples to
25

---

Published as a conference paper at ICLR 2026
-1000
-1000
-27/0
-2000
-20000
-3000
+000
-19900
+0000
MLMIC
MUNC
MLYI
MCRC
20930
Figure 9: MLCFRG (hidden dimension size is 15) Figure 10: MICFRG: (hidden dimension size is with square-exponential kernels.
20) with square-exponential kernels.
Extra Experiments: Effect of Sample Sizes We also study the effect of sample sizes: we test the effect of using the same strategy of sample sizes as in (Fujisawa & Sato, 2021). We used preconditioned squared-exponentaal kernels this example. The hyperparameters of kernels is
chosen by median heuristics. For both Figure IT and Figure T2 MLMCRG use the sample size strategy in (Fujisawa Sato (02T) with starting sample size no 5. Figure [T] MLCFRG use the same sample size strategy as MLMCRG. In Figure 12 sample size is fixed to be 1 for MLCFRG for iterations.
-20.0
-2000
115000
-0000
-0000
20030
200300
-0000
-0000
MEVERGG
MUNDRG MLCTIG
HLMERG
MUNDAG
MLLFFG
MILFAGG
MUCIEG
326 250 300 350
408
23E 356
300
53 300 353 250
Figure 11: MICFRG (Hidden Dimension Size is 5). Same strategy of sample sizes.
Figure 12: MLCFRG (Hidden Dimension Size 5). Sample size is fixed to be for MLCFRG.
26


Reviews:
[{"review_id": "FCHsqH0vCc", "reviewer": "Reviewer_8axS", "summary": "his paper proposes Multilevel Control Functionals (MLCFs), a novel variance reduction technique for efficiently estimating computationally expensive integrals. The core idea is an elegant combination of two powerful variance reduction strategies: Multilevel Monte Carlo (MLMC) and non-parametric, Stein-based Control Functionals (CFs).\n\nMLMC methods accelerate estimation by using a telescoping sum of a hierarchy of low-fidelity, cheap approximations to the expensive, high-fidelity integrand. The variance is reduced by estimating the small differences between levels.\n\nThis paper's key insight is to treat the MLMC difference terms themselves as integrands that can be further variance-reduced. The authors propose applying a non-parametric Stein-based Control Functional to each level of the MLMC estimator. This results in an estimator for the sum of these variance-reduced differences, which has significantly lower variance.\n\nThe authors provide a solid theoretical analysis, including a variance bound that demonstrates a faster convergence rate than standard MLMC, particularly for smooth integrands and densities in low-to-moderate dimensions. They also derive the optimal sample allocation across levels to minimize this variance bound under a fixed computational budget.\n\nFurthermore, the paper extends this framework to variational inference (VI) by proposing the Multilevel Control Functional Re-parameterized Gradient (MLCFRG) estimator. This estimator applies the MLCF idea to the multilevel gradient estimator for the ELBO (MLRG), and a practical, efficient recursive update is provided.\n\nThe method's effectiveness is demonstrated empirically on a synthetic example, a boundary-value ODE, a Bayesian inference problem (Lotka-Volterra system) using MCMC samples, and a variational inference problem (Bayesian Neural Network). In all cases, the proposed MLCF and MLCFRG methods outperform their respective baselines (MLMC, CF, MLMCRG).", "questions": "Hyperparameter Sensitivity: How sensitive is MLCF's performance to the choice of kernel and its hyperparameters? The paper uses several different kernels. How much of the performance gain is attributable to careful, and potentially expensive, hyperparameter tuning via marginal likelihood maximization at each level?\n\nRelation to Multi-Fidelity BQ: How does MLCF compare to other methods that combine multilevel/multi-fidelity approaches with kernel-based methods, such as the multilevel Bayesian Quadrature (Li et al., 2023) you cite? A direct empirical comparison, even on the ODE problem, would be very insightful to position MLCF in the literature.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Novelty and Significance: The paper's main contribution—hybridizing MLMC with Stein-based Control Functionals—is both novel and highly intuitive. It directly addresses a practical and important problem: the high computational cost of integration in many scientific and machine learning applications. This combination is powerful, as it leverages the \"divide-and-conquer\" strength of MLMC and the non-parametric variance-reducing power of CFs.\n\nTheoretical Soundness: The method is supported by a strong theoretical foundation. The paper provides a clear variance bound that explains why and when MLCF should be effective (dependence on smoothness and dimensionality). The derivation of the optimal sample allocation adds to the method's practical utility.", "weaknesses": "Computational Cost of Control Functionals: The primary weakness, acknowledged by the authors, is the cubic computational cost of inverting the kernel Gram matrix to construct the control functional, where the cost scales with the number of design points. This cost is incurred per level for the standard MLCF estimator and per iteration for the MLCFRG estimator. While the authors argue this is negligible if the integrand is sufficiently expensive, this scaling severely limits the number of points that can be used to build the CF, which in turn limits the achievable variance reduction.\n\nPractical Complexity: The method adds a new layer of complexity compared to the relatively simple MLMC. A user must now select an appropriate kernel (e.g., Mateŕn, SE) and its hyperparameters (e.g., length-scale) for each level. The paper suggests maximizing the marginal likelihood (in the Appendix), but this is a non-trivial, costly optimization problem in itself, adding to the overall computational burden.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "6gRYxWOMcQ", "reviewer": "Reviewer_u3Vy", "summary": "This paper introduces Multilevel Control Functionals (MLCFs), a novel variance reduction technique for Monte Carlo estimators. MLCFs extend traditional control variates by combining two powerful ideas: non-parametric Stein-based control variates and multi-fidelity methods . The authors provide theoretical analysis showing that MLCFs achieve faster convergence rates when integrands and densities are smooth and dimensionality is moderate. They validate MLCFs empirically on differential equation (DE) tasks (e.g., Bayesian inference for ecological models) and extend the framework to variational inference (VI), demonstrating improved performance on Bayesian neural network (BNN) benchmarks.", "questions": "see weaknesses", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1.The combination of Stein-based control variates with multi-level methods is innovative.\n2. The convergence rate analysis is a key strength, even though I did not validate the details of mathematical theory.\n3. The empirical evaluations are well-chosen and cover diverse use cases:", "weaknesses": "1. The abstract and theory section reference MLCFs working when “dimensionality is not very high,” but this is underspecified. What is the upper bound of dimensionality for MLCFs to remain competitive? \n2. Maybe more baseline methods (mentioned in related works) could be included in the experiments.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "qOZi0BBv9g", "reviewer": "Reviewer_NiZU", "summary": "This paper introduces a novel variance-reduction method for computing expectations using Markov chains and variational inference. The proposed approach, Multi-Level Control Functionals (MLCF), integrates multi-fidelity modeling with non-parametric control variates based on the kernelized Stein discrepancy. The authors prove that the estimator is unbiased, establish theoretical bounds on its variance, and derive the level-wise optimal sample allocation. The method is empirically validated on inference tasks involving dynamical systems and Bayesian neural networks.", "questions": "- Please provide the converged training ELBO and test log-likelihood values corresponding to Figures 5 and 6.  \n- At what dimensionality does the integrand cease to be considered “moderate”?  \n- What can be said about the relationship between the magnitude of $L$ and the extent of variance reduction?  \n- How does $L$ affect computational cost?  \n- How should practitioners choose or tune $L$ in practice?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The paper presents a novel and well-motivated approach to variance reduction, a key challenge in Bayesian inference.  \n- The method is theoretically rigorous, providing explicit variance bounds and optimal sample allocation formulas.  \n- The empirical results are convincing, demonstrating that MLCF consistently yields lower-variance estimators across the tested scenarios.", "weaknesses": "- The impact of varying the fidelity level $L$ is not explored, making it difficult to assess its practical importance.  \n  *(See related questions below.)*\n\n### Minor\n- In the BNN example, it is not immediately clear from the main text that variance reduction is applied to  \n  the gradient estimator of the ELBO. This is mentioned in the appendix; consider moving or referencing it earlier for ease of reading.  \n  (A single clarifying sentence would suffice.)", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)