On Average-Case Error Bounds for Kernel-Based Bayesian Quadrature
Authors: Xu Cai
Thanh Lam
Jonathan Scarlett
Venue: Accepted by TMLR ()

Abstract: In this paper, we study error bounds for Bayesian quadrature (BQ), with an emphasis on noisy settings, randomized algorithms, and average-case performance measures.  We seek to approximate the integral of functions in a Reproducing Kernel Hilbert Space (RKHS), particularly focusing on the Mat\'ern-$\nu$ and squared exponential (SE) kernels, with samples from the function potentially being corrupted by Gaussian noise.  We provide a two-step meta-algorithm that serves as a general tool for relating the average-case quadrature error with the $L^2$-function approximation error. When specialized to the Mat\'ern kernel, we recover an existing near-optimal error rate while avoiding the existing method of repeatedly sampling points.  When specialized to other settings, we obtain new average-case results for settings including the SE kernel with noise and the Mat\'ern kernel with misspecification.  Finally, we present algorithm-independent lower bounds that have greater generality and/or give distinct proofs compared to existing ones.

Published in Transactions on Machine Learning Research (06/2023)
On Avssage-Case Error Bounds for Kernel-Based Bayesian Quadrature
Xu Cai
Cair Miirounnu.cee edu
Deportment of Computer Science National University of Singapore
Chi Thanh Lam
lamchithanh1199 1997 19770mmmill..l...mmmmmmmmm com
Department Computer Science
National University of Singapore
Jonathan Scarlett
scarlett@ccopp.nusccom edu.
Department of Computer Science, Department of Mathematics, Institute of Data Science National University of Singapore
Reviewed on OpenReview: https: oppenreveee, net/ forum? id- JJrKbg33355
Abstract
In this paper, we study error bounds for Bayesian quadrature (BQ), with an emphasis on noisy settings, randomized algorithms, and average-case performance measures. We seek to approximate the integral of functions in a reproducing kernel Hilbert space (RKHS), particularly focusing on the Matérn-v and squared exponential (SE) kernels, with samples from the function potentially being corrupted by Gaussian noise. We provide two-step meta-algorithm that serves as general tool for relating the average-case quadrature error with the f2-function approximation error. When specialized to the Matérn kernel, recover existing near-optimal error rate while avoiding the existing method of repeatedly sampling points. When specialized to other settings, we obtain new averago-case results for settings including the SE kernel with noise and the Matérn kernel with misspecification ion. lower bounds via distinct proof techniques to the
Finally, we present existing ones.
1 Introduction
The inlegration of black-box functions a fundamental problem with numerous applications, with Bayesian inference being a prominent example. The method Bayesian quadrature (BQ) O'Haaga 1991 Rasmussen
& Ghahramani 2003) has particularly gained popularity, adopting Bayesian modeling techniques to model
the unknown function and reduce the required number of function evaluations. Mathematically, the goal is to approxinate the quantity
T()) - f(x)p(x)dx,
(1)
where p(x) is a known weighting function, but we only have black-box access to f(x). It is common to study this problem for worst-case functions in given class, and with noiseless function queries and deterministic Kanagawa & Hennig 2019)). However, this is also substantial motivation to understand average-case performance measures, noisy queries, and randomized
algorithms (c.g., see (Novak 1988
Kanagawa et al. 2016
algorithms (e.g., see (Novak
1988
Plaskota 1996:
Novak 2016
As an example motivating noisy sottings, when performing intograls in Bayesian inference, function evaluations themsclves may be implemented using a randomized subroutine whose variations can be modeled introducing noise terms. As rather different example, in the same way that noisy Bayesian optimization

---

Published in Transactions on Machine Learning Research (06/2023)
(BO) can be used to find maximal sensor readings in sensor network (e.g., see Krause & Ong 2011
noisy BQ methods could be used to find average (or weighted average) readings.
1.1 Overview of Contributions
In this paper, our focus is on kernel-based BQ methods, corresponding functions lying a reproducing kernel Hilbert space (RKHS). We focus on the widely-used Matérn and Squared Exponential (SE) kernels, given by (Rasmussen & Williams 2006)
xx( 21-v [(v)
22 x - x'1l
V2v 1x
x'Il 1
Ju
(2)
1
x x'12
xsE(x,X)) = exp
(3)
2/2
where Ju is the modified Bessel function, r is the gamma function, I is the length-scale, V is the Matérn smoothness parameter, and ||-I| denotes the Euclidean norm. Our main contributions are outlined as follows.
Meta-Algorithm. Our main algorithm a "meta--lgorithh" that builds on an idea originally used for Sobolev (and other) functions noiscless scttings, and appears be traced back (Bakhvalov 1959). The mela-algorithm runs any algorithm for estimating the function in 12-norm using the first hall of the samples, runs Monte Carlo (MC) method using the remaining samples to estimate the residual, and then combines the two. We give general theorem relating the L of the first step with the overall integration error, and show that this has broad implications beyond the settings for which the idea was used previously, as outlined below.
Similar ideas have also been proposed for variance reduction as control variates (Ripley 2009), and recently cast into a general functional approximation problem as control functionals (CF) (Oates al. 2017). A notable distinction between
algorithm and CF is that consider the effect of the noise, different assumptions are placed on p(x), and the specific function classes that we consider different (Matérn and SE kernels). See Section 1.2 for further discussion.
Resulting Upper Bounds. Some specific applications of our general theorem are as follows:
By using a maximun-variant algorithm from (Vakili al. 2021a] in the first part, we establish that for the Matérn-y kernel, one attains order-optimal error bound of (TT --1 + oT- 4) (with time horizon T, noise variance o², Matérn smoothness v, and dimension d). This scaling was previously derived in (Plaskota 1996), but used method based on repeatedly re-sampling a carefully-chosen set of points many times to reduce noise, which may be undesirable in practice.
We additionally provide corollaries of our main theorem providing results that appear to be new (though related results with other settings or criteria do exist), including for the Matérn kernel with misspecified smoothness, and average-case performance for the SE kernel in both the noiseless and noisy settings.
Lower Bounds. While our main contributions are those outlined above, we also provide some results regarding algorithm-independeet lower bounds. For the Matérn kernel with noise, we provide an alternative proof of the S + -T lower bound from (Plaskota) 1996) that establishes a single source of difficulty for both terms in the bound (which were previously handled separately). Moreover, establish that the Q(oT-=) ) lower bound holds much more generally in noisy settings, including for the SE kernel.
Our results are summarized in Table 1 along with various existing results that we discuss in the following subsection. Following the above outline, the paper structure is briefly described as follows: We formally introduce the problem Scction [2] state lower bounds Scction 3l provide the algorithm and its resulting upper bounds in Section and provide experimental results in Section to the appendices.
Most technical proofs are deferred
2

---

Published in Transactions on Machine Learning Research (06/2023)
Table 1: Summary of some of the most related existing results and our results. Results with an underline were partially or fully existing but are re-derived in our work (e.g., upper bounds based on distinct algorithms, or lower bounds via alternative proof techniques). Results with double underlines are new to this paper, to our knowledge.
Summarized Results
Noiseleas Worst
Noiseless Average
Noisy Average
Q(T-8-1)
2
Matérn Lower
Novak
2016
Thm. 1)
Novak
2016
Thm. 3)
Plaskota
1996
Thm. 4)
0
T
8
b)
O
O(T
oT 8)
Matérn Upper
Novak
2016
Thm. 1)
Novak
2016 Thm. 3)
Plaskota
1996 Thm. 4)
Q(T
Cr3)
SE. Lower
Kuo et. al. 2017 Thm. 1.1)
N/A
9(oT-3)
(Linest algs. only)
O(e-Cra)
SE Upper
Kanagawa & Hemnig 2019
o ( b-b
0 T +)
Cor, 1.4)
1.2 Related Work
Numerical Integration and Bayesian Quadrature. Extensive early work on numerical integration
appeared in the literature on information based complexity, e.g., see (akhvvllov] 1959 Novak 1988
Traub
2003 Novak & Wožniakowski] (008)) and the references therein. Function classes considered included Sobolev, Hölder, and others, with the Sobolev class being particularly related to the Matérn RKHS (sce
Appendix A). More explicit. use of kernel-based methods appeared in (Narcowich & Ward 2002 Wendland
2004 Wendland & Rieger] 20055 Rieger & Zwicknagl, (2010), and Bayesian quadrature from probabilistic
perspective (Rasmussen & Williams 2006 Kanagawa et al. 2016 2018
Wynne et al. 2021 has gained
particular popularity due to its role in statistical machine learning.
Upper Bounds Under the Matérn Kernel. For d-dimensional Matérn-v RKHS functions with query budget T, early literature such as (Bakhvalov] 1959 Novak 1988) proved that in the noiseless setting, the 2), whereas by considering the average-case error of randomized algorithm. this can be reduced to O(T-5- 5-1). An extensive survey of the noiscless setting can be found in (Novak) 1988), where it is also noted that basic MC sampling attains O(T-²) error; this observation also extends immediately to the noisy setting with o = O(1) (Plaskota 1996).
best possible worst-case (deterministic) error is O(T
As we hinted above, the L bounds for the Matérn kernel with noisc in (Wynne et al. 2021 will be particularly useful for our purposes, including in misspecified settings. Some implications of these results for noisy BQ were also noted in (Wynne et al. 2021), but they resulted in worse scaling (namely, O(T z with s v + 4) than that shown in Table In fairness, however, several of the results in (Wynne et al, 2021) can also be applied to settings with non-stochastic noise, whereas we only handle i.i.d. random noise.
Upper Bounds Under the SE Kernel. For functions in the SE RKHS, we observe in Table 1 that there are non-minor gaps between lower and upper bounds. For a different setting in Rd with p(x) in being a Gaussian distribution, (Kuo et al. 2017) gives a worst-case lower bound for linear algorithms, by decomposing SE functions using an orthogonal basis the 122 space. This lower bound has been further improved in KKrvvnen et al., 2021) using another orthogonal basis the native RKHS space, but it only applics to (linear) Gauss-Hermitte rules. Overall, research on the SE kernel appears have mainly focused on noiseless worst-case upper bounds, and exponential convergence has been developed via scattered data approximation Wendland 2004) and sampling inequalities for infinite smooth functions (Rieger & Zwicknagl 2010
A notable work related to ours is the adaptive BQ (ABQ) algorithm from (Kanagawa & Hennig 2019). While their algorithm is adaptive, we see in Table [ that our non-adaptive algorithm yields a slight improvement over theirs, alboit with a weaker averago-case guarantce. Similar obscrvations also apply to Matérn functions, where the ABQ algorithm only achieves O a) scaling (Kanagawa & Hennig 2019 short of the optimal scaling.
Cor. C.4), which falls
3

---

Published in Transactions on Machine Learning Research (06/2023)
Differences from control functionals. As mentioned earlier, while the control functionals (CF) literature (e.g., Oates et al. (0077] and our work share a similar approach in exploiting a tractable surrogate function with a bias-correction term in estimating the intractable integral, there are several differences between them:
Our main focus is on the noisy setting, whereas CF algorithms typically designed for noiso-free scenarios.
CF algorithms are typically developed under the restriction that p(x) is intractable, so that point- wise access to the gradient information of logp(x) is necessary, for instance, by minimizing kernel Stein discrepancy using Stein operators (South et al. 2022). However, our setting assumes access to p(x) and explores regret bounds accordingly.
The focus of the CF literature has been accelerating the 0 ((}}) convergence rate of MC or Quasi-MC (Oates & Girolami 2016 Portier & Segers 2019) under suitable assumptions, but to our knowledge less attention has been paid on attaining optimal or near-optimal convergence rates with matching lower bounds. We focus on this goal with an emphasis on the Matórn and SE kernels.
To highlight a representative example of the last dot point, we note that an O(T-7/6) rate is established in
(Oales et al. 2017 Thm. 3) under suitable assumplions on a class of gradient-based kernels, but this cannot encompass our noiseless results for the Matérn kernel in which, depending on the smoothness, it may be impossible to attain a rate of O(T-)/)), or it may be possible to attain a much faster rate.
Equivalence between Random Features and Bayesian Quadrature. An interesting equivalence between worst-case kernel-based quadrature and random feature based function approximation in I norm has once established in [achh 2017). While the main focus in (Bach) 2017) is the worst-case noiseless setting, it is noticed in Sec. 3.1 therein that their methods have a certain tolerance to noise. Noisy error bounds can be obtained via this fact, but they are higher than ours, with the first of the two terms (e.g.. soe 9 below) typically being a factor VT smaller in our work, and being order-optimal in several cases of interest. We crucially rely randomization to achieve this, whereas the proposed algorithm in (Bach 2017 is deterministic. We note that for the SE kernel, one of our results (Corollary 4 uses a result from Bach 2017 as an intermediate step.
Relationship with Bayesian Optimization. The extensive literature on Bayesian optimization (BO) is also related to our work in the sense of iterating acquisition functions. However, BO turns out to be a strictly harder problem. As surveyed in detail in the noiseless setting in (Novak 1988 Sec. 1.3), noiseless BO is closely related to estimating the function in 7 norm, and noiseless BQ is closely related to estimating it in L norm, with the former being strictly harder. Viewed differently, key difficulty in BO is identifying a single short and narrow "bump" (Bull 2011 Scarlett et al. 2017), whereas in BQ the same bump contributes negligible amount to the integral. We elaborate on this connection in Appendix F showing that BO-type techniques yield suboptimal BQ bounds for the Matérn kernel. On the other hand, in Section 4 we will see that this approach gives a fairly good result for the SE kernel in the noiseless (but not noisy) setting.
Other Variants of BQ Finally, various works on BQ have explored more sophisticated techniques and variations such adaptive sampling (Kanagawa & Hennig. 2019), active area search (Ma et al., 2014), and settings with multiple related functions [Gessner et al. 2019 (among many others). We are not aware of any (beyond those outlined above) that are directly related to our study of theoretical error bounds in relatively standard settings.
2 Problem Setup
Let f D - R be real-valued function on the compact domain D = [0, 1jd By shifting and scaling, our results readily extend to arbitrary rectangular domains. In addition, our upper bounds easily extend to general compact domains.
We consider the class (B) of functions whose RKHS norm ||x is upper bounded by some constant
B > 0. We focus in particular the Matérn-y kernel (see ② ) whose function class Hm is norm-
4

---

Published in Transactions on Machine Learning Research (06/2023)
equivalent to the Sobolev class (see Appendix A for details). We also derive results for the SE kernel (see 3 where the function class is denoted as HSE B).
Let p(x) be a known and bounded density function, i.e., p(x) E [0, (0ax] for some Pmax > 0, and p(x))xx = 1. We define P(Pmax) to be the set of all functions satisfying these conditions. Our goal is to estimate the integral of an RKHS function f D -> weighted by p(x), defined in (1) Before forming an approximation of this integral, the algorithm takes T observations: At time step 6, select X E D, and observe yt f(xx) Et: where 6 (N0 o is Gaussian noise. The final approximate integral is denoted by î.
As is commonly done for RKHS functions, we will use Bayesian methods based on GP prior GP(O, k) and a Gaussian noise model. After observing t noisy samples, the posterior distribution is also a GP with the following posterior mean and variance:
[e(x) = kx(x)) + ApI1) ye:
(4)
(x ((xxx))k((x) (K, + daIe) (x(),,
(5)
K. = [k(x,, x xc) [k(x),xe)) e.r is the kernel matrix, I, is the identity matrix of dimension and Ay and A are hyperparameters. Despite the Bayesian GP posterior needs A = Ay = o², as commonly accepted in recent RKHS-basee BO literature, it is often useful to consider "fictitious" GP posterior that the parameter X or A may differ from o². In this case, as revealed in Kanagawa et al. 2018), the GP posterior mean is shown to be equivalent to the approximating function in kernel ridge regression (e.g., see (Caponnetto parameter. On the other hand, when A = O, the posterior variance is also known as the power function (Santin & Haasdonk 2017), and will be frequently used in this paper. Nonetheless, we will continue to use the terms GP posterior mean and variance to refer to (4) and (5), respectively. Unless explicitly stated otherwise, when mention X, we consider common choice A =
where = [y),...., ".......,..'" Mr]",, k,(x) = [k(x;, x)] [x(x;,))]
De Vito 2007)), where A treated as the regularization
We consider both adaptive algorithms (i.e., the algorithm observes y1,.... Ut-1 before choosing xt) and non- adaptive algorithms (i.c., all XI...., Xy are chosen in advance). In fact, we will prove our lower bound for adaptive algorithms and upper bounds for non-adaptive algorithms, thus establishing the stronger type of result in both cases.
Using the shorthands H4 (B) (as well as HM HsE when the kernel k is specialized as Matérn SE) and P = P(pmaxx) the settings in Table 1 are sunmarized as follows, the last of which is the one we primarily focus on:
Noiseless Worst-Case Error: E (T) = SUPpEP, |T 7/.
Noiseless Average-Case Error: Elg((T) (T) = sup suppEP,fEH. E[[ Î1], where the expectation is over the randomized algorithm.
Ek (T,o) = suppep, [ Î1]. where the expectation is over the randomized algorithm and the noisc.
Noisy Average- Case Error:
We have omitted a notion of worst-case error for the noisy setting, since there are subtle issues in posing such setting in mcaningful manner (e.g., see (Plaskota, 1996)). For instance, deterministic given the noisy observations, it can still obtain randomness by taking digits after the 1000th decimal point (say) of the observed values yt- That is, the randomness from the noise alone could still give the same effect as using a randomized algorithm.
if the algorithm is
3 Lower Bounds
Our lower bounds for the noisy setting are stated as follows; note that we allow a to vary with T.
Thcorem 1. (Averago-Case Noisy Lower Bounds) Consider our problem setup with constant parameters (B,vd, I), noise variance o², and query budget T. Then, for any algorithm (possibly adaptive and/or mndomized) for estimaling I, we have the following lower bounds on the average-case error:
5

---

Published in Transactions on Machine Learning Research (06/2023)
Algorithm 1 Two-batch integral estimation meta-algorithn
1: Input: Function f, domain D, time horizon T. function estimation algorithm ESTIMATEF UNC
2: Use ESTIMATEF UNC with 2 samples to produce function estimate
for t _ + I,..., do
4: Sample X2 ~ p(x) independently
5: Observe V/+ = fxx) +
6: end for
7: Compute the approximate integral í, - Jo (())((xx
Compute the residual R 7 E1-//2+1 f(xi))
Output i = i1
1. For Matérn kernel with v + d 2 1, we have
Ea (T, ( oT o)-
(6)
2. For any kernel such that maxx k(x, x) and fo.uj k(x, xx)xx # 0 for some x² (e.g., the SE Matérn kernel with x² = =), we have:
g = = O)
Moreover, these lower bounds hold even under the fixed weight function p(x)
3.1 Proof of Theorem
The first term in (6) follows directly from noiseless average-case error bounds (see Appendix B The second term can also be established by considering constant-valued functions (Plaskota| 1996), or as a special case of (7) (which is proved below). On the other hand, is also of interest to establish a single hard subset of functions that yields both terms in (6) in unified manner, thus establishing a single source of difficulty for both terms. We provide such approach in Appendix D considering functions composed of several small "bumps". The idea is that with too few samples the algorithm cannot reliably determine which bumps are positive and which are negative, and if too many of these are uncertain, then a certain level of error is unavoidable.
It remains to prove 7 To do so, we consider the function f(x) = A(x, x)), which is bounded and has a non-zero integral by assumption. Moreover, since k(, k(,x*) x²) is trivially the RKHS class, we have that the scaled function fe(x) cx(x) has RKHS norm at most B for all 0 below a suitably-chosen threshold.
We consider the sub-class of functions fo with RKHS norm most B, and show that f(oT- à) queries are needed even for this sub-class. To see this, note that sampling any point x corresponds to observing c, scaled by fxx), with (0, noise. Without loss of optimality, we can assume that the algorithm only samples at the psint(s) where [x(] is largest (so that e is maximally scaled while the noise level is fixed), and by assumption we have [ffxx) < 2o. Then, the noisy integration problem on this sub-class reduces to the noisy estimation of e with Gaussian noise, which is well known incur ( error (even when e is known be below an arbitrarily small constant threshold), e.g., see (Plaskota 1996).
4 Upper Bounds
In this section, we introduce our main mela-algorithm and derive upper bounds on the average-case error. We follow the high-level idea of combining function estimation methods with Monte Carlo estimation on the residual, previously proposed for the noiscless sctting (c.g., (Bakhvalov) 1959 Soc. 2) and (Novak] 2016 Thm. 3)), but with different details account for the noise. The meta-algorithm is shown Algorithm and is described follows. The samples are performed in two batches, but still in non-adaptive manner

---

Published in Transactions on Machine Learning Research (06/2023)
Algorithm 2 Maximum variance sampling; can be used as ESTIMATEPUNC in Algorithm
Input: Function J, domain D, GP prior GP(O, k), GP noise parameters Ap and Aa time horizon 2: for t= -..... do
3: Select x, = arg maxxeD E 0x-1(x)
Receive yt = f(x,)) 4
Update or using X1,..,,,;
end for
Update (ryz(xx) using X1, ,T,///, 01..... Mr/2
8: Output MTy2
(i.c., the second batch can be chosen without knowing the first batch) The first batch uses T/2 samples to construct approximation f of f, which forms initial estimate 4 So (x))xx The subroutine ASTIMATEFUNC for doing so is kept general for now, but by default, we will use maximum variance sampling as shown in Algorithm 2 with the estimated f being given by the GP posterior mean Pryz- Note that these T/2 samples are non- adaptive because the posterior variance of a GP does not depend on any observations. The remaining T/2 samples estimates the residual R between the difference of the true integral and this value I:
R R=I- - -
Having fixed and hence f, we can refine our overall estimate of I by estimating the residual R. To do so; we use the last T/2 samples to construct a Monte Carlo estimator f of R, i.e.,
R= 2
T
Σ (yt - f(xx)),
T=T/2+1
in which each x, is sampled from p(x). Hence, the approximated integral is simply i= 1 + ft.
In the following theorem, we establish the upper bound of the average-case integration error of Algorithm ! in terms of the noisy I.2 distance between f and f, with extra noise term. This theorem general and can be applied to variety of kernel choices, noise settings, misspecified settings, etc., we will see below. In the following, all omitted proofs are deferred to Appendix E
Theorem 2. (Quadrature-L² Relationship) For any ftaed J Hr and E P(Pmaxx, consider running Algorithu 1 with T observations being corrupted by Gaussian independent noises N0, and obtain the estimates f and i. Then, the average error satisfies
B - Í1 <
أ V20T
(8)
where the expectation is averaged over the randomness of the points queried and the noise.
For the Matérn kernel, building on the results of (Wynne al. 2021 Thmn 4), we obtain the following two corollaries.
Corollary 1. (Average- Caserge-Case Noisy Matérn Upper Bound) For HM, consider the setup of Theorem where the subroutine ESTIMATEFUNG for producing I is mazimun- sampling (Algorithm ? with porameters, >a = and À = (TT Then, the avernge-case integration error is upper bounded by
EM (T,
'Our main theorem remains true when the first batch consists of adaptive sampling, but all of our corollaries will be based non-adaptive sampling.
2e focus maximum variance sampling for concreteness, but analysis and results unchanged when, at time arbitrary point satisfying (x) 2 7 1 is found for some C (0,1) (c.g., 7 4). This requirement is potentially much casier to attain in practice instead of insisting the global maximum.
We consider of points. Alternatively, any other quasi-uniform set could be pre-specified, such uniform grid.
0 order exploit an existing result that the resulting maximum maximum-varianes variance rule gives quasi-uniform

---

Published in Transactions on Machine Learning Research (06/2023)
As discussed in Section 1.2 Corollary 1 matches the lower bound (Theorem up to constant factors, and shows that the overall difficulty of noisy BQ is as roughly hard either noiseless BQ or univariate Gaussian mean estimation with variance o² (whichever is harder). Since the algorithm used is non-adaptive but the lower bound applics even to adaptive algorithms, we conclude that there no adaptivity gap in terms of scaling laws in this case.
Similarly to prior works such (Kanagawa et al. 2016 Teckentrup 2020 Wynne et al. 2021 we now
consider the setting in which the smoothness hyperparameter v is unknown, and Algorithm 2 is accordingly modified follows:
Instead of assuming that is given (via the prior ko), we assume that the algorithm is given a sequence of estimates {d=) These kept generic, but could correspond to estimates that are iteratively updated more data is collected.
When forming the GP posterior at time t, the parameter in is used.
Following [Wynne et al. 2021 p given the time horizon T, we define -ー =nffsr de and +t sup,et De and we assume that the number of distinct values in {Dr] is upper bounded by fixed constant.
Corollary 2. Awcrage-Cass Noisy Matérn Upper Bound Misspocified Sctting) For f E Hus, consider the setup of Theorem & where the subroutine ESTIMATEFUNCE for producing f is the above-described modification of marimumm-urrinne sampling with parameters {Dr)t=i satisfying ++ > v, and with Xa = and X
(TT- Then, the average-case integration error is upper bounded by
eMg (T, -T, - minte md
+ 3)
(10)
As a special case of this result, when -ー = ++ (and both are greater than in accordance with the assumption y+ > v); the scaling precisely reduces to that of Corollary [ Thus, interestingly over-estimating the smoothness paramcter is not harmful in terms of scaling laws; similar obscrvations were made in prior works such as [Kanagawa et al.) 2016) (Remark 2 therein). More broadly, Theorem 2 can be used convert other T,2 guarantees from (Wynne et al. 2021 (or from other works) to BQ guarantees, but for brevity only formally state Corollary
We now turn to the SE kernel without noise, by utilizing results from kernel-based interpolation on cubes [Wendland] 2004 Rieger & Zwicknag] 2010), and obtain the following corollary.
Corollary 3. (Average-Case Noiseless SE Upper Bound) For f E HSE, consider the setup of Theorem where ESTIMATEFUNC chosen to take samples on uniformby-spucedd grid and return 1T7/2 with X 0. Then, the average-case error is upper bounded by
SSE (T) O O O +
(11)
for some consland C, >> 0.
The same result can also be obtained with exponentially small but non-zero X (see Appendix E.3 but we provide the formal statement with X 0 for simplicity.
We specifically build on (Bach] 2017), which analyzes weighted random Fourier features from function approximation point of view. The approximation requires sampling from a leverage function distribution based on an infiiitedimmensional integral operator. With the help of their analysis, we obtain the following corollary.
For SE kernel with noise, we follow the use of random features (Rahimi & Recht 2007
Corollary 4. (Average-Case Noisy SE Upper Bound) For HSE, under the setup of Theorem 2 there erists a random fealure sampling algorithm that produces an approrimation I Algoritlmm 1 such that the average-cas case error is
ave (T,O) C +
4This can be generalized so that these definitions restrict 2 N some finite but simply take N reduce notation.
8

---

Published in Transactions on Machine Learning Research (06/2023)
for some constant C, 0.
In Section D highlighted that BQ is related to RKHS-based optimization problems. We briefly note that the techniques from such works also be used to obtain BQ bounds via Algorithm 2 (using it for T obscrvations rather than only the first half) or sampling on uniform grid, but the resulting scaling is lypically highly suboptimal. Specifically, in Appendix F
we show that this approach leads to the following:
cM STgg(To) = o 0 z=za (log T) Totza + oT wia (log T)
(12)
eM (T) = 0 0
(13)
4+2 (T, σ) 0 (log 7) d-1 2 T o(log T) 8 r $). eSE
avg
(14)
(T) 0(e TT²)
(15)
These four results correspond to the Matérn noisy poisy/noisee noiseless and SE noisy /noiseless settings respectively, Algorithm 2 can be used directly to obtain the first three results, whereas (15) is based on sampling
uniform grid. While these bounds above are weaker than our earlier results, the noiseless SE case 115 gives
almost as fast decay to zero as
5 Experiments
In this section, we conduct simulation studic to compare our two-batch algorithm with its component parts, maximum variance sampling (MVS) and Monte Carlo sampling (MC). These experiments serve to verify that Algorithm can be effective in practice, but we will also discuss some potential gaps between the theory and practice. In particular, our goal is not to establish state-of-the-art practical performance.
We refer to Algorithm as MVS-MAT or MVS-SE when the kernel is Matérn or SE, and similarly for Algorithm with MVS-MC-MAT or MVS-MC-SE. To attain better understanding of how MVS-MC performs with respect to time, modify it to alternate between MVS samples and MC samples, instead of doing all of one followed by all of the other. Mathematically, this does not change the behavior at the final time step.
5.1 Setup
GP model. We adopt the common choice v = 332 for Matérn-w kernel, and the GP mcan noise hyperparameter Ay is fixed as 10-5 for both kernels [ The lengthscale is left as a free parameter, as is an additional scale parameter that we introduce (multiplying (2) to permit functions with varying ranges. Except where stated otherwise, these two parameters are learned by maximizing the data log-likelihoo [Rasmussen & Williams 2006) using the built-in SciPy optimizer based on L-BFGS-B, which is also used for finding the maximum variance point. We seek to solve the BQ problem with constant weight. function, i.., p(x) 1.
Evaluation. We compare MVS-MC against its two components, MVS and Monte Carlo. For MVS and MVS-MC, we perform 100 trials, with each trial using distinct set of 3 random initial points (but these are common to both mothods). For MC, since every round is already randomized, we simply run 100 trials without initial points. For all functions, we consider time horizon of = 250, and evaluate the performance using the mean absolute error, with the ground truth value (and also î Line of Algorithm I being determined by trapezoidal rule with 105 uniformly-spaced grid points (without noise). Errors plotted with logarithmic scale for the purpose of more casily distinguishing values close to Error bars in plots indicate +0.55 standard deviation with respect to the 100 trials.
Thee code can be found at htttpa://githan 60ur theory choices such A 0(7 2) with unspecified implied constants, which makes it difficult choose the parameter in way that exactly matches the theory. However, since the theory dictates choosing An be small, we set be small here accordingly.
httm//aatree//eeeeeeee wwweelieeeedddiiiiaaaae Bayesian--usdrrauurr
9

---

Published in Transactions on Machine Learning Research (06/2023)
Synthetic3d
Synthelk3d MAS-NC-M0
SE Syntheticdd
Synthelic3d
0.31
1976 MCSE
His Hol
His 58
0.32
0.33
Emor
0.08
0.66
201
0.68
30-1
200
ala a's ah d' de NVS MC Splits
55
100 Querics
260
0.2
02
e'a
Oucrics
MVS MC Spits
Figure 1: Results for synthotic functions.
Addleyld 0.5
Npinezd
Energy TS o- 0.1
Morta-Cart
10a
Maraa Care MoS--ATT NOS-MGHHUT
Manta-Cara HSS MAT
331
Mento-Cats
HiS
MY5-4E MYS-MCSE
His
Exrol
H4%
Enrrrr
10
2/0
350
116
100
2.0
50
200
100
250
Queries
Querfes
Queries
Quertes
Ackloyld-a
Keanned-
Alpine7d
Energy TS -i 0.1
6.12
MYS-MC MAT
MiS 4C44
Mes HC-HHT
Nos-
625
MYS-CCWA
Mis MC-42
3.60
6:2
Entrr
0.00
6:00
24
0.04
3.82
6:0
32
25
8.0
87
2/ MVS-MC Spils
MVS-MC Spits
NVS-MC Spiiss
MVS-MC Splits
Figure2 2: Results for benchmark and real data functions.
Varying the split size. In our theoretical analysis, we let MVS-MC use half of the samples for each component (MVS and MC). However, in practice, it may be beneficial to allow different splits. In our experiments, additionally explore the effect of the split fraction, ie., the fraction allocated to the MVS part, and plot the final error with respect to that fraction. A split of 0 corresponds to MC, and a split of corresponds to MVS. For these plots we use a linear y-axis scale, since the values are mostly further from zero and casicr to compare on such a scale.
5.2 Results
We provide selection of our results here, and give additional results in Appendix G
5.2.1
Synthetic Kernel-Based Functions
We first simulate on a set of synthetic functions following the method of (Janz et al. 2020 where each function constructed by sampling m = 30d points, x1 ..... &m, uniformly on [0, 1J,, and aj... am uniformly [-1, 1]. The function is then defined as f(x) > 0,k(,, x). The length-scale and (in the case of the Matérn kernel) are set to be fixed 0.2 and 3/2 respectively.
We let MVS and MVS-MC know the kernel hyperparameters exactly (i.c., they match the ones used to produce the functions). Some results for d = 3 and a E (0.1, are shown in Figure [ and further (d, pairs are shown in Appendix G We observe that MC performs well at the higher noise level, but performs poorly at the lower noise level, which aligns with the theory. At the final time step, MVS in fact performs well both cases, though it can perform poorly in the low-query regime (c.g.. Matérn 4D). While MVS-MC is not universally best, it is generally able to capture the benefits of both MVS and MC in this experiment.
10

---

Published in Transactions on Machine Learning Research (06/2023)
5.2.2 Benchmark Functions and Sensor Measurement Function
We consider a varicty of well-known benchmark functions, namely, Ackley, Alpine, Gramacy-Lee, Griewank and Keane; see (Bingham 2013) for the descriptions. We also consider an experiment for sensor measurement data, which is described in Appendix G
A small selection of the results are shown in Figure and the rest
are shown in Appendix G
These results again indicate that MC is better suited to higher noise levels, but they provide somewhat less clear picture for MVS and the ideal MVS-MC split. In general, the error as a function of the split fraction can increase, decrease, or exhibit "U-shaped" behavior, though we again found the choice of 0.5 to usually a good one (even if not always optimal).
Notably, in both the syuthetic and benchmark experiments, the error can drop suddenly at 1.0, this is where the MC component theoretical benefits of MVS-MC over MVS may not always be observed in practice (e.g., because the theory hides important constant factors, only captures the behavior for very large T). It may be of interest to seek refined theory addressing these findings in future work.
longer used and MVS-MC simply becomes MVS. This indicates that the
6 Conclusion
We have developed framework for relating average-case quadrature error with T.2-function approximation error, allowing us to derive a number of both existing and new results. In addition, we explored algorithm- independent lower bounds with greater gencrality and/or distinct proofs compared to cxisting oncs. In future work, it may be interest to address some the remaining gaps, such as those between the upper and lower bounds for the SE kernel (see Table 1
Acknowledgments
This work was supported by the Singapore National Research Foundation (NRF) under grant number R- 252-000-A74-281. C. T. Tam was supported by the Singapore-MIT Alliance for Research and Technology (SMART) PhD Fellowship.
References
Francis Bach. On the equivalence between kernel quadrature rules and random feature expansions. J. Mach. Learn. Research (JMLR), 18(1):714-75)) 2017.
Nikolai Sergeevich Bakhvalov. On the approxinate calculation of multiple integrals. Vestmik MGU, Ser. Math. Mech. Astron. Phys. Chem, 4:3- 18, 1959.
Derek Bingham. Virtual library of simulation experiments: Test functions and datasets. https://ww/ sfu,
ca--ssuraano/ wwerr/nnnonnnnnnieee Lande..xhmml 2013.
Adam D Bull. Convergence rates of efficient global optimization algorithms. J. Mach. Learn. Research
(JMLR), 12(10), 2011.
Xu Cai and Jonathan Scarlett. On lower bounds for standard and robust Gaussian process bandit optimization. In Inl. Conf. Mach. Learn. (ICML). PMLR, 2021.
Andrea Caponnetto and Ernesto De Vito. Optimal rates for the regularized least-squares algorithm. Foundations of Compulational Mathemalics, 7:331-366 2007.
Thomas M Cover and Joy A Thomas. Elements of Information Theory. John Wiley & Sons, Inc., 2006.
William Feller. An introduction lo probabilily theory and ils applications. Wiley, 1957.
Alexandra Gessner, Javier Gonzalez, and Maren Mahsereci. Active multi-information source Bayesian quadrature. In Tnl. Conf. Art. Intel. Stats. (AISTATS)) 2019.
11

---

Published in Transactions on Machine Learning Research (06/2023)
Kerstin Hesse. A lower bound for the worst-case cubature error spheres of arbitrary dimension. Numerische Mathemalikk, (03(3):413-433) 2006.
David Janz, David R. Burt, and Javier González. Bandit optimisation of functions the Matérn kernel RKHS. In Int. Conf. Art. Intel. Stats. (AISTATS), 2020.
Motonobu Kanagawa and Philipp Hennig, Convergence guarantees for adaptive Bayesian quadrature methods. In Conf. Neur. Inf. Proc. Sys. (NeurIPS), 2019.
Motonobu Kanagawa, Bharath Sriperumbudur, and Kenji Fukumizu. Convergence guarantees for kernel- based quadrature rules misspecified settings. In Conf. Neur. Inf. Proc. Sys. (NeurIPS), 2016.
Motonobu Kanagawa, Philipp Hennig, Dino Sejdinovic, and Bharath K Sriperumbudur Gaussian processes and kernel methods: A review on connections and equivalences. https://arxvvoor//aaaa 0001110078299 2018.
Toni Karvonen, Chris Oates, and Mark Girolami. Integration in reproducing kernel Hilbert spaces of Gaussian kernels. Mathematics of Computation, (0(331):200-22233 2233, 2021.
Andreas Krause and Cheng Soon Ong. Contextual Gaussian process bandit optimization. In Conf. Newr. Inf. Proc. Sys. (NeurIPS), 2011.
Frances Kno, Ian Sloan, and Henryk Wožniakowsli. Multivariate integration for analytic functions with Gaussian kernels. Mathemalics of Compulation, 86(304):829 853, 2017.
Madison Lee, Shubhanshu Shekhar, and Tara Javidi. Multi-scale zero-order optimization of smooth functions in an RKHS. In IEEE Inf. Symp. Inf. Theory (ISIT), 2022.
Yifei Ma, Roman Garnett, and Jeff Schneider. Active area search via Bayesian quadrature. In Int. Conf. Art. Intel. Stats. (AISTATS), 2014.
Francis J Narcowich and Joseph D Ward. Scattered data interpolation on spheres: Error estimates and locally supported basis functions. SIAM Journal on Mathematical Analysis, 33(6):1333311 1410, 2002.
Erich Novak. Deterministic and stochastic error bounds in numerical analysis, volume 1349. Springer, 1988.
Erich Novak. Some results on the complexity of numerical integration. In Monte Carlo and Quasi- Monte Carlo Methods, pp. 161-183. Springer, 2016.
Frich Novak and Henryk Wožniakowski. Tractability of malliiaariate problems: Standard information for functionals, volume European Mathematical Society, 2008.
Chris Oates and Mark Girolami. Control functionals for quasi-monte carlo integration. In Artificial Intelligence and Statistics, pp. 56 65. PMLR, 2016.
Chris J Oates, Mark Girolami, and Nicolas Chopin. Control functionals for monte carlo integration. Journal of the Royal Statistical Society. Series B (Statistical Methodology), pp. 695-718, 2017.
Anthony O'Hagan. Bayes hermite quadrature. J. Stat. Plan. Inference, 29(3)2245 260, 1991.
Bayes-hermite
Leszek Plaskota. Worst case complexity of problems with random information noise. Journal of Complerity, 12(4):416-439 1996.
François Portier and Johan Segers. Monte carlo integration with a growing number of control variates. Journal of Applied Probability, 56(4):1168-1186, 2019.
Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. In Conf. Neur. Inf. Proc. Sys. (NeurIPS), 2007.
Carl Edward Rasmussen and Zoubin Ghahramani. Bayesian Monte Carlo. Conf. Neur. Inf. Proc. Sys. (NeurIPS), pp. 505-512. 2003.
12

---

Published in Transactions on Machine Learning Research (06/2023)
Carl Edward Rasmussen and Christopher K. I. Williams. Goussian processes for machine learning. Adaptive computation and machine learning. MIT Press, 2006.
Christian Rieger and Barbara. Zwiiknngl. Sampling inequalities for infinitely smooth functions, with applications to interpolation and machine learning. Advances in Computational Mathematics, 32(1):103 129, 2010.
Brian D Ripley. Stochastic simulation. John Wiley & Sons, 2009.
Klaus Ritter. Average-case analysis of numerical problems. Number 1733 in Lecture Notes Mathematics. Springer Science Business Modia, 2000.
Gabricle Santin and Bernard Haasdonk. Convergence rate of the data-independent p-groedy algorithm kernel-based approximation. Dolomites Research Notes on Approzimation, 2017.
Gabriele Santin and Robert Schaback. Approximation of eigenfunctions in kernel-based spaces. Advances in Computational Mathematics, 42(4):973-99) 1909,,93--1 2016.
Jonathan Scarlett and Volkan Cevher. An introductory guide to Fano's incquality with applications in statistical estimation. htlps://arxiv.org// wabpiiiwwww//a//////a/aa/b/ss/ httpg///b//1110..0000 000555 2019.
Jonathan Scarlett, Tija. Bogunovic, and Volkan Cevher. Lower bounds on regret for noisy Gaussian process bandit optimization. In Conf. Lean. Theory (COLT), 2017.
Leah F South, Toni Karvonen, Chris Nemeth, Mark Girolami, and Chris J Oates. Semi-exact control functionals from sard's method. Biomelrika, 109(2):351-2) 16657, 2022.
Niranjan Srinivas, Andreas Krause, Sham Kakade, and Matthias Seeger. Gaussian process optimization in the bandit setting: No regret and experimental design. In Int. Conf. Mach. Learn. (ICML), 2010.
Yitong Sun, Anna Gilbert, and Ambuj Tewari. But how does work in theory? linear SVM with random features. In Conf. Neur. Inf. Proc. Sys. (NeurIPS), 2018.
Aretha L. Teckentrup. Convergence of Gaussian process regression with estimated hyper-parameted and applications in Bayesian inverse problems. STAM/ASA Journal Uncertainty Quantification, 8(4):1310- 1337, 2020. doi:
Joseph F. Traub, Information-Batt Based Complezity, pp. 850- 854. John Wiley and Sons Ltd., GBR, 2003. ISBN 0470864125.
Sattar Vakili, Nacime Bouziani, Sepehr Jalali, Alberto Bernacchia, and Da-Shan Shiu. Optimal order simple regret for Gaussian process bandits. In Conf. Neur. Inf. Proc. Sys. (NeurIPS), 2021a.
Saltar Vakili, Kia Khezeli, and Victor Picheny. On information gain and regret bounds in gaussian process bandits. In Int. Conf. Art. Intel. Stats. (AISTATS), pp. 82 90. PMLR, 2021b.
Holger Wendland. Scattered data approzimation, volume 17. Cambridge university press, 2004.
Holger Wendland and Christian Rieger. Approximate interpolation with applications to selecting smoothing parameters. Numerische Mathematik, (4144)72277788 2005.
George Wynne, Srançois-Xavies Briol, and Mark Girolami. Convergence guarantees for Gaussian process means with misspecified likelihoods and smoothness. J. Mach. Learn. Research (JMLR), 22(123):1 40, 2021.
Wenjie Xu, Yuning Jiang Emilio T Maddalena, and Colin N Jones. Lower bounds the worst-case complexity of efficient global optimization. Xrxiv preprint ar arXiv:2209.000 2022.
13

---

Published in Transactions on Machine Learning Research (06/2023)
A Connections Between Matérn RKHS and Sobolev Class
We briefty summarize some definitions and propertics of the Sobolev function class, as stated and used in the work (Wynne et al. 2021 that we will later build on.
For integer-valued the Sobolev function class of order on domain D is defined as
WD(D)= E
Σ
10 fll < oo).
(16)
and the associated norm is given by
-0 EΣ 100h 11 1/2
1/2
MDWWIDD) =
(17)
la|se
where each a a d-dimensional multi-index with |a] = =i= dg, and 0 f(x) is the weak derivative order d a. For a fractional (i.c., non-integer) order s E with 2 (i.c., v > 0), the Sobolev function class of order s is defined on Rd as
WS(R ( LR())
/ gd (1+ <o)
(18)
where f(E) is the Fourier transform of f(x). The associated norm is given by
IVllwz(R= = (1-x
1/2
For functions defined over D S Rd, the fractional Sobolev norm is taken as the infimum over all functions on Rdd that match the original function when restricted to D:
lfllg = inf [I'lllweee)) f E (R") and f'(x) f(x) VX D).
(19)
We note that are slightly abusing notation here, since in the integer case 19) only matches 117 up to constant factor (rather than exactly), but this will not impact any of our results in which the constant factors are kept implicit.
In the following, we will use W as a shorthand for (D) Here we state a known result that characterizes the equivalence between Sobolev space and RKHS of Matérn kernels.
Lemma 1. (Sobolev Space& & RKHS of Matérn kernels (Teckentrup) 2020 Prop. 3.3)) Let HM be the RKHS of the Matérn-v kernel on D = [o, 1]d (with any fized positive length scale), and let s = v + d/2. Then, the the Sobolew space Wz. That is, Hm = Wi, and there exist constants C1 C2 > 0 such that for any f E Hm (or equivalently, f E W2) we have
Malérn RKHS is norm-equivalent
cillfllwz I Il I callfllwg.
(20)
B Existing Results for the Noiseless Setting
The lower bound for the noiseless setting (i.c., o² = 0) presented below restates the result of (Novak
d 2016, Thm. 1 & 3), with a slight generalization to fractional Sobolev spaces where & = v + 2 may be non- integer. Here and subsequently, the kernel parameters v, I, dimension d, and RKHS norm B are all treated constants, and consider the limit. T x0.
Thcorem (Noiscless Lower Bound) Consider the noiscless problem setup with constant parameters (B, v. d) satisfying + d 1, and time horizon T. Then, we have the following lower bounds for the Matérn kernel in the noiseless selting:
14

---

Published in Transactions on Machine Learning Research (06/2023)
1. For any algorithm (possibly adaptive), the worst-case enror satisfies
e wat S2(T $-)).
2. For any algorithm (possibly adaptive), the average-case error satisfies
Cavg(() -2(1
Moreover, these lower bounds hold even under the fixed weight function p(x)
While Theorem is already known for Sobolev spaces with integer-valued 8, fractional bounds are less commonly stated in early literature, though interpolation techniques exist for generalizing integer results to fractional ones (see, e.g. (Ritter) 2000, Prop. 8)). Overall, we believe it useful to present, a self-contained proof for the purpose of (i) completeness in handling the non-integer case, and (ii) introducing tools that will be re-used in the noisy setting (Appendix D
These scaling laws are known to be order-optimal, since matching upper bounds (using non-adaptive algorithms) have been cstablished, c.g., NNoaak 1988 Sec. 1 & 3) for detailed summary. For the average-case criterion, see also Corollary with O.
C Proof of Theorem 3 (Noiseless Lower Bounds)
We generally follow the analysis of (Novak 2016 Thm. 1 & 3), while adapting the notation to match ours, and making some minor adjustments in the analysis.
C.1 Function Class Construction for the Matérn Kernel
We first describe bounded-support class of functions consisting of multiple bumps. The following lemma introduces the associated function and some useful properties.
Lemma 2. (Bounded-Support Function Construction Bull 2011 Lem. 5), (Cai & Scarlett
2021
Lem. 4))
Let lxx) = exp (1-īxx]= {[|x||2 < 1} be the d-dimensional bump function, and let go(x)
KO h
2x ) i.e.
rescaled version of for some
and > 0. Then, go satisfics the following properties:
(x(x)= for all outside the l2-ball of radius centered at the origin;
[o(x) [0, E] for all x, and go(0) =
l19o||ии I ci noo (1)" when is the Matérr-v kernel on Rd, where e₁ is constant.
In accordance with this lemma, let 9x(x) h( wo be the rescaled bump function with values in e0, and radius (iee, diameter w). For suitably-chosen M, we consider M such bumps with disjoint supports by shifting, and accordingly consider 2M functions, one for each possible sign pattern of the bumps. Mathematically, all functions of the form f(x) := Σ (x).(x) are contained in Hm, where ò; E {+1, -1} and (x) the shifted version of gu(x). Since the bumps form d-dimensional grid of step size in each dimension and the domain is D [0, 1]d, have
1
M
(21)
In fact, since the bumps have spherical support instead of rectangular, we could fit more of them into [0, 1]
but doing so would only impact the constant factors, which do not seek to optimize.
A ID example of f(x) is illustrated in Figure 3
Letting lo = Ji go(x)ds denote the integral of we have
1o M)
(22)
15

---

Published in Transactions on Machine Learning Research (06/2023)
Figure 3: 1D illustration of function in the restricted Matérn function class.
which follows from the fact that the integral of h(x) is constant, the vertical scaling shrinks the integral by G, and the horizontal scaling shrinks the integral by ((w)).
To evaluate the RKHS norm of each f, we use the equivalence of the RKHS norm of f and the Sobolev norm, as stated in Lemma ! Since each f is a sum of M disjoint-support functions, we can define these disjoint regions as Di,..., Das First considering the case of integer-valued (sce (17) for the relevant definition),
have:
1/2
Iffllw3
Σ 1100 fxx))(1)
M
2
1/2
ja
3igi (x)
12)
cx S
M
2 )
1/2
8,00 gi (x)
12)
ox S
i=1
M
2
2/2
8,00 (x(x)
dx)
co S
M
1/2
100 (xx) dxx)
(23)
ox S
D.
M
1/2
(1gg(x)) (1
M
I1g0 Iw3.
(24)
where 123
holds since the regions Di,..., DM are disjoint.
For non-integer values of s v 2. there always exists N+ such that < n and proceed
analogously to the analysis (Hesse 2006). Note also that our assumption + 3 ensures that n 2 1. From the definition of fractional Soboley norm in 19 and using Hölder's inequality, we can "interpolate" the
16

---

Published in Transactions on Machine Learning Research (06/2023)
fractional Sobolev norm between two integer Sobolev norms (we specify the constants X, n, P. and q below):
U/2 a1+ (1x1))M))))
inf
/ XEE (1+ 11x11))))
(25)
(1+ ||x (2) |-E)) [(E) (1 F(E))
XXED
(1 (1+ x ( (--
(5~A) Alp (2 nix) $ (f (1 (1+1 (1+
(26)
I
|| ||
24-2n
(27)
wh+1
M
I1ool I| 2 |1900
(28)
O
M
2n+2225 |1
|1900
12x-2
(29)
Hn-
n-1
22+2222
2 2n
. 0
(30)
(w-
32n+2
(w+++1 $ ) 220-20
Mc²
(31)
w
Me²
0
(32)
122
where:
25 holds since we are considering functions with support only within the subset D, and extending them to IRd would simply set f'(x to zero outside of D.
have chosen x = (n + 1)(s n), " = 2(s n) (0,2), atl-s > and q 8- 0, so that the pre-condition of Hölder's incquality holds ( + = 1).
In
26
In addition, (7 uses the identitics (s Aip n, (2 1)pp and Ag + which follow from direct substitutions and simplifications.
28 directly applics
uses Lemma 1 which establishes the equivalence of norms between the Sobolev space W and Matérn RKHS H. with parameter v = (i.e., the notation Hy equivalent to HM, but shows the explicit dependence on v/).
29
30
substitutes the RKHS norm bound for go, i.e. || HMc = 0 10 ), as demonstrated in Lemma
31
cancels (w with ( 222022 so that the exponent to w becomes (2n (2n-d) +2+ 2m) =
2d
8.
uses v d
Again using the norm equivalence between Sobolev function and the Matérn RKHS, can bound the RKHS norm of ff as follows in view of 32 and
21:
||
(33)
M
M
Equating the right-hand side of satisfying
33 with B and rearranging, it follows that || S B with choice of
B
0
(34)
M3+1
17

---

Published in Transactions on Machine Learning Research (06/2023)
C.2 Completion of the Proof of Theorem
3
The finite function (sub-)class described in Section C.1 is denoted by HM HM- We first consider the worst-case criterion. Suppose there is algorithm estimating T for function f HM, and the algorithm has budget of T M The best this algorithm can do is determine the signs of M out of the M bumps. For the unexplored regions, being in the worst-case setting, we can make an adversarial argument: The adversary either makes all these bumps negative or all them positive, whichever leads to a higher error. This leads two feasible values of that differ by ((M))) which turn implics that adversarially choosing the worse of the two must give f ((M Hence, when T
M 2: we have
EM west
(T) 2 2 2
To)
(35)
=MMM
(36)
= =Q({),
(37)
where
366
uses
22
and
37 uses
21
Substituting
34
and T = M into
2
37
we oblain
B
wst (T)
Ω
(27)
T
which establishes the desired worst-case lower bound.
For the averago-case criterion, we procced slightly differently. We first note that the supremum in E (T) supseNa -[[T-P] fl] can be lower bounded by the average with respect to f drawn uniformly from HM (and also still averaged over any randomness in the algorithm). Again considering budget of T 2 there M
must at least 2 regions with no samples, and for each such region, the associated integral is +10 or - lo with equal probability. Thus, conditioned on the observed samples, the posterior distribution of can be expressed a sum of two independent random variables, one of which of the form
M/2 S 2
8, Io.
(38)
where for notational convenience, we assume (without loss of generality) that it is the first 2 regions that have samples. That is, Su is random variable expressing the posterior unccrtainty in the M non- sampled regions. There may also be additional uncertainty due more than M regions being non-sampled, but for proving a lower bound, it sulfices to consider the case that there is no additional uncertainty beyond S Sa.
Since 8, is equiprobable on {+1, -1}, we have E[S = 0. Moreover, the variance of 38 is given by
M/2 Var(5; 1o) 2 12
03
(39)
4
since 0.. .dx l to are ii.d. random variables each having variance 1. In the limit of large M, if we consider the normalized sum
M
E
Zs
05 1
then by the central limil theorem (CLT) Feller 1957
Ch. VIII), as M 0o, Zsy convergences in
distribution to the standard normal distribution (O,, 1). In other words Sy is asymplotically distributed as o' When estimating Gaussian (or an asymptotically Gaussian quantity), an average absolute error proportional to the standard deviation is unavoidable; for instance, if we know that Z N/(0
18

---

Published in Transactions on Machine Learning Research (06/2023)
then our best guess of Z is 2 = 0, but this still gives E|| - 21] = - =
Thus, the lower bound for
the average case regret is ENg (T') 2 D(os ) 22 Mo)), and by substituting 21
22
34
obtain
(T) 2 ( VM Ω (T
avg
D Unified Derivation of Both Terms in Theorem
1 (Lower Bound)
In this appendix, we show that analyzing a suitably chosen hard subset of functions provides an alternative derivation of both terms in Theorem ! in a unified manner. Specifically, we consider the same class HM that was used in the noiseless case in Appendix [ but the analysis requires substantial modifications and is much more technical. We proceed in several steps.
Step 1: Reduction to simpler problem. Similarly to the noiseless setting above, we start by lower bounding EM ( oT, = supseиx Biir - A]] by B[|r- where now the average is taken over three sources of randomness: the uniform distribution over 2 functions in Hm (with parameters M and e), the randomization in the algorithm, and the noise. Moreover, once I is randomized, letting the algorithm be deterministic is without loss of optimality, so we assume this is the case (this is simply an instantiation of Yao's minimax principle).
We claim that in order to attain a lower bound with the preceding prior, it sulfices to attain a lower bound in the following simplified setup:
There exists an unknown collection of signs S. E [1] for i i= ......... M, each taking either value independently with probability 1
The goal of the algorithm is to estimate To EM Si, where To is the integral of single (positive) bump in the original BQ problem.
At time i, the algorithm may select an index i, (possibly in an adaptive manner) and observe Yt = eS,, + e with independent noise E ~ N(0,0²²).
We observe that this problem is exactly equivalent to our original problem in the case that the BQ algorithm is constrained to select midpoints of the bumps, where the bump takes its highest absolute value (i.e., e). Intuitively, this is without loss of optimality because such points have the highest signal, and are thus the most informative.
To make this more formal, we note that sampling a point with absolute value [f(x) = c E (0,e) gives y/ =S. +e, with S, being the associated bump sign, and this informatiion-hheoreecaalll equivalent to observing f eSt, + fe.. Since & > 1, this simply amounts still observing eSt,, but with more noise, and this extra noise could always be artificially added anyway. Hence, sampling at the midpoint is without loss of generality or optimality]
We proceed by studying this simplified problem.
Step 2: Establish hardness of estimating most sign values. As a stepping stone to characterizing the difficulty of estimating To MMM S;, we provide an auxiliary result on the hardness of estimating $ (S,,...,Svv) to within a certain Hamming distance. Although estimating each individual sign is a harder problem than estimating their sum (which may appear concerning from the perspective of proving lower bound), this will turn out to be useful intermediate step. We let S = (S.,..,, denote an estimate of S based on the queries.
Lemma 3. In the simplified setup with discrete queries is E (1,..., M3 (rather than x E D), consider any (possibly adaptive) deterministic algorithm that produces an estimate S of S. Then, there crists sufficiently
The locations with [f()] = carry information, so we can assume without loss of generality that they are never sampled.
19

---

Published in Transactions on Machine Learning Research (06/2023)
small constant such that we require a time horizon of
0²
c M max 1, 2
(40)
in order to obtain E[dn(s. respect to S uniform on the 2d possibilities, as well as the random noise.
$)] S M. Here du denotes the Hamming distance, and the expectation is with
Proof. By standard variant of Fano's inequality with approximate recovery (e.g., see (Sarlett & Cevher 2019 Thm. 2)), we have
M P (H
T(S:S) 21 log(2M) log ( N mmxx
S) 2
(41)
4
where .(:-) denotes the mutual information Cover & Thomas 2006), and Nmax is the number of vectors in {-1, 1 within Hamming distance A of any fixed vector (e.g., the all-ones vector). We have Nmax = Em (?) & M (M/44) from which a standard bound on the binomial coefficient gives Nmax I (MHz(/U((+o() with H2) = glog + (1 log T- being the binary entropy function. Since H2(1/4) is strictly smaller than log 2, substitution into (41 gives
M 4
I(S;S)+1 +1 ((M)
P (S 2
21
(42)
Moreover, following standard steps, we can upper bound the numerator as follows (Scarlett & Cevher 2019 Sec. 3):
Use the data processing inequality to write S( < I( Y), with (I, Y) being the length-T collection of sampled inputs and observed outputs by the algorithm.
Use the chain rule for mutual information to upper bound I( by corresponding sum over time indices: I(S Y) < I I(S:Y/1))
To simplify the last expression, we note that given 4, the only entry of that impacts Y, is St,, 80 we can further write I(S; 4,) < ((S,,;
Recall that when St, equals some value s E (-,, the corresponding observation is y N(s o²). Hence, by the relation between mutual information and KL divergence (Scarlett & Cevher 2019 Sec. 3.3), the preceding mutual information is further upper bounded by the KT, divergence between (O,, o²) (this is the same regardless of whether has a +1 -1 coefficient), which is 22
and
20
Substituting the preceding findings back into the preceding inequality 1(;;1, Y) S ET (S:Y it follows that (SS I, Y) S Uzs Hence, T c. Mgg with a small enough constant c, then the right-hand side of (42) exceeds z. The fact that SS 2 M with probability exceeding 2 then implies that IE [S,n(S,))) 3
a
The preceding argument proves the lemma when g 2 1. On the other hand, if o < 1, then the requirement in 44) simply becomes T 2 cM, which claim to be trivially necessary for attaining E[du(S, S)] 88 long as c < J. To see this, note that with any smaller number of samples, a quarter (or more) of the indices cannot even be sampled once. When this is the case, the algorithm cannot do any better than guessing the corresponding S; values, getting each one correct with probability 3.
The contrapositive statement of Lemma 3 that if T < c M max {1, 2 then it must hold that [ddu(S, S)] > M the preceding inequality can be written as
8. Furthermore, by writing the Hamming distance as sum of indicator function [ + S.).
M P(S +
Σ
S.1 8"
(43)
M T=1
20

---

Published in Transactions on Machine Learning Research (06/2023)
Step 3: Characterize the posterior uncertainty. In the argument that follows, are not directly interested in IP(S: * Si], but instead P(S; * S: D). where D (I, Y) contains the T pains of the form (it, Dt) collected throughout the course of the algorithm. This conditional probability can be viewed as representing the posterior uncertainty of S, with a value of 1 mcaning complete uncertainty, and a value of or mcaning complete certainly.
2
The probability in 43 is taken with respect the joint randomness of and the noise (which enters via D). While the decomposition P[S]P(D[S] is most natural, it is useful to consider the opposite form P[D]P(() D). that can analyze the posterior distribution P[
1 with respect to D (i.e., the With this in mind, we claim that the following holds with probability at least To equation to follow depends on D but still contains randommess via P[S)
M P
D
(44)
M
16"
To see this, assume by contradiction that this were only to hold with probability less than 18- Then, letting A denote the event that (44p holds, we would have
P 3
PlS. + S.DD
M
M
M Σ P(S. + 1 & A1]
M
E
P(S. + S.DD E A)
+E
M
where the first line uses the tower property of expectation. Then, the two terms are bounded as follows:
Using what we assumed by contradiction and upper bounding # DEEEPE + 3 [D] S the first term is at most 16.
M
Using the opposite incquality also at most
(44 for D & A, and upper bounding 10} S 1, the second term is
Thus, obtain 1 EE probability at least (with respect to D).
P(S, + 3.] S 1 which contradicts 43 and we conclude that 44 must hold with
We now apply a similar argument to the preceding one, but considering the uniform distribution over M implicit in 444 (as opposed to the distribution of D). Omitting the details to avoid repetition, it follows
1. that at least M of the indices in {1.....M} have P(S, + St(D] > 32: any smaller number than M 32 would
32 contradict (44h.
The above findings are summarized in the following lemma.
Lemma 4. Under the setup of Lemma 3 if
M max {1, {1
(45)
for sufficiently small 0, then with probabilily at least (with respect to D), there exist al least M indices such that P|S; # S, [D] 32
Step 4: Central limit theorem. We now return to the problem formed in Step 1, where the goal is to estimate lo EM S,, and the algorithm does not necessarily form any entry-by-entry estimate S. We continue to let D denote the samples collected, and we let D; denote the subset of D corresponding to times when is =
Lemma Under the uniform prior on S (Si...., SM,, conditioned on any collection of samples D, we have that the signs (S,...,SMs) remain conditionally independent.
21

---

Published in Transactions on Machine Learning Research (06/2023)
Proof. This is immediate from the fact that consider independent prior (namely, the uniform prior over all 2M sign patterns) and assume that the noise terms between times are independent. Thus, whenever some index in is selected, the resulting observation bears information about S; but bears no information about any the other Sj.
By Lemma conditioned on D, the posterior distribution of > M random variables. Moreover, by Lemma 4 when T satisfies (45
iml S is sum of independent =1-valued it holds with probability at least 16 that 32 fraction of these indices have strictly positive posterior variance. This, in turn, implies that
3
at least To EE S, has a posterior variance of (2M1)))
Having constant fraction of strictly positive-variant variance terms is sufficient for applying the central limit theorem for independent but non-identical random variables (Feller 1957 Ch. VIII). By doing so, we find .M that To Di=1 S is asymptotically Gaussian; the mean is inconsequential for our purposes, and the variance scales as (HM)), Le., the standard deviation R(VM Io). As highlighted in Appendix C when we have posterior standard deviation of Q(VM- Io), we incur S(VM Io) error. Since we have shown that this is the case with constant probability, follows that the average error ((VM 1o).
Step 5: Simplification. Recall from Appendix C that in our function class, we have M = 1 + Lo @(w ( and e O , Hence, the scaling M( To) can be expressed Q VM now complete the proof by considering two cases:
We
the maximum in 455 is achioved by the first term, then we have M = ((T) (supposing that is as high as possible subject. to 45)). Moreover, the above-established fact gives an Q xx=)) lower bound.
If the maximum in 445 is achieved by the second term, then we get M = =0 0 or equivalently
VM 0 (T) Hence, the lower bound is Ω (T) for both kernels.
Combining these two cases, obtain final lower bound of EM (T) ?(max x This is equivalent -(T5-1 of ))), and the proof of Theorem is complete.
E Proofs of Average-Cass Case Upper Bounds
E.1 Proof of Theorem 2 (General Guarantee for Algorithm 1
Our first result establishes that â is an unbiased estimator of R.
Lemma 6. Condition on arbilrary fixed values of y1,-.., (and hence, fixed i and R), and consider the resulting distribution of it due to the randomness in X1/241 XT and e1/241, ET- We have
B(R] R.
(46)
22

---

Published in Transactions on Machine Learning Research (06/2023)
Proof. Recalling that f is the initial estimate of based on the first samples, we have
E[
E
(ye - f(xx))
f=1/2-1
T
2
Elyt - f(xx)]
T
(=T/2+1
2
B(f(xx) f(x.) + te
T
(=1/2+1
2
Σ
B(f(x,) -(x))) +
Σ
[|(e]
T
T=//241
t=T/2+1
2
Σ
p(x) f(x)]dx +
D
(=//2+1
p(x)[f(x) f(x))xx= R.
where the second and fourth equalities are due to the lincarity of expectation, and the fifth equality holds since X2 ~ p(x) and E[c]] =.
In addition, the variance of the residual estimator Å is bounded according to the following.
Lemma 7. Under the selup of Lemma O we have
Var R < 2pmax
1ff
fl1E. +
20
(47)
T
where the variance is with respect to the randomness in XT/2+1..... and e ET.
Proof. We have
T
Var[ Var
2
Σ
(y - f(xx))
T
(=T/2+1
4 Var
(f(x,) - f(xx)) + Ei)
7
t=T/241
4
Σ
Var f(x)) f(xi) tc
T
1=1/2+1
4
Σ
Var(f(x,) f(xx,)
Σ
Var(e]]
T
T
T/2+1
T/2+1
4
20
Σ
Varls(x:)
f(xx))
(48)
R2
T
T/2+1
where on the third line we use the independence of x, and e, across t, on the fourth line we use the fact that 6 is independent of f(xx) fxxx; and the last line is due to Varles] = g². Moreover, by the definition of
23

---

Published in Transactions on Machine Learning Research (06/2023)
variance, we have
Var (xxi -B((())((((() f(x))) Elf(x)) f(x1))
I {((((x)) f(x)))
= S p(x)(((x))) (xx))
S Pmaxlls - f11221
(49)
where the last inequality is due our assumption that p(x) [O, Pmax] for all x. Substituting 49 into (8
obtain [47 as desired.
We can now analyze the error of our algorithm averaged over all samples, including the first T/2. Let E, (respectively, Ezl-) denote averaging with respect to the randomness from the first (respectively, second) batch. We first note that conditioned on the first T/2 samples, we have
E2 | i R1]
V E2 (I i [
E2 (R - [8)²|
(Ba [B2[(t] (B(BB
ft) "1
Ba
V Var[R]
(50)
2pmax T=
flf - fll 2OT
51
where the first. inequality follows from Jensen's inequality, the third line holds due to Lemma and the last two steps are due to Lemma and the elementary inequality v S v vT.
Then, incorporating the randomness from the first T/2 samples, obtain
E
| i
Al]
Eg
11 -- -fl]
/ 22max T 1 [1[1 - -f11z=] -20T-2
(52)
where we applied the tower property of expectation, followed by (51 (note that E, f[]:a] = Bill--fl1zz since no quantities from the second batch are present). This concludes the proof of Theorem
fllza]
E.2 Proofs of Corollaries 1 and (Noisy Matérn Upper Bounds)
It remains to upper bound the average-case T.--rrorr B [| between the true function and the estimate f. To do so, we will use results from WWyne et al. 2021 on the L2 estimation error for functions in Sobolev spaces, and adapt them to our problem for functions in the Matérn RKHS. These results depend on various technical assumptions made (Wynne et al.] 2021), some of which are trivially satisfied in our setting: our domain [0, 1]d implies their Assumption I, our Matérn-v functions with finite RKHS norm implies their Assumptions 2, 4 and 5. Moreover, their Assumption 3 only pertains to the misspecified setting, imposing the requirement that (2t has finitely many values (as we also assume for Corollary 2
We consider the commonly-used notions of fill distance hx and separation radius 4x, which are widely used (e.g., see (Wynne et al. 2021 Rieger & Zwicknag] (010)) and the associated convergence rates for kernel-based interpolation methods. For a point set X C D, the two distances are defined as
1 min ||x yll.
hx sup inf ||x yll. xeX yen
Qx 2 xxyEXxx/y
(53)
Intuitively, sufficiently small fill distance implies that covers the whole domain D. It is known that quasi- uniform point sets achieve the optimal order of hx and gx, at (TT (see, e.g. (Novak & Wošniakowsii
24

---

Published in Transactions on Machine Learning Research (06/2023)
2008 Thm. 4.17)). Morcover, (Santin & Haasdonk 2017
Cor. 11) shows that greedily minimizing the GP posterior variance with X = 0 leads to asymptotically uniform points, which in turn achieves small fill distance. Considering
sampling (Algorithm 2) in the first batch in our meta-algorithm (Algorithm 1. we are able to directly make use of the following results from (Wynne et al. 2021
Lemma 8. (L²-Error Result from (Wynne et al., 2021 Thm. 4)) Let f be any function in HM. and X = with Xa = 0, which outputs por with parameter Au- with hx I ho, and when v known, the average
{x,]E=1 be the sequence of points selected by Algorithm Then, there erists a constant ha > 0 such that VX noisy -eror
- -TT1 O AWllflla (1 B)
(54)
where E = ET) is the vector of noise terms
Moreover, under the modified version of Algorithm 2 for the misspecified setting (described just above Corollary where w and wt are also defined), we have
Elll p1]l 0 (bxin(uw )au -u Ox ux-"" +xu Affllme + 11 (hx X + 1) (1x
(55)
To obtain this lemma from (Wynne et al.
2021 Thm. 4), we substitute 0, Ty + -
+ 214 Th + μウ
+ d
2
Tk -ー + d. - 2 and m(.) in the notation therein. Setting A O and lx O(T b), 54 can be simplified as
Elll - pr 1 $ B + + [|][=11)
(56)
Similarly, by setting A,. = (T- hx = OTT and qx OTT k) 55
reduces to
Bill - [T]| = 0 miniy -it 1
B + T bB|lel))
(57)
Note that since the noises are i.idd Gaussian with zero mean, we have
T
||||||] <
B Blill
Ble2y2l V Varlei] + Varler/2] = =V 2 200
-V20
(58)
For E Hx and known v, substituting 58 and 56 into 52
we obtain
E 1 - f, - 1a|] = 0 VPmax'' $ (T-
8 T +E][e[]]) + oT
= -0 0 -1 + t)
which yields Corollary
For the misspecified setting, substitute
57
58
into
52 and obtain
E [1 1, - [R] = ( V/max T (T minto mint=-.w)+ wett
$B + T (Ellelll)) + tof
- (T- min(y [0) T + of o
which yields Corollary 2
E.3 Proof of Corollary (Noiseless SE Upper Bound)
Related to the fill distance discussed in the previous subsoction, has been shown in Wendland 2004) that the convergence orders for infinitely smooth functions turn out to depend exponentially on the fill distance hx. In particular for SE kernel, the 122 guarantee is illustrated in the following lemma.
25

---

Published in Transactions on Machine Learning Research (06/2023)
Lemma 9. (Noiscless L Guarantee (Wendland] 2004 Thm. 11.12)) Let f be any function in HSE, there exist constants ho, C. C. such that VX S D with h S ha, and when per calculated with X and X s exp [2(C, )o)) the I?² approzimation error bounded by
bx
Il - -rrlll = ( (ogghbx)))
(59)
As opposed to the Matérn kernel, the scaling of hx induced via maximum variance sampling for SE kernel functions not clear; see the discussion in (Santin Haasdonk 2017 Soc. 4.1). To overcome this, we adopt
the simpler approach of sampling on uniformly-spaced grid, i.e., X {( k1
I....... ...... & }){A; {o,..., 1}}
with N = T to ensure that xx = T (e.g., see Kanagawa & Hennig 2019 Thm. 4.3)). The substitution of hx = T into 55 then yields
|f - prllia = =
By applying Theorem 2 with 0 and À = 0 (Lemmn 9 still holds), the proof of Corollary is complete.
E.4 Proof of Corollary (Noisy SE Upper Bound)
The bulk of this subsection is devoted to introducing definitions and results from (Bach 2017). Although we are working lowards the noisy setting, all results stated are for the noiseless setting until staled otherwise.
It is often useful to study an RKHS through an integral operator M, which leads to an isometry with 22 space with mcasure dp:
(EN() f [(x)((x, )dp(x))
For the moment, consider any kernel (including SE) that can be written in the form
x(x,x)= = / 6(x, w)olx w)dp(w)
for some o(x,)) (dp) + (dp).. Let {x+TT=1 be iid samples drawn from density g with respect to measure dp, and pick weights B such that the approximation of
T
3
f()
o(x,)))
v((x))
which belongs to the Hilbert space Hh formed by the approximated kernel through T random features {wi) with the same density q:
I 1 i=1 g(w.) o(x, wi)
(x. y) T Σ
The weights 3 are chosen to solve the following minimization problem with Lagrange multiplier À > 0:
min llf flle2(00) + TAIBB113
B
with the solution
1 8
+ +11) -1 S. +11)"
T
FOOT
|| < T
(60)
T
where t. R.R._ + (dp)) is an operator:
T Σ
3
o(wi)))
V V((w))
1=1
26

---

Published in Transactions on Machine Learning Research (06/2023)
The quantity " in turn defines th following empirical integral operator Σ: (d - (d(p):
1 E g(wi) ($N()
(S, (i. 1) (1100)0
T
This allows us write f
f= >=>س>>>(>> -
(61)
where E1/2 is the unique positive self-adjoint square root of E, and forms bijection from (dpp) to the RKHS H4 Bach 2017).
With the above definitions, it was shown in Bach 2017) that the noiseless L2-error between f and f corresponds to the eigenvalue decay of the integral operator Σ if is the optimized distribution:
q(x)
M4 (x)²
14
(62)
where ,4 is the i-th largest cigenvalue of Σ, and (x) the corresponding cigenfunction. Specifically, Ilf 1ff-flleat (llp has geometric error (e.g., exp(-iz)) if μn decays coometrically/expoeeniie The guarantee that we make use of is formally stated as follows.
Lemma 10. Bachh 2017 Prop. 2) For the optimized distribution defined in 62 and the estimate f defined in 66 let > 0 and d, = AII and ossume that T 2 5d, log( 16da
Then, it holds with probability
least that
inf Il 31 1
sup - fl-lle - $2V..
Note that dy represents notion of effective dimension effective degrees of freedom. For the SE kernel integral operator, the eigenvalue decay satisfies the following (e.g.. see (Santin & Schaback 2016 Thm. 15)):
14 exxp(((C.))))
(63)
for some constant Ce 0. Moreover, Sun et al. 2018 Lem. 6) shows that, for #4 decaying according to 63 it holds that dx = O(log ))) Rearranging T T- (odx)og(>>) in Lemma 10 gives us that
d, = @ log T
T
Then, by equating 0 (1oTT) with (log ) and rearranging, we obtain
vz= o exp
T (1C
(64)
log
for some constant C, > 0. This determines the - between f and f in Lemma 10 in the absence of noise.
To account for the effect of noise, use the following lemma to extend Lemma 10
Lemma 11. (Bach) 2017 Sec. 5) Under the preceding setup, if each function query is corrupted by independent noise with variance not ezceeding gxx,)0² in the i-th entry, then we have
inf Ilf
sup $1 1
Ellif fllea(a)) <2vā+ o||l|||
While the presence of q(x;) in the precoding statement seems complicated, is fortunately considerably simplified in our case due to the following.
Lemma 12. Bach 2017] Sec. 4.4) For ahift-inwariant kernels in [0, 1)at (including SE), the optimized distribution q is the uniform distribution when the wniform measure dp is used.
27

---

Published in Transactions on Machine Learning Research (06/2023)
Therefore. have q(.) = 1 (i.c., the uniform distribution over [0, 1|9), and the additional noisy term in Lemma 11 is at most r||/8|| x 2oT due to the upper bound on || given in 600. Combining with (64)
the noisy L guarantee becomes
E [[I - fllmaaal] = o exp --
) + 3)
T
inf 1l $1
sup
log log T
Iflla S2
and substituting into our general BQ guarantee (Theorem obtain Corollary
F
Alternative Analysis Based on Confidence Bounds
In this section, we present analysis based constructing confidence intervals of the true function each time step, which is a popular strategy in the analysis Bayesian optimization (BO) algorithms. Our analysis is most closely related to that of the BO simple regret in (Vakili et al. 2021a
F.1 Noisy Setting
In this subsection, present our analysis in a slightly different way than our earlier results. Specifically, we first state the results with high probability (i.e., 1- 1-6), 8), and then convert them to expectations in the same manner as done in (Vakili et al.
2021a App. F), by setting o= V
We first provide several existing results that we use to prove results in 12 and
Lemma 13. Leee et al. 2022 Prop. & Remark 5) Let Hr. be the set of functions whose RKHS norm is
upper bounded by a constant B > 0. Then f is L-Lipschitz continuous with some constant L depending only on the kernel paramelers.
Tn the following, we use the shorthand notation = (xi, XT), and similarly for other quantities indexed by t. Recall also the posterior mean and variance defined in 44] (5), with parameter X >> 0.
Lemma 14. (Confidence Intervals Vakiii et al. 2021a Thm. 1)) Fir function f satisfying Ilfll S B, and assume Gaussian noises with variance o². Assume further that x1-1 are independent of EI4T, i.e., the points are chosen non-adaptively. For a fived x E D, for any E (7), define the upper and lower confidence bounds as
Uf(x) = (x(x) + (B+) (x( = - (B+
with B(8) = 22g and 8 (0, 1). Then, we have for any x E that
f(x) < (x)
w.p. at least 1 -
(65)
f(x) > (x)
w.p. at least - 6.
(66)
Lemma 15. (Adaption of (Srinivas et al.
2010 Lem. 5.4)) Letting yr sup Tyner fir)) where f.r (f(xir) R denotes the function values at the points N1...., Xr, we have
T
Σ
(xi) S log (1 + x²
2
(67)
1=1
Let D be finite subdomain of D j0, 1jd with d/2 points, with equal spacing of width 1 in cach
vT dimension. For any x E D, let [x] = arg min x ||x - xx112. By construction, we have, for any x e D that
< Vd 0
-
(68)
28

---

dlog + 220g 1 σ Published in Transactions on Machine Learning Research (06/2023)
By Lemma 13
the function f is L-Lipschitz. Thus we have for any D that
lf(x) [((x) L||x - (x]ا|22 0 V
(69)
For any fixed x E D, applying Lemma 14 gives the following with probability at least 1 - _
2|
f(x) 2 pT(x) B + 3
oT OT(xx.
(70)
D
By union bound over all x D, we have, for all x D simultaneously that
f(x) px - B + 3
[D] (x()x.
(71)
with probability at least 1 - 6 2. Similarly, we have, for all x E D that
f(x) S pT(() + B + 8 D )
or(xx),
(72)
with probability at least 1 - -.2 Combining 71 and 72 and again applying the union bound, we have, for
all D that
[f(x) NT (x)(x) s B + 3
or(x),
(73)
with probability at least 1 1- 8. We now extend the above upper bound to any point x in the domain D:
lf(x) pT(x)] S [(x) 5((x) + f(xx - pr ((x) + po (xx)g) pr(())
< [(x) - {(x1)) + ([((x)))--mm (T(());) + ( - mT(x))
s 0 I, (/4 + B 3 5 1 (T(()))))) 0
(74)
where the second inequality is by applying the triangle incquality, and the third incquality is due to and 69).
73
We will now show an upper bound on the posterior variance ox)).. Due to the decreasing property of posterior variance, we know that (x) S (x) for all x and t. Furthermore, due to the maximum variance sampling strategy, we have (xx) 2 o--1 (x). Thus, we have
or ([x))) S 0. ((x)) 1xx)
(75)
for all x and / 5 T. Squaring and averaging over / E [T] gives
( < T Σ o7 (x1).
(76)
t=1
and applying Lemma 15
wC obtain
2nr
or (x((1)p)) v
(77)
Toog
(1+
Substituting 77 and D = Td/2 into 74 and recalling the definition of 3(-) in Lemma 14
we obtain
(s(x - (T((x) s 0
(TT
+
2
(78)
Tlog (1+
=0
aoy (dlog T log 1))
B
(79)
T
29

---

Published in Transactions on Machine Learning Research (06/2023)
Finally, the absolute error of the above algorithm can be upper bounded by
/r S, - -pr(x)d & max f(x) - pr pr((x) / p(xx)xx - o( (VTT xEn
(dlog + log 2))
B (BBoo//
T
(80)
with probability at least 1 8. For the Matérn-v kernel, "r = O (log T) क-ं ) and for the SE kernel. YT (llg((T))d)) (see (Vakili et al. 2021b]). Thus, we obtain the noisy upper bounds in 12 and upon setting = /T as outlined at the start of this subsection.
F.2 Noiseless Setting
For the noiseless setting, we first state two useful lemmas, the first giving standard deterministic confidence bound, and the second relating the posterior variance and the fill distance hx (see (3))
Lemma 16. (Kanagawa et al. 2018 Cor. 3.11) For any f E Hy with llfl B, it holds that Lx s fxx & (x( for any and x E D, where
Ux(x) M4-1 (x) + Bo,-1(x))
[,(x) = (x) (x) - Bo,-1
and where M--1(-) and 0-)1(-) are given in (4 5 and calculated with X == 0.
Lemma 17. (Adaption of (Santin & Haasdonk 2017 Thm. 3)) Consider the set of points X : fxi]T=1
oblained by Algorilim 2 Then:
For functions in the Matérn RKHS, we have
max or = (h xED
(81)
2. For functions in the SE RKHS, there erists Cp > such that
max oT = 0
(82)
xED
For the Matérn kernel, as discussed in Appendix E.2 the fill distance being hx #), which is the best possible. Hence, the desired upper bound (13 for
is known that maximun-variance sampling leads to
the Matérn kernel follows directly from (81) along with Lemma 16 and similar steps to 78 (0)).
As mentioned in Appendix E.3 the fill distance of MVS on the SE kernel is not clear. To facilitate a similar
treatment, one may adopt the same grid points in 82 to derive an upper bound of O(e-C,TT However, a more recent result offers a more explicit dependence on the constant using grid points, as demonstrated below.
Lemma 18. Xu et al. 2022
Eq. (32))) Consider the domain D [0, 1]d and grid- brsdd subset X C D. Then, sampling each point once yields the following upper bound the noiseless posterior standard deviation (with x ===) 0):
max (x 0 (00--(115)
(83)
xED
Then, combining with Lemma 16 and proceeding similarly to 80), we obtain the noiseless upper bound in
15
G Additional Experimental Results
G.1 Further Synthetic Experiments
In Figures and we present additional results on synthetic kernel-based functions and benchmark functions. Overall, while there is no definitive ordering between the methods in general, observe similar
30

---

Published in Transactions on Machine Learning Research (06/2023)
findings to those discussed in Section 5.2 In particular, for the plots of error vs. split fraction, the trend can be decreasing (particularly at low noise), increasing (particularly at high noise), or "U-shaped", but 0.5 is generally reasonable choice. As we already discussed in Section 5.2, there is often a sudden drop at 1.0.
We note that some the MVS curves exhibit non-monotone behavior (e.g., for Alpine-2D), We believe that this is because the only randomness in MVS the noise and the initial points, whereas MC and MVS-MC have much more randomness due being randomized algorithms. When there is limited randomness and few queries have been made, the algorithm is essentially outputting uncertain guess, and it can happen that this guess luckily has a low error, but then this luck diminishes more samples are taken. In contrast, MVS-MC and MC have enough internal randomness to "average out" the lucky and unlucky scenarios.
G.2 Sensor Measurement Data
We consider the problem of estimating average sensor reading from limited queries, which each query consists of reading the value at given time instant. Note that since the algorithms consider are non- adaptive, the query times can be pre-computed. The data set consists of energy consumption readings for London Houscholds that took part in the UK Power Networks led Low Carbon London Project, between November 2011 and February 2014 8
We construct time-series signal (shown in Figure 7 of length 19,548 by sampling the data at intervals of one hour. Our goal is to estimate the average energy consumption during this period. Although the domain is now discrete, we let MVS and MVS-MC work on the contimous space, and round the selected decimal value to the nearest point in the data set. To create noisy BQ problem, we artificially add Gaussian (0, o²) noise to each query, with a E (0,0.1, 0.1, 0.5). The results are shown in Figure in this case, we found the various methods to perform relatively similarly to each other.
8Tee data be downloaded Bata..oonoo.ggo .o..
31

---

Published in Transactions on Machine Learning Research (06/2023)
Ackleyld
Ackleyld
Ackleyld D.5
102
20*
10
MVSING-SE
Ernor
WVS MC-S
Enror
MC-SE
13-4
10-4
10-
10a
100
250
100
150
200
Queries
Quaries
(a) Ac
Alpine2d
Alpine2d
Morte
Nomte
10
Error
Error
Eror
30:2
10-1
13-2
Queries
(b) Alpine 2D 2D
MAT_ Syntheiiteei
MAT_ Synt_Syceeii a=
MAT_Synthetic4d
Nome Carb
MYS-CCAAT
Ertr
Ener
1.
150
Queries
Queries
(c) Matérn-Syuthetien
SE Syntheticld - a=
Syntheticcd- a= 0.1
SE SE_Syytheeiccd 0.5
Monts Carde
Monte His-ss MUS-C--S
Monta Curoo
NVS se
wy se
NVS-NC-S2
WSSNC-5:
Error
Error
20-4
10-1
50
103 Queries
100 Queries
150
203
Queries
(d) SF-Synthetic-4D
Figure 4: Comparison of algorithms and the effect of noise.
32

---

Published in Transactions on Machine Learning Research (06/2023)
MAT Synthetic3d d=
MAT Synthetic3d d.1
MAT Synthetic3d o D.5
MISAIC-AT
+ NVS-MCMAT
+
NVS-ME-MTT
0.09
3.08
0.09
0.09
3.04
Froor
0.05
0.04
Erop
0.07
0.05
3.03
0.03
0.02
0.04
0.01
0.03
0.2 6.3 0.4 05 06 0.7
088 09 1.0
0.1
0.2
0.3
044
06 0.7
69 3.0
6.1 0.
0.3 0.4 05 D.6 0.7 08 09 1.0 MVS-MC Splits
MVS-MC Splits
MVS-MC Splits
Matérn Synthetic 3D
MAT Synthetic4d
MAT Synthetichd
MAT Synthetic4d - am 0.5
W/S-MC-HHI Hol
+
MVS-MCNAT
0.14
NUS-HCMAT
3.025-
0.12
0,600
3.035
0.08
0.05
0.230
0.055
0.000
0.02
2.1 022
0. 04 2. D. 07 0. 0. 1.0
0.2 0.2 5.3 0. 0.5 06 na
0.1
0.2
0.3 0.4 05 0.6 07 08 09 10 MVS-MC Splits
MVS-MC Splits
MVS-MC Splits
(b) Matérn Synthetic 1D
SE Synthetic3d am
SE Synthetic3d 0.1
SE Synthet p.5.
9.14
0.14
5.10-
0.10
0.10
0.08
0.05
0.88
0.05
0.00
0.02
MVS-MC Spiils
Splils
Splits
(c) SE Synthetic 3D
SE_Syythetic4d
Synthetic4d
SE Synthetic4d o- 0.5
MYS MC.
MVS-NCSSE
WVG IC.
0.250
0.04
3.125
0.33
0.03
grar
3:375
0.02
3.350
0.02
3.025
3.00
0.7 88 09 10
MVS-MC Splits
MVS-MC Splits
MVS-MC Splits
(d) SE Synthetic 4D
Figure 5: Comparison of different MVS-MC splits, i.e., the fraction of rounds for which MVS is used.
33

---

Published in Transactions on Machine Learning Research (06/2023)
Ackleyld-
Ackleyld a=0.1
Adkleyld o= 0.5
0.125
0.1
0.12
NVS-HCMAT MVS-MC-SE
0.12
+-
NYS-ME-MA
NVS-MC-SE
0.100
0.10
0.10
0.08
Erre
0.09
0.00
0.050
0.81
0.08
0.025
0.02
0.04
0.00
0.02
53
52
0.5
0.
3.5
5.6
5.7
0.5 0.0 1.0
0.1
0.2
0.3
0.4 0.5 0.6 0.7 0. 59 3.0
0.1
0.
0.3
04 05 D.6 0J 0.0 09 10
MVS-MC Spilts
MVS-MC Spiits
MVS-MC Spliss
Ackley-1D
Ackley2d a= 0
Ackley 2d 0.1
Ackley2d 07 0.5
0.30
WWSMC-WA
0.05
HUS-MEMTT MVS-MC-SE
0.08
MISNC-SE
0.250
SC--
3.125
3.06
0.84-
3.300
3.035
0.02
0.02
0.025
0.1. 0.2
0.3 0.6 0. 0.5 0.7
0.0 09 100
0.1 0.2 0.3 0.4
D.
05 10
MVS-MC Splits
MVS-MC Splits
MVS-MC SpiRt
(b) Ackley-2D
Griewank2d
Griewank2d 0.1
Griewank2d a 0.5
0.16
3.04-
0.18
0.03
0.12
0.10
0.02
0.08
0.05
0.03
0.04
3.50
0.02
100
0.5
MVS-MC Spiils
Splits
Splits
(c) Griewank-2D
Keane2d
Keane2d
Keane2d
3.54
CYSHC-WAT MVS-C-SE
0.05
HUS-MCMAT MVS-NC-SE
NUS-MC-MT NVS-MC-SE
3.56
0.06
Eror
Froor
0.10
0.64
0.08
0.05
0.01
0.02
033 0.7
68 09 10 MVS.MC Splits
0.1
02 03 04 03 DE MVS.MC Spiits
MVS-MC Spiits
(d) KD
Figure 6: Comparison of different MVS-MC splits, i.e., the fraction of rounds for which MVS is used.
34

---

Published in Transactions on Machine Learning Research (06/2023)
of pluN
(40001)
20-
mill
10
ApT
Ort
Jon
Apt
Ju
Det
Jen
Figure 7: Time-scries function for energy measurements.
Energy TS o=
Energy TS o= 0.1
Energy 0.5
Morte-Cardo MVS-MC-MAT MUS-SE
10a
Morte-Cario
Monte-Caric
MVS MC MAT MVS.SE
MVS MC. MAT MVS-SE
Error
MVS-MC-SE
Error
MVS-MC-SE
MVS-MC-SE
100.
10*
107
50
1000
150
200
250
50
100
150
200
250
50
300
150
200
250
Querles
Queries
Querles
Figure 8: Results for time series energy data.
35


Reviews:
[{"review_id": "Tj402r39to", "reviewer": "Reviewer_s6oR", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Both content and the presentation of this article fail to meet the standards of the journal. Indeed, the article lacks a paragraph, that follows the introduction, where the content of each section is described. In my opinion, the presentation may be challenging . For instance, the word ‘quadrature’ appears on the first page, while the first example of a quadrature rule appears on the sixth page. Morevover, it seems that the authors are not aware of the literature on control functionals. See, for instance, [1], where a variant of the meta-algorithm was proposed and analysed theoretically. \n\n\n[1] Oates, Chris J., Mark Girolami, and Nicolas Chopin. \"Control functionals for Monte Carlo integration.\" Journal of the Royal Statistical Society. Series B (Statistical Methodology) (2017): 695-718.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This article studies error bounds for kernel-based quadrature rules. These algorithms are widely used in the literature of numerical integration. The focus was put on the noisy settings where the algorithm is randomised, or the integrand is noisy. In particular, a meta-algorithm was proposed that build a quadrature rule using an L2-function approximation of the integrand. ", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "I believe that a major revision is necessary in order to include a comparison of the contributions of this work to the existing work on the literature of control functionals.", "broader_impact_concerns": "No", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "R8Qebn1JPr", "reviewer": "Reviewer_jwFj", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n+ The paper is well-written; all the formulas and the theorems are well-described. Insights of the theorems are also well-explained in the paper.\n+ I think the paper targets an interesting problem. It seems to me that the theoretical analyses are meaningful as understanding the error bounds of BQ is important. Although I can’t dig into detail to confirm the correctness of the theoretical analyses as I’m not very familiar with the theoretical analyses of BQ.\n+ The paper also tries to evaluate their algorithm to see if it is effective in practice.\n\nWeaknesses:\n+ I have one concern regarding the experiments. I don’t understand why to obtain better understanding of the method, in the experiments, the paper modifies to alternate between maximum variance sampling and Monte Carlo sampling? Isn’t this result in a different method? I see the justification of “mathematically, this does not change the behavior of the final step” is not convincing. Why don’t we run experiments with the proposed algorithm? And why alternating will not change the behavior of the proposed method? And is this invariant behaviour only for the final step or for all the iterations? I think more rigorous justifications and explanations need to be provided. Ideally, as a reader I would expect to see the performance of the original proposed algorithm – the algorithm that is used to conduct the theoretical analysis.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper provides several theoretical analyses for the error bounds of Bayesian Quadrature (BQ) in the noisy settings, randomized algorithms, and average-case performance measures. In particular, the paper introduces a two-step meta-algorithm that acts as a means to derive the average-case quadrature error with the L2-function approximation error. The theoretical analyses can be applied to different types of kernels including the Square Exponential and the Matern kernels. The paper also conducts some experiments in order to verify the behaviours of their algorithm.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "+ Please answer my concerns in the Weaknesses section.\n+ [minor] The plots in Figs. 1 and 2 are a bit small and the font sizes are also a bit small. Maybe it’s better to increase the sizes to ensure readability.\n", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "Xh7xBu8HAI", "reviewer": "Reviewer_towK", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "As detailed in the paper, little theoretical work on Bayesian quadrature (BQ) in a noisy setting exists, so any contribution to this direction would obviously be welcome. It is therefore a shame that the connection of the paper to Bayesian quadrature is tenuous at best. The paper is really about fundamental limits of computation in certain function spaces when the integrand evaluations are noisy than rather than error bounds for some commonly used Bayesian quadrature method in the noisy setting (the lower bounds obviously apply to such \"standard\" BQ methods as well). That is, the upper bounds apply to a certain two-step \"meta-algorithm\" that combines standard BQ with a Monte Carlo residual, but I have never seen this method used in the BQ literature. It seems that the numerical experiments in Section 5 do not show this algorithm outperforms standard BQ (I take it that MVS-SE and MVS-MAT refer to standard BQ, i.e. integration of Eq (4) with the integration points obtained from maximum variance sampling).\n\nIt is good that the authors are explicit that the lower bounds for Matérns have been derived by Plaskota a while ago and that they merely provide new proofs. Although there are statements here and there about how the proofs differ from those by Plaskota, one is left wondering how significant these new proofs really are: Can the technique yield new results in some cases? Is it simpler or more general?\n\nThe paper is generally reasonably well written. However, there are many issues with mathematics, both real and presentational.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proves upper and lower bounds on the average-case error of randomized numerical integration in RKHSs of Matérn kernels (i.e., Sobolev spaces) and in the RKHS of the squared exponential (SE) kernel, both for noisy and noiseless integrand evaluations. Upper bounds are proved by the means of a two-step \"meta-algorithm\" that uses maximum variance sampling to select half the samples and random sampling to select the remaining half.\n\nThe bounds for Matérns are not new, having been proven by L. Plaskota in 1996, but the authors claim to supply new proofs that may have some advantages over those by Plaskota.\n\nFor the SE kernel only upper bounds are proved. These appear to be new, though, as I will detail later, there are better bounds in the literature.\n", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Critical:\n\n1. Better bounds for the SE kernel are available in the literature. From Eq (11.11) in Theorem 11.22 of H. Wendland's 2005 book \"Scattered Data Approximation\" (a fairly standard source of BQ upper bounds) one easily obtains a noiseless worst-case upper bound of the form $\\exp(-C T^{1/d} \\log T)$, which is better than any bound in this paper. Wendland's Theorem 11.22 can also be used to improve Lemma 17 in Appendix F and the results that exploit this lemma. The only occurrence of the aforementioned bound in the BQ literature that I am aware of is as Theorem 2.20 in the 2019 thesis \"Kernel-Based and Bayesian Methods for Numerical Integration\" by T. Karvonen. A more explicit upper bound is probably obtainable by combining Eq (1.5) in Theorem 1.2 of https://arxiv.org/abs/2209.12473 with a basic tensor-product construction.\n\n2. In appendices, there are substantial issues in how the Fourier definition of the Sobolev norm is used. In Eqs. (18) and (19), the norm is defined by integrating over $D$ (an issue that is also present in the integrals on p16). This is not the correct way to proceed. Instead, one has to consider restrictions on $D$ of functions defined on $\\mathbb{R}^d$ (to define the Fourier transform you need a function on the whole real space rather than only on its subset). See e.g. Section 3.1.1 in Teckentrup (2020). \n\n3. Perhaps I have missed something here, but I think that Wenzel et al. (2021) only consider the case $\\lambda = 0$, which would mean that their results cannot be used (at least without modification) to prove Corollaries 1 and 2 in Appendix E.2. [Moreover, p24 refers to Theorem 11 in this paper but there is no such theorem in the paper - there is Lemma 11 though but that is not about the fill-distance.]\n\n4. From Appendix B one gets the impression that the authors are claiming a generalisation to fractional Sobolev spaces of the standard lower bounds for integration in Sobolev spaces. Such fractional versions are well known (at least for the worst-case error, and I would be very surprised if things were different for the average-case error). I do not recall what a standard reference for this is, but see at least Proposition 8 on p143 of K. Ritter's 2000 book \"Average-Case Analysis of Numerical Problems\". For approximation (as opposed to integration) the most general treatment is given in E. Novak & H. Triebel \"Function Spaces in Lipschitz Domains and Optimal Rates of Convergence for Sampling\", Constructive Approximation, 2006.\n\n5. The paper is full of references to Novak (1988). Most of these references are for specific facts about convergence rates etc. and should therefore provide more specific pointers than the entire book (e.g., theorem number).\n\n\nOther:\n- p2: \"maximxum\"\n- p3: \"...where the ABQ algorithm only achieves $O(T^{-\\nu/d})$ scaling\" I believe this is only because the proof technique that Kanagawa and Hennig use is sub-optimal.\n- p3: \"are have\"\n- p6: What is $\\mu_0$ in the input for Algorithm 2?\n- p7: \". we specifically...\"\n- p8: \"...note that that the techniques...\"\n- p13: In Appendix A it is stated that the norms in Eqs (17) and (19) are equal. In reality these define norms that are equivalent, which one sees by applying the binomial theorem to the polynomial in the Fourier transform and using the properties of the Fourier transform; after this it is seen that each $L^2$-norm of a weak derivative in Eq (19) is multiplied by a binomial coefficient not present in Eq (17).\n- p16: What are $\\mathcal{H}^{n - \\frac{d}{2}}$ etc. in Eqs. (25)-(30)?\n- p28: \"usefull\"\n- p28: I think there should not be \"with probability one\" in the statement of Lemma 15.\n", "broader_impact_concerns": "-", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: Two of the three reviewers identified strengths of this paper, with one being more positive.

The key concern of Reviewer jwFj was in the way the experiments were run (by alternating between a maximum variance sampling (MVS) phase and a Monte Carlo (MC) sampling phase). From my look at the paper, I believe the authors’ choice here was justified, as both sampling methods (as mentioned by the authors in the paper) are non-adaptive. Hence, for any time point corresponding to the last time within a MC sampling phase, the distribution of the sample collected thus far will be the same as long as an equal number of rounds were allocated to each of MVS and MC. Therefore, I believe this concern is resolved.

The authors give both upper bounds and lower bounds, which I discuss in turn.

**Upper bounds**

As far as I am aware, the results for the squared exponential (SE) kernel in the noisy setting are new. For the Matérn kernel, the authors show that a two-phase sampling algorithm can recover the rate that Plaskota (1996) previously showed for another method. Plaskota's method involves repeatedly resampling at the same points in a specific set of points. The authors' method instead is unlikely to ever sample the same point twice. To consider why the author's method may be advantageous, let us consider active learning. In active learning methods, it is sometimes considered unrealistic to be able to request labels at the same point multiple times (this may be akin to asking a person to repeatedly label the same image multiple times). I do think it is a nice advantage that the authors' method therefore avoids repeated sampling, while still avoiding adaptive sampling.

**Lower bounds**

Next, I briefly discuss the lower bound. While the lower bound rate the authors' achieve is not new, the way in which they achieve this lower bound, via a single hard set of functions to get both terms in the lower bound, appears to be new. This adds sufficient interest to the lower bound and may be useful in future work.

In total, I believe the revised version is suitable for publication at TMLR.
Decision: Accept (poster)