Nested Expectations with Kernel Quadrature
Authors: Zonghao Chen
Masha Naslidnyk
Francois-Xavier Briol
Venue: ICML 2025 poster (2025)
Keywords: Kernel Quadrature
Monte Carlo
Kernel Ridge Regression
Primary area: general_machine_learning->kernel_methods
TL;DR: This paper proposes a new method to estimate nested expectation based on kernel quadrature.

Abstract: This paper considers the challenging computational task of estimating nested expectations. Existing algorithms, such as nested Monte Carlo or multilevel Monte Carlo, are known to be consistent but require a large number of samples at both inner and outer levels to converge. Instead, we propose a novel estimator consisting of nested kernel quadrature estimators and we prove that it has a faster convergence rate than all baseline methods when the integrands have sufficient smoothness. We then demonstrate empirically that our proposed method does indeed require the fewest number of samples to estimate nested expectations over a range of real-world application areas from Bayesian optimisation to option pricing and health economics.

Nested Expectations with Kernel Quadrature
Zonghao Chen
Masha Naslidnyk'
François-Xavier Briol
Abstract
et al., 2024), of acquisition functions in active learning and Bayesian optimisation (Ginshourger and Le Riche, 2010; Yang et al., 2024), of objectives in distributionalll--rbuus optimisation (Shapiro et al., 2023; Bariletto and Ho, 2024; Dellaporta et al., 2024), and of statistical divergences (Song et al., 2020; Kanagawa et al., 2023). Computing nested expectations is also key task beyond machine learning, including in fields ranging from value of information for decision making (Giles and Goda, 2019; Mala, 2024) fi- nance and insurance (Gordy and Juneja, 2010; Giles and Haji-Ali, 2019), manufacturing (Andradóttir and Glynn, 2016) and geology (Goda et al., 2018).
This paper considers the challenging computa- tional task of estimating nested expectations. Ex- isting algorithms, such as nested Monte Carlo or multilevel Monte Carlo, are known to be consis- tent but require large number of samples at both inner and outer levels to converge. Instead, we propose a novel estimator consisting of nested kernel quadrature estimators and we prove that it has faster convergence rate than all base- line methods when the integrands have sufficient smoothness. We then demonstrate empirically that our proposed method does indeed require fewer samples to estimate nested expectations on real-world applications including Bayesian opti- misation, option pricing, and health economics.
The estimation of nested expectations is particularly chal- lenging since there are two levels of intractability: the in- ner conditional expectation, and the outer expectation, both of which must be approximated accurately in order to ap- proximate the nested expectation / accurately. The most widely used algorithm for this problem is nested Monte Carlo (NMC) (Lee and Glynn, 2003; Hong and Juneja, 2009; Rainforth et al., 2018). It approximates the inner and outer expectations using Monte Carlo estimators with N and T samples respectively. NMC is consistent under mild conditions, but has relatively slow rate of convergence. Depending the regularity of the problem, existing results indicate that we require either ((A-3) or O(A evalu- ations of g to obtain a root mean squared error smaller or equal to A. This tends to be prohibitively expensive; for example, we would expect in the order of either 1 or 100 million observations to obtain an error of A 0.01. This is infeasible for many applications where obtaining samples evaluating g is expensive.
1. Introduction
We consider the computational task of estimating nested expectation, which is the expectation of a function that itself depends on another unknown conditional expecta- tion. More precisely, let Q be Borel probability measure with density on è and I's a Borel probability measure with density Po on X' S Rdxx which is parameterized by E C Rde, Given integrable functions f R + and g: X θ + I. we are interested in estimating:
I = Ea-Q If (EX~2a ()(X
(1)
1 / J g(x, ())(())d (0))d0.
inner conditional expectation
This issue has led to the development of number of meth- ods aiming to reduce the cost. Bartuska et al. (2023) pro- posed replacing the Monte Carlo estimators with quasi- Monte Carlo (QMC) (Dick et al., 2013). This algorithm, called nested QMC (NQMC), requires only O(A 2.5) function evaluations obtain error of size A (so that we only need in the order of 100, 000 observations for an error of A = 0.01). However, NQMC requires strong regu- larity assumptions which may not hold in practice (a mono- tone second and third derivative for f). Separately, Bujok et al. (2015); Giles and Haji-Ali (2019); Giles and Goda (2019) proposed to use multi-level Monte Carlo (MLMC) and showed that this can further reduce the number of func-
outer expectation
Nested expectations arise within a wide range of tasks, such as the computation of objectives in Bayesian experimen- tal design (Beck et al., 2020; Goda et al., 2020; Rainforth
'Department of Computer Science, University College London *Department of Statistical Science, University Col- lege London. Correspondence to: Zonghao Chen <zong- mao.chen.22@@ml uck.
Proceedings of the 1nnd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Nested Expectations with Kernel Quadrature
tion evaluations 2) (so that we only need in the order of 10, 000 observations for error of A 0.01). The algorithm has relatively mild assumptions on f and g. which makes it broadly applicable but sub- sub-optimal for ap- plications where f and g are smooth and where we might therefore expect further reductions in cost.
with h E (T)r,, and N independent and identically dis- tributed (i.i.d.) realisations æ1:N from n, standard MC ap- proximates the expectation of h under as follows:
N
1
Ex-= [b(x)]
~
Σ n=1
h(ra)-
N
To fill this gap in the literature. we propose a novel algo- rithm called nested kernel quadrature (NKQ), which is pre- sented in Section 3. NKQ replaces the inner and outer MC estimators of NMC with kernel quadrature (KQ) estimators (Sommariva and Vianello, 2006). We show in Section 4 that NKQ requires only Ō(A dxx' dia
For the nested expectation / in (I), the use of a MC esti- mator for both the inner and outer expectation leads to the nested Monte Carlo (NMC) estimator (Hong and Juneja, 2009; Rainforth et al., 2018) given by
1 T
N Σ
* function evalua- *x
ÍNMC = T T
Σ
f
1 N
g(=(9, 12
(2)
tions to guarantee an error smaller or equal to A. Here Ô denotes O up to logarithmic terms, 8x, se are constants relating to the smoothness of f and g in X' and θ and we have sx > dx/2 and se > de/2. In the least fa- vorable case, therefore recover the a) of NMC, but when the integrand is smooth and the dimension is not too large, we are able to have cost which scales better than O(A -2) and the method significantly outperforms all competitors. In those cases, we may only need in the order of a few hundred or thousands observations for an error of A = 0.01. This fast rate demonstrated numerically in Section 5, where we show that NKQ can provide signifi- cant accuracy gains in problems from Bayesian optimisa- tion to option pricing and health economics. Moreover, we show that NKQ can be combined with QMC and MLMC, providing an avenue to further accelcrate convergence.
where 01:r are T' i..d. realisations from Q and (t) i.i.d. realisations from Po, for each t E (1.....7) .T). The root mean-squared error of this estimator goes to zero at rate O(N ; +TT 5) when f is Lipschitz continuous (Rain- forth al., 2018). Hence, taking N T 0(A -) leads to algorithm which requires N x T O function evaluations to obtain error smaller or equal to A. When f has bounded second order derivatives, the root anean-squared error converges at the improved rate of O - - T- (Rainforth et al., 2018). Taking N = VT = 0(A+)) therefore leads to an algorithm requiring N x N 0(A-3) function evaluations get an error of A (Gordy and Juneja, 2010; Rainforth et al., 2018). De- spite its simplicity, NMC therefore requires a large number of evaluations to reach given A.
FNNN are N
2. Background
As a result, two extensions have been proposed. Firstly, Bartuska et al. (2023) proposed to use (2), but to replace the i.i.d. samples with QMC points. QMC points are points which aim to fill X' in somewhat uniform fashion (Dick et al., 2013), with well-known examples including Sobol or Halton sequences. Bartuska et al. (2023) used randomized QMC points, which removes the bias of standard QMC by using a randomized low discrepancy sequence (Owen, 2003). For nested expectations, they showed that nesting randomized QMC estimators can lead to a faster conver- gence rate and hence smaller cost of O(A -2.5). How- ever, the approach is only applicable when Po and Q Lebesgue measures on unit cubes (or smooth transforma- tions thereof), and the rate only holds when f has mono- tone second and third order derivatives.
Notation Let N+ denote the positive integers and N N, U (O}. For h X' S Rd + R, FI.N and h(TEN) are vectorized notation for INN] E R Nxd and (M<Z)), M(xx)] E R respectively. For a vector = (ai...., aa] T E Rd, define ||al|: = (E= 1/)2/0 For a distribution π supported on X' and 0 < P S cc.
(T) is the space of functions h- X' - R such that [| == -====x[((XX)]
< do and L (T) is the space of functions that are bounded r-almost everywhere. When is the Lebesgue measure Lx over X, we write (X)) := Lx(xx). For B E N, (2)) denotes the space of functions whose partial derivatives of up to and including order , are continuous. For two positive sequences {an IREN, and (bn JoEN,, an x bn means that lim, aa is a positive constant, On = Obb.) means that lim bn
be en < oo and
Alternatively, Bujok et al. (2015); Giles (2015); Giles and Haji-Ali (2019): Giles and Goda (2019) proposed to use multi-level Monte Carlo (MLMC), which decomposes the nested expectation using a telescoping sum the outer in- tegral, then approximates each term with MC. The inte- grand with the l'th fidelity level is constructed as the com- position of f with an inner MC estimator based on N, sam- ples. More precisely, the MLMC treatment of nested ex- pectations consist of using:
a, Öbr,) means that a, positive ro
b")"") for some
Existing Methods for Nested Expectations Standard Monte Carlo (MC) is an estimator which can be used to ap- proximate expeccationns/nterrals through samples (Robert and Casella, 2000). Given an arbitrary function h. X' - RR
2

---

Nested Expectations with Kernel Quadrature
IMLMC
Σ
1
Te
(f(Jez) (()e-1a)) + (Jo.,)
1
To
where (,)(.) Ex (x.)) -) is the kernel mean embed- ding (KME) of in the RKHS He (Smola et al., 2007). Minimizing the right hand side with additive regulariser term | over the choice of weights leads to the fol- lowing KQ estimator:
Tt.
To
1=1
N₁
where Je, =
1
Σ ( 0.) for l E (O.....I), L): (3)
NE
ikQ (x(INN) + +NNNN h(ZIN), M(ILN))(((.N)),
(5)
Under some regularity conditions, Theorem 1 from Giles (2015) shows that taking Nt- O(2") and Te = 0 leads an estimator requiring O(A-2) function evaluations to obtain root mean squared error smaller or equal to A. Although MLMC has the best known efficiency for nested expectations, Ni and Ti need to grow exponentially with 1, and we therefore need a very large sample size for its theoretical convergence rate to be- come evident in practice (Giles and Haji-Ali, 2019; Giles and Goda, 2019). MLMC also requires making several challenging design choices, including the coarsest level to use, and the number of samples per level. Most impor- tantly, MLMC as well as all existing methods fail to ac- count for the smoothness of the functions f and g.
where IN is the N N identity matrix, K ((NI:N,II:)) E RRNNN is the Gram matrix and X 2 0 is regularisation parameter ensuring the matrix is numer- ically invertible. The KQ weights are given by W1:N Pe (F((N) + NAIx) and are optimal when X = 0.
KQ takes into account the structural information that h E H4 so the absolute error |/ - I goes to as N - co. Specifically, when the RKHS H4 is norm- equivalent to the Sobolev space (x') (s > 2), KQ achieves the rate O(N- (Kanagawa and Hennig, 2019; Kanagawa et al., 2020). This is known to be minimax opti- mal (Novak, 2006; 2016), and significantly faster than the ((N rate of standard MC. Interestingly, existing proof techniques that obtain this rate take X = 0 in (5) and re- quire the Gram matrix K to be invertible, whilst the new proof technique based kernel ridge regression in this pa- per obtains the same optimal rate while allowing a positive regularization
fast rate
Kernel Quadrature Kernel quadrature (KQ) (Som- mariva and Vianello, 2006; Rasmussen and Ghahramani, 2002; Briol et al., 2019) provides an alternative to standard MC for (non-nested) expectations. Consider an arbitrary function : X' + R and distribution π on X, and sup- pose we would like to approximate Ex~- [h(X]] KQ is an estimator which can be used when h is sufficiently reg- ular, in the sense that it belongs to reproducing kernel Hilbert space (RKHS) (Berlinet and Thomas- 2004) Hy with kernel k. We recall that for a positive semi-definite kernel k: : X' x X + R. the RKHS H4 is a Hilbert space with inner product HH,, and norm || | (Aronszajn, 1950) such that: (i) k(z.-) E Hg for all x E X, and (ii) X,
N- _ (log N) 2:12 which improves nu- merical stability when inverting K. (See Remark B.2)
Despite the optimality of the KQ convergence rate, the rate constant can be reduced by selecting points I1:N other than through i.i.d. sampling. Strategics include impor- tance sampling (Bach, 2017: Briol et al., 2017), QMC point sets (Briol et al., 2019; R. Jagadeeswaran, 2019; Bharti et al., 2023; Kaarnioja et al., 2025), realisations from deter- minental point processes (Belhadji et al., 2019), point sets with symmetry properties (Karvonen and Särkkä, 2018; Karvonen et al., 2019) and adaptive designs (Osborne et al., 2012; Gunter et al., 2014 4; Briol et al., 2015; Gessner et al., 2020). Most relevant to our work is the combination of KQ with MLMC to improve accuracy in multifidelity settings (Li et al., 2023).
the reproducing property holds, i.e. for all h E Ha. h(z) = k(x, -)) An important example of RKHS is the Sobolev space W (2)) (s >> 5), which consists of func- tions of certain smoothness encoded through the square in- tegrability of their weak partial derivatives up to order s,
(x()) ( E L2(X): : Dd E L2(x) for all /3 E Nd with |8| S s}, E (4)
Two main drawbacks of KQ compared to MC are the worst- case computational cost of ((N3) (due to computation of the inverse of the Gram matrix), and the need for a closed- form expression of the KME j/x Fortunately, numerous approaches can mitigate these drawbacks. To reduce the cost, one can use geometric properties of the point set (Kar- vonen and Särkkä, 2018; Karvonen et al., 2019; Kuo et al., 2024), Nyström approximations (Hayakawa al., 2022; 2023), randomly pivoted Cholesky (Espeerly and Moreno, 2023), or the fast Fourier transform (Zeng et al., 2009), To obtain a closed-form KME, KQ users typically refer to existing derivations (see Table in Briol et al. (2019) or Wenger et al. (2021)), or use Stein reproducing kernels
where f denotes the B-th (weak) partial derivative of f.
Assuming h E H4 and Ex~rIVk(XX k(X, x) < co, the KQ
estimator Iko = ΣN=1 Wn ((,,)) uses weights obtained by minimizing an upper bound on the absolute error:
1 1--xe ikol
Ex [b(x]] -
Wn wnh ((
$ |1lllbcce
Σ w (zn--) Ilaa
3

---

Nested Expectations with Kernel Quadrature
(Oates et al., 2017: 2019; Si al., 2021; Sun al., 2023).
glus By)
0r
In this paper, tackle both drawbacks through change of variable trick. Suppose we can find a continuous trans- formation map & such that 21:N = ((u1:N) where u1:N are samples from a simpler distribution U of our choice. A direct application of change of variables theorem (Sec- tion 8.2 of Stirzaker (2003)) proves that Ex-= h(X)
F(O)
Fko(O0)
(2)
22
X
g .80
+ R to ho o U + R and the kernel quadrature estimator be- comes lko = pn (u(ui:)) (Ku + N (h 4) (uI:N), where Ku = kz (uI.N, U1:N). ). The measure U is typically chosen such that the KME is known in closed-form, and the KQ weights (u1:N) ) (Ku + NAIN) can be pre- computed and stored so that KQ becomes weighted aver- age of function evaluations with (IV) computational com- plexity. The main technical challenge of using the change of variable trick is to find such transform map $. See Ap- pendix F. for further details.
Ju ((((w)) (u), so the integrand changes from
0,
02 02
0y
(3)
Stage I
Stage II
Figure 1. Illustratioon of NKQ. In stage I, we estimate J(Ot) using (k) = > g",, 0e) for all t fl...., T'}. In stage II, we estimate I with ÍNKQ = Σ'= w. Fxq(Ot) where (0r) == (O))) The shaded areas depict Po (for stage and Q (for stage II).
(0) 0,) using KQ estimator: tions Fi.N and g(
(c)
Before concluding, we note that the KQ estimator is of- ten called Bayesian quadrature (BQ) (Diaconis, 1988; O'Hagan, 1991; Rasmussen and Ghahramani, 2002; Briol et al., 2019; Hennig al., 2022) since can be derived as the mean of the pushforward of a Gaussian measure h conditioned on ((IINN) (Kanagawa et al., 2018). The advantage of the Bayesian interpretation is that it provides finite-sample uncertainty quantification, and it also allows for efficient hyperparameter selection via empirical Bayes.
(k. == upe. (1 (K + NAxIN) "a 0.).
(7)
Here hx is a reproducing kernel on X. pro. ()) Ex-Pa, [Xx((X))] is the KME of Po. and K (t)
(x x ((xx) is an N x N Gram matrix. Using the same kernel hx for each t E fl...., T) is not essential, but we assume this to be the case for simplicity. Given these KQ estimates, we then we apply the function f to get (0)) = (()oo(0)).
3. Nested Kernel Quadrature
Stage We use a KQ estimator to approximate the outer expectation using the output of Stage
We can now present our novel algorithm: nested kernel quadrature (NKQ). To simplify the formulas, we write
İNKQ :00 (2(0 + T>elT) (ixo0)r) (8)
J(0) = Ex-E, lx(X, 0). F(0) = f(J(O), (6)
Here ke is a reproducing kernel on O, pg
so that the nested expectation in (1) can be written as E [F)].. We will assume that we have access
Ea-alke )] the embedding of Q and K (01:, 01T) ) a T x T Gram matrix.
OuT = 10....,
.i....,.0 0-] E
Combining stage I and II, NKQ can be expressed in single equation as a nesting of two quadrature rules:
(1) = 1 ./
[00] EX
Ti:N
İNKO Σ
T
Σ w (e B.)).
(9)
g ( g( g0 (Ii:N 0.) = (0(110 ..)......gl .......99(% 0.)] E R
(8)
(1)
for all t E (1...., T}. and f is a function that can be eval- uated. We do not specify how the point sets are generated, although further (mild) assumptions will be imposed for our theory in Section 4. Using the same number of fune- tion evaluations N per O, is not essential, but assume this as it significantly simplifies our notation. Given the above, we are now ready to define NKQ as the following two-stage algorithm, which is illustrated in Figure 1.
where wi,..... WN. are the KQ weights used in stage I for jko(O)) and w w are the KQ weights used in stage II. Although these weights are stage-wise optimal when Xx = Ny thanks the optimality of KQ weights, it is unclear whether they are globally optimal due the non- linearity of f. Note that NMC can be recovered by taking all stage I weights to be 1/N and all stage II weights to be 1/7', which is sub-optimal. In addition to the algorithmic simplicity of our proposed estimator NKQ, we demonstrate its superior performance in terms of both rate of conver- gence (Section 4) and empirical performances (Section 5).
Stage I For each E T),, we estimate the inner conditional expectation J evaluated at 0 with N observa-
4

---

Nested Expectations with Kernel Quadrature
NKQ inherits the two main drawbacks of KQ. Firstly, solv- ing the linear systems obtain the stage I and II weights has worst-case computational complexity of O(TN3 + T Secondly, NKQ requires closed-form KMEs at both stages: MFs, for all E (1,...,7) 7 in stage I, and don in stage II. Fortunately, we can often use the approaches dis-
kernel ridge regression. Interestingly, the proof in this pa- per leads to much better rate for CKQ than the best known rate for CBQ (see Remark 4.2).
4. Theoretical Results
cussed in the previous section to reduce the complexity to
In this section, we derive a convergence rate for the abso- lute error [λ - I] as the number of samples N, T + co. Before doing so, we recall the connection between RKHSs and Sobolev spaces. A kernel k on R is said to be trans- lation invariant if <xxxx') = Va - x²) for some positive definite function v whose Fourier transform (w) is a fi- nite non-negative measure on Rd (Wendland, 2004, The- orem 6.6). Suppose X' has a Lipschitz boundary, if k is translation invariant and its Fourier transform ((w) de- cays as O(1 + ||=1)) when w + oo for then its RKHS H4 is norm equivalent to the Sobolev space (x( (Wendland, 2004, Corollary 10.48). More specifically, it means that their set of functions coincide and there are constants cı, C2 > 0 such that ei [1 S || (g(()) I f2]lhl/H4 holds for all h E Hy. In this pa- per, we call such kernel a Sobolev reproducing kernel of smoothness 8. An important example of Sobolev kernel is the Matérn kernel- the RKHS of Matérn-u kernel is norm-equivalent to W5(x) with = 1 d/2. All Sobolev kemelss are bounded, i.e. supzex k(z, x) I s for some positive constant k. When the context is clear, we use |[flls,2 = [lflllz(x) to denote the Sobolev space norm.
O(TN- + T) and obtain closed-form kernel embeddings.
NKQ requires the selection of hyperparameters, including for the kernels in both stage I and II. We typically take kx and ke to be Matérn kernels whose orders are deter- mined by the smoothness of f and (as justified by The- orem I; see Section 4 for details). This leaves us with a choice of kernel hyperparameters which include length- scales 7x, T 0 and amplitudes Ax, Ae > 0. The lengthscales are selected via the median heuristic. The 25x ax (log N) ==ax""
d/2,
regularizers are set to Xx No #1g+2 "a (logT) da following Theo-
and Xe = To 28 rem 1, where A λo,e are selected with grid search over {0.01, 0.1, 1.0}. Finally, we standardise our function val- ues (by subtracting the empirical mean then dividing by the empirical standard deviation), and then set the amplitudes to Ax = Ae- = 1. This last choice could further be im- proved using a grid search, but we do not do this as we do not notice significant improvements when doing so in experiments and this tends to increase the cost.
Before presenting our theoretical results, we briefly com- ment on the connection with existing KQ methods. If we could evaluate the exact expression for the inner con- ditional expectation (0) pointwise, then (following (5)) the KQ estimator for would be IkQ = [o(01rr + IT) F (f1:r)) Comparing with (8), NKQ can thus be seen as KQ with noisy function values (ko((1eT) (re- placing the exact values F(f)r) in (8)). Although it is proved in Cai et al. (2023) that noisy observations make KQ converge at a slower rate, we prove that the stage II observation noise is of the same order as the stage I error, and consequently, we can still treat stage II KQ as noiseless kernel ridge regression and the additional error caused by the stage observation noise would be subsumed by the stage error (See Remark 4.1). NKQ is also closely related to family of regression- methods for estimating con- ditional expectations (Longstaff and Schwartz, 2001; Chen et al., 2024b). Indeed, with a slight modification of Stage II in (8), we can obtain an estimator of J(0) that we call conditional kernel quadrature (CKQ)
Theorem 1. Let X = [0, 1]dx and Q = [0, 1jdd, Suppose OLT are i..dd samples from Q and z are Li.d samples from Po, for all t E {1,.... .T}. Suppose furthrr that kx and ke are Sobolev kernels of smoothness 8x >> dx/2 and se > /2/22, and that the following conditions hold
(1) There exist Go.e, Gi,e.Goxx. GI,x > 0 such that Go I 9(0) I G1e and Go.x S po(x) S G1xx for any O E 2 and x E X.
(2) There exists S1 0 such that for any @ E O and any
B E Ndo with |3| S se, (D9g(=0)lsss S St.
(3) There exist S2, S3 > 0 such that for any x E X,
119(2, sll S S2 and ||0 + p=())llae.2 S S33 S 1.
(4) There exists S. > 0 such that derivatives of f up to and including order se + 1 are bounded by Sa.
Then, there exists No. To E N+ such that for N> No, T and Xo x to obtain the following bound
To we can take Ax x 2 (log N) 2xxt2 2xa+2
T 2 2 (log T)
(0k == (ee ( + TAgIT) (ko(O
Twke
(10)
CKQ highly resembles conditional BQ (CBQ) (Chen et al., 2024b); the difference is in stage II, where CBQ uses het- eroskedastic Gaussian process regression whilst CKQ uses
S T (GNN #x
(log
1 "x
c0 + C.T Te (log T)
1nt1
which holds with probability at least -T, C1, Cz are
5

---

Nested Expectations with Kernel Quadrature
two constants independent of T.,,
Method
Cost
Corollary 1. Suppose all assumptions in Theorem 1 hold.
If w[we set N O
and T O - ")),
de ), then N x
NMC
O 0(A
T OlA <x=
de x =e]
samples are sufficient to guarantee that I INKel= S A holds with high probability.
NQMC
O(A-2.6) 2..)
MLMC
O(A -)
To prove these results, we can decompose - INKol into the sum of stage I and stage II errors, which can be bounded by terms of order N *x (log.V) - and T (log T)
δ A dx de
NKQ (Corollary 1)
2011 11
Table 1. Cost of methods for nested expectations, measured through the number of function evaluations required to ensure | -- f < A. NMC rate taken from Theorem of Rain- forth et al. (2018), NQMC rate is taken from Proposition 4 of Bar- tuska et al. (2023), MLMC rate taken from Section 3.1 Giles (2018). Smaller exponents r in A -" indicate cheaper method.
respectively; see Appendix C. Interest- ingly, note that the stage II error does not suffer from the fact that we are using noisy observations (ko(0 and maintain the standard KQ rate up to logarithm terms (see Remark 4.1). We emphasize that our bound indicates that the tail of 1 - ÍNKQ] is sub-exponential. This contrasts with existing work on Monte Carlo methods, which only provides upper bounds on the expectation of error with no constraints on its tails (Giles, 2015; Bartuska et al., 2023).
same order as the stage I error Ō(N subsumed by it.
2x xx
) and is therefore
Remark 4,2 (Convergence rate for CKQ). For the CKQ es- timator JcxQ (defined in (10) that approximates the para- metric / conditional expectation J(]) uniformly over all O O the ervor can be upper bounded in the same way as 'x NKQ. The stage I error and can be shown to be O(N dx) using the same analysis from (C.13) to (C.16); and the stage error and can be shown to be Ö(T de using the same analysis from (C.18) to (C.34). Combining the two er- ror terms, we have || holds with probability at least 1 - ST- The rate is better than the best known rate ONN 4x +T of CBQ proved in Theorem 1 of (Chen et al., 2024b) since s0 The intuition behind the faster rate is that CKQ benefits from the extra flexibility of choose regularization parame- ters Ax, Xe; while CBQ, as two stage Gaussian Process based approach, is limited to choose Xe equal to the het- eroskedastic noise from the first stage. may be possible to modify the proof of (Chen et al., 2024b) to improve the rate further but this has not been explored date.
We now briefly discuss our assumptions. Assumption (1) is mild and allows (2(P)) (resp. (2()) to be norm equiva- lent to (x') (resp. L2(e)), which is widely used in statis- tical learning theory that involves Sobolev spaces (Fischer and Steinwart, 2020; Suzuki and Nitanda, 2021). Since our proof essentially translates quadrature error into gen- eralization error of kernel ridge regression, Assumptions (2), (3), (4) ensure that functions f.g and the density p have enough regularity so that the regression targets in both stage I and stage II belong to the correct Sobolev spaces. These are more restrictive, but are essential to obtain our fast rate and are common assumptions in the KQ literature. Assumptions (2), (3), (4) can be relaxed if mis-specification allowed; see e.g. Fischer and Steinwart (2020); Kana- gawa et al. (2023); Zhang et al. (2023). Theorem shows that for NKQ to have a fast convergence rate, one ought to use Sobolev kernels which are as smooth as possible in both stages. Furthermore, when sx = se = oo (c.g. when the integrand and kernels belong to Gaussian RKHSs), our proof could be modified to show an exponential rate of con- vergence in a similar fashion as Briol (2018, Theorem 10) or Karvonen al. (2020).
| - (g) = Ö(N
2x dx +T de) "M
7--
de
In Table 1, we compare the cost of all methods evaluated by the number of evaluations required to ensure /ī - /| < A. We can see that NKQ is the only method that explicitly exploits the smoothness of g. P. f in the problem so that it outperforms all other methods when dx de < 2.
Remark 4.1 (Noisy observations in Stage II of NKQ). Note that NKQ employs noisy observations (0r, Fx in stage II KQ rather than the ground truth observations {0r, F(O)) Although Cai et al. (2023) establishes that KQ with noisy observations converges at a slower rate than KQ with noiseless obser- vations, a key distinction in our setting is that, as shown in (C.30), the observation noise in stage I1 KQ is of order Ö(N
Remark 4.3 (Combine NKQ with MLMC and QMC). We have previously mentioned that KQ could potentially be combined with other algorithms to further improve effi- ciency, and studied this for both MIMC and QMC. For the former (i.e. NKQ+MLMC), we derived a new method called multi-level NKQ (MLKQ), which closely related to multilevel BQ algorithm of Li et al. (2023) and for which we were able to prove a rate of Ō(A ilarly to NKọ, when dx sx + de
"x dx ). whereas the noise in Cai et al. (2023) remains at a constant level. As a result, we can still use KQ in stage II as if the observations (O) are noiseless, and the additional error it introduces happens to be of the
2ne Sim- < 2, the rate for MIKQ a
6

---

Nested Expectations with Kernel Quadrature
is faster than that of NMC, NQMC and MLMC. However, the rate we managed prove is slower than that for NKQ. and a slower convergence was also observed empirically (see Figure 6). We speculate that the worse performance is caused by the accumulation of bias from the KQ estimators at each level. See Appendix D.2 for details.
NMC (N - 27
NMC (N = vT)
NKQ (N 7)
NKQ (N vT)
Synthetio Experiment
10
10-
10
We also consider combining NKQ and QMC. In this case, we expect the same rate as in Theorem can be recovered by resorting to the fill distance technique in scatter data approximation (Wendland, 2004). This is confirmed em- pirically in Section 5, where we observe that using QMC points can achieve similar or even better performance than NKQ with i.i.d. samples.
10
10
10
10a
10
Cost-
N x
d=1
d = 10
d 20
10a
10-
5. Experiments
10-
We now illustrate NKQ over range of applications, in- cluding some where the theory does not hold but where we still observe significant gains in accuracy. The code to reproduce all experiments is available at https:// gi htwpa/////uucccccccenness.t
10
10
10'
10" 10"
Cost
Figure 2. Synthetic experiment. Top: Verification of theoretical results. The thin grey lines are theoretical rates of A Cost"* and A Cost 1/3 Bottom: Comparison of NKQ and NMC as dimension increases. Results are averaged over 1000 inddpnn- dent runs, while shaded regions give the 25%-75% quantiles.
Synthetic Experiment We start by verifying the bound in Theorem 1 using the following synthetic example: Q = U[0 1], P. = U[0 1], g(x, 0) = अ² + 0 and f()) = 2 in which case = 0.4115 can be computed analytically. We estimate 7 with i.i.d. samples 01. U|0, 1] and i.i.d. samples x=(2) ~ U[0, 1] for t E f1....,TY T}. The assump- tions from Theorem 1 are satisfied with sx = 8e = 2 and dx= = de- =1 1 (see Appendix F.2). Therefore, to reach the absolute error threshold A, we choose N T A -0.5 for NKQ following Corollary 1. On the other hand, based on Theorem 3 of Rainforth et al. (2018), the optimal way of assigning samples for NMC is to choose N = VT A-1
the kernel lengthscale using the median heuristic provides very good performance. In the right-most plot, we see that NKQ with Matérn kernels outperforms Matérn- 2 kernel, indicating practitioners should use Sobolev kernels with the highest order of smoothness permissible by Theorem 1.
Risk Management in Finance We now move beyond synthetic examples, starting in finance. Financial in- stitutions often face the challenge of estimating the ex- pected loss of their portfolios in the presence of poten- tial economic shocks, which amounts to numerically solv- ing stochastic differential equations (SDEs) over long time horizons (Achdou and Pironneau, 2005). Given the high cost of such simulations, data-efficient methods like NKQ are particularly desirable.
In Figure 2 Top, we see that the optimal choice of N and T suggested by the theory indeed results in a faster rate of convergence for both NMC and NKQ. For this synthetic problem, we confirm that both the theoretical rates of NKQ (Cost- 1 and NMC (Cost- 3- from Theorem and Rainforth et al. (2018, Theorem 3) are indeed realized. We also adapt the synthetic problem to higher dimensions (dx = de= = d) in (F.42) and observe in Figure 2 Bottom that the performance gap between NKQ and NMC closes down as dimension grows. Such behaviour is expected be- cause the cost of NKQ Ö(A 4x * - and therefore de- grades as the dimensions dx and de increase; whilst the cost of NMC remains the same.
Suppose shock occurs time 1 and alters the price of an asset by factor of + for some s 2 0. Conditioned on the asset price S(n) = 0 at the time of shock, the loss of an option associated with that asset at maturity < with price S(C) = x can be expressed as [(6) = Ex-F, (x))),, where g(x) = (xx) - (1 + 8)2) measures the shortfall in option payoff and the distribution P. is induced by the price of the asset which is described by the Black-Scholes formula. The payoff function we consider is that of but- terfly call option: (xx) = max - K1, 0) + max K2, 0) - 2 max(c (K1 + 2,, 0) for K1, K2 2 O. Since we incur a loss only if the final shortfall is positive, the expected loss of the option at maturity can be expressed as = Eo-q [max(Ex [x() 0)]. Under this setting,
We also conduct ablation studies, which are reserved for Figure 5 in the appendix. In the left plot, we see that the result are not too sensitive to Ag. although very large values decrease accuracy whilst very small values cause numerical issues. In the middle plot, we see that selecting
7

---

Nested Expectations with Kernel Quadrature
NMC NKQ
NMC (QMC)
MLMC
adjusted life-years) under treatment in set of possible treatments C, and @ represent the additional variables that may be measured. Then, J0) Ex-P, [ge (X, 0)] repre- sents the expected patient outcome given the measurement of 0. The EVPPI is defined as / = 1 max (72,c), where 4 = L [maxxec Jc(0)] and 1z.c = Eeا [J) and therefore consists of [C] + 1 nested expectations.
NKQ (QMC)
Risk Management in Finance
10
20-1
We follow Section 4.2 of Giles and Goda (2019), where both Po and Q are Gaussians, and 9ı (2 0) 10' + x72834)) - (œ1 + x223344) and ga (2,, 0) (0g + x15316617) - (I10 + x11 The ex- act practical meanings of each dimension of x and 0 can be found in Appendix F.4, but includes quantities such as 'cost of treatment' and 'duration of side effects'. Here we have dx 17 and de = 2, the former being relatively high dimensional. The ground truth EVPPI under this setting is I == 538 provided in Giles and Goda (2019).
10
Cost
Health Eomomins
10
Cost
Figure 3. Top: Financial risk management. Bottom: Health eco- nomics. Results are averaged over 100 independent runs, while shaded regions give the 25%-75% quantiles.
For estimating both L and 22... Assumptions (2)(3) are sat- isfied with infinite smoothness s = se = oo, but the max function in / is only in (RR which violates Assump- tion (4). As a result, for estimating I we take hxx to be Gaussian kernel and ke to be Matérn- kernel (so as to be conservative about the smoothness in 0). For estimating 2e. we select both kx and lig to be Gaussian kernels. For NKQ, we choose N = T = A-1 whereas for NMC, we choose N = VT = A 1 For MLMC, we use I. = 5 levels and allocate the samples at each level following Giles and Goda (2019). We run NKQ and NMC with both i.i.d. sam- ples and QMC samples. In Figure Bottom, we present the mean absolute error of NKQ, NMC and MLMC with increasing cost. We can see that NKQ consistently outper- forms other baselines.
de dx = and 3.077 can be computed analytically.
In this experiment, Assumptions (2)(3) are satisfied with se 5x 1, but the max function is not in (IR) which violates Assumption (4) (see Appendix F.3). Nevertheless, we still run NKQ with xix and ke being Matérn- 1
2 kernels and choose N= = T = A1 for NKQ following Corol- lary 1. For NMC, we follow Gordy and Juneja (2010) and choose N vT = A -1 For MLMC, we use L lev- els and allocate samples at each level following Giles and Goda (2019).
In Figure 3 Top, we present the mean absolute error of NKQ. NMC and MLMC with increasing cost. We see that NKQ outperforms both NMC and MLMC as expected. For each method, we obtain the empirical rate r by linear re- gression in log-log space, and compare this against the the- oretical rate in Table 1. For NMC, our estimate of f 2.97 matches theory (r- = 3), but when using QMC samples in- stead, our estimate of = 2.74 shows we under-perform compared to the theoretical rate (r = 2.5). This is likely because the domains are unbounded and the measures are not uniform, breaking key assumptions. Finally, for NKQ, we obtain
Bayesian Optimization We conclude with an applica- tion in Bayesian optimization. Typical acquisition func- tions are greedy approaches that maximize the immediate reward, while look-ahead acquisition functions optimize accumulated reward over a planning horizon, which results in reduced number of required function evaluations (Gins- bourger and Le Riche, 2010; González et al., 2016; Wu and Frazier, 2019; Yang et al., 2024). The utility of a two-step look ahead acquisition functions can be written as the fol- lowing nested expectation.
1.90 for i.i.d samples and f 1.91 for QMC samples which match (and even slightly outperform) the theoretical rate (r = 2).
(=;; D) := Es,p [9(fr..) + max
[g (frDrs ]]]]
2'
where f.p. are the posterior distributions given data D and D' : D u (z, f(D(2)). In this experiment, the prior is Gaussian process with zero mean and Matérn-0.5 co- variance so the posterior SID remain a Gaussian process. The initial starting data Do consists of 2 points sampled uniformly from a prespecified interval. Here, g is the re- ward function and we use q q-expected improvement (Wang et al., 2020) with q = 2 so z = (21, 22) and g( 2) =
Value of Information for Healthcare Decision Making In medical decision-making, a key metric to evaluate the cost-benefit trade-off of conducting additional tests on pa- tients is the expected value of partial perfect information (EVPPI) (Brennan et al., 2007; Heath et al., 2017). For- mally, let Ie denote the patient outcome (such as quality-
8

---

Nested Expectations with Kernel Quadrature
Dropwace
Adday
Cosine8
10
10a
NKQ NMC
NMSE
10-
MLMC
10-1
10
0
20
50
50
60
Suwcends
Scoonds
Secoads
Figure 4. Bayesian optimization with look ahead funcciontoonuuncionnn are averaged over 100 independent runs, while shaded regions give the 25%-75% quantiles.
function. The plots are NMSE against accumulative wall clock time. Results
(apx=12(ffp(zs) - m, ,, 0). The constant maax is the maximum reward obtained from previous queries. Al- though JID (resp. JID²) is a Gaussian process, we only ever consider its evaluation on z (resp. z'), and we there- fore only have integrate against two-dimensional Gaus- sians. Notationally speaking, fipr(23 22) correspond to xr and (11, #2) correspond to 0 in (1) (i.e. dx = de 2), but we use the notation of fyzr, Syp to stay consistent with the GP literature. As a result of the max operation, 8xx = but we do not have sufficient smoothness in θ.
future work.
6. Conclusion
This paper introduces a novel estimator for nested expec- tations based on kernel quadrature. We prove in Theo- rem that our method has a faster rate of convergence than existing methods provided that the problem has sufficient smoothness. This theoretical result is consistent with the empirical evidence in several numerical experiments. Ad- ditionally, even when the problem not as smooth as the theory requires, NKQ can still outperform baseline meth- ods potentially due the use of non-equal weights.
We benchmark NKQ, NMC and MLMC on three synthetic tasks from BoTorch (Balandat al., 2020). For NKQ, both xx and ke are M kernels since we want to be conservative about the smoothness. Although both Q and Pa are Gaussian so closed-form KMES are available, we use the "change of variable trick" which maps Gaus- sian distributions to two uniform distributions over [0, 1]d (see Appendix F.5) to reduce the computational complexity of NKQ to (T N). To reach specific error threshold A 0.01, following Table 1, we choose N T A for NMC and N = T = A for NKQ. For MLMC, we use the same code as Yang et al. (2024). The normalized mean squared error (NMSE) max.eDD (us max.ee2 Jun(() Juu is used as performance metrie, where Do (resp. Ds) is queried data at initialization (resp. after S iterations), Jan is the black box function to be optimized and (a=( is the maximum reward.
Following our work, there remain number of interest- ing future problems and we now highlight two main ones. Firstly, we propose a combination of KQ and MLMC that we call MLKQ in Appendix D.2. However, we believe our current theoretical rate for MILQQ is sub-optimal due to the sub-optimal allocation of samples at each level. Further work will therefore be needed to determine whether this is viable approach in some cases. Secondly, for applica- tions where function evaluations are extremely expensive, NKQ could be extended to its Bayesian counterpart. This would allow us to use the finite sample uncertainty quan- tification for adaptive selection of samples, which could further improve performance. Finally, establishing mini- max lower bounds for nested expectation remains an open and compelling problem. The main difficulty lies in its two-stage structure. To the best of our knowledge, exist- ing minimax lower bounds for two-stage problem typically reduce the problem to a one-stage problem before deriving the bound; see, for example, Chen and Reiss (2011, Chap- ter 3), Meunier et al. (2024, Appendix F.1). and Zhang al. (2025). However, remains unclear how to directly estab- lish meaningful minimax lower bounds of genuinely two- stage problems, such as our nested expectation.
In Figure 4, we compare the efficiency of each method by plotting their NMSE against cumulative computational time in wall clock. We can see that NKQ achieves the lowest NMSE among all methods under a fixed amount of computational time in all three datasets, even though the assumptions of Theorem are not all satisfied. Since the Dropwave, Ackley, and Cosine8 functions are syn- thetic and computationally cheap (see Appendix F.5), we expect the advantages of NKQ be more pronounced for Bayesian optimization on real-world expensive problems. Furthermore, many other utility functions in Bayesian optimization-suu as predictive entropy cearch-involve nested expectations (Balandat et al., 2020). We leave the empirical evaluation of NKQ on these utility functions to
Impact Statement
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here.

---

Nested Expectations with Kernel Quadrature
Acknowledgments
stochastic collocation methods with importance sam- pling for Bayesian optimal experimental design. Inter- national Journal for Numerical Methods in Engineering, 121(15))3482-335) 2020.
The authors acknowledge useful discussions with Philipp Hennig and support from the Engineering and Physi- cal Sciences Research Council (ESPRC) through grants [EP/S021566/1] (for ZC and MN) and [EP/Y022300/1] (for FXB).
Ali Behzadan and Michael Holst. Multiplication in Sobolev spaces, revisited. Arkiv för Matematik, 59(2):275-32) -366, 2021.
References
Ayoub Belhadji, Rémi Bardenet, and Pierre Chainais. Ker- nel quadrature with DPPs. Advances in Neural Informa- tion Processing Systems, 32, 2019.
Yves Achdou and Olivier Pironneau. Computational meth- ods for option pricing. SIAM, 2005.
Robert A Adams and John JF Fournier. Sobolev spaces. Elsevier, 2003.
Colin Bennett and Robert C Sharpley. Interpolation of op- erators. Academic press, 1988.
Aurélien Alfonsi, Adel Cherchali, and Jose Arturo Infante Acevedo. Multilevel Monte Carlo for computing the SCR with the standard formula and other stress tests. Insurance: Mathematics and Economics, 100:234-260, 2021.
Alain Berlinet and Christine Thomas-Agnan. Reproduc- ing Kernel Hilbert Spaces in Probability and Statistics. Springer Science+Business Media, New York, 2004.
Ayush Bharti, Masha Naslidnyk, Oscar Key, Samuel Kaski, and François- Xavier Briol. Optimally-weightted estima- tors of the maximum mean discrepancy for likelihood- free inference. In International Conference on Machine Learning, pages
Aurélien Alfonsi, Bernard Lapeyre, and Jérôme Lelong. How many inner simulations to compute conditional expectations with least-square Monte Carlo? arXiv preprint wrXhv:2209.041553 2022.
2023.
Alan Brennan, Samer Kharroubi, Anthony O'Hagan, and Jim Chilcott. Calculating partial expected value of per- feet information via Monte Carlo sampling algorithms. Medical Decision Making, 27(():448 470, 2007.
Sigrun Andradóttir and Peter W. Glynn. Computing Bayesian means using simulation. ACM Transactions on Modeling and Computer Simulation, 26(2):1-26 2016.
Nachman Aronszajn. Theory of reproducing kernels. Transactions of the American mathematical society, 68 (3):337-404, 404, 1950.
François-Xavéas Briol. Statistical computation with ker- nels. PhD thesis, University of Warwick, 2018.
François-Xavier Briol, Chris Oates, Mark Girolami, and Michael Osborne. Frank- Wolfe Bayesian quadrature: Probabilistic integration with theoretical guarantees. In Advances in Neural Information Processing Systems, pages 1162-1170, 2015.
Francis Bach. On the equivalence between kernel quadra- ture rules and random feature expansions. Journal of Machine Learning Research, 18(19):1-38, 2017.
Maximilian Balandat, Brian Karrer, Daniel Jiang, Samuel Daulton, Ben Letham, Andrew G Wilson, and Eytan Bakshy. Botorch: A framework for efficient Monte-carlo Bayesian optimization. Advances in Neural Information Processing Systems, 33:21524 525321,38, 2020.
François-Xave -avier Briol, Chris Oates, Jon Cockayne, Ye Chen, and Mark Girolami. On the sampling prob- lem for kernel quadrature. Proceedings of the Inter- national Conference on Machine Learning, pages 586- 595, 2017.
Nicola Bariletto and Nhat Ho. Bayesian Nonparametrics meets Data-Driven Robust Optimization. Advances in Neural Information Processing Systems, 2024.
François-Xaa Xavier Briol, Chris Oates, Mark Girolami, M. A. Osborne, and D. Sejdinovic. Probabilistic integration: A role in statistical computation? (with discussion). Statis- tical Science, 34(1):1-22, 2019.
Arved Bartuska, Andre Gustavo Carlon, Luis Espath, Se- bastian Krumscheid, and Raul Tempone. Double-loop randomized Quasi-Monte Carlo estimator for nested in- tegration. 19X1v.2302.14119, 2023.
Karolina Bujok, Ben M Hambly, and Christoph Reisinger. Multilevel simulation of functionals of bernoulli random variables with application to basket credit derivatives. Methodology and Computing in Applied Probability, 17: 579-604, 2015.
Joakim Beck, Ben Mansour Dia, Luis Espath, and Raúl Tempone. Multilevel double loop Monte Carlo and
10

---

Nested Expectations with Kernel Quadrature
X. Cai, C. T. Lam, and J. Scarletr. On average-case error bounds for kernel-based Bayesian quadrature. Transac- tions Machine Learning Research, 2023.
Michael B. Giles. MLMC for nested expectations. Con- temporary computational mathematics-4 celebration of the 80th birthday of ian sloan, pages 425-442, 2018.
Xiaohong Chen and Markus Reiss. On rate optimality for ill-posed inverse problems in econometrics. Econo- metric Theory, 27(3):497-521, 2011. ISSN 0266- 4666, 1469-4360. URL http:/ /ttp://tttstto.g////////// stable/27975490.
Michael B. Giles and Takashi Goda. Decision Decision-making der uncertainty: using MLMC for efficient estimation of EVPPI. Statistics and computing, 29:739-751, 2019.
Michael B. Giles and A-L Haji-Ali. Multilevel nested sim- ulation for efficient risk estimation. SIAM/ASA Journal on Uncertainty Quantification, 7(2):497-525,
Zonghao Chen, Aratrika Mustafi, Pierre Glaser, Anna Ko- rba, Arthur Gretton, and Bharath K Sriperumbudur. (de)- regularized maximum mean discrepancy gradient flow. arXiv preprint arXiv:2409. 14980, 2024a.
David Ginsbourger and Rodolphe Le Riche. Towards Gaussian process-based optimization with finite time horizon. In mODa 9 Advances in Model-Oriented De- sign and Analysis, pages 89-96, 2010.
Zonghao Chen, Masha Naslidnyk, Arthur Gretton, and François-Xavier Briol. Conditional Bayesian quadra- ture. Conference on Uncertainty in Artificial Intelli- gence, 2024b.
Takashi Goda, Daisuke Murakami, Kei Tanaka, and Kozo Sato. Decision-thcoretic sensitivity analysis for reservoir development under uncertainty using multilevel Quasi- Monte Carlo methods. Computational Geosciences, 22 (4):1009-1020, 2018.
Charita Dellaporta, Patrick O'Hara, and Theodoros Damoulas. Decision making under the exponential fam- ily: Distributionally robust optimisation with Bayesian ambiguity sets. arXiv:2411.168291 2024.
Takashi Goda, Tomohiko Hironaka, and Takeru Iwamoto. Multilevel Monte Carlo estimation of expected informa- tion gains. Stochastic Analysis and Applications, 38(4): 581-600, 2020.
Persi Diaconis. Bayesian numerical analysis. Statistical Decision Theory and Related Topies IV, pages 163-175,
1988.
Josef Dick, Frances Kuo, and Ian H. Sloan. High- dimensional integration: The Quasi-Monte Carlo way. Acta Numerica, 22(April 2013);133-288, 2013.
Javier González, Michael Osbore,, and Neil Lawrence. Glasses: Relieving the myopia of Bayesian optimisation. In Artificial Intelligence and Statistics, pages 790-799.
PMLR, 2016.
Mona Eberts and Ingo Steinwart. Optimal regression rates for SVMs using Gaussian kernels. 2013.
Michael B. Gordy and Sandeep Juneja. Nested simulation in portfolio risk measurement. Management Science, 56 (182: 2010.
David Eric Edmunds and Hans Triebel. Function spaces, entropy numbers, differential operators. (No Title), 1996.
Arthur Gretton. Introduction to rkhs, and some simple ker- nel algorithmss. Advanced Topics Machine Learning. Lecture Conducted from University College London, 16 (5-3):2, 2013.
Ethan N. Epperly and Elvira Moreno. Kernel quadrature with randomly pivoted cholesky. In Advances in Neural Information Processing Systems, pages 65850-65868. 2023.
Tom Gunter, Michael A Osborne, Roman Garnett, Philipp Hennig, and Stephen J Roberts. Sampling for inference in probabilistic models with fast Bayesian quadrature. Advances in Neural Information Processing Systems, 27, 2014.
Lawrence C Evans. Partial differential equations, vol- ume 19, American Mathematical Society, 2022.
Simon Fischer and Ingo Steinwart. Sobolev norm learning rates for regularized least-squares algorithms. Journal of Machine Learning Research, 21(205): 2020.
Hanyuan Hang and Ingo Steinwart. Optimal learning with anisotropic Gaussian SVMs. Applied and Computa- tional Harmonic Analysis, 55:337-367, 2021.
Alexandra Gessner, Javier Gonzalez, and Maren Mah- sereci. Active multi-information source Bayesian quadrature. In Uncertainty in Artificial Intelligence, pages 712-721. PMLR, 2020.
Satoshi Hayakawa, Harald Oberhauser, and Terry Lyons Positively weighted kernel quadrature via subsampling. Advances in Neural Information Processing Systems, pages 6886 6900, 2022.
Michael B. Giles. Multilevel Monte Carlo methods. Acta Numerica, 24:259-324 2015.

---

Nested Expectations with Kernel Quadrature
Satoshi Hayakawa, Harald Oberhauser, and Terry Lyons. Sampling-based Nyström approximation and kernel quadrature. In International Conference Machine Learning. volume 202, pages 17778-12699, 2023.
Toni Karvonen, Chris Oates, and Mark Girolami. Integra- tion in reproducing kernel Hilbert spaces of Gaussian kernels. Mathematics of Computation, 90(331):2009- 2233, 2020.
Anna Heath, Ioanna Manolopoulou, and Gianluca Baio. A Review of Methods for Analysis of the Expected Value of Information. Medical Decision Making, 37(7):747- 758, 2017.
Frances Y Kuo, Weiwen Mo, and Dirk Nuyens. Construct- ing Embedded Lattice-Based Algorithms for Multivari- ate Function Approximation with a Composite Number Points. Constructive Approximation, 2024.
Shing Hoi Lee and Peter W. Glynn. Computing the distri- bution function of a conditional expectation via Monte Carlo: Discrete conditioning spaces. ACM Transactions on Modeling and Computer Simulation, 13(3):238-255 2003.
Philipp Hennig, Michael A. Osborne, and Hans Kerst- ing. Probabilistic Numerics: Computation as Machine Learning. Cambridge University Press, 2022.
Jeff L.. Hong and Sandeep Juneja. Estimating the mean of a non-linear function conditional expectation. In Proceedings of the 2009 Winter Simulation Conference, pages 1223-1236, 2009.
Kaiyu Li, Daniel Giles, Toni Karvonen, Serge Guillas, and François- Xavier Briol. Multilevel Bayesian quadra- ture. In International Conference on Artificial Intelli- gence and Statistics, pages 1845-1868. PMLR, 2023.
Tuomas Hytonen, Jan Van Neerven, Mark Veraar, and Lutz Weis. Analysis in Banach spaces, volume 12. Springer, 2016.
Jihao Long, Xiaojun Peng, and Lei Wu. duality analysis of kernel ridge regression in the noiscless regime. arXiv preprint arXiv:2402.157/1. 2024.
Vesa Kaarnioja, Пia Klebanov, Claudia Schillings, and Yuya Suzuki. Lattice rules meet kernel cubature. arXiv:2501.09500 2025.
Francis A Longstaff and Eduardo Schwartz. Valuing american options by simulation: simple least-squares approach. The review of financial studies, 14(1):113- 147, 2001.
Heishiro Kanagawa, Wittawat Jitkrittum, Lester Mackey, Kenji Fukumizu, and Arthur Gretton. A kernel stein test for comparing latent variable models. Journal of the Royal Statistical Society Series B: Statistical Methodol- ogy, 85:986-1011, 2023.
Firdous Ahmad Mala. Value of information for healthcare decision- making. Technometrics, 66(4):677-678, 2024. doi: 10. 00080/0040170222222220022
Dimitri Meunier, Zhu Li, Tim Christensen, and Arthur Gretton. Nonparametric instrumental regression via ker- nel methods is minimax optimal, 2024. URL https:
Monotobu Kanagawa and Philipp Hennig. Convergence guarantees for adaptive Bayesian quadrature methods. In Advances Neural Information Processing Systems, pages 6237-6248, 2019.
/ / //arxiv.org////// h2tpr/////////1/11a1/
Erich Novak. Deterministic and stochastic error bounds in numerical analysis, volume 1349. Springer, 2006.
Monotobu Kanagawa, Bharath K. Sriperumbudur, and Kenji Fukumizu. Convergence analysis of determin- istic kernel-based quadrature rules in misspecified set- tings. Foundations of Computational Mathematics, 20: 155-194, 2020.
Erich Novak. Some results on the complexity of numer- ical integration. Monte Carlo and Quasi-Monte Carlo Methods: MCQMC, Leuven, Belgium, April 2014, pages 161-183, 2016
Motonobu Kanagawa, Philipp Hennig, Dino Sejdinovic, and Bharath K Sriperumbudur. Gaussian processes and kereel methods: A review on connections and equiva- lences. arXiv preprint 187.00288, 2018.
Chris Oates, Mark Girolami, and Nicolas Chopin. Control functionals for Monte Carlo integration. Journal of the Royal Statistical Society B: Statistical Methodology. 79 (3):695-718, 2017.
Toni Karvonen and Simo Särkkä. Fully symmetric kernel quadrature. SIAM Journal on Scientific Computing, 40 (2):697-720, 2018
Chris J. Oates, Jon Cockayne, François-Xavits Briol, and Mark Girolami. Convergence rates for a class esti- mators based on Stein's identity. Bernoulli, 25(2):1141- 1159, 2019.
Toni Karvonen, Simo Särkkä, and Chris Oates. Symmetry exploits for Bayesian cubature methods. Statistics and Computing, 29:1231-122 2019.
Anthony O O'Hagan. Moyes-Hermite quadrature. Journal of Statistical Planning and Inference, 29:245-2600 1991.
12

---

Nested Expectations with Kernel Quadrature
Michael Osborne, Roman Garnett, Zoubin Ghahramani, David K Duvenaud, Stephen Roberts, and Carl Ras- mussen. Active learning of model evidence using Bayesian quadrature. Advances in Neural Information Processing Systems, 25, 2012.
Ingo Steinwart. Support Vector Machines. Springer, 2008.
Ingo Steinwart, Don Hush, Clint Scovel, et al. Optimal rates for regularized least squares regression. In COLT, pages 79-93, 2009.
Art B Owen. Quasi-Monte carlo sampling. Monte Carlo Ray Tracing: Siggraph, 1:69-88, 2003.
David Stirzaker. Elementary Probability. Cambridge Uni- versity Press, 2003.
Fred J. Hickernell R. Jagadeeswaran. Fast automatic
Zhuo Sun, Alessandro Barp, and François-Xavee Xavier Briol. Vector-valued control variates. In International Confer- ence Machine Learning, pages 32819-32846, 2023.
Bayesian cubature using lattice sampling. Statistics and Computing, 29(6):1255-122) 2019.
Tom Rainforth, Rob Cornish, Hongseok Yang, Andrew Warrington, and Frank Wood. On nesting Monte Carlo estimators. In International Conference Machine Learning. pages 4267-4276. PMLR, 2018.
Taiji Suzuki and Atsushi Nitanda. Deep learning is adap- tive to intrinsic dimensionality of model smoothness in anisotropic Besov space. Advances in Neural Informa- tion Processing Systems, 34:3609-3621, 2021.
Tom Rainforth, Adam Foster, Desi R. Ivanova, and Fred- die B. Smith. Modern Bayesian experimental design. Statistical Science, 19(1):100-114 2024.
Roman Vershynin. High-dimensional probability: An in- troduction with applications in data science, volume 47. Cambridge university press, 2018.
Carl Rasmussen and Zoubin Ghahramani. Bayesian Monte Carlo. In Advances in Neural Information Processing Systems, pages 489-496, 2002.
Jialei Wang. Scott C Clark, Eric Liu, and Peter I Frazier. Parallel Bayesian global optimization of expensive func- tions. Operations Research, 68(6): 850-1865, 2020.
Klaus Ritter. Average-case analysis of numerical prob- lems. Number 1733. Springer Science & Business Me- dia, 2000.
Holger Wendland. Scattered data approximation, vol- ume 17. Cambridge university press, 2004.
Christian Robert and George Casella, Monte Carlo Sta- tistical Methods. Springer, 2000.
Jonathan Wenger, Nicholas Krämer, Marvin Pförtner, Jonathan Schmidt, Nathanael Bosch, Nina Effenberger, Johannes Zenn, Alexandra Gessner, Toni Karvonen, François-Xavier Briol, et al. ProbNum: Probabilistic Numerics in Python. arXiv preprint 12.00021 2021.
Walter Rudin. Principles of mathematical analysis, vol- ume 3. McGraw-hill New York, 1964.
Alexander Shapiro, Enlu Zhou, and Yifan Lin. Bayesian distributionally robust optimization. SIAM Journal Optimization, 3379):127- 2023.
Jian Wu and Peter Frazier. Practical two-step lookahead Bayesian optimization. Advances in Neural Information Processing Systems, 32, 2019.
Shijing Si, Chris J. Oates, Andrew B. Duncan, Lawrence Carin, and François-Xavier Briol. Scalable con- trol variates for Monte Carlo methods via stochas- tic optimization. Proceedings of the 14th Confer- ence on Monte Carlo and Quasi-Monte Carlo Methods. arXiv:2006.074877 2021.
George Wynne, François-Xarin Xavier Briol, and Mark Girolami. Convergence guarantees for Gaussian process means with misspecified likelihoods and smoothness. The Jour- nal of Machine Learning Research, 22(1):5468-555 55077, 2021.
Alex Smola, Arthur Gretton, Le Song, and Bernhard Schölkopf. A hilbert space embedding for distributions. In International conference algorithmic learning the- ory, pages 13-31. Springer, 2007.
Shangda Yang, Vitaly Zankin, Maximilian Balandat, Ste- fan Scherer, Kevin Thomas Carlberg, Neil Walton, and Kody JH Law. Accelerating look-ahead in Bayesian op- timization: Multilevel Monte Carlo is all you need. In Forty-first International Conference on Machine Learn- ing. 2024.
Alvise Sommariva and Marco Vianello. Numerical cuba- ture on scattered data by radial basis functions. Comput- ing. 76:295-310, 2006.
Yang Song, Sahaj Garg, Jiaxin Shi, and Stefano Ermon. Sliced Sscore matching: A scalable approach to density and score estimation. Uncertainty in Artificial Intelli- gence Conference, pages 574-584, 2020.
Xiaoyan Zeng, Peter Kritzer, and Fred J. Hickernell. Spline methods using integration lattices and digital nets. pages 529-555, 2009.
13

---

Nested Expectations with Kernel Quadrature
Haobo Zhang, Yicheng Li, Weihao Lu, and Qian Lin. On the optimality of misspecified kernel ridge regression. In International Conference Machine Learning, pages 11331-41353, PMLR, 2023.
Zhen Zhang, Xin Liu, Shaoli Wang, and Jiaye Teng. Minimax optimal two-stage algorithm for moment es- timation under covariate shift. In The Thirteenth International Conference on Learning Representa- lions, 2025. URL httpss///ppereeee...m..
forum?id=ock Eorum?2droc4y772X31
14

---

Nested Expectations with Kernel Quadrature
Supplementary Material
Table of Contents
15

---

Nested Expectations with Kernel Quadrature
Additional notations For two normed vector spaces A, B, A = B means that A and B are norm equivalent, i.e. their sets coincide and the corresponding norms are equivalent. In other words, there are constants c, C2 such that ei ||Å||A A S | I c2]lA||a holds for all h E A, written as || ||A = 1-Ila. For A C B, A is said to be continuously embedded in B if the inclusion map between them is continuous, written as A + B. ||T)| denotes the norm of an operator T A + B. For function f X' C IRd - R and a E N4, we use 0g f to denote the standard derivative Oxs ... -r f and D2 f to denote the weak derivative. For f E Wz(X), we use ||fll,,2 = (l to denote its Sobolev space norm. X means < up to some positive multiplicative constants.
A. Existing Results on Kernel Ridge Regression
In this section, we present Proposition to 3 which are adaptation of theorems from Fischer and Steinwart (2020) applied to Sobolev spaces. These propositions are foundations of the proof of Theorem I in Appendix C.
In the standard regression setting, we are given N observations ( 3 which are i.d sampled from an unknown joint distribution on X' R. Here, X' C Rdd is compact domain. The marginal distribution of P on X is m, and the conditional distribution P( x) satisfies the Bernstein moment condition (Fischer and Steinwart, 2020). In other words,
there exists constants o, L > 0 independent of such that
[ ly h-
x) S mm m- 2
1
(A.1)
is satisfied for r-almost all x X' and all m 2. For example, (A.1) is satisfied with L oo when P(- a)) is Gaussian distribution with bounded variance do- Additionally, (A.1) is also satisfied when there is noise in the observation s = L- = 0, which will be discussed in Appendix B.
In a regression problem, the target of interest is the Bayes predictor h" : X' + R. - E[y x = z]. One way of estimating h is through kernel ridge regression (Fischer and Steinwart, 2020): given a reproducing kernel : x' + R, the kernel ridge regression estimator hy + I is defined as the solution to the following optimization problem (A > 0):
N
1
hu = arg min
+ N
[v (y--((a)) h(:)))
(A.2)
H. is the reproducing kernel Hilbert space (RKHS) associated with kernel k. Fortunately, it has the following closed-form expression (Gretton, 2013, Section 7)
'y = k(.,I1:N) (k(ZIN, #1:N) + NAIN) y1:N-
We also introduce an auxiliary function Ix : X - 3' which is the solution to another optimization problem:
AMIFFIRA + ST- (y (w- ((=)(()) dy)
hy = arg min JEHL
(A.3)
In regression setting, it is of interest to study the generalization error between the estimator hy and the Bayes optimal and particularly its asymptotic rate of convergence towards 0 as the number of samples N tend to infinity. The generalization error can be decomposed into two terms, through triangular inequality,
predictor h", I|A h- [l)))()) -h"]ll
ll - h²" (z)(()) I !l - hxllnggm) + Il - h" -h lI))))))
(A.4)
where the first term I! - (x) is known as the estimation error and the second term // - 11 (1''l22(r) is known as the approximation error. Next, we are going to present propositions that study these two terms separately under the following list of conditions.
(S1) is a Sobolev reproducing kernel of smoothness s > 2-
(S2) is probability measure X' with density p: X R. There exist positive constants Go, G1 such that Go <
p(x)= < G1 for any X.
(S3) The Bayes predictor h* E W (2').
(S4) There exists universal constants a, L > 0 such that (A.1) holds.
16

---

Nested Expectations with Kernel Quadrature
Proposition 1 (Approximation error). Under Assumptions (S1)-(S4),
112 - 1 (=1)) S 114" |l, Ab.
Proof. This is direct application of Lemma 14 of (Fischer and Steinwart, 2020) with B = and 7 = 0. Proposition 2 (Estimation error). Suppose Assumptions (S1)-(54) hold. Let (A)) be the effective dimension defined in Lemma I, and ka be defined in Lemma 2. N > Ax, then with probability at least 4e"
5767² N (12DA- de + -21- 1 + M
L
l1k2 - )xll22(m) S
A 1
(A.5)
N
where and M are constants independent of N, and gx, Axto Lox are defined as follows
9A = log (2)
Ax,- Sk TgA # La = max { (1, th (11 + (al1/"/l..2))
Proof. This proposition is special case of Theorem 16 in Fischer and Steinwart (2020) under the following adaptations towards our Sobolev space setting: 1) Lemma 1 proves that (V(A) S DA & and Lemma 2 proves that kis S M for 2) ||1 h ll (T) upper bounded by A 4x (11h"l (r) + ka [||h" proved in Corollary 15 of Fischer and Steinwart (2020). ||,,] is the norm of the covariance operator defined in (E.40).
Proposition 3. Suppose Assumptions (SI)-(S4) hold. For Ax, and L, defined above in Proposition 2, if N A then with probability at least 1
57672 T.2DA 24 + 12 24 ||'11 + 2M
L3 N
h1 h- La(ar
A
26
+ 111 1111~12222 A.
N
Proof. By the triangle incquality in (A.4), combining Proposition and Proposition 2 finishes the proof.
B. Noiseless Kernel Ridge Regression (Kernel Quadrature)
In this section, we present the upper bound on the generalization error ||h hx (x)) in Proposition 3 adapted to the noiseless regression setting, which will be employed the proof of Theorem in the next section. Our proof follows the outline of the proof for Theorem 1 in (Fischer and Steinwart, 2020), modified for our choice of regularization parameter X. Note that this section is of independent interest to some readers as it presents the first standalone proof on the convergence rate of kernel quadrature that I): it allows positive regularization parameter À >> 0 and 2): it provides convergence in high probability rather than in expectation. The closely-related work is Bach (2017) which requires i.i.d samples from an intractable distribution; and Long et al. (2024) which provides a more general analysis on noiseless kernel ridge regression in both seel-spee and mis-specified setting.
Suppose we have N observations FINN which are i.i.d sampled from an unknown distribution π on X' along with N noiseless function evaluations h" (I1:N) where h² : X C Rd >> R. The setting appears for instance when the measurement of the output values is very accurate, or when the output values are obtained as result of computer experiments.
Proposition 4. Let X C Idd be compact, and IaNN be N i..d. samples from ा. Define hxx()) (-) = (-,, #1:N) (M(IINN FIIN) + NN) -M (x and suppose conditions (S1)-(S4) are satisfied. Then, if AN X
N 2 (log N) N) 2422 there exists an No such that for all N > No.
11" hix (la(()) S C (log N)
11.2
(B.6)
holds with probability at least 4e for constant C(X, Go, G1) that only depends on X, Gu. G1.
Proof. Notice that hAv is precisely the solution to the optimization problem defined in (A.2) only with yi replaced by h" (x.). Similarly, we define hxy as the solution to the optimization problem defined in (A.3) only with y replaced by
h" (x). Note that Assumption (S4) instantly satisfied with T. 0.
Similar the proof of Proposition 3, we decompose the generalization error into an estimation error term || h ll L2(r)) and an approximation error term | h² [2(())
17

---

Nested Expectations with Kernel Quadrature
Approximation error
Take AN = N 2 (log N) 2n12 tten from Proposition 1, we have
| h (z) S A = N (log N) a | 1l,,2-
Estimation error
Recall all the constants 9Ax= AxN. and LAN defined Proposition 2. Since I. 0, we know the + (ell/"/ls.)) In order to apply Proposition 2, we need to check there indeed exists
constant LAx = N No such that 2 AxN,r is satisfied for all N 2 No. To this end, we are going to verify that 1im AxN.T/A /N O.
Notice that
AANIT N
8h T9AN AN 8(log ~ kār log (2e
| AN
lim
N
N
|>
where (AN) and ka are defined in Lemma 1 and Lemma 2. Since lima-+so AN = lim N 3 (log N) 2 0, there exists N such that Au & |n for all N 2 N'. Therefore, as N tends to infinity,
AANTT
lim
I lim 8(log - + k log (4cN(Aw))
N N-000
N
N-000
s lim 8(log N) + k log 4eDA
N
= lim 8(log N) - - k log (4eD) +
8(log N) kar log (N( N
lim
(B.7)
N
< lim 16(1og N) - k22 log (N)
N
where M and D are constants defined in Lemma and Lemma 2. So there exists N such that N 2 Ax for all
N 2 Nº. Taking No max(N', ]",, then we have N 2 Axyi for all N 2 No. From Proposition 2, we know that with probability at least 4e-
576- N
Ihay - -ux (2
20 |
M (114" + M11h"l11.. (lx.2) N 1
57672
(M2 (log N) =11 2nr || |1:011 + 2² (11h" (em)))) + Mllh" 11,,2) N 3 (log N) N) $ 2-4
(glog
< 5 N 2 (log N) 2 MM² 11 11 + 2 (11""|(=~(>]) + M| M11h""aa22)
So have, with probability least 1 TT
|hxN NAw (1 S 24TN * (logN) N) 2ي = (IM + M²) ||," + MI ll
Combine approximation and estimation error Combining the above two inequalities on approximation error || hn (l)())))))))))) and estimation error |lhx, hxw (l we have with probability least 1 4T
| - huwlle S 24 (log ~ (a1+ M M²) |1 + M (11h111l (12).(())
Finally, following the arguments of Lemma 2 that the operator norm of W (2') + (o) is bounded, we have ||* ()) I R I1A""l,, where R is constant that depends on X, Go, G,. With probability at least 4e
- 1 (IL S 2 á (log N) * (1+ (1) + M²) ||* = CTN (log (logN)= + |1.
for C = 24(1 + (1 + R)M+ + M2), which concludes the proof.
Corollary 2. Let X be a compact domain in Bd and Z1.N are N i.id samples from n. Ex~rlk(X, ((,NI (K(1IN, NN - h² (œ1:) ) is the KQ estimator defined in (5). Suppose conditions (Al)- (A3) are satisfied. Take AN x N (log N) 2يي2 then there exists No >> O such that for N No.
Ixo
ixe / h (c))dT(a) S C (log
(B.8)
holds with probability at least 4e Here C is constant that is independent of N.
18

---

Nested Expectations with Kernel Quadrature
The proof of Corollary is direct application of Proposition after observing the following.
Ixo - h- f
- h" h=(x) dx(x) = 11h- (hh"-hwwliaaee) hhw hAw (=(x))d(()) = $ (x)
-
Remark B.1. We prove in Proposition that the generalization error of hax in noiseless regression setting is Ó( *). which is faster than the minimax optimal rate O(N wia) in standard regression setting. The fast rate expected because we are the noiseless regime so "overfittin"" is not a problem hence our choice of regularization parameter AN x N- (1og N) - decays to O at a faster rate than As x N- in standard kernel ridge regression (Fischer and Steinwart, 2020, Corollary 5). The Ó(N-2) rate also optimal (up to logarithm terms) and cannot be further improved because matches the lower bound of interpolation (Sections 1.3.11 and 1.3.1 of Novak (2006), Section 1.2, Chapter of Ritter (2000)).
Remark B.2 (Comparison existing upper bound of kernel (Bayesian) quadrature). The upper bound in Corollary 2 matches existing analysis based on scattered data approximation in the literature of both kernel quadrature and Bayesian quadrature (Sommariva and Vianell., 2006; Briol et al.. 2019; Wynne et al., 2021) and is known to be minimax opti- mal (Novak, 2016; 2006). Existing analysis takes X = 0 and requires the Gram matrix K( T1:N) to be invertible, in contrast, our result allows a positive regularization parameter AN z N 3 (log N) 2=2 which improves numerical stability of matrix inversion in practice. One closely-related work is Bach (2017), but it requires i.dd samples from an intractable distribution.
C. Proof of Theorem 1
Remark C.I. In this section, we use p(3;; 0) to denote the density po so that we can use p(a:;)) -) to denote the mapping 9 + (z).). Although introduce shorthand notation of kernel mean embedding in the main text, Pt, Ex [xx,rI(XX, .).. in this section we going to write it out with its explicit formulation.
For any 0 E O, Fio θ + IR and JkQ O + R are two functions that generalize the definition of Fk)(0,) and (k)
in (7) to all 0 E @. To be more specific, for any E O, given samples x : [a] .....IN 200 samples from Pa.
consisting N i.d.
(0; (0)) ( = ( a (kx(= (0) NAxIN)
g(
(C.9)
Fko (0)) 5(Jko(0; -1))))
(C.10)
where we explicitly specify the dependence of samples Vi:N on 9 in the above two equations. Next, we define
N p(:;: 0)
[k)(0) = E 10 ( = 1 jko(0:
II
(C.11)
N II 1=1
(k) = E Lig NB [Foo0:2200) (000:11000 / (xo(0:IINN)
p(z;; O)
(0) which marginalize out the dependence on samples TINN- We can see that JKQ E L2( since g(x,-) E W (0) (2) z L₂(Q) from Assumption (1) and (3); and p(za-) Lz(Q). Also FRQ E Lz(() because f is Lipschitz continuous from Assumption (4). Therefore, the absolute error [ - isKol can decomposed as follows:
INKO|
F(O))(0)d0 ( ko((0, (
(KT()
F(0))(0yd0 1 (ko(0))(0)(0))
+ 1/ Fko(O)g(O)d0) (1 (a(O, (1)))(()) (ko 01r)) + TheIr) Fko(O1:T)
19

---

Nested Expectations with Kernel Quadrature
S Bena
- - (0]
+ Fko(-) - (ka 01r) + (helT) Fko F
(C.12)
12(0)
Stage
Stage IT erroe
The last inequality holds because | (Q) I | (lla(Q)) Next, we analyze Stage error and Stage error separately,
Stage Error From Assumption (4), Lipschitz continuous and the Lipschitz constant bounded by S4.
[kk(() F(0)
E #I.N~~ ~o
FRQ (0: (0::)=1)) F(0)
S (o) LNN ~Pa
Fko (0: x ((xx) F(0)
- S4 1. ~P
| (0; = J(O}
(C.13)
where the first inequality holds by Jensen inequality and the last inequality holds by Lipschitz continuity of f. Define
g(x, 0: 2: kx(2, ((xx (0)
(N NAx (
(C.14)
Here 0(:,0==1:N) (Po) because the Sobolev reproducing kernel kx is bounded and measurable; and g(: E (2(P))
by Assumption (3). Thus,
g(x, 0))p(r; O)dr s [0, - s(.,0) Lg(P))) [kQ(0; - - 0:
(C.15)
Based on Assumption (2), g(, 0) E W2* for any 0 E O. Therefore, based on Proposition 4, if one takes Ax,N x
N - (logN) 2xx then there exists No such that for N No.
100.0:2:: - g(, L2(Ba)
I ETN
(log N) "x (((-0))lexx..
(C.16)
holds with probability at least -T- The probability is taken over the distribution of 2 ie Po. Here e is constant independent of N. Hence, with Lemma 4, have
E (o) Po 0(1, 0, g( 0) (2(F0)
S SeN 3x (log N) 2x4
11 01lxx.2
(C.17)
By plugging the above inequality back into (C.15), obtain
E
LP [ko((f: x (M) J(0)
S seN=
2x *x (log "x
2x (9((,6)]le
Therefore, the Stage I error be upper bounded by
E F(I) (x) S4 (0)
[ko(0:= (0)
J(0) (1o(;;0) 0))sx.2
~Po
SS,iee 38
(log N)
N) 2x42
C (log N)
5x
ext1
(C.18)
where C3 == S4S,C E is constant independent of N.
Stage II Error The upper bound on the stage II error is done in five steps. In step one, we prove that 3ko(
(0) (0) given fixed samples #1.N- In step two, show that J E W (O). In step three, we upper bound [kko(; x (. )l through the triangular inequality that [| x (0) Il & || + || - - llse. In [ (Jxo(<s=10)) through marginalizing out the samples x In the last step, we use kernel ridge regression bound proved in Proposition to upper bound the stage error.
step four, we upper bound (O) E Egg
20

---

Nested Expectations with Kernel Quadrature
Step One. In this step, we are going to show that Jkq lies in the Sobolev space W2- given fixed samples (0)
T--N- Notice
that the dependence of JkQ(0) on is through two mappings: 0 + kx(xx "i )p(;; 0)dr and + g( 0). We (0)
(0)
going to show that 0 + Sx (x, (0) )p(a; d) lies in the Sobolev space n (0) for any i E (l...., N). To
this end, we are going demonstrate it possesses weak derivatives up to and including order se that lie in (+). Take
p ④ + R to be any infinitely differentiable function with compact support in θ (commonly denoted 4 E C (+)), with its standard, non-weak derivative of order 3 denoted by &hep. Since 8 + p(x;0) 0) E "2" (Ø), for any |3| S se it has weak derivative 0 + p(x; 0) E (0). Then,
p(p(2;))) (2; 0)dx @ J 1x(e, 200) / ((0) ((; Is ((0) / kx 2 (0) ( (it) (-1)181 / kx((, (00) / 0d (00p(x; Oyd0de (101 - (-1))88
(C.19)
/ 0 03 / La ! >(<; Oydzd0.
In the above chain of derivations, we are allowed to swap the integration order in (i) by the Fubini theorem (Rudin, 1964) because kx is bounded and the fact that 0 + ((0) p(x;;)) E (i))) since Dzp((;;)) E L2(@) (Assumption (3)) and &p E (+): (ii) holds by definition of weak derivatives for D6 (z; 0); and (iii) holds again by the Fubini theorem. By definition of weak derivatives, (C.19) shows that fx (2, x; (0) (x(:; f) has a weak derivative of order 3 of the form
Dy Uxx kx(x, (1z;0)dx = =/s kx(x, x Di 0)dr
Also, since xx bounded and - p)p(2;0) 0) E (O), the weak derivative above is in (O). Consequently, we have
Σ
D f kx )p(;; 0)dz 2 d0
Σ
1 1/
( de kix x (1) (o)p(o;))de( do
e
181500
Sala
Vol(X)
Σ
1815
(())
1///2/ D8p(1;; ? drdd
Vol(X)
Σ x²
B|5
(136)
Vol(A*)
/ ||p(1; dr.
x
In the above chain of derivations, (i) holds because 1xx (xx)de S Vol(2*) (zx) Fdr for compact X, (ii) holds because kx is upper bounded by x and Se p(g 0) < oo from Assumption (3), (iii) holds because p(z; W (0) for any x E X based on Assumption (3). Also, one can interchange the order of integration in (iii by the Fubini's is theorem (Rudin, 1964).
( As a result, for any i,j E (1,..., ..). we have fi, 00 + Sx 1x (z. ((;;; 0)dr E W5 and f2g : 0 + g(x)) 0) W (9) from Assumption (3). Therefore, we know from Lemma 3 that their product fú,; fai E W (0) hence kkq as linear combination of fi. f2. is in (0).
Step Two. In this step, we are going to show that 0 fx g(x, (),(2; 0)de is also in the Sobolev space W (0). Since both .)x, E (7)(() and p)=; E (O). we know from Lemma 3 that A + g(,, 0) p(x, 0) E W (0). By following the same steps as in (C.19), we obtain that for any [8] < se.
D& / g(2, 0)p(; 0)dx= / Dg D2( (((x0))((;;0))) dr.
(C.20)
We are now ready to study the Sobolev norm of J,
2 / Dg / p(x; p(::0)9(o, ()g(() 0)dr
|
Σ 18
do
21

---

Nested Expectations with Kernel Quadrature
Σ
/
1
D& (p(o; dr do
2
18/5.
(18)
Vol(22)
/
Σ
//108 (p(x; D0se0, d0da
101500
(w)) vol(xt) / lIp(2; )9(( >)1 dr
x vo1(2))2253
(C.21)
(d Here, (i) holds by (C.20), (ii) holds since (x))dx(2= < Vol(x) fx (f(x) drr for compact X, (iii) follows from the
definition of Sobolev norm, and (iv) holds by Lemma and Assumption (3) that 119(z, oll.s. S S2, [pp(;; el I Sg.
Step Three In this step, study the Sobolev norm of SkQ for some fixed TINN' by upper bounding it with
|J Jko(; 11:))l sls.. + || [sl.. Since gl .-) E W (+) by Assumption (3), it holds that g(x, (0) *x(2,=1NN (kx (kx "IN.N. )= + N
g(x,;=1[))
) is in W (0) for any fixed (10) Therefore,
|p(e; (9(2,)) - 0(0,, = (0)
11
< |,([;; (((x;-)llan 1l.a.2 (o(c,)) - ((t,,-=1111 (0)
20.2
16.2
S S3 o(z,)) ((2,;;=11.) (0)
(C.22)
where the first inequality holds by Lemma 3 and the second inequality holds by Assumption (3) that [|p(2;; (1(xx)))loos >1a0, S Sg. Now, we consider the Sobolev norm of JxQ J,
|| | (2) - 1 Σ
S Do / p(z:0) (ola, (o(,,0))((x,))) ((2, 0:2(1h)) dr
2
d0
19.2
1815
1
Dg (plo;)) (g(x, 0) g(2, 0: 0:2(.2)) )dz d0
1815
(1))
Vol(2))
Σ
0(x, 0; x )²² 1 d0dx 1 D8 18 (p(x: (o(z, -(x,)))90000000000))
x
1,8
vol(x) /, |ple:) (gla,)) g(z, xx(()))))))) (0)
1 dr
s0,2
so / - (0) 2
Hxx, ;;
dæ,
(C.23)
where the above chain of derivations (i) (iv) follow the exact same reasoning as (C.19) and (C.21), Next, notice that
1 [llee - 002, - 1 (0) dr =
1 Σ
/
D (gla, 0) - g(x, 0; (; I' d0dr
s0.2
13 See
Σ
12 0) DE 0: (1 Lx(x)
d0.
(C.24)
By Assumption (2), De 0) E (xx(x) for any |3| < 8e. Therefore, by applying Proposition 4 with h*(+)) = + NAxIN) - -D2 B; ( we get that
Dgg(. 6), and hxx-)) = 0/00(,0;=1 kx(. (0)
(kx a (9)
(0)
0; 11xx))))
S CTN
2 (log N) 2x |pas(:0)1 1xx2
(C.25)
(Po)
holds with probability at least 1 -T, for a that only depends on X. Go,x, GI, From Assumption (1), we know that Ta (Po) 2 (xx) (they are norm equivalent) and [[] I Vol(At) - (llx(())) for any f (z()))). Therefore,
22

---

Nested Expectations with Kernel Quadrature
for any 0 E ④ and any |3| S se, with probability at least 4e
| 0) - (000-106600 0; riN
l2((2)) - e-Vol(2)) G xx (log N) 03)
xx
<E>Vol(2))
#x (log N)
2 St.
(C.26)
By plugging (C.26) into (C.24), and then plugging the result into (C.23), we get that with probability at least
4e
|ko((=2)))) (0) J
S
Σ
GT
C
V # (log N) 3x
se.2
|
-(re
$80 + de - 1)
2012 (C.27)
de -1 GTG Vol(2)) " (log N)
"x
By combining this result with the bound
| Vol(X))22g proven in (C.21), we get that with probability least - 4e and any N No it holds that
|| |kolss0
S |ko((s - 3 -J + +|JJ1l
10.2
18.2
880 + de -1
1x11 (log N) xx Vol( #
er Cr S1
de 1
$ 2Vol(2)
(C.28)
where Ng is defined the smallest integer for which the first term is subsumed by the second term.
Step Four: In this step, are going upper bound the Sobolev norm of FRo- From Chapter 5, Exercise 16 of (Evans, 2022), we have Fko = f o jkQ is in w (9) because f has bounded derivatives up to including se + 1 and [: 11th ( lls.2 I Vol(X))253 with probability least 4e- proved in (C.28). Hence, || )lea. S Co
(()
holds with probability at least 4e- Next, recall the definition of (k) in (C.11),
(0) = 100 ~Pa ~Ba [Fko (0: (0:=12) (0) )
... Jx.-.J/ (xo(0:=11))(aa plx (0) ; 0)p(122)))1))) "000 p(x ; dz dx dx
(0)
we know that lip(at (0); ( )lou, S S3 from Assumption (3) and | (0) llea,2 S Co proved above. Therefore, from Lemma 3 we have [pp(z)") ;)) = )lee) is bounded, so x!" + p(x!") :-) Fio(- Bochner integrable with respect the Lebesgue measure Lx. From Lemma 5, we have,
For any i = I,...,N, N.
is
N FFOllle $ Sx... Jx. Prool: ( II p(=10)))) a=l
dx ( day ...dx
dx
$0.2
1 | "x [ [
II (0)
: )1so.
dx{ dx! (0)
10.2
s CoSA Vol(x))
s Co.
(C.29)
The last inequality holds by Sy < from Assumption (3) and X 10, 1]xx so Vol(2') = 1.
Step Five. We are now ready to upper bound the stage II error, which was defined as
Stage II error = Fio(.) ()) (-,61.) (ke (0.r,:1:)) TheIr)
Fro (fE:)
The idea is to treat the stage II error the generalization enror of kernel ridge whhr can be bounded via Proposition 3. Given i..dd observations (01, (k, (00)) (..)), (fr, (ko (0r, (0T) )), the target of interest in the context
23

---

Nested Expectations with Kernel Quadrature
of regression is the conditional mean, which in our case is precisely (k) E =./NFF PkQ 1.2 defined in (C.11).
Alternatively, (fo =I.N) can be treated as noisy observation of the target function Fkc where the observation noise
is defined as r: Q + R with r(0: ((0:=1i. = Fk)(0:=):N) (0) - (O). So we automatically have E. uig [r(0]] = 0. For
any positive integer m 2 2,
E _0000
Po Pa [[][] = E 10 ~o
Fro(0; 10) =ko(0)
(8)
- 2 1
E
FROl0: , F(0) m
+2m-1 |Ko(() F(0)
= 2m 1
E
(6)
~Nw
FRol0: (1)) x - F(0)
+ 2
E B_
FiN
~F, F. F - F(0)
< 2
E
(6)
Fxo(0:=100) (0) - F(0)
+ 2
E (0) ~po
Pko x0 (0)
- F(0)
FRo x111 F(0) m
2 E (0) ~
Jko x0 ((0 m
< 2 ST E (0)
IN
~P'
(t) x 2 mSST eT N ax (log N)
22
m 2x11
xx
(C.30)
In the above chain of derivations, (i) holds because (a + b)md S 22--1(a" + 6Tt). (ii) holds because we know from (C.15) and (C.16) that [ =ffN)) - ((0) I CTN के (log N) xx S₁ holds with probability at least - -e- and so E B,00r |kkq (0; (0) I( can be bounded via Lemma 4. Therefore, by comparing (C.30) with (A.1), we can see that the observation noise r indeed satisfy the Bernstein noise moment condition with
o o= = 2S; 5% (log.) +x+1 = Cy.N"
5x (log (log.N) N) 'x+1
for C7 := 2 constant independent of N. Before we employ Proposition 3, we need to check the Assumptions (S4). Assumption (S1) is satisfied for our choice of kernel ke. Assumption (S2) is satisfied due Assumption (1). Assumption (S3) is satisfied due to (C.29). Assumption (S4) is satisfied for the Bernstein noise moment condition verified above. Next, we compute all the constants in Proposition 3 in the current context. N/ (Aa) is the effective dimension defined 2se defined in Lemma 2 upper bounded by a constant Me.
in Lemma upper bounded by Da>a d /2s0 ka with a de
| |I is the norm of the covariance operator defined in (E.40). Hence
[Ee ll + ) [Eel
Dh. == log (Ne)
Axo.r = sk Ao
2
doa t L,
(1 + kr (kelFFell)))) }.
Lxe == max
L
Applying Proposition 3 shows that, for T Ase.7.
|| |iko- - 1s (::0r:)) (Kie (01.T, 01:) + IT) -1 Akro (0I:T)
(zz((1)
57672 T
di J
1- dus
1.2
des F de
+ 13100
|| + 2M3
A
20
+
(C.31)
T
holds with probability at least 1 - TT, We take do x 2 (log 1)
2,0
20++
di then similar to the derivations from (B.7),
AxeiT
lim T
lim 16(log T) T000
k2 r log (T) 0.
(C.32)
T
It means there exists a finite To > 0 such that T > Axour holds for any T > To. Notice that, with probability at least
1 4e-
( = =|| (0) S Roll Fkollaao S ReCs
(C.33)
24

---

Nested Expectations with Kernel Quadrature
based on (C.29) and the fact that W2 (0) Lee (0) with (0) + L.00 (+)] S RS. we have
Lae I maxxLL,T " + (log T)
1a+1 Th
20g 2ag-da
(Re + MMM C6}
+x max(C;N xx (log N)
1x+1 xx
1a+1 2ia xig
T
(log T)
(Re + Me) Cs}.
So the above (C.31) can be further upper bounded by
576r²
CGN
2 2x (log N)
2xg DeT(log 7)
dx
11
27e M
+1
"log T)
21 70
2ng de
d0
C
T
max CGN 2 (log N)
12
dx
T
2 50
+1
(log T)
22
20g
de 40
(Re + Me)² C/
5767
29
+
2M
T(log T)
-0
T
T
+ CCT
2d0 (log T"
+2 4g
-77672
CGN
2 (log N)
22 dx
Dlog -> 1g + MMT
1et1
2-a
de (log T)
1e-
200 de
48
C/
+ 57672
2M max CGN 2 (log N)
xx
xx
T 1-1
T
2.2
"
(log T)
ve 11
-0
20 de
(Re Ma)² 08 ) (log (ogg?)=o 22
"8
+ CCT
CCT
- 22 (log T)
12 "
I 5767 C - (log N)
xx xx
20 De + 5767
M
(log T)
22 dg
C²
+ +57672 2. 2M C2 - N (log N)
12 xx
+ 57672 2. T
2002
(log T)
200 12 "a
(Re + Ma)² CB + CZT (log T)
20g 2
xxx+ = (C x (log N) xx + C
24
2+a+2 40
Jo
(log T)
Cs, Cg are two constants independent of N,T. In (i), we use max(d1, (2) S a1 + 02 we also use the following
20-1 2ng -de (log
20
d0 S (logT)
"a
(log2) S 1.
Therefore, we have that,
Stage II error := Fike - k(, (1.1
(a..:r))(ee) (kk (k 01.)) + TAgIr) - -=
(2(
S CsN 2x (log N) N) xx + CoT as
(log T) 22/1
(C.34)
holds with probability at least Se
Combine stage and stage error
Combining the stage error of (C.18) and the stage II error of (C.34), we obtain
1 INKo
< Stage I error Stage II error
< CN # (log N) Tلل + (CON 2x (log N) N Jxt
+ CGT t (log T) 2لll
(log -0
S r (Css + Cs)N Xx (log N) 4x
2x+1
+ COT de (logT)
"8
ett! de
=: (G #x (log N) 2 + C2T 2a (log T) 2911
holds with probability at least 1 - ST Here C₁ C2 are two constants independent of N, T so the proof concludes here.
D. Multi-Level Nested Kernel Quadrature
In this section, we are going to introduce a novel method that combines nested kernel quadrature (NKQ) with multi-level construction as mentioned in Section 4.
25

---

Nested Expectations with Kernel Quadrature
D.I. Multi-Level Monte Carlo for Nested Expectation
First, we briefly review multi-level Monte Carlo (MLMC) applied nested expectations Ee-glf(Ex-> (9(X, 0))] introduced in Section 9 of Giles (2015) and Giles and Goda (2019). At each level E, we given T2 samples θ1:T, sampled Li.d from Q and we have Ne samples Ii:N, sampled i..d from Po. for each t = 1...., Tg. The MLMC implementation is to construct an estimator Pe at each level such that can be decomposed into the sum of Pe-
(8.)
L e PP Po:=
Is E [Pt] -Ea-e [Po] +
Σ
1 N2
Σ
.0)
(8)
9
The estimator Ye for Clie - Pt-1] is
T2
N₂
N₂ 1
N₂
Y
T:
1
N₂ n=1
Σ
9
-
N--1
1
Σ n=1
9
x 0.)
Ne-1
Σ
g
(x
n=Nz-1+1
To
No
Yo = To
Σ n=1
01)
No
Compared with (3) in the main text, notice that here the 'antithetic' approach which further improves the perfor- mance of MLMC (Giles, 2015, Section 9). The MLMC estimator for nested expectation can be written as
L.
MMMM
Σne Yr. 2
(D.35)
At each level E, the cost of Yd is x;N2xT;) and the expected squared error - - provided that f has bounded second order derivative (Giles, 2015, Section 9)² Here the expectation is taken over the randomness of samples. So the total cost and expected absolute error of MLMC for nested expectation can be written as
Cost O
Σ No x Te.
E - |MLMC| = O No x T
Σ 8=
(D.36)
Theorem 1 of (Giles, 2015) shows that, in order to reach error threshold A, one take Ne 2 and Tç x 2- A -2 Therefore, one has IE [ IMLMC| = O(A) along with Cost O(A-3)
D.2. Multi-Level Kernel Quadrature for Nested Expectation (MLKQ)
In this section, present multi-level kernel quadrature applied to nested expectation (MLKQ). Note that MIKQ is different from the multi-level Bayesian quadrature proposed in Li et al. (2023) because our MLKQ is designed specifically for nested expectations. At each level f, we have T7 samples 01.,, sampled i.i.d from Q and we have N₂ samples xi sampled i.i.d from Pe, for each & ,....,II, Te. Different from MLMC above, we define
(0.)
L Σ 1=1
N ~Ba (Jko (B:2fiiv))
I n E/ L
Ee [PNkq.o] +
Egng [PrQQ PNKQ.e-1), PsKQ,
E 10
The estimator Y for defined in (8).
- P when f 1 is the difference of two nested kernel quadrature estimator
YNKQ,F := Eone live (0, 01 (Ke.T, + TT) (Eico (1 (0..z)) - lko (DITAiZ ( ))
(01rz)) where Exo(O i (oNr))) is a vectorized notation for [Fko (01: x ) ... (ko
(0rz) ) E RTx and similarly for 1INe
(9122)) At level 0, YNKQ.O Ea-q [ke (0, 01: )] K K.rr Tode,o Fro (01:1.). The multi-level (ko(01::;:111N2
SSection 9 of (Giles, 2015) uses variance EE[?], which is equivalent the expected square error since Y; unbiased estimate of Eanq.Pe P
26

---

Nested Expectations with Kernel Quadrature
nested kernel quadrature estimator is constructed
L MMLKQ = Σ
YNKO,E-
Same as MLMC above, the cost of YNKQ, is (NNxxi)) The following theorem studies the error Pikq.e-1]l-
PPNKQ.
Theorem 2. Let X [0, 1jdx and O [0, 1, At E 2 1, 015.................... Sr, are Te i.i.d. samples from Q and x No are Ne i.i.d. samples from Fo, for all t E {1,... ,Te}. Both kemelss kx and ka are Sobolev reproducing kernlls of smoothness > dx/2 and se > d2.2. Suppose the Assumptions (1), (2), (3), (4) in Theorem I hold. Suppose
2 Ne-1 > No > Ne-1. Then, for sufficiently large Ne 2 and To 2 1, with Xx,t x No -2ax
4x
20x (log Ne) and
2re
d0, x T
1x 2x (log (log.N)
YNKQ, Ee [PNKe.. Pukq,e-1] S (Ne
xx N)xx
Te
20 dag
holds with probability at least 12e
The proof of the theorem is relegated to Appendix D.3.
2x
P = Under Theorem 2, the expected error El -aeiieeee---ellll F Ō(N,
xT, 2n 16
based Lemma 4,
up to logarithm terms. Here, the expectation is taken over the randomness of samples. Therefore, similarly to (D.36), the total cost and expected absolute error of multi-level nested kernel quadrature can be written as
L 1
L
2x ux
Cost = O
N₂
x
E M = Ô
Σ
N₂
x
T2
2 40
(D.37)
1=
If we take N₂ x 2x
Te x
20ea dia -
2ogs
2
then the error E
|
MMLKOl = Ö(A)) and the cost is
Σ
L
2
dx'
20g -re-daa
A
1
dx
de T
Σ 2( (xxx-2) -2))
2)0
A
dx 1
-rg--
Tag
de
O(A
-1
dx xx
de xe
Equivalently, to reach error O(A), the cost Ö(A -1
de
Remark D.1 (Comparison of MLKQ and MLMC). To reach given threshold A, the cost of MIKQ is Ö(A which is smaller than the cost of MLMC ((A when the problem has sufficient smoothness, i.e. when dx ax
-1
20x
dx de
2
+ de 2.
so
If we compare (D.36) and (D.37), the superior performance of MLKQ can be explained by the faster rate of convergence in terms of N₂ al each level when dx
< 1. Nevertheless, we can see in (D.37) that the MLKQ rate at each level in terms xx
Te is OTT, Tettes ) which is slower than the MLMC rate O(T,²) +) in (D.36). An empirical study of MLKQ is included in Figure 6 which shows that MLKQ is better than MLMC in some settings but both are outperformed by NKQ by huge margin. A more refined analysis of MLKQ is reserved for future work.
D.3. Proof of Theorem 2
The proof uses essentially the same analysis as in Step Five of Appendix C which translates |Y; Ee (Pr - Pr-1]] into the generalization eтror of kernel ridge regression. First, we know that by following the same derivations as in (C.29) that
Fxxo (0; (0: Pro. E W (0) and | S Co.
(0) = E [00 Pe
[iko (0:21.0) (0) ))}, FkQ.e-1 E (Q) and | S Co.
Fkq.e-1 (0) = co
~a
27

---

Nested Expectations with Kernel Quadrature
Given i.i.d. observations (01, Pko (01, (01) target of interest in the context of regression is the conditional mean, which in our case is precisely
1 ) - ( (0)
(0rz) 1 I:Ne=1 )),................ ")....,((0) (OT, ( & "1.Ne Fko ( F L1:Ne-1 )), the
0 + (0) - F (0) - E
"ifi ~Bo [ (0: - ziine ~ee
[ (0; (0:
can be viewed as noisy observation of the true function FKQ. - F where the noise satisfied the following condition. For each E and positive integer m 2 2, similar to (C.30) we have,
Alternatively, (f, x (0, wi:v ))
B [iko (0:211h) - lio (0:21.Ne.)) )]
I "LNN
Eko (0:I(:Ne) - 1E, "1:88-1
~Ea Fiee 00 (0:z(.Ne
~Po
S 2m E Fkoe (0:u(:N)) (0) - E Boye "ii ~Po
Pkoe (0:u(.Ne)
(0)
m
FRo 00:0 (.: E [000 "INNe- 2 E
(0:21.)..) Fko (0: (8)
~oo
N, +x (log N₂)
2xt1
-m 2x NoT1
(log No-1) +x
xx
+
1x+1 "x
S N (log No)
where the second last inequality follows by replicating the same steps in (C.30), and the last inequality is true because 2dx/sx N > N₂ > N As result, by replicating the steps for (C.31), we have
YNKQ. Ee P1
(FROO. Fx0,e-1) - ke (.01:1,) (Ke,r + -Tdd,elt.)) (Fzo (Orar (11,))) - Pkoo (Or.yi 1) L2(1)
(0::7,)
de
2x ax
22x12 1- de dx
de
20e dp
S 7² T TT NO.7
tee N₂
(log No)
+ N0,
20
| PrQu-1 1120 + A0,e"
2n
|kqq PKQ. 11
+ +|| - P Ae,e.
(D.38)
holds with probability at least 1 - TT Next, we are going to upper bound || |FRo.e - FKQ
To this end, notice that
FkQ. - Fl s F - 111 + +2|| - 11
Using the same steps in (C.29) and (C.27) subsequently, we have
||FK.. - -Fllc S |Fx (:=((f)))
F
S... Vol(AX) $ Jkke (zx((NNe) - - S N,
5x ax
(log No)
holds with probability at least 1 TT Similarly, we have || F'lls0 S NN
2x
(log Ne-1) "x
holds with
ex
ex+1 "x
probability at least 1 - 4T Consequently, we have || - F S Ne at least 1 Se Therefore, plugging it back to (D.38), we obtain
xx
(log No)
holds with probability
|NKOO - Eawe EEaee.NKo. - 1]12
?
To A
20
dg
N₂
28x
dx (log No)
22x12
+
1-
20
do
To N₂
2xx dx
(log Ne)
2x
4x
de De
a
26x (log No)
t2
2x xx
2ox12
+\g,
Te
de
N₂
dx
+ N₂ xx
(log Ne)
dx
holds with probability at least 17 Therefore, by taking d x T
2o
we obtain with probability at least
1 Se
YNKQ. Eang[PPKO Pikq S TTę 20 1d
N₂ (log No)
xx
The proof is concluded.
28

---

Nested Expectations with Kernel Quadrature
E. Further Background and Auxiliary Lemmas
All the results in this section are existing results in the literature. We provide them here and prove some of them in the specific context of Sobolev spaces explicitly for the convenience of the reader.
More technical notions of Sobolev spaces and the Sobolev embedding theorem In the main text, we provide in (4) the standard definition of Sobolev spaces W (2) when E N. Actually, Sobolev spaces Wx()) can be extended to that are positive real numbers. Such extension could be realized through real interpolation spaces (see (Bennett and Sharpley, 1988, Definition 1.7)). (x* := [W3 (2), W wk++1 (2)]n. where k E N, E (k: k + 1),r 8 [s]. Actually, such interpolation relations hold for any 0 I s, t and 0 < r < 1 (Adams and Fournier, 2003, Section 7.32),
W (x') (2) (2) n221
k= (1 (1-11117 rr--rsste +tt
(E.39)
A special case of the above relation is (x())= = [2)x)) (2)) 11/1122
The Sobolev embedding theorem (Adams and Fournier, 2003), when applied to (2')),, states that if $ (where d is the dimension of X¹). then (x) can be continuously embedded into (X)), the space of continuous and bounded functions. In other words, for every equivalence class [ʃ] E (A)), there exists a unique continuous and bounded representative f E (Z)), and the embedding map :: (X) + (2)), defined by [((f)) = f. is continuous. This contimuous embedding I can be written as (2)) + (x)). Since every continuous linear operator is bounded, we have ([) + (X)) bounded by a constant that only depends on S.
More technical notions of reproducing kernel Hilbert spaces (RKHSs) For bounded kernels, supzex x(x,x) S K, its associated RKHS H can be canonically injected into (1( using the operator i H - TT(TT), I + with its adjoint 6 (T) - H. given by f-) = k(zz,) (x) la and its adjoint can be composed to form L2(x) endomorphism T == L called the integral operator, and H endomorphism
Σx = G, / k(,2x) @ k(,x))d((z))
(E.40)
(where & denotes the tensor product such that (a 29 b)c = (b, c) for a, b, E H) called the covariance operator. Both 2i and Ta are compact, positive, self-adjoint, and they have the same eigenvalues 91 Section 2 of Chen et al. (2024a) for more details.
2..2. 2 0.
Please refer to
Lemma 1 (Effective dimension (()). Let X C Rdd be a compact domain, a be probability measure on X with density p: X + R. k. : X x X + R is a Sobolev reproducing kernel of order >> d {emm 1m are the eigenvalues of the integral 2.
operator Ta. Define the effective dimension N : (0, oo) + [0, oo) as (A) == Em
pm em-x If p(x)) 2 G >> for any
x. E X, then N(A) < DX-# # with constant D that only depends on G and X.
Proof. First, we study the asymptotic behavior of the eigenvalues (em) m21 of the integral operator T.. Theorem 15 of (Steinwart et al., 2009) shows that the eigenvalues 2m share the same asymptotic decay rate as the squares of the entropy number e?, (Iz) of the embedding I= : (x') + (T). Denote Lx as the Lebesgue measure on X. Since p(x) 2 G for any x E X, know dix
dr < G G- Vol so (TTT) + L2(X))l S Vo-V) -1, and consequently we have from Equation (A.38) of Steinwart (2008) that
Cm (I+) S em (Icx) (2( L2(x)) < Vo1 em (Ixx).
Moreover, (Edmunds and Triebel, 1996, Equation 4 on p. 119) shows that the entropy number em (Iex) S êm s/d for some constant é, so we have em (Ix) S G- cm -s/d and consequently we have Om x e (Tx) S G (=) 222 -2s/d_ =: 02m 22/d
Next, we have
Pm
/
00
dt = X X C2 M22/dd
C2
/ 00
C2 22/d dr
1 Ae, 1 2s/d
m21
2m
m21
J0
Strictly speaking, the definition of (4) extended to real numbers 8 actually corresponds to the complex interpolation space of Sobolev spaces. Fortunately, complex interpolation spaces and real interpolation spaces coincide under Hilbert spaces (Hytonen et al., 2016, Corollary C.4.2), which is precisely our setting since p = 2.
29

---

Nested Expectations with Kernel Quadrature
1 dr A
1
xd
A
C2
du
x 00
D)
1 1 2
sin
xd
1+ TC2
2x
where D is constant that depends on the domain X' and G.
Lemma 2. Let X Rd be a compact domain, x be probability measure on X with density p p: k Sobolev reproducing kernel of order operator Tx. If there exists Go, Gi > 0 such that Go p(x) S &Gifor any X, then
{em, Cm 1220 the eigenvalues and eigenfiunctions of the integral
ka
sup Σ ZEX
2m cm (x) M,
(E.41)
holds for any 2s < a. Here, M is constant that depends on and G1, Go-
Proof. If 4. (x)) can be continuously embedded into (X') the space of bounded functions (Adams and Fournier, 2003, Case A, Theorem 4.1 Hence, the operator (x') > Lo (2t) is continuous linear operator between two normed vector spaces, hence a bounded operator. And (T(T) is norm equivalent to L2(X)) because Go S p(x) < G1 for any æ E X. Notice that & defined here is exactly | lloo defined in Equation 16 (Fischer and Steinwart, 2020), so we know from Theorem of Fischer and Steinwart (2020) that
sup Σ (x) [L2(T), W2()) + (x)
TEA
m21
Notice that [L2( (3)))),, (z) (x))). (X), and notice the fact that W)'(X) + Loo (2) for any z. the right hand side of the above equation is bounded. Therefore, have (E.41) holds for any 2x a.
d
Lemma 3. Let X C Rdd be a bounded domain with Lipschitz continuous boundary and W2(() be a Sobolev space with 2- If functions X - R and g X TR. lie in (X), then their product 9 also lies in W W2(X) and satisfies Ilf glla I ||glle
Proof. This is Theorem 7.4 of Behzadan and Holst (2021) with 81 82 and pi pp 2
Lemma 4. For a positive valued random variable R, and c. such that P(R S cT) 21 exp(-r) ( for any positive т. holds that E Rm S coml for all integers m 2 1. Co is some constant that only depends on c, m.
Proof. Notice that R is essentially sub-exponential random variable. Since a sub-exponential random variable equivalent to the square root of a sub-Gaussian random variable, from Proposition 2.5.2 of Vershynin (2018), we have
E[R"M E 2m, < (m + 1) 2c.m! Here I' denotes the gamma function and co is some constant that only depends on
Lemma 5. For mapping F from compact domain X C Bd Hilbert space H, given measure a X, if F is w-Bochner integrable, then F(x)dp(x) H and additionally J F(z)du(z)lla < ([(l)))))dd)))
C
Proof. This is Definition A.5.20 of Steinwart (2008).
F. Additional Experimental Details
F1. "Change of Variable" Trick for Kernel Quadrature
In the main text, we have shown that the two major bottlenecks of KQ/NKQ are:
The closed-form KME Ex-E [x((,,))) x)]-
The (NN) computational cost of inverting the Gram matrix k(2 R1:N).
30

---

Nested Expectations with Kernel Quadrature
Ay = 0.001
Ag 1.0
0.01
Ap 10.0
D.I
10.0 median
Matern-1.5
Matern0.3
0.1
7 1.0
Gaussian
10
10
10
16-4
10
10
10
10"
10-3
10
10
10
10
10 10
10
10 Cost
30
10
10"
10
10
10
10
107
10
Cost
Cret
Figure 5. Further ablation studies in the synthetic experiment. Left: NKQ with different proportionality coefficients A for regularization parameter Ax, Àa- Middle: NKQ with different kernel lengthscales y both stages. Right: NKQ with different kernels in both stages. The nested Monte Carlo (NMC) in blue is presented as benchmark in all figures.
Fortunately, both two challenges be solved with the "change of variable" trick. Here, we only present the idea for KQ but the same holds for NKQ in both stages.
The integral of interest is 1 = fx h(z) Suppose we can find continuous transformation & such that X = (UU), where U ~ U is another random variable which is easy to sample from. Then the integral can be equivalently expressed as 1 = Su h (((()))U(()) by direct application of change of variables theorem (Section 8.2 of (Stirzaker, 2003). Now the integrand changes from X' + R to h 4b 21 +
and the kernel quadrature estimator becomes
IxQ- = Bu-uilul (U, uINN)] (uu(NN ) + NAIx) (h 1 (uEN))
Here ku is reproducing kernel on U. Since U is simple probability distribution, can find its closed-form KME in Table in Briol et al. (2019) or the ProbNum package (Wenger et al., 2021), which addresses the first challenge. Additionally, notice that both the Gram matrix k(uINN, uI.N) and the KME Eu--r[K(U,, wi:)) are independent of l and (k(U, M1:N)] (k((1,, "1.N), + NAI x) can be pre-computed and stored. As
KQ T so the KQ weights w WI:N- result, KQ becomes a simple weighted average of function evaluations Σn= w h(x,). Therefore, the computational cost reduces linear cost (NN) and hence the second challenge is addressed. The downside of the "change of variable" trick is that the Sobolev smoothness of h & L - IR is unclear when + is not smooth, lose the theoretical convergence rate from Theorem 1.
E.2. Synthetic Experiment
Assumptions from Theorem 1 We would like to check whether the assumptions made in Theorem hold in this syn- thetic experiment. Recall that we use both kx and kis to be Matérn-3/2 kernels so we need to verify Assumptions (1) (4) with se sxx = 2.
1. Both distributions Po and Q are uniform distributions over [o, 1]. so Assumption (1) is satisfied.
2. D 0) E (2') and p) 0) E (22) for 3 = 0, 1, 2 so Assumption (2) is satisfied.
3. Both g(x, -), p(x,)) E (e) so Assumption (3) is satisfied.
4. f (RR) so Assumption (4) is satisfied.
The synthetic problem can be modified to have higher dimensions d. this synthetic experiment, we set both dx = dg
d. For a = jai....,.a. oa] E Rd, define ||a]|: = (Ei a?) l/b.
U(O, 114 0 ~ 1J,, g(x,0) 0) = || + |0112-5 f(z) =2
(F.4)
The true value of the nested expectation can be computed in closed-form: mean absolute error NMC and NKQ as dimension d grows. We see that NKQ outperforms NMC by a huge margin in low dimensions, but the performance gap closes down higher dimensions, which is expected because the rate proved in
16 2
25 d. In Figure 2, we study the
49 291
31

---

Nested Expectations with Kernel Quadrature
Corollary 1 is O(A Xx the same.
which becomes larger when dimension increases yet the smoothness of the problem remains
In Figure 5, we conduct a series of ablation studies on the hyperparameter of NKQ in the synthetic experiment. Although Theorem suggests choosing the regularization parameters Ax, λe that are proportionate to N- - and T respec- tively, it is unelear in practice how to pick the exact proportionality coefficients Ag. Figure 5 Left shows that À = 1.0 and Au = 0.1 give the best performances, while using Àg too big (Au 10.0) suffers from slower convergence rate and using An too small (Ao = 0.01, 0.001) causes numerical issues when N, T become large. Figure 5 Middle shows that kernel lengthscale, if too big (7 = 10.0) or too small (7 = 1.0), would result in worse performance for NKQ and that the widely- used median heuristic is good enough to select satisfying lengthscale. Figure 5 Right shows that NKQ with Matérn-3/2 kernels has better performance than with Matérn- 1/2 kernels, which agrees with Theorem I indicating that it is preferable to use Sobolev kernels with the highest permissible orders of smoothness. Interestingly, we see that NKQ with Gaussian kernels has similar performance as with Matémn-3/2 kernels. Similar phenomenon have been shown both thcoretically and empirically that kernel ridge regression with Gaussian kernels are optimal in learning Sobolev space functions when the lengthscales are chosen appropriately (Hang and Steinwart, 2021; Eberts and Steinwart, 2013).
F3. Risk Management in Finance
In this experiment, we consider specifically an asset whose price S(T) at time r follows the Black-Scholes formula S(() So exp (r) d²r/2) for 2 0, where a is the underlying volatility, So is the initial price and W is the standard Brownian motion. The financial derivative we are interested in is butterfly call option whose payoff at time T can be expressed as (()()) = max (S()) K1, 0) + max(S(T) K2, 0) 2 max(S(T) - (K1 + K2)/2, 0). We follow the setting in (Alfonsi et al., 2021; 2022; Chen et al., 2024b) assuming that a shock occur at time 7, at which time the option price is S1)) = 0, and this shock multiplies the option price by + 8. The option price maturity time Ç is denoted as S(C) = z. To summarize, the expected loss caused by the shock can be expressed as the following nested expectation:
((((())) S(J()) = max (J(0).0), J(O) = f g(x))P((d)), g(x) = <(x) (5((1 + s)x).
Following the setting in (Alfonsi al.. 2021: 2022: Chen et al., 2024b), we consider the initial price So 100, the volatility 0.3, the strikes K₃ 50, K₂ 150, the option maturity Ç 2 and the shock happens 1) 1 with strength 0.2. The option price at which the shock occurs are OuT sampled from the log normal distribution deduced from the a"
Black- Scholes formula 01. Q = Tognormal(log So 2 1), an)) Then (!()) are sampled from another log normal
distribution also deduced from the Black-Scholes formula 11:N Po. Lognormal(log 0. o² (C n), 02 (P((())) n)) for t= I.
(1)
In this experimental setting, although both only depends on and it is combination of piece-wise linear functions so E W3. (2)). The probability density function of Pa is infinitely times differentiable
Notice that log normal distribution LogNormal o²) can be expressed as the following transformation from uniform distribution over [0, 1].
u ~ U[,, 1], cxp( (u)o + m) ~ LogNormal(m,
Here, -p-1 is the inverse cumulative distribution function of standard normal distribution. Therefore, we can use the "change of variables" trick from Appendix F.I such that we have closed-form KME against uniform distribution from Probnum (Wenger et al., 2021), and also the computational complexity of NKQ becomes O(N x T). Although 1p-1 is infinitely times differentiable, we still use Matém-1/2 kernels in both stages to be conservative of the smoothness of the integrand after applying the "change of variables" trick.
F.4. Health Economics
In the medical world, it is important to compare the cost and the relative advantages of conducting extra medical exper- iments. The expected value of partial perfect information (EVPPI) quantifies the expected gain from conducting extra experiments to obtain precise knowledge of some unknown variables (Brennan et al., 2007):
EVPPI E max [0] max EE [..(O)] Jz(0) / ge(z, ()F
32

---

Nested Expectations with Kernel Quadrature
NMC
NMC (QMC)
MLMC MLKQ
NKQ
NKQ (QMC)
Syuthotic Experiment
Risk Management Finorer
Health Ecomcenics
10"
10-4
10P
10
10:
10-4
101
10
107
10
Cost
Cost
Coet
Figure 6. Comparison of all the methods including MLKQ on the synthetic experiment (Left), risk management in finance (Middle) and health economics (Right).
where e is set of potential treatments and g. measures the potential outcome of treatment c. EVPPI consists of [C] + nested expectations.
We adopt the same experimental setup delineated in (Giles and Haji-Ali, 2019), wherein and 0 have joint 19- dimensional Gaussian distribution, meaning that the conditional distribution I& is also Gaussian. The specific meanings of all x and 0 are outlined in Table 2. All these variables are independent except that 01, 02, Ta, $14 are pairwise correlated with a correlation coefficient 0.6. We are interested in estimating the EVPPI of a binary decision-making problem (C = (1, with 91 (2r, 0) (0 a0g.2g + (77 - (a1 + x223344) and (2 0) = 10 102 + 2154166717 (x'0 + #11W11244). The ground truth EVPPI under this setting is 538 provided in (Giles and Goda, 2019).
For estimating I with NKQ, select xx be Gaussian kernel and ke to be Matérn 1/2 kernel, because I contains max function which breaks the smoothness so we use Matérn-1/2 1/2 kernel to be conservative. For estimating 22,- with NKQ, we select both to be Gaussian kernels because both 91 and the probability densities are all infinitely times continuously differentiable. We have access to the closed-form KME for both Matémn-122 and Gaussian kernels under a Gaussian distribution from Probnum (Wenger et al., 2021).
F.5. Bayesian Optimization
For NKQ, we use the change of variable trick such that the Gaussian distribution of JIDs (21, 22) after S iterations can be expressed as the pushforward a fixed uniform distribution U over [0, 1J2 through continuous mapping &s. As a result, the KQ weights Eu-u(le Mi (U, W1:x)] )] (hu ("),N "1.N) + NAIx) -I become independent of S, and can therefore be precomputed and stored in advance. We apply this Shange-of--vrialte trick to both Stage and Stage II KQ steps in our NKQ algorithm, resulting an overall ONN x T) computational cost at each iteration, matching that of NMC.
The formulas of the synthetic Dropwave, Ackley, and Cosine8 functions are provided below:
1 cos 12 12v/>² +
firopwave (z,y)= = 0.5 (x² + y²) +
(x, E (-5.12,.5112],
2
fAcMley (x) = 20 exp ( (-0.2 ||) exp
cos (2TTa)) - 20 + exp(1), x (-32.768, 32.768)²
2
fCosine ((x)= Σ
cos (JTTz;) I E 1-1.1]
1=1
33

---

Nested Expectations with Kernel Quadrature
Variables
Mean
Std
Meaning
Xi
1000
1.0
Cost of treatment
X2
0.1
0.02
Probability admissions
X3
5.2
1.0
Days of hospital
X1
400
200
Cost per day
Xs
0.3
0.1
Utility change if response
X6
3.0
0.5
Duration of response
X7
0.25
0.1
Probability of side effects
Xs
-0.1
0.02
Change in utility if side effect
X,
0.5
0.2
Duration of side effects
X10
1500
1.0
Cost of treatment
X1
0.08
0.02
Probability of admissions
X12
6.1
1.0
Days of hospital
X13
0.3
0.05
Utility change if response
X14
3.0
1.0
Duration of response
X15
0.2
0.05
Probability of side effects
X16
-0.1
0.02
Change in utility side effect
X17
0.5
0.2
Duration of side effects
0.7
0.1
Probability of responding
0.8
0.1
Probability of responding
Table Variables in health economics experiment.
34


Reviews:
[{"review_id": "2WfOdTACiK", "reviewer": "Reviewer_4w5G", "summary": "The paper introduces Nested Kernel Quadrature (NKQ), which a novel method for estimating nested expectations (i.e., integrals where the outer expectation involves a function of an inner expectation). Compared to existing methods such as Nested Monte Carlo (NMC) and Multilevel Monte Carlo (MLMC), which are suitable for functions with relatively more generic assumptions on the smoothness, NKQ considers the case of stronger regularity assumptions and exploits smoothness in the integrand to achieve faster convergence rates by applying existing results from the kernel quadrature (KQ) literature to both the inner and outer expectations. Both theoretical proofs and numerical experiments on synthetic problems, Bayesian optimization, option pricing and health economics are provided to justify the findings.\n\n## Update After Rebuttal\nThe reviewer agrees that the feasibility of the proposed method needs to be discussed from a practical perspective in the revised version of the manuscript. However, the reviewer is inclined to keep the original score to accept this paper, as most of the theoretical contribution are solid and sound from the reviewer's perspective.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you very much for your careful consideration of our paper and for checking the proof of the main theorem. We are very happy that you agree that our convergence rate results (Theorem 1) are technically sound and that our experiments are well-designed.\n \nThank you for also bringing up the issue of minimax rates and for highlighting papers we were unaware of —this is indeed a fundamental and important question. To the best of our knowledge, most existing approaches on minimax rates for two-stage algorithms first reduce the problem to a one-stage regression by assuming that the target of the second regression is already known. For instance, in the suggested Zhang et al (2025 ) paper, the minimax lower bound provided requires the assumption that the oracle density ratio between the source and target distributions is known (Theorem 1). Similarly, in the literature on nonparametric instrumental variable regression— see for example Theorem 4 and Appendix F.1 in https://arxiv.org/pdf/2411.19653, and Chapter 3 of https://arxiv.org/pdf/0709.2003 — minimax lower rates are provided under the assumption that the conditional expectation operator is known, even though in practice their proposed algorithms would learn this operator from data.\n \nTranslating this approach to our setting of nested expectations, we could assume access to the true value of the inner expectation. This effectively reduces the problem to a standard (single-stage) expectation, for which the minimax optimal rate is well-known: $\\mathcal{O}(T^{-s_\\Theta/ d_\\Theta})$ [Sommariva and Vianello, 2006a, Section 2 of https://arxiv.org/pdf/2307.06787]. We agree with you that this is an interesting point to discuss and will do so in the camera-ready version of our paper.\n \nThat said, it clearly remains an open challenge to establish minimax lower bounds for genuinely two-stage problems without reducing them to one-stage settings. Developing such a theory would be highly nontrivial and beyond the scope of our paper, but also a very exciting direction for future research.", "claims_and_evidence": "Here is the main claim made in this paper: NKQ achieves faster convergence rates compared to existing methods like NMC and MLMC when the input functions are sufficiently smooth. Both theoretical and empirical results are included to support the claim. Theoretically, the authors proved Theorem 1 and Corollary 1 to establish the faster convergence rates under general regularity assumptions. Empirically, the authors have conducted extensive experiments on synthetic problems, Bayesian optimization, option pricing and health economics to verify that NKQ performs better on sufficiently smooth functions compared to standard methods like NMC and MLMC.", "methods_and_evaluation_criteria": "Yes. Please refer to the \"Experimental Designs or Analyses\" section below for a detailed discussion.", "theoretical_claims": "Yes. The convergence rate result (Theorem 1) and its derivation are technically sound under the stated assumptions. The main idea is based on decomposing the error into inner and outer quadrature errors and applying existing bounds from the kernel quadrature literature. The assumptions on smoothness, boundedness, and properties of the kernel are well-justified.", "experimental_designs_or_analyses": "The experiments are well-designed and cover a range of scenarios, such as synthetic experiments that validate theoretical rates, applications in finance, health economics, and Bayesian optimization, ablation studies on choices of kernels and regularization parameters, as well as careful comparison with baseline methods like NMC, MLMC and NQMC.", "relation_to_broader_scientific_literature": "This paper, which proposes the NKQ method and briefly shows that NKQ can also be combined with other advanced Monte Carlo methods like MLMC, should be mainly situated as a combination of Monte Carlo methods and kernel quadratures.", "essential_references_not_discussed": "To the best of the reviewer's knowledge, all essential references are properly discussed in the manuscript. For some other references that are loosely related, the authors are encouraged to take a look at the series of work listed in the \"Other Comments or Suggestions\" section below.", "other_strengths_and_weaknesses": "The main strength of this paper is that it provides not only solid theory but also extensive experiments to justify the findings. One potential weakness is that in terms of theoretical contribution, it might be beneficial to discuss whether the proposed estimator is minmax optimal or not by deriving an information theoretic lower bound. For instance, one may refer to the bounds derived in [1,2]\n\nReferences:\n\n[1] Zhang, Zhen, Xin Liu, Shaoli Wang, and Jiaye Teng. \"Minimax Optimal Two-Stage Algorithm For Moment Estimation Under Covariate Shift.\" In The Thirteenth International Conference on Learning Representations.\n\n[2] Blanchet, Jose, Haoxuan Chen, Yiping Lu, and Lexing Ying. \"When can regression-adjusted control variate help? rare events, sobolev embedding and minimax optimality.\" Advances in Neural Information Processing Systems 36 (2023): 36566-36578.", "other_comments_or_suggestions": "N/A", "questions_for_authors": "N/A", "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": {"ethical_review_concerns": "N/A"}}, {"review_id": "CQwoiMJ0d8", "reviewer": "Reviewer_uEeW", "summary": "This paper presents a method to estimate nested expectations using (essentially) bilevel application of kernel quadrature. The central result is that this gives a convergence rate that depends on the smoothness constants for each of the two functions in the expectations. This recovers the best known rates in worst-case scenarios, but can be more data-efficient. Separate to that, ideas are proposed to make *computationally* efficient, essentially by changing variables in the integral so it is possible to pre-compute the weights needed for the estimator.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you very much for your consideration and strong support of our paper (including the very kind comment about the effort we put in the writing). You mention leaning towards a ‘strong accept (5)’ and we remain at your disposal in case there is anything we can do to help in this respect.\n \n## 1. Acronyms\nThank you for pointing this out. We will do our best to limit the use of acronyms to concepts which are repeatedly used, and will check that all acronyms have been properly defined.\n \n## 2. Knowledge of constants\nThank you for this interesting and important question – we will certainly add a discussion of this point in the camera-ready version of the paper. Firstly, most constants in Theorem 1 actually have no impact on algorithm design. The constants named $S$ or $G$ only impact $C_1$ and $C_2$, and $C_1$ and $C_2$ do not impact the use of the algorithm in any way.\n \nThe only important constants for the algorithm are the smoothness parameters $s_\\mathcal{X}$ and $s_{\\Theta}$, which help us identify the appropriate level of smoothness to use for kernels $k_{\\mathcal{X}}$ and $k_{\\Theta}$. This is very common in this literature [Sommariva and Vianello, 2006a, Briol et al., 2019], and even poor choices of these smoothness parameters does not necessarily mean the rate will be impacted [Kanagawa et al., 2020, Wynne et al., 2021]. Although we initially mentioned this very briefly below Theorem 1, we will now use the extra space available to us in the camera-ready version to discuss this more extensively.", "claims_and_evidence": "Yes, both clear and convincing.", "methods_and_evaluation_criteria": "Yes, the evaluations methods are appropriate. The empirical evidence isn't huge, but it's very well done and I think more than adequate given the theoretical contribution.", "theoretical_claims": "I only briefly perused the proofs in the appendix. They look fine, but I didn't really check.", "experimental_designs_or_analyses": "Not sure what I would check here. The experiments/analyses seem fine.", "relation_to_broader_scientific_literature": "Nested expectations are an important variant of monte carlo estimation, a classical problem of immense interest. Nested expectations are important enough that even small improvements/insights are worthwhile. I believe this paper easily clears that bar.", "essential_references_not_discussed": "I am not aware of any works that are not cited.", "other_strengths_and_weaknesses": "This is an extremely minor point that should have no impact on anything, but the paper is pretty heavy with acronyms. In a few cases (e.g. MLMC) I couldn't actually find the original definition. Consider reducing these to the few that are standard and/really used repeatedly (QMC, RKHS, MC, KQ, KME) and/or providing a key in the appendix.\n\nOverall, I think this is a good paper and recommend acceptance. The only reason I shy away from strong accept is the fact that I'm not an expert in this area.", "other_comments_or_suggestions": "I'd like to remark that I think this paper is very well-written. It's quite technical, but the authors have done a truly excellent job and I think it's nearly as clear and easy to read as is possible given the content. I really appreciate this!", "questions_for_authors": "One question I had was to what degree is \"knowledge of constants\" needed? Which of the constants mentioned in Thm. 1 are actually used in the *construction* of the estimator, and not in just the analysis? I'd love to see some analysis of how this estimator compares to previous work in different \"constant knowledge regimes\". If no smoothness constants are known, does the guarantee reduce to previous bounds?\n\n## Update\n\nAfter looking at the other reviews, I still recommend accepting this paper. That said, I do agree with the concerns of other reviewers that the scope of the method isn't as large as the authors claim, due to the need to reparameterize. (I have read the authors' response to this, but I do not find it very convincing.) To be clear, I think this is *fine*. I think that the method applies to many problems of real interest, so I do not see this as a problem that needs to be solved. But I agree that the paper should be modified to make this limitation more clear, to avoid over-claiming. This is my only concern.\n\n(To be clear, I still recommend acceptance. My rating of 4 should be considered to be already taking this into account—I would probably recommend 5 if only the paper had made the limitation in scope more clear.)", "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": "MavDOBLcqo", "reviewer": "Reviewer_duxH", "summary": "This paper proposes a novel kernel quadrature estimator for nested integrals involving two integration levels: the inner integral $g(X, \\theta)$, which computes the conditional expectation over $X$, and the outer integral, which evaluates the expectation of this conditional expectation over $\\theta$. We establish an improved convergence rate of the mean squared error for this nested integral estimator. The effectiveness of the proposed method is demonstrated through synthetic experiments and three real-world applications.\n\n## update after rebuttal\nI share Reviewer fCUr’s concerns. The computational aspects remain unclear to me. I maintain original score 3.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your consideration of the paper. We are very glad to see that you think the paper is very clearly written and appreciate our effort to include real-world applications. We have carefully addressed all your comments below, but please do let us know if there is anything we can do to help.\n \n## 1. The core idea appears incremental\nWhile the core idea may appear simple from an algorithmic standpoint, we view this simplicity as a key strength rather than weakness. In practice, we believe that methods that are conceptually clear and easy to implement are  much more likely to be adopted. In this specific instance, we show that our proposed method can significantly outperform existing approaches despite this algorithmic simplicity, and believe there is therefore clear evidence of the value of our proposal.\n\nThe other very strong element of novelty is in terms of the theoretical guarantees we are able to obtain (and the associated proof techniques). Although we appreciate this is something that not all users of our method might be interested in, it is still highly relevant to practitioners. Specifically, we prove a significantly faster convergence rate in $T$ than what is suggested by the existing literature (e.g., [Chen et al., 2024b]). Our rate of $T^{-s_\\Theta/d_\\Theta}$ notably improves upon the previously established rate of $T^{-1/4}$. This means that one does not need to take $T$ as large as existing results might suggest, which is important for practitioners who might otherwise have wasted significant computation on trying to approximate the outer expectation.\n\n## 2. Application to other acquisition functions in Bayesian optimization\nThank you very much for raising this as we were actually not aware that these other acquisition functions were also nested expectations. The reason we decided to focus on the setting of look-ahead optimisation is that there is an entire paper on this topic (in fact written by some of the authors of the BoTorch package) [Yang et al., 2024]. We therefore assumed that this example would be particularly relevant to the BO community. We will add a broader discussion of the applicability of our method to these additional acquisition functions to the camera-ready version.  \n\n## 3. Change of variable trick\nThank you very much for your thoughts on this. We agree that extending our presentation and discussion of this trick could be very helpful and will do this in the camera-ready version. In the meantime, see below for an answer to each of your sub-comments relating to this.\n \n1.   ### Applicability of the change-of-variable trick \nWe want to clarify that the change-of-variable trick is in fact applicable to all experiments in our paper. We do however agree with you that it is less helpful when the number of function evaluations is already inherently limited since the cost is relatively low in that case anyway. In this paper, the finance and health economics experiments were selected specifically because they are expensive, and so having a better estimator such as NKQ is particularly important. There are however numerous related settings where this is not the case and the change-of-variable trick could still be useful.\n \n2. ### Cubic complexity\nBO with Gaussian processes indeed requires cubic cost at each iteration, however this cost is cubic in the number of queries of the BO algorithm, not the number of samples $N$ or $T$ for estimating the nested expectation (i.e. the acquisition function). Therefore, precomputation reduces the cost in $N$ and $T$ for optimising the acquisition function at each iteration of the BO algorithm, while the cost of GP posterior update is still cubic and is the same across all methods regardless of how the acquisition function is being approximated. We will make this distinction clear in the revised manuscript.\n \n3. ### Popular libraries already cache quadrature weights, and I presume the authors intend to leverage such caching to avoid repeated computations\nWe do indeed use the idea of caching for the weights of NKQ in estimating the acquisition function so that we do not need to re-compute the weights at each iteration. We are not, however, using NKQ when computing the marginal likelihood in GP hyperparameter selection because it is a standard expectation not a nested expectation.\n \n4. ### Comparison to QMC\nWe conduct an additional experiment with Quasi Monte Carlo ($N=T=\\Delta^{-2}$), as shown here https://github.com/Anonymous65536/nest_kq/blob/master/bo.pdf,  demonstrating that while QMC outperforms standard Monte Carlo in estimating the nested expectation within the utility function, it is less efficient than our proposed NKQ method across all three Bayesian optimization datasets. \n\n## 4. Equating kernel quadrature directly with Bayesian quadrature\nThank you for pointing this out. We of course agree and will further highlight the distinction between frequentist and Bayesian approaches in the camera-ready version of the paper.", "claims_and_evidence": "Their primary theoretical contribution regarding the improved convergence rate is established in Theorem 1. The empirical advantages of the proposed estimator are demonstrated in Figures 2–4.", "methods_and_evaluation_criteria": "The absolute error of the integral estimate, defined as $\\Delta = \\lvert I - \\hat{I} \\rvert$, is a standard metric in the kernel quadrature literature, and evaluating its convergence rate is a common criterion in theoretical analyses.", "theoretical_claims": "Although I did not verify the correctness of the proof in detail, the results appear reasonable given the known optimality properties of kernel quadrature estimators.", "experimental_designs_or_analyses": "The experiments are conducted systematically, and the analysis appears reasonable.", "relation_to_broader_scientific_literature": "Numerical integration forms the foundation of probabilistic machine learning and plays a critical role in numerous scientific and engineering applications.", "essential_references_not_discussed": "One of the strengths of this paper is its comprehensive reference list. Although the references are somewhat biased toward statistical approaches, I find the current coverage satisfactory.", "other_strengths_and_weaknesses": "Strengths\n- The paper is very clearly written. Due to the inherent complexity of nested integration tasks, however, a summary table of notations in the appendix would be beneficial. Occasionally, it was challenging to recall the precise meaning of each notation.\n- The effort to include real-world applications is commendable, especially given that many papers in this community typically limit themselves to theoretical results without real-world experimental validation. Our community would greatly benefit from increased attention to practical applications beyond purely statistical interests.\n\nWeakness\n- The core idea appears somewhat incremental. Although the improved convergence rate is valuable, the idea essentially is applying kernel quadrature twice. Given that kernel quadrature is already known to outperform Monte Carlo and Quasi-Monte Carlo methods in terms of convergence rate, it is unsurprising that its nested version would exhibit superior performance without needing empirical validation. Thus, this work feels more like an “exploitation” of known results rather than an exploration of novel concepts. Readers with a statistical background might enjoy the proof technique, but the fundamental idea itself does not seem particularly innovative.\n- Additionally, the application to Bayesian optimization was somewhat disappointing. When reading the introduction, I expected applications to prominent information-theoretic acquisition functions, such as Predictive Entropy Search (PES), Max-value Entropy Search (MES), or Sample Average Approximation (SAA)-based acquisition functions proposed in the original BoTorch paper [1] and Wilson et al. [2]. These functions involve nested integrals, are widely used in the Bayesian optimization community, and improvements here would have significant practical impact. Instead, the authors chose a two-step look-ahead acquisition function, which, although intriguing, lacks popularity in both practical settings and theoretical studies within the Bayesian optimization community.\n- I found the change-of-variable trick interesting, but it requires further elaboration. Specifically, the claimed computational efficiency does not seem particularly relevant for the first two examples in finance and healthcare with expensive querying, making it important to clarify its effectiveness in the iterative Bayesian optimization context. The authors state that the computational complexity is linear, $\\mathcal{O}(N \\times T)$, assuming pre-computation of quadrature weights with the Gram matrix inversion. However, Bayesian optimization involves iterative updates of Gaussian process surrogate models, typically incurring cubic complexity per iteration. It remains unclear how pre-computation mitigates this inherent complexity. Given that popular libraries like GPyTorch and BoTorch already cache quadrature weights during hyperparameter optimization via marginal likelihood, I presume the authors intend to leverage such caching to avoid repeated computations in the acquisition function estimation—am I correct in this interpretation? If the proposed method indeed offers computational advantages comparable to simpler approaches like quasi-Monte Carlo, this should be explicitly demonstrated through experiments for wider audience.\n\n- [1] https://arxiv.org/abs/1910.06403\n- [2] https://arxiv.org/abs/1805.10196", "other_comments_or_suggestions": "- From a Bayesian perspective, I am slightly hesitant about equating kernel quadrature directly with Bayesian quadrature. Although these terms are often used interchangeably in the literature, their foundational assumptions and underlying interpretations differ significantly. I believe the authors are likely aware of this nuance, but highlighting this distinction explicitly would be beneficial to readers who, like myself, appreciate the philosophical differences between these approaches. \n- The references should be refined for consistency and accuracy; for instance, \"Conditional bayesian quadrature\" should be corrected to \"Conditional Bayesian quadrature,\" capitalizing \"Bayesian\" appropriately.", "questions_for_authors": "See weakness.", "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": "rIEoIK5Kji", "reviewer": "Reviewer_fCUr", "summary": "The authors study the problem of estimating nested expectations, i.e.,\nexpectations over two entities where we can condition on one of the entities\n$E\\[g(X, \\Theta)\\] = E\\[E\\[g(X, \\theta) | \\Theta = \\theta\\]\\]$. These\nexpectations are common in Bayesian statistics. While there are existing\ntechniques to estimate such expectations via Monte Carlo (e.g. nested QMC or\nMLMC), the authors propose their own method nested kernel quadrature (NKQ)\nas an alternative approach. The authors demonstrate that when the integrands\nare sufficiently smooth, their proposed method can have a faster convergence\nrate than aforementioned methods. The authors then employ their method in a\nfew different applied settings, comparing the accuracy of their method to\nalternatives. They demonstrate that their method does indeed perform well,\nmatching the rates found in their theoretical analysis.\n\n## update after rebuttal\nBased on the authors feedback and other reviews, I'm inclined to keep my score of a weak reject. I think the answers have not accurately explained the applicability and computational resources required to apply this method, and without a more moderate explanation of these assumptions, I am not inclined to accept this paper.", "questions": "", "limitations": "", "rating": 2, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your consideration of the paper. We are very glad to see that you think the idea of using NKQ to estimate these expectations is interesting and the paper is well-written and easy to follow.\n \n## 1. The pushforward map for change of variable trick\nThank you for your question about this. We agree that we could have fleshed out a bit more the details of this trick (see also our response to reviewer duxH) and now plan on providing an extensive discussion using the additional space available for the camera-ready version.\n \nThis approach is actually general and widely available. For example, for one dimensional distributions, a simple pushforward map is the inverse CDF, in which case the base measure is the uniform distribution. This is typically called inverse transform sampling. Such pushforward maps can be obtained through a multitude of simulation tricks which are widely popular in the Monte Carlo literature; see [Chapter 2 of https://www.cs.fsu.edu/~mascagni/Devroye.pdf] for a very detailed discussion. Additionally, the idea of using such pushforward transformations is the foundational idea behind the field of simulation-based inference, in which case the pushforward map is simply called a ‘simulator’ or ‘generator’ [https://arxiv.org/pdf/1911.01429]. In this particular case, we note that the literature considers exactly the setting where the distribution $p$ is intractable, although it might not cover all cases that you mention.\n \nIf the distribution $p$ is known up to a normalization constant, instead of using the change of variable trick, we have a nicer alternative than the change of variable trick by using the Stein reproducing kernel which constructs a new kernel with a standard base kernel and the score function of $p$; see the third bullet point of Section B.1 in https://arxiv.org/pdf/2406.16530 of how it can be applied for kernel / Bayesian quadrature.  \n\nWhen the pushforward map does not exist, our method would incur the computational burden of inverting the Gram matrices. Fortunately, a wide range of techniques already exist in the literature that aim to reduce this cost. See the detailed discussion between lines 138–152 in the right column of page 3 of the paper. We will add a clear sentence in the paper emphasizing the necessity of using these approximate methods in this case.\n\n## 2. Minors: typos and big O notations\nThank you for catching these, we will fix these in the revised version of the paper.", "claims_and_evidence": "The authors provide a thorough theoretical analysis demonstrating the\nsuperior convergence rate of their method. This is also accompanied by three\ndifferent empirical studies that corroborate the efficacy of their approach.", "methods_and_evaluation_criteria": "Yes, the datasets and benchmarks used in the paper are reasonable\nexamples to demonstrate the utility of their method.", "theoretical_claims": "There is a section that would merit more discussion. While the authors\ndescribe their \"change of variables trick\" at the end of page 3 and in\nappendix F.1, they do not discuss how difficult it is to determine the\npushforward map $\\Phi$ for a given problem. This could be flushed out more,\nand is the main weakness of the paper.", "experimental_designs_or_analyses": "Yes, the empirical studies appear reasonable. The demonstrate the soundness\nof their method on some examples which do not perfectly satisfy their\nassumptions.", "relation_to_broader_scientific_literature": "The method described here could provide a more efficient approach for\npractitioners of applied statistics or machine learning when one needs to\nestimate an expectation over two entities, assuming all their assumptions\nhold (more on this below).", "essential_references_not_discussed": "To the best of the reviewer's knowledge, there aren't any missing references.", "other_strengths_and_weaknesses": "The idea of using a nested quadrature to estimate these expectations is\ninteresting and worth the analysis demonstrated in this paper.\n\nThe paper is very well written and easy to follow. \n\nHowever, the biggest issue with this paper is the significance of this\nmethod. The authors do point out that using quadrature methods can have two\nbottlenecks for large $N$: they must compute the kernel mean embeddings and\nalso invert the $N\\times N$ Gram matrix, which can be $\\Theta(N^3)$. While\nthe authors do propose the \"change of kernels\" trick as a workaround, they\ndo not discuss the limitations this imposes on the distributions $p$ and\n$q$.\n\nWhile it is quite easy to ascertain the pushforward map $\\Phi$ when $p$ and\n$q$ are univariate or multivariate Gaussians, it does not appear very simple\nto find this mapping for any general multivariate $p$ or $q$. When $p$ is a\nmultivariate distribution with an unknown normalization constant, one\ntypically employs MCMC to estimate expectations. If it were so easy to\nacquire $\\Phi$, MCMC would not be necessary.\n\nThe authors could do a better job explaining in what circumstances it is\npossible to obtain $\\Phi$, or at least warn that in most cases, one must\nactually use the approximate methods (e.g. Nystrom methods) already\nmentioned in the paper.", "other_comments_or_suggestions": "There are some minor typos (search for the word \"considere\") and also\nsome missing periods in Appendix F. Moreover, the authors commonly\nuse big-$O$ notation when they actually should be using $\\Omega$ or\n$\\Theta$ to represent \"at least\" asymptotically.", "questions_for_authors": "[Q1] Under what conditions on $p$ and $q$ can one obtain the pushforward\nfunction $\\Phi$? If one cannot acquire this function, what are the best-case\ncomputational resources needed to employ these methods?", "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:

Meta-review: This is overall a solid paper that should be accepted to the paper. That said, all reviewers highlight that there are computational issues that are not sufficiently addressed by the 'change of variables trick' and these issues should be clearly stated and discussed in the final version.
Decision: Accept (poster)