When can Regression-Adjusted Control Variate Help? Rare Events, Sobolev Embedding and Minimax Optimality
Authors: Jose Blanchet
Haoxuan Chen
Yiping Lu
Lexing Ying
Venue: NeurIPS 2023 poster (2023)
Keywords: Information-theoretic Lower Bounds
Sobolev Embedding Theorem
Quadrature Rule
TL;DR: ML-based control variates improve Monte-Carlo rate with sufficient smoothness. Truncated Monte-Carlo is optimal without sufficient smoothness, while ML-based control variates can't boost the rate due to rare and extreme events.

Abstract: This paper studies the use of a machine learning-based estimator as a control variate for mitigating the variance of Monte Carlo sampling. Specifically, we seek to uncover the key factors that influence the efficiency of control variates in reducing variance. We examine a prototype estimation problem that involves simulating the moments of a Sobolev function based on observations obtained from (random) quadrature nodes. Firstly, we establish an information-theoretic lower bound for the problem. We then study a specific quadrature rule that employs a nonparametric regression-adjusted control variate to reduce the variance of the Monte Carlo simulation. We demonstrate that this kind of quadrature rule can improve the Monte Carlo rate and achieve the minimax optimal rate under a sufficient smoothness assumption. Due to the Sobolev Embedding Theorem, the sufficient smoothness assumption eliminates the existence of rare and extreme events. Finally, we show that, in the presence of rare and extreme events, a truncated version of the Monte Carlo algorithm can achieve the minimax optimal rate while the control variate cannot improve the convergence rate.

When can Regression-Aduused Control Variates Help? Rare Events, Sobolev Embedding and Minimax Optimality
Jose Blanchet
Haoxuan Chen ICME
Department of MS&E and ICME
Stanford University
Stanford University Stanford, CA 94305
Stanford, CA 94305
j jose. wwchet@@tanffor..c.
haoxuancdstanf ord. edu
Yiping Lu
Lexing Ying Department of Mathematics and ICME Stanford University Stanford, CA 94305 Lexing@stanford. edu
Courant Institute of Mathematical Sciences
New York University New York, NY 10012 yiping. lu@nyu.uuc.... edu
Abstract
This paper studies the use of a machine learning-based estimator control variate for mitigating the variance of Monte Carlo sampling. Specifically, we seek to uncover the key factors that influence the efficiency of control variates in reducing variance. We examine a prototype estimation problem that involves simulating the moments of a Sobolev function based on observations obtained from (random) quadrature nodes. Firstly, we establish an ieformation-theoreee lower bound for the problem. We then study specific quadrature rule that employs a nonparametric regression-adjusted control variate to reduce the variance the Monte Carlo simulation. We demonstrate that this kind of quadrature rule improve the Monte Carlo rate and achieve the minimax optimal rate under a sufficient smoothness assumption. Due the Sobolev Embedding Theorem, the sufficient smoothness assumption eliminates the existence of rare and extreme events. Finally, show that, in the presence of rare and extreme events, truncated version of the Monte Carlo algorithm can achieve the minimax optimal rate while the control variate cannot improve the convergence rate.
Introduction
In this paper, we consider a nonparametric quadrature rule on (random) quadrature points based on regression-adjuster control variate 1 2 4 To construct the quadrature rule, we partition our available data into two halves. The first half used to construct nonparametric estimator, which then utilized control variate to reduce the variance of the Monte Carlo algorithm implemented over the second half of our data. Traditional and well-known results [1 Chapter 5.2] show that the optimal linear control variate can be obtained via Ordinary Least Squares regression. In this paper, 9 10), which uses a non-parametric machine learning- based estimator regression-adjustr control variate. We aim to answer the following two questions:
we investigate a similar idea for constructing quadrature rule 3
5
7
Is using optimal nonparametric machine learning algorithms to construct control variates an optimal way to improve Monte Carlo methods? What are the factors that determine the effectiveness of the control variate?
37th Conference on Neural Information Processing Systems (NeurIPS 2023).

---

p' we can directly embed Rare and extreme event
(a)
(b)
Figure 1: According to the Sobolev Embedding Theorem [T. the Sobolev space W can be embedded in L where 5F = á. When is large enough, as shown in (a), the smoothness assumption can rule out the existence of rare and extreme events. When is not sufficiently large, specifically < 2dg-dp, there may exist peak (a.k.a rare and extreme event) that makes the Monte Carlo simulation hard. Under such circumstances, the function's 2q-th moment is unbounded.
2pg
To understand the two questions, we consider basic but fundamental prototype problem of estimating moments of Sobolev function from its values observed on (random) quadrature nodes, which has wide range of applications in Bayesian inference, the study of complex systems, computational physics, and financial risk management []. Specifically, we estimate the q-th moment f (x of f based on values f(11). ... fzz,) observed on n (random) quadrature nodes 21. Iy E Q for function f in the Sobolev space Ws.(SS),, where f Rdd The parameter @ here is introduced to characterize the rare events' extremeness for estimation. To verify the effectiveness of the non- parametric regression adjusted quadrature rule, we first study the statistical limit of the problem by providing minimax information-theoreet lower bound of magnitude n max((( 2)9-1,-4 5 --
We also provide matching upper bounds for different levels of function smoothness. Under the sufficient smoothness assumption that d(2g-p) find that the non-parametric regression
2p4
adjusted control variate f can improve the rate of classical Monte Carlo algorithm and help attain minimax optimal upper bound. In (3.4) below, we bound variance (f - f of the Monte Carlo target by the sum of the semi parametric influence part Sq f2g f)² and the propagated estimation error [f f)22 Although the optimal algorithm in this regime remains the same, we need to consider three different cases to derive an upper bound on the semi- -arammertetri influence part, which is the main contribution of our proof. We propose new proof technique that embeds the square of the influence function (q and estimation error (f - f)² the Sobolev Embedding Theorem [1]. The two norms used for evaluating (/r-1)² and (S- f)²
2 in appropriate spaces via f}²
should be dual norms of each other. Also, we should select the norm for evaluating (f f)- f)² in way that it's easy to estimate f under the selected norm, which helps us control the error induced by (f 1--ft f)?. .. A detailed explanation how to select the proper norms in different cases via the Sobolev Embedding Theorem is exhibited in Figure the first regime when 8 > d.
in 0 (52) and attain final convergence rate of magnitude n - - - For the second regime when d(2g-p) s < d. the smoothness parameter is not large enough to ensure that f E L (Ω). Thus, 2[24-2)
p-
we evaluate the estimation error (j - î)² under the I norm and embed the square of the influence function (1 in the dual space of _ (S2). Here the validity of such embedding is ensured by the lower bound d(2g=p)
p(2q-2) 2) on Moreover, the pamimearaettric influence part still dominant in the second regime, so the final convergence rate the same as that of the first case. In the third regime, when d(2g-p) < d(2g-r) the aemi-paraietric influence no longer dominates and the p(2g-2)- 2j
2pg
final converge rate transits from n : to (((1-))-1
When the sufficient smoothness assumption breaks, i.e. d(2g-2) according to the Sobolev 2pg
dp < 2q. This Embedding Theorem [T], the Sobolev space W embedded in L La-op and d-sp indicates that rare and extreme events might be present, and they are not even guaranteed to have bounded 2 norm, which makes the Monte Carlo estimate of the q-th moment have infinite variance. Under this scenario, we consider a truncated version of the Monte Carlo algorithm, which can be proved to attain the minimax optimal rate of magnitude o - --1 In contrast, the usage of

---

regression-adjusted control variates does not improve the convergence rate under this scenario. Our results reveal how the existence of rare events will change answers the questions raised at the beginning of the section.
We also use the estimation of a linear functional as example investigate the algorithm's adaptivity to the noise level. In this paper, we provide minimax lower bounds for estimating the integral of fixed function with general assumption on the noise level. Specifically,, we consider all estimators that have access to observations { fIz; + 6i) of some function f that is s-Hölder smooth, 1id ~ n "O,, for some y > 0. Based on the method of two
Lid Uniform(O, ([0, 199) and 6 where 2; fuzzy hypotheses, we present a lower bound of magnitude n __ -.------ 1--2))))) 2), which exhibits smooth transition from the Monte Carlo rate the Quasi-Monte Carlo rate. At the same time, our informatiionthhoreee lower bound also matches the upper bound built for quadrature rules taking use of non- con-maramttric gggression-adjustee control variates.
1.1 Related Work
Regression-Adjuste Control Variate The control variate method a technique used for variance reduction in Monte-Carlo simulation. Consider the task of estimating the expectation EX for some random variable X. The idea of control variate method to introduce another random variable Y correlated with the random variable X, such that the random variable X - Y has smaller variance than X. Since E.X = E[ - Y| E{{] and Y is deterministic, one may obtain a variance reduced estimator of IE [X] by summing up E and empirical estimate of E[X YJ. Such random variable Y is called a control variate. Regression-adjusee control variate, in particular, refers to the case when Y is obtained by applying regression methods to observed data samples of X.
Regression- -adjustee control variates have shown both theoretical and empirical improvements in a wide range of applications, including the construction of confidence intervals [1[33, randomized trace-estimation [5], dimension reduction [16), causal inference [17], light transport simulation
18
MCMC simulation
[19], estimation of the normalizing factor [10) and gradient estimation
20
21.
It is also used as a technique for proving the approximation bounds on two-layer neural
networks in the Barron space [22]
Regarding literature most related our work, we mention 3 7 8 10), which also study the theoretical properties of nonparametric control variate estimator. However, the theoretical analysis in 3 7 does not provide a specific convergence rate in the Reproducing Kernel Hilbert Space, which requires a high level of smoothness for the underlying function. In contrast to prior work, our research delves into the effectiveness of a non- paramataeeec regression-adjusee control variate in boosting convergence rates across various degrees of smoothness assumptions and identifies the key factor that determines the efficacy these control variates.
Quadrature Rule
There is a long literature on building quadrature rules in the Reproducing Kernel Hilbert Space, including Bayes-Hermin quadrature [23]
24
25| 26 27]. determinantal point
processes 128 29 30 31], Nyström approximation [32 33 kernel herding 34 35 36) and kernel
thinning 37 38 39] Nevertheless, the quadrature points chosen in these studies all have the ability to
reconstruct the function's 's information, which results in suboptimal rate for estimating the moments.
Functional Estimation
There are also lines of research that investigated the optimal rates of 47
estimating both linear 88
40
41
42
43
44 45 46
48 49 [50] and nonlinear [51 52 53 54 55
[56 57 (58 59
functionals, such as integrals and the L norm. However, as far as the authors know, previous works on this topic have assumed sufficient smoothness, which rules out the existence of rare and extreme events that are hard to simulate. Additionally, existing proof techniques are only applicable in scenarios where there is either noise constant level of noise present. We have developed a novel and unified proof technique that leverages the method of two fuzzy hypotheses, which allows account for not only rare and extreme events but also different levels of noise.
60
61
62
63
64
1.2 Contribution
We determine all the regimes when quadrature rule utilizing a nonparametric estimator as control variate to reduce the Monte Carlo estimate's variance can boost the convergence rate of estimating the moments of a Sobolev function. Under sufficient smoothness assumption,
3

---

Smoothness
Embed fin L
Embed L
Fstimation in I.2
Estimation
Choose embedding good for both evaluating the semi-parametric hardness and function estimation
1/2
--
Truncate Monte Carlg
Regression-adiissee Control Variate
Minimax rate
We summarize the minimax optimal rates and the corresponding optimal algorithms with respect to the function smoothness here. When the function is smooth enough, gggression-adjusee control variates can improve the Monte Carlo rate. However, when there exist rare and extreme events that are hard to simulate, truncating the Monte Carlo estimate directly yields a minimax optimal algorithm. Above the transition point of algorithm selection
Figure 2:
d(2g-P) while the transition 2pg
point of the optimal convergence rate s d(24-p) build the optimal convergence guarantee for p(2q-2)"
any algorithm that utilizes a pepressioonaadjuseee control variate f. we need to embed the square of the influence function 1=² in appropriate space via the Sobolev Embedding Theorem and evaluate the estimation error (ʃ f)² under the dual norm of the associated with the chosen space, which allows us to achieve optimal emiipparametric efficiency. Our selections of the metrics in different regimes shown this figure.
which rules out the existence of rare and extreme events due to the Sobolev Embedding Theorem, the regression-adjusted control variate improves the convergence rate and achieves the minimax optimal rate. Without the sufficient smoothness assumption, however, there may exist and extreme events that hard to simulate. In this circumstance, discover that truncated version of the Monte Carlo method is minimax optimal, while ongression-adjuseee control variate can't 't improve the convergence rate.
As far the authors know, our paper is the first work considering this problem without assuming that the underlying function f is uniformly bounded. All previous work assumed that > d which implies f E L L (St) and neglects the possibility of spike functions. As a p. result, they were unable to discover the transition between the two regimes described above.
Under the assumption that d(2q-p) the main difficulty in establishing the convergence guarantee lies in determining the right evaluation metric for function estimation. To select suitable metric, we introduce new proof technique by embedding the influence function into an appropriate space via the Sobolev Embedding Theorem and evaluating the function estimation in the corresponding dual norm to achieve optimal semi- semi-paramettr efficiency. Our selection of the proper embedding metrics is shown in Figure 2
2pg
To study how the regression adjusted control variate adapts to the noise level, we examine the linear functionals, ie. the definite integral. We prove that this method minimax optimal regardless of the level of noise present in the observed data.
1.3 Notations
Let || be the standard Fuclidean norm and = [0, 1]d the unit cube in Rdd for any fixed d e N. Also, let [(-} denote the indicator function, i.e, for any event A we have 1 1{A} 1 if A is

---

true and 1{A} = 0 otherwise. For any region R S Ω2 we use (( :- Sq 1fx E Rjdz to denote
the volume of R. Let ((2) denote the space of all continuous functions f: f Ω2 + R and [-] be the rounding function. For any s and f e ((2)) we define the Hölder norm || (c))(0) by
|D [() D* [())
II == max [
| fll (0) + max |M= NER,
sup
(1.1)
|| - y1|--1a)
The corresponding Hölder space is defined as C* == { E C(2) : II < { When
8 = 0, we have that the two norms || (l) and || |1,,) are equivalent and (²) = (1)).
Let No == (0) be the set of all non-negative integers. For any & E No and S S oc, we define the Sobolev space W (S2) by
W*P(2)) = {s ( D (1)) N6 satisfying |a| s $
(1.2)
Let (c)+ denote max(c,0) 0} for any
R. Fix any two non-negative sequences fan Jn= and (bn 1x0 1n=1- We write an S bas or an = O(b)). to denote that ag < Cb, for some constant C write an 2 bn or a, = w(b)), to denote that a, 2 cb, for some
independent of n. Similarly, constant e independent of n. We use an = (b,.) to denote that On = O(h)) and an = whh)).
2
Information-Theerrti Lower Bound on Moment Estimation
Problem Setup To understand how the non-parametric regression-adjusted control variate improves the Monte Carlo estimator's convergence rate, we consider prototype problem that estimates function's 's q-th moment. For any fixed
N and f E (S),()), we want to estimate the q-th moment 1 := Jn f with n random quadrature points { [z, Ji=1 с Q. On each quadrature point x; (i = 1, ..,,,,,,, ,.), we can observe the function value yi := fx;).
In this section, we study the information-theeorttee limit for the problem above via the method of two fuzzy hypotheses (65). We have the following information-theoreti lower bound on the class H444 that contains all estimators H S x R - IR of the q-th moment r7.
Theorem 2.1 (Lower Bound on Estimating the Moment) When p > 2 and q < < 2g, let Hi denote the class of all the estimators that use quadrature points in and observed function values fyr (r to estimate the q-th moment of f. where {x,)= are independently and
identically sampled from the uniform distribution on Q. Then we have
inf sup AIEHt! (e)
E(
[ (1: ( r3
2 n max - (-1(5-3)-1,-3311
(2.1)
Proof Sketch Here we give a sketch for our proof of Theorem 2.1 Our proof is based on the method of two fuzzy hypotheses, which is a generalization of the traditional Le Cam's 's two-point method. In fact, each hypothesis in the generalized method is constructed via a prior distribution. In order attain lower bound of magnitude A via the method of two fuzzy hypotheses, one needs to pick two prior distributions Po: Pl on the Sobolev space (S) such that the following two conditions hold. Firstly, the estimators Iÿ differ by A with constant probability under the two priors. Secondly, the TV distance between the two corresponding distributions P. and P, of data generated by Ao and μu is of constant magnitude. In order to prove the two lower bounds given in (2.1). we pick two different pairs of prior distributions as follows:
Below we set mn ( and divide the domain Ω into md small cubes Ωn Ω2 Pm,,, each of which has side length m -1 For any p E (0, 1), we use Up Wp to denote the discrete random variables satisfying Pe(,
(wu - 1) p and IP(,p- 1) = Puu, = 1) 1 p.
(I) For the first lower bound in (.. we construct some bump function g WSP(SS) satisfying supp(g) C Ω and Ig Sa, g(x)ds O(m ((m)(=++))0 Now let's take some sufficiently small
constant E (0, 1) and pick no, μ to discrete measures supported on the two finite sets (11= 9
and ( On the one hand, the difference between the q-th moments under f4 and μr can be
5

---

lower bounded by 0(1x(3-))11) with constant probability. On the other hand, K ( can be upper bounded by the KI. divergence between V and L which of constant magnitude.
(II) For the second lower bound in (2.I). we set M > 0 to be some sufficiently large constant and A For any 1 j m", construct bump functions f, W** satisfying supp(f,) S2 and / = Ja, fz(r)dr = @mm=ks=- for any S j S m and 1 I k S s. Now let's pick fo: M1 to m 3 Sj] and M [ (1) S3 be discrete measures supported on the two finite sets M
where {w, mm Jj=1 and {wo ( Jj=1 are independent and identical copies of w and W- respectively, On the one hand, applying Hoeffding's 's inequality yields that the g-th moments under /An and Ja differ by O(n -- 2) with constant probability. On the other hand, note that K o can be bounded by the KL divergence between two multivariate discrete distributions (u)), (0))) and (wD) (1) .... : ). where {w, Ji=1 and {w), (1) ym are independent and identical copies of w and Wing respectively. Hence, K ((Poo|(P)) ) is of constant magnitude.
Combining the two cases above gives us the minimax lower bound in 2.1 We defer a complete proof of Theorem 2.T to Appendix B.2
3 Minimax Optimal Estimators for Moment Estimation
This section is devoted to constructing minimax optimal estimators of the q-th moment. We show that under the sufficient smoothness assumption, a regression-adjustee building minimax optimal estimators. However, when the given function is not sufficiently smooth, we demonstrate that a truncated version of the Monte Carlo algorithm is minimax optimal, and control variates cannot give any improvement.
control variate is essential for
3.1 Sufficient Smoothness Regime: Non-parametric Regression-Adjuster Control Variate
This subsection is devoted to building a minimax optimal estimator of the q-th moment under the 1 which guarantees that functions in the space W* are sufficiently smooth. assumption that å > 20 From the Sobolev Embedding theorem, we know that the sufficient smoothness assumption implies WS)(SS) C L (S2) c L (S2), where p- uniformly sampled quadrature points (a: N and corresponding observations {u+ [(=:))=_ of f, the key idea behind the construction of our estimator H&. is to build nonparametric estimation f of f based on a sub-dataset and use f as a control variate for Monte Carlo simulation. Consequently, it takes three steps to compute the numerical estimation of 1 for any estimator HG S2" xIR" + R. The first step is to divide the observed data into two subsets S {( )=))==n:> S2 ((E,, w)]]
á. Given any function f E W>*P(S) along with
of equal size and use machine learning algorithm to compute nonparametric estimation fi. of f based on Sj. Without loss of generality, we may assume that the number of data points is even. Secondly, we treat fi.y as a control variate and compute the q-th moment 1 Using the other dataset
S2. we may obtain a Monte Carlo estimate of 13 f1.4
as follows: 17 7
R " 2 (v)
11.4
fly (a:)) Finally, combining the estimation of the q-th moment 7 = Ja fi.y (x)dx with the
h.
estimation of 17 1", h gives us the numerical estimation returned by HC:
2 frg (miHi (m))))))))))) = M.g( ) +
Σ
12
v/
Mi. (1i)).
(3.1)
n +=
We assume that our function estimation f is obtained from an 2 oracle K y satisfying Assumption 3.1 For example, there are lines of research 49
Ω x R
W
50
59
60
considering
how the moving least squares method [66] 67 can achieve the convergence rate in

---

Assumption 3.1 (Optimal Function Estimator as an Oracle) Given any function f E (9)
and n N, let {aili= be n data points sampled independently and identically from the uniform
distribution on Qt Assume that for s > dd there exists an oracle Kn : 52" x R" + (
2pg
that estimates f based on the n points { {x,);- Vi=1 along with the n observed function values {f(x;) [==
and satisfies the following bound for any satisfying 1 E |max{ $ .0), max p: 1 1 > }1
7 In -
[B (BAD [1Kn (a( - - -1
(3.2)
A construction of the desired oracle and a complete proof of the upper bound above (up to logarithm factors) is deferred Appendix [E] bound that matches the information- -hteoreticctc lower bound in Theorem
Based on the oracle above, we can obtain the following upper
Theorem 3.1 (Upper Bound on Moment Estimation with Sufficient Smoothness) Assume that 2dg-dp, Let {ri)== be n quadrature points independently and 2pg
p > 2. q < 2q and >
identically sampled from the uniform distribution on St and {y: == f(x;) R= be the corresponding
n observations of I E WS). Then the estimator HGG constructed in (31 above satisfies
E (ac) | (
S n max (5 1331,-5-1
(3.3)
Proof Sketch Given a non- non-parametric estimator f of the function f. we may bound the variance of the Monte Carlo process by (fa- - f0)² and further upper bound it by the sum of the following two terms:
1f9 fa x (f AP + (f (f-
(3.4)
semti parametricc influnce
estimation ettor propagation
The first term above represents the parametarectticc influence part of the problem, as g is the influence function for the estimation of the q-th moment f". The second term characterizes how function estimation affects functional estimation. If consider the special case of estimating the mean instead of a general q-th moment, i.e, q = 1, the semi parametric influence term will disappear. Consequently, the convergence rate won't 't transit from n 3 à to -9 - = in the special case.
Although the algorithm remains unchanged in the sufficient smooth regime, we need to consider three separate cases to obtain an upper bound on the integral of the semi-parametric influence term [ - 12 in (.4). An illustration of the three cases is given in Figure 2
From Hölder's 's inequality, we know that Jn ( (x)( - 2 dx can be upper bounded by || (ay]1(f |l( f) where || II cay and ||]| (9) are dual norms. Therefore, the main difficulty here is to embed the function I different spaces via the Sobolev Embedding Theorem under different assumptions on the smoothness parameter 8. When the function is smooth enough, i.e. >>
d. we embed the function f in (??()) and evaluate the estimation error f under the L norm. Then our assumption on the oracle (3.2) gives us an upper bound of magnitude n 3 on Ilf - flli: which helps us further upper bound the semi- cami-parametric influence part p-
So (=) - f(x))22z dx by n up to constants. When d(2g-n) p(24-2) < <
d we embed
the function f in T. 24 (5) S I. - (S2) and evaluate the estimation error f- under the T norm. Applying our assumption on the oracle (3.2) again implies that the semi-parametric influence part Sa (24-
f( can be upper bounded by - - up to constants. When d(2g-p) we embed the function f in Le and evaluate the error of the oracle in p(2--2j-
d124-p) <
2pg
L where = 3. Similarly, we can use (.. to upper bound the semi-parametric
influence part Jaen (z) - 2 by n 20(
The upper bound on the propagated estimation error Saen (f(x) f(x)) in 3.4) can be derived by evaluating the error of the oracle under the 2 norm. i.e by picking r = 24 in 3.2 above, which yields upper bound of magnitude n: 2 -1
7

---

The obtained upper bounds on the semi- pamimiareettric influence part and the propagated estimation eтror above provide us with clear view of the upper bound on the variance of f fr, which is the random variable we aim to simulate via Monte-Carlo in the second stage. Using the standard Monte-Carke algorithm to simulate the expectation of fr then gives us an extra n factor for the convergence rate, which helps attain the final upper bounds given in 333 A complete proof of Theorem 3.1 is given in Appendix C.I
3.2 Beyond the Sufficient Smoothness Regime: Truncated Monte Carlo
In this subsection, we study the case when the sufficient smoothness assumption breaks, i.e. 1 According to the Sobolev Embedding theorem, we have that WP is embedded in L "p. dp- < 2g, the underlying function f is not guaranteed to have bounded 2 norm, which indicates the existence of rare and extreme events. Consequently, the Monte Carlo estimate of f's q-th moment must have infinite variance. which makes it hard to simulate. Here we present a truncated version of the Monte Carlo algorithm that can achieve the minimax optimal convergence rate. For any fixed parameter M > 0, our estimator is designed as follows:
Since -
implies T-sp
H (a (u.2 =
max {
min M). -M]
(3.5)
In Theorem 3.2 provide the convergence rate of the estimator (3.5) by choosing the truncation parameter M in an optimal way.
Theorem 3.2 (Upper Bound on Moment Estimation without Sufficient Smoothness) Assuming that p 2, q < P < 24 and s <
2p4 we pick M O(nZ-á) .). Let {si)i=z be quadrature points independently and identically sampled from the uniform distribution on and {v. = f(x))]= be the corresponding n observations of f E WS.P(S).. Then we have that the estimator H constructed in (3.5) above satisfies
2dg-dp,
E Ex: Lo
Hy (0,)1 (wHH=1 -17 17
S n a
(3.6)
Proof Sketch The error can decomposed into bias and variance parts. The bias part caused by the truncation in our algorithm, which is controlled by the parameter M and can be bounded by ((:l))e)))> my lf According to the Sobolev Embedding Theorem, W**(s)) can be embedded in the space LP where 1 p"
.. As [s()) M implies [](xx)] < MM-P"1 (=))) the bias can be upper bounded by M9-2", Similarly, the variance is controlled by M and can be upper bounded by
M ". Combining the bias and variance bound, we can bound the final error as M F² + M 4
Vn
By selecting M O(nF* O(n= 2), obtain the final convergence rate 2( }) A complete proof of Theorem 3.2 is given in Appendix C.2
Remark 3.1 164 has shown that the convergence rate of the optimal non-paranne ann-parametric regression-
based estimation is n - which is slower than the convergence rate of the truncated Monte Carlo estimator that we show above.
4 Adapting to the Noise Level: Case Study for Linear Functional
In this section, we study how the regression-adjusted control variate adapts to different noise levels. Here consider the linear functional, i.e. estimating function's definite integral via low-noise observations at random points.
Problem Setup We consider estimating I, f (z)dr, the integral of over Ω for fixed
function f E (+2) with uniformly sampled quadrature points (xx) C Ω. On each quadrature point v; (i = 1,................ n), we have noisy observation yı = f(x,) [[x;)+6, +.. Here the e, are independently and identically distributed Gaussian noises sampled from (/ n -27), where [o).
8

---

4.1 Information-Theown Lower Bound on Mean Estimation
In this subsection, we present minimax lower bound (Theorem 4.1) for all estimators H 12" xR= R of the integral y of function f E (*(2) when one can only access noisy observations.
Theorem 4.1 (Lower Bound for Integral Estimation) Let H denote the class of all the estima- tors that use n quadrature points |i= and noisy observations {y= f(x)) + ai to estimate the integral of f. where {xi) b and {c))== are independently and identically sampled from the
uniform distribution on St and the normal distribution N (0,, ((,,n22) 22) respectively. Assuming that ry [0] and > O, we have
inf sup (x. HEHA (EC
(1 1/1
2 n maax(-22-5-22
(4.1)
Remark 4.1 Functional estimation is a well-studied problem in the literature of nonparametric
statistics. However, current coormation-theorer lower bounds for functional estimation [5T 52 53
56 [5 58
65 [68] assume a constant level of noise on the observed function values. One essential idea for proving these lower bounds is to leverage the existence of the observational noise, which enables to upper bound the amount of information required to distinguish between two reduced hypotheses. In contrast, we provide a minimax lower bound that is applicable for noises at any level by constructing two priors with overlapping support and assigning distinct probabilities to the corresponding Bernoulli random variables, which separates the two hypotheses. A comprehensive proof of Theorem 4.1 is given in Appendix D.2
4.2 Optimal Nonparametric Regression-Adjusted Quadrature Rule
In the discussion below, we use the nearest-neighbor method as an example. For any k the k-nearest neighbor estimator fu-NN f given by fu-N((=) = 1 ビラー W1(0), where { is a permutation of the quadrature points { such that 1|0 - 2|| S
21,2,--
b =11 ... S |1 =|| holds for any E Ω. Moreover, we Tki := t 1س=1
to denote the collection of the h nearest neighbors of z among { for any 2 E Ω. For any 1 1 i < 12 2 we take D, C Ω to be the region formed by all the points whose nearest neighbors
contain Ty Le, D; = { 2 E Ω: :0:4 E Tr. Our estimator ARNNN can be formally represented as
HN-NN (1:))- (ud
V(D.)
2
Σ
12
2
n
Σ j=1
1z E D2lvi)
y1
1/4
T1
n
1=
i=
++1
Jq fRoN(
Ei=g+1 3+1
In the following theorem, we present an upper bound on the expected risk of the estimator H
Theorem 4.2 (Matching Upper Bound for Integral Estimation) Let {xi)i be n quadrature points independently and identically sampled from the uniform distribution on Ω and {y = f(x)) + 4 be the corresponding a noisy observations of f E (*)), where {e) are independently and identically sampled from the normal distribution N(,, (200, Assuming that E [0, co) and E (0, 1), we have that there exists l E N such that the estimator HENN constructed above satisfies
] maxi - -
E (1. 1a (u
11
(4.2)
Remark 4.2 Our upper bound in Theorem 4.2 matches our minimax lower bound in Theorem 4.1 which indicates that the adljstedd quadrature rule associated with the nearest neighbor estimator is minimax optimal. When the noise level is high (7 4). the control variate helps

---

to improve the rate from n (the Monte Carlo rate) to n. 2 via eliminating all the effects of simulating the smooth function. When the noise level low (7 > z). we show that our estimator NN can achieve the optimal rate of quadrature rules [46]. We defer complete proof of Theorem Appendix D.3
5 Discussion and Conclusion
In this paper, we have investigated whether a non-parametric regression-adjusted control variate can improve the rate of estimating functionals and its minimax optimality. Using the Sobolev Embedding Theorem, we discover that the existence of infinite variance rare and extreme events will change the answer to this question. We show that when infinite variance rare and extreme events are present, using con-parametric machine learning algorithm control variate does not help to improve the convergence rate, and truncated Monte Carlo minimax optimal. When the variance of the simulation problem is finite, using a regression-adjuster control variate via optimal non-paramttric estimator minimax optimal.
The assumptions we made in this paper, such as boundedness of the domain l and constraints on the parameters p, q. might be too restrictive for some application scenarios. We left relaxations of these assumptions as future work. One other potential direction is to investigate how to combine importance sampling with regression-adjje adjjusted control variates. Also, the study of how gegresssoonadjjusee control variates adapt to the noise level for non-linear functionals [62 63 may be of interest. Moreover, another intriguing project analyze how the data distribution's information achieve both better computational trackability and convergence rate g
can be used to
Acknowledgments and Disclosure of Funding
Jose Blanchet supported part by the Air Force Office of Scientific Research (AFOSR) under award number FA9550-20-1B 20-1-0397 and the National Science Foundation (NSF) under award number DMS-1915967. Haoxuan Chen is supported by the T. Lo Graduate Fellowship Fund. Yiping Lu is supported by the Stanford Interdisciplinary Graduate Fellowship (SIGF). Lexing Ying is supported by the National Science Foundation (NSF) under award number DMS-201 1699999999999 and DMSS-200 16..
References
[1] Søren Asmussen and Peter W Glynn. Stochastic simulation: algorithms and analysis, volume 57. Springer, 2007.
[2] Russell Davidson and James G MacKinnon. Regression methods for using control variates in monte carlo experiments. Journal of Econometrics, 54(1-3):203-222, 1992.
[3] Chris Oates and Mark Girolami. Control functionals for quasi-monte carlo integration. In Artificial Intelligence and Statistics, pages 56-65. PMLR, 2016.
[4] Fred Hickernell, Christiane Lemieux, and Art B Owen. Control variates for quasi-monte carlo. Statistical Science, (0(1): - 31,, 2005.
[5] Roland Assaraf and Michel Caffarel. Zero-variance principle for monte carlo algorithms. Physical review letters, 83(23):4682, 1999.
[6] Antonietta Mira, Reza Solgi, and Daniele Imparato. Zero variance markov chain monte carlo for bayesian estimators. Statisties and Computing, 23:653-666 6625, 2013.
[7] Chris Oates, Mark Girolami, and Nicolas Chopin. Control functionals for monte carlo integration. Journal of the Royal Statistical Society. Series B (Statistical Methodology), pages 695-718, 2017.
[8] Chris J Oates, Jon Cockayne, François- Xavier Briol, and Mark Girolami. Convergence rates for a class of estimators based on stein's 's method. Bernoulli, 25(2): 141 - 1159, 2019.
[9] Leah South, CJ Oates, A Mira, and Drovandi. Regularised zero-variance control variates. arXiv preprint
2018.
10

---

[10] David Holzmüller and Francis Bach. Convergence rates for non-log-conese sampling and log-partition estimation. arXiv preprint arXiv:2303.03227, 2023.
[11] Robert A Adams and John JF Fournier. Sobolev spaces. Elsevier, 2003.
[12] Anastasios N Angelopoulos, Stephen Bates, Clara Fannjiang, Michael Jordan, and Tijana Zrnic. Prediction-powered inference. arXiv preprini arXiv:2301.09633 2023.
[13] Yaniv Romano, Evan Patterson, and Emmanuel Candes. Conformalized quantile regression. Advances in neural information processing systems, 32, 2019.
[14] Raphael A Meyer, Cameron Musco, Christopher Musco, and David P Woodruff. Hutch++: Optimal stochastic trace estimation. In Symposium on Simplicity Algorithms (SOSA), pages 142-155. SIAM, 2021.
[15] Lin Lin. Randomized estimation spectral densities of large matrices made accurate. Nu- merische Mathematik, 136: 13:1183-213, 2017.
[16] Aleksandros Sobczyk and Mathieu Luisier. Approximate euclidean lengths and distances beyond johnson- paansoon-indenstrauss, arXiv preprint arXiv:2205.12307. 2022.
[17] Hanzhong Liu and Yuehan Yang. Regression-adjusted average treatment effect estimates in stratified randomized experiments. Biometrika, 107(4))935-94)) 2020.
[18] Thomas Müllrr, Fabrice Rousselle, Alexander Keller, and Jan Novák. Neural control variates. ACM Transactions on Graphics (TOG), 39(6): 2020.
[19] Denis Belomestny, Artur Goldman, Alexey Naumov, and Sergey Samsonov. Theoretical guarantees for neural control variates in meme. arXiv preprint 2rXiv:2304.01111, 2023.
[20] Jiaxin Shi, Yuhao Zhou, Jessica Hwang, Michalis Titsias, and Lester Mackey. Gradient estimation with discrete stein operators. Advances in Neural Information Processing Systems, 15:25829-28841 2022.
[21] Hao Liu, Yihao Feng, Yi Mao, Dengyong Zhou, Jian Peng, and Qiang Liu. Action-depedent control variates for policy optimization via stein's identity. arXiv preprint arXiv: 2017.
[22] Jonathan W Siegel and Jinchao Xu. High-order approximation rates for shallow neural networks with cosine and relu activation functions. Applied and Computational Harmonic Analysis, 58:1-26, 2022.
[23] Anthony O'Hagan.
hermite quadrature. Journal of statistical planning and inference, 29(3):245-260, 1991.
[24] Motonobu Kanagawa, Bharath K Sriperumbudur, and Kenji Fukumizu. Convergence guarantees for kernel-based quadrature rules in misspecified settings. Advances in Neural Information Processing Systems, 29, 2016.
[25] Francis Bach. On the equivalence between kernel quadrature rules and random feature expan- sions. The Journal of Machine Learning Research, 18(1):714-751, 2017.
[26] Toni Karvonen and Simo Sarkka. Fully symmetric kernel quadrature. SIAM Journal on Scientific Computing, 4002):667-AA200 2018.
[27] Motonobu Kanagawa and Philipp Hennig. Convergence guarantees for adaptive bayesian quadrature methods. Advances in Neural Information Processing Systems, 32, 2019.
[28] Ayoub Belhadji, Rémi Bardenet, and Pierre Chainais. Kernel quadrature with dpps. Advances in Neural Information Processing Systems, 32, 2019.
[29] Ayoub Belhadji. An analysis of ermakov-zolotukhin quadrature using kernels. Advances in Neural Information Processing Systems, 34:27278-27225 2021.

---

[30] Rémi Bardenet and Adrien Hardy. Monte carlo with determinantal point processes. Annals of Applied Probability, 2020.
[31] Guillaume Gautier, Rémi Bardenet, and Michal Valko. On two ways to use determinantal point processes for monte carlo integration. Advances in Neural Information Processing Systems, 32, 2019.
[32] Satoshi Hayakawa, Harald Oberhauser, and Terry Lyons. Positively weighted kernel quadrature via subsampling. arXiv preprint arXiv:2/07.09597, 2021.
[33] Satoshi Hayakawa, Harald Oberhauser, and Terry Lyons. Sampling-basee nyström approxima- tion and kernel quadrature. arXiv preprint (tt::2301.09517 2023.
[34] Yutian Chen, Max Welling, and Alex Smola. Super-samples from kernel herding, arXiv preprint 3203.34.2 2012.
[35] Simon Lacoste-Julien, Fredrik Lindsten, and Francis Bach. Sequential kernel herding: Frank- wolfe optimization for particle filtering. In Arlificial Intelligence and Statisticx, pages 544-552.
PMLR, 2015.
[36] Ferene Huszár and David Duvenaud. Optimally-weighted herding is bayesian quadrature. arXiv preprint arXiv: 1204. 1664, 2012.
[37] Wilson Ye Chen, Lester Mackey, Jackson Gorham, François-Xavier Briol, and Chris Oates. Stein points. In International Conference on Machine Learning, pages 853. PMLR, 2018.
[38] Raaz, Dwivedi and Lester Mackey. Kernel thinning. arXiv preprint arXiv:2 2021.
[39] Raaz Dwivedi and Lester Mackey. Generalized kernel thinning. 10.05:2100.0159,, 2021.
arXiv preprint
[40] Erich Novak. Deterministic and stochastic error bounds in numerical analysis, volume 1349. Springer, 2006.
[41] Joseph F Traub, GW Wasilkowski, H Wozniakowski, and Erich Novak. Information- complexity. SIAM Review, 36(3)5514-514, 1994.
[42] Novak and H Wozniakowski. Tractability of multivariate problems, volume i: Linear information, european math. Soc., Zürich, 2(3), 2008
[43] Erich Novak and Henryk Wožniakowski. Tractability of Multivariate Problems: Standard Information for Functionals, volume 2. European Mathematical Society, 2008.
[44] Nikolai Sergeevich Bakhvalov. On the approximate calculation of multiple integrals. Journal of Complexity, (4):55022 516, 2015.
[45] Aicke Hinrichs, Erich Novak, Mario Ullrich, and H Wožniakowski. The curse of dimensionality for numerical integration of smooth functions. Mathematics of Computation, 83(290):2853-
2863, 2014.
[46] Erich Novak. Some results on the complexity of numerical integration. Monte Carlo and Quasi-Monte Carlo Methods: MCQMC, Leuven, Belgium, April 2014, pages 161-183, 2016.
[47] Aicke Hinrichs, David Krieg, Erich Novak, Joscha Prochno, and Mario Ulirich. On the power of random information. Multivariate Algorithms and information-based complexity, 27:43-64, 2020.
[48] Aicke Hinrichs, David Krieg, Erich Novak, and Jan Vybíral. Lower bounds for integration and recovery in l2. Journal of Complexity, 72:101662, 2022.
[49] David Krieg and Mathias Sonnleitner. Random points are optimal for the approximation of sobolev functions. arXiv preprint arXiv:2009. 1275, 2020.
[50] David Krieg, Erich Novak, and Mathias Sonnleitner. Recovery of sobolev functions restricted to iid sampling, Mathematics of Computation, 91(338):2715-2733 2022.
12

---

[51] Lucien Birgé and Pascal Massart. Estimation of integral functionals of a density. The Annals of Statistics, 23(1):11-29, 1995.
[52] David L Donoho and Michael Nussbaum. Minimax quadratic estimation of quadratic func- tional. Journal of Complexity, 6(3):290-323, 1990.
[53] David I. Donoho. One-sided inference about functionals of a density. The Annals of Statistics, pages 1390-1420, 1988.
[54] David L Donoho and Richard C Liu. Geometrizing rates of convergence, ii. The Annals of Statistics, pages 533-667, 1991.
[55] David L Donoho and Richard C Liu. Geometrizing rates of convergence, iii. The Annals of Statistics, pages 668-701, 1991.
[56] James Robins, Lingling Li, Eric Tchergen, Aad van der Vaart, et al. Higher order influence functions and minimax estimation of nonlinear functionals. Probubility and statistics: essays in honor of Duvid A. Freedman, 2:335-421, 2008
[57] Jiantao Jiao, Kartik Venkat, Yanjun Han, and Tsachy Weissman. Minimax estimation of functionals of discrete distributions. IEEE Transactions on Information Theory, 61(5):2835- 2885, 2015.
[58] Akshay Krishnamurthy, Kirthevasan Kandasamy, Barnabas Poczos, and Larry Wasserman. Nonparametric estimation of renyi divergence and friends. In International Conference on Machine Learning, pages 919-927. PMLR, 2014.
[59] Peter Mathé. Random approximation of sobolev embeddings. Journal of Complexity, 7(3):261- 281, 1991.
[60] Stefan Heinrich. Randomized approximation of sobolev embeddings, ii Journal of Complexity, 25(5):455-472, 2009.
[61] Stefan Heinrich. Randomized approximation of sobolev embeddings, iii. Journal of Complexity, 25(5):473-507 2009.
[62] Yanjun Han, Jiantao Jiao, and Rajarshi Mukherjee. On estimation of I, 1,-norms in gaussian white noise models. Probability Theory and Related Fields, 177(3-4): 1243 1003)1223-99 1294, 2020.
[63] Oleg Lepski, Arkady Nemirovski, and Vladimir Spokoiny: On estimation of the I, norm of regression function. Probability theory and related fielss, 11::221-253, 1999.
[64] Stefan Heinrich. On the complexity of computing the l norm. Journal of Complexity, 49:1-26, 2018.
[65] Alexandre B Tsybakov. Introduction to nonparametrie estimation, 2009. URL https://doi. wwg/10.1007/b13799 Revised and extended from the, 9(10), 2004.
[66] Holger Wendland. Local polynomial reproduction and moving least squares approximation. IMA Journal of Numerical Analysis, 21(1):285-30) 300. 2001.
[67] Holger Wendland. Scattered data approximation, volume 17. Cambridge university press, 2004.
[68] Yanjun Han, Jiantao Jiao, Tsachy Weissman, and Yihong Wu. Optimal rates of entropy estimation over lipschitz balls. The Annals of Statistics, 48(6):3228 - 3250, 2020.
[69] Gérard Biau and Lue Devroye. Lectures on the nearest neighbor method, volume 246. Springer, 2015.
[70] Aleksandr Reznikov and Edward B Saff. The covering radius of randomly distributed points on a manifold. International Mathematics Research Notices, 2016(19):6065--6994 2016.
[71] Rafat Latala. Estimation of moments of sums of independent real random variables. The Annals of Probability, 25(3):1502-151: 1997.
[72] losif Pinelis. Optimum bounds on moments of sums of independent random vectors. Siberian Adv. Math, 5(3):141-150, 1995.
13


Reviews:
[{"review_id": "HAhjY0G0xA", "reviewer": "Reviewer_pkMo", "summary": "In this papers the authors consider estimating the $q^{\\rm th}$-moment of a function $f$ (i.e. $\\int f^q(x)dx$) by observing samples of $x_i,f(x_i)$ when $x_i$-essentially follows a uniform distribution over (a compact) domain of integration. The paper first develops an information theoretic lower bound (using standard testing arguments) for estimating $\\int f^q(x)dx$ for $f\\in W^{s,p}$ (Sobolev Space of smoothness $s\\in \\mathbb{N}$ and order $p\\geq 1$). In order to match the lower bounds the authors subsequently specialize the results in two domains -- high and low smoothness respectively. In the high smoothness regime, a suitably bias corrected estimator based on an initial ML based estimator of $f$ (that satisfies some desirable initial properties) is shown to be optimal whereas in low smoothness regimes, owing to existence of rare and extreme events due to unboundedness of the $2q^{\\rm th}$-moment of $f$, the problem needs to addressed differently and the authors show that a truncated version of the classical Monte Carlo method can provide optimality guarantees. Finally the paper also provides some details on the noisy version of the problem (i.e. when one observes $x_i,f(x_i)+\\varepsilon_i$, $\\varepsilon_i\\sim N(0,\\sigma_n^2)$) to obtain connections to the nonparametric regression literature.", "questions": "If one is allowed to choose the design points carefully from a suitable class of distributions, can the authors discuss how the minimax rates change based on the class of distributions of $x$ -- indeed for uniform distribution one understands the problem and the problem does not change much for compactly supported (known/sufficiently regular) densities of $x$ which are bounded away from $0$ on its compact support. However, when one goes beyond this class, some properties might/might not change (e.g. for entropy estimation in density models it makes a fundamental density when the density is bounded away from $0$ versus when its not) since not having uniform distribution like samples in all sub-domains renders the integration complexity different. ", "limitations": "None noted.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "(i) A complete paper with minimax upper and lower bounds for estimating $q^{\\rm th}$-moment of a function $f$ (i.e. $\\int f^q(x)dx$) by observing samples of $x_i,f(x_i)$.\n\n(ii) Considers both high and low smoothness regimes of the problem.", "weaknesses": "Not much discussion on a a very well developed non-parametric theory of estimating $\\int f^q(x)dx$) by observing samples of $x_i,f(x_i)+\\epsilon_i,\\epsilon_i\\sim N(0,\\sigma^2)$ where both $1/\\sqrt{n}$ and (slower than)non-$1/\\sqrt{n}$-rates along with efficiency bounds  $1/\\sqrt{n}$ are presented. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely thank the reviewer for the insightful feedback on our work. Below are our responses to the questions raised in the review:\n\n1. Relation between our work and classical non-parametric theory\n\nRegarding a missing discussion on the classical non-parametric theory of estimating functionals via noisy data samples, we would like to mention that previous work on non-parametric estimation of functionals have been cited and discussed in line 97-105 of the paper. Below we also provide a more detailed discussion on the relation between our work and previous work on the estimation of nonlinear functionals (moments). \n\nWe remark that our work is the first to consider the case when the underlying function $f$ is in the Sobolev space $W^{s,p}(\\Omega)$, where $s < \\frac{d}{p}$. This implies that $f$ is not embedded in $L^{\\infty}(\\Omega)$ and may have a spike. In contrast, previous work all assumed that the underlying function is sufficiently smooth. \n\n \nSecondly, the technique used in our proof of the information-theoretic lower bounds is a bit different from that of previous work [1,2]. In one case of our proof, we pick the two priors to be two discrete distributions of functions such that any function’s sign on each grid in the divided domain is a discrete random variable whose distribution depends on $n$ (the number of data points). This enables us to calculate the amount of information required to distinguish between the two hypotheses even when there is no observational noise. We think that the technique we used here might be useful for the nonparametric statistics community. \n\nIn addition, our proof of the upper bound illustrates how the idea of doubly robust estimation can help us design estimators to achieve minimax optimality when there’s no observational noise. Under the other regime when observational noise exists, however, the behavior is completely different. As the influence function is not exactly zero in our estimator, we use Sobolev Embedding Theorem to embed the influence function $f^{q-1}$ and the error function $f-\\hat{f}$ into different $L^p$ spaces, which yields the desired upper bound. For the sake of completeness, we list a simplified version of the Sobolev Embedding Theorem below: \n\nFor any non-negative integer $s$ and $1 \\leq p < r \\leq \\infty$ satisfying $\\frac{1}{p}-\\frac{s}{d} = \\frac{1}{r}$, we have $W^{s,p}(\\Omega) \\subset L^r(\\Omega)$, $\\textit{i.e,}$ every sobolev function in $W^{s,p}$ is also a $L^r$ function.\n \n2. Design of points\n\nAs a response to the question on the design of points, we would like to point out that if the given density supported on the unit cube $\\Omega$ is both upper bounded and lower bounded by fixed constants, we can claim that the convergence rate remains unchanged and minimax optimal. However, the difficult part is that we don't have any prior information about the distribution from the given data points.\n\nMoreover, the design of points is in fact related to one important take home message emphasized in our work. Specifically, the minimax optimal convergence rate derived in our paper matches that of the Quasi-Monte Carlo method. From the perspective of information theory, our result reveals that a set of uniformly sampled points doesn’t lead to any loss in information compared to a set of carefully chosen points. For the Quasi-Monte Carlo method, most of the computational resources are used to locate the positions of desired points. After the points are picked, we may simply take the sample average of the function’s values on those points as our estimation. However, when we choose to use a set of uniformly sampled points, we no longer need to spend resources locating the desired points. The tradeoff is that our algorithm for estimation is more complicated, which leads to an extra computational cost. \n\nAs a side-note, we also provide some intuitive explanations that may help people understand the tradeoff between algorithms and design of points in Reproducing Kernel Hilbert Space (RKHS). Authors of [3,4] showed that one can compress the data points to present functions in RKHS with convergence rate $n^{-1}$. Once the data points are selected, one only need to take an expectation over the function values at the selected points to estimate the function's integral. In our algorithm, we first use half of the data to run kernel regression, which has a convergence rate of $n^{-\\frac{1}{2}}$. Then we use the Monte Carlo method to estimate the error between the original function and the result given by kernel regression. Thus, the convergence rate of our algorithm is also $n^{-1} =n^{-\\frac{1}{2}} * n^{-\\frac{1}{2}}$, where the first $n^{-\\frac{1}{2}}$ is the magnitude of the error induced by kernel regression and the second corresponds to the convergence rate of the Monte-Carlo Method. \n\nWe will take all the suggestions given by the reviewer into account while revising our manuscript. Finally, we would like to thank the reviewers again for their time and efforts, which help us improve the quality of our manuscript.\n\n[1] Han, Y., Jiao, J. and Mukherjee, R., 2020. On estimation of $l_r$-norms in Gaussian white noise models. Probability Theory and Related Fields, 177(3-4), pp.1243-1294.\n\n[2] Tsybakov, A.B., 2004. Introduction to nonparametric estimation, 2009. URL https://doi.org/10.1007/b13794.\n\n[3] Dwivedi, R. and Mackey, L., 2021. Kernel thinning. arXiv preprint arXiv:2105.05842.\n\n[4] Dwivedi, R. and Mackey, L., 2021. Generalized kernel thinning. arXiv preprint arXiv:2110.01593.\n"}, {"review_id": "v4wVZRI6Fk", "reviewer": "Reviewer_Qqae", "summary": "The authors study the theoretical properties of Monte Carlo estimators of the moments of a Sobolev function with nonparametric regression-adjusted control variates. In particular, they show that when a certain smoothness assumption is satisfied, then the regression-adjusted rule achieves the minimax optimal rate (where the minimum is over estimators, the maximum over test functions/integrands). When the smoothness assumption considered is not satisfied, the Monte Carlo estimate of the moments has infinite variance, and the authors show that to again achieve the minimax optimal rate one needs to resort to a truncated estimator. Further, the authors study the performance of the regression-adjusted estimator for integral estimation with *noisy* observations, a somewhat nonstandard setting, and provide some lower and upper bounds for the estimator's absolute error.", "questions": "- What are other valid choices beyond KNN for Section 4? \n- It would improve the paper to give example classes of $f$'s of interest", "limitations": "The authors should add a couple of sentences about limitations in the conclusions.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "- A good theoretical contribution to the literature on minimax-like results for Monte Carlo estimators \n- Good effort in making the notation and results accessible to a wider audience \n- Interesting setting with noisy function evaluations", "weaknesses": "- The paper reads a bit like a laundry list where results are presented in a sequence without much of a \"story\" connecting them; in particular, sections 3 and 4 feel disconnected, the former being about minimax results with noiseless observations specifically about moment estimation, the latter being about integral estimation (q=1 ? ) with noisy observations. Why does the setting change, what is the motivation? \n- Regarding your contribution, you mention \"existing proof techniques are only applicable in scenarios where there is either no noise or a constant level of noise\" - here are you talking about results for moment estimation or integral estimation ? Can you clarify in both settings which aspects are missing from the existing literature and how your results filll that gap ? \n- While I understand that the contribution is of theoretical nature, could you make more effort motivating why the setting of noisy observations is interesting / relevant to applications of Monte Carlo ? Further on this point, could you give intuitions for when the main smoothness assumption is expected to hold or not hold (beyond rare events) ?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We are really grateful for the reviewer's valuable feedback on our work. For questions raised in the review, we list our answers below:\n\n1. A coherent story connecting all the results\n \nWe already built a “story” connecting the results presented in the paper as the reviewer suggested. The story is to investigate whether a machine learning based control variate can help improve the Monte Carlo methods in terms of convergence rate.  We find out that if the random variable simulated via the Monte Carlo method is of finite variance, then we can always use a non-parametric control variate to improve the convergence rate. However, if rare events of infinite variance exist, a non-parametric control variate no longer helps to improve the convergence rate. This \"story\" is repeated several times in the title, abstract, introduction, and contribution sections. \n\n2. Our contributions in the paper\n\nWe've already listed all our contributions in section 1.2 (which is the \"story\" introduced before) in an intuitive way. We discuss our contributions in proof techniques here.\n\nFirstly, we emphasize that our paper is the first to consider the case when the underlying function $f$ is in the Sobolev space $W^{s,p}(\\Omega)$, where $s < \\frac{d}{p}$. In this regime, $f \\in W^{s,p}(\\Omega)$ may not be embedded in $L^{\\infty}(\\Omega)$. Different from literatures of doubly robust estimators, here we are considering a nonlinear functional, so the influence function is not exactly zero in our estimator. Therefore, we need to use Sobolev Embedding Theorem to embed the influence function $f^{q-1}$ and the error function $f-\\hat{f}$ into \"dual\" $L^p$ spaces to obtain the desired upper bound. For the sake of completeness, we list a simplified version of the Sobolev Embedding Theorem below: \n\nFor any non-negative integer $s$ and $1 \\leq p < r \\leq \\infty$ satisfying $\\frac{1}{p}-\\frac{s}{d} = \\frac{1}{r}$, we have $W^{s,p}(\\Omega) \\subset L^r(\\Omega)$, $\\textit{i.e,}$ every sobolev function in $W^{s,p}$ is also a $L^r$ function.\n \nSecondly, the technique used in our proof of the information-theoretic lower bounds is a bit different from that of previous work [1]. In one case of our proof, we pick the two priors to be two discrete distributions of functions such that any function’s sign on each grid in the divided domain is a discrete random variable whose distribution depends on $n$ (the number of data points). This enables us to calculate the amount of information required to distinguish between the two hypotheses even when there is no observational noise. We think that the technique we used here might be useful for the nonparametric statistics community. \n\nIn addition, for the problem of estimating integral based on noisy observations, our contribution is to characterize the minimax optimal convergence rate for any noise level ranging from zero to $O(1)$. \n\n3. Applications of our theoretical analysis\n\nWe already discussed related application in subsection 1.1 related work \"Regression-Adjusted Control Variate (RACV)\" of the paper. Here we will list the most important applications of RACV along with their references. These include gradient estimation [2], statistical inference in biology [3], causal inference [4], estimation of the normalization factor and sampling [5], MCMC simulation [6]. \n\nWe admit that our theoretical results are derived under a few assumptions (The given function is in a sobolev space; The data points are uniformly sampled from a unit cube; etc). We will discuss how these assumptions limit the application of our theory while revising the conclusion section. \n\nFinally, we would like to thank the reviewer again for all the helpful comments, which help us improve the quality of our manuscript.\n\n[1] Tsybakov, A.B., 2004. Introduction to nonparametric estimation, 2009. URL https://doi.org/10.1007/b13794.\n\n[2] Shi, J., Zhou, Y., Hwang, J., Titsias, M. and Mackey, L., 2022. Gradient estimation with discrete Stein operators. Advances in neural information processing systems, 35, pp.25829-25841.\n\n[3]  Angelopoulos, A.N., Bates, S., Fannjiang, C., Jordan, M.I. and Zrnic, T., 2023. Prediction-powered inference. arXiv preprint arXiv:2301.09633.\n\n[4] Liu, H. and Yang, Y., 2020. Regression-adjusted average treatment effect estimates in stratified randomized experiments. Biometrika, 107(4), pp.935-948.\n\n[5] Holzmüller, D. and Bach, F., 2023. Convergence rates for non-log-concave sampling and log-partition estimation. arXiv preprint arXiv:2303.03237.\n\n[6] Belomestny, D., Goldman, A., Naumov, A. and Samsonov, S., 2023. Theoretical guarantees for neural control variates in MCMC. arXiv preprint arXiv:2304.01111."}, {"review_id": "9Ohi8t88Tv", "reviewer": "Reviewer_gXt5", "summary": "This paper studies whether we can learn a control variate to reduce variance in Monte Carlo sampling. ", "questions": "Is there a simple and intuitive toy example to better understand problem and proposed solutions?", "limitations": "To someone who's not familiar with this area of research, the paper lacks introduction and is very confusing. ", "rating": 3, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "This paper studies whether we can learn a control variate to reduce variance in Monte Carlo sampling. ", "weaknesses": "To someone who's not familiar with this area of research, the paper lacks introduction and is very confusing. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": ""}, {"review_id": "14KZvxnyCF", "reviewer": "Reviewer_Z34k", "summary": "This paper presents theoretical results on estimates of the integral of f^q based, and when regression adjusted control variates can lead to better Monte Carlo estimators.", "questions": "The main thing that would improve my opinion of the work would be if the authors could add some details linking the theory to practice.\n\nThese were questions/comments I had when reading the paper that relate to the current presentation, and are suggestions that may make the paper more accessible to a wider audience: \n\nThe paper would be more accessible if an informal definition and intuition for a Sobolev space is given. \n\nI like the idea of something like Figure 1 to try and get some intuition about the dichotomy between situations where rare/extreme events are and are not possible. But at the moment I think this is of limited use — it just shows a curve with and without one. Could more be said, I.e. how does the Sobolev spaces in (a) and (b) differ (i.e. a has larger s perhaps?) and what does this mean about the functions (e.g. a is smoother than b). It may not be clear why smoothness as defined by the Sobolev space immediately means that you do not get extreme events as in the figure, as compared to assumptions on the derivatives of the function. Also the definition of extreme event seems to depend on q. Again some intuition here would help.\n\nPresumably there is some link in that if f is in a Sobolev space, then f^q is in a different space — and really what matters are the properties of f^q?\n\nWhat if Omega is unbounded (as you use a uniform distribution for the points in Theorem 2.1)? Or equivalently how does the distribution on the points affect the results in (2.1). At one point in the introduction you mention x_i uniform on the unit cube, but initially you just say that Omega is in R^d. I think things could be clearer.\n\nFigure 2 mentions truncated Monte Carlo — has this been defined at this stage (even briefly/informally)? \n\nSection 4: I think the introduction here could be clearer. That is stating that you now consider situations where you observe f(x_i) with error, and evaluate how the noise level impacts performance.", "limitations": "NA", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The results look to be novel and improve existing ones in this area. Furthermore, this is an area that is important in practice and for which strong theory can be impactful.\n\nI am not familiar with the historical work on this problem, and thus find it hard to say just how impressive the theoretical developments are -- thus the main focus of my review has been on their practical usefulness. It seemed that there is potential for these to be important, but (see Weaknesses) the presentation currently limits this.", "weaknesses": "The paper's presentation limits its accessibility to those who are very familiar with the theoretical area, which means that it is unlikely that the work would have impact on practice. As just one example, the introduction is written under the assumption that the reader knows what a Sobolev space is.\n\nSimilarly, there is little attempt to relate the theory to practice. E.g. how would a user know what Sobolev space their function is in? This may be straightforward in some situations (but is not explained). In general settings where we want to estimate E_pi(f(X)) and area using MCMC to simulate from pi -- things are complicated if the domain of X is unbounded, as I believe the theoretical results assume require to transform X to e.g. lie on the unit hyper-cube. The expectation is still an integral of a function over the hyper-cube, but perhaps it is less clear what the function is/what its properties are.\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We would like to thank the reviewer for the feedback on our work. Below are our responses to the questions raised in the review:\n\n1. Connection between application and our theoretical work\n\nWe already discussed related application in subsection 1.1 \"Regression-Adjusted Control Variate (RACV)\" of the paper. Here we will list the most important applications of RACV along with their references. These include gradient estimation [1], statistical inference in biology [2], causal inference [3], estimation of the normalization factor and sampling [4], MCMC simulation [5]. \n\nMoreover, we would like to emphasize that our main contribution is to provide a theoretical understanding of RACV. We examine all possible cases and pinpoint the regimes when RACV can help us obtain a minimax optimal estimator. We will make sure to clarify our main contributions in our next version. \n\nWe admit that our theoretical results are derived under a few assumptions (The given function is in a sobolev space; The data points are uniformly sampled from a unit cube; etc). We will discuss how these assumptions limit the application of our theory while revising the conclusion section. \n\n2. Intuitive Explanation of mathematical concepts (Sobolev spaces and embedding)\n\nWe have given the definition of Sobolev space in section 1.3. Intuitively, one may interpret $W^{s,p}(\\Omega)$ as the space formed by all the functions $f$ satisfies the following two conditions: (1) $f$ can be differentiated up to $s$ times; (2) The $i$-th derivative of $f$ has finite $L^p$ norm for any $0 \\leq i \\leq s$. \n\nMoreover, we also list a simplified version of the Sobolev Embedding Theorem below: \n\nFor any non-negative integer $s$ and $1 \\leq p < r \\leq \\infty$ satisfying $\\frac{1}{p}-\\frac{s}{d} = \\frac{1}{r}$, we have $W^{s,p}(\\Omega) \\subset L^r(\\Omega)$, $\\textit{i.e,}$ every sobolev function in $W^{s,p}$ is also a $L^r$ function.\n\nNow let's consider dividing the unit cube $\\Omega$ into $n$ grids, each of which has side length $n^{-1/d}$. We will then use a bump function supported on one grid  $\\Omega'$ as an example to provide some intuition on the connection between existence of rare events and the smoothness parameter $s$. For simplicity, here we just take $s=\\frac{d}{p}$ to be the threshold. Our bump function is given by\n\n$$\ng(x) = n^{(-\\frac{s}{d}+\\frac{1}{p})}K(n^{\\frac{1}{d}}(x-c))\n$$\n\nwhere $c$ denotes center of the grid and $K(y) :=  \\prod_{i=1}^{d}\\exp(-(1-y_{i}^2)^{-1})$ is compactly supported on $[-1,1]^d$. We will then proceed to verify that  $g \\in W^{s,p}(\\Omega)$. In fact, for any $|t| \\leq s$, we may use change of variable to directly compute the $L^p$ norm $\\|D^{t}g\\|_{L^p}$ of the $t$-th derivative $D^{t}g$, which is given by:\n\n$$n^{\\frac{|t|-s}{d}} \\|D^{t}K\\|_{L^p}$$\n\nNow let's consider different values of $s$. On the one hand, when $s >\\frac{d}{p}$, we have that $f \\in L^{\\infty}(\\Omega)$ for any Sobolev function $f \\in W^{s,p}(\\Omega)$ via the Sobolev Embedding Theorem above. This corresponds to Part (a) of Figure 1 in our paper, where the function is uniformly bounded without any extreme event. On the other hand, when $s <\\frac{d}{p}$, from the bump function above we have that the power $-\\frac{s}{d}+\\frac{1}{p}$ in the leading coefficient $n^{-\\frac{s}{d}+\\frac{1}{p}}$ is positive. As the bump function $g$ constructed above lies in the Sobolev space $W^{s,p}(\\Omega)$ for any $n$, we may take the limit $n \\rightarrow \\infty$ to deduce that a peak (rare event) can exist in this bump function, which corresponds to Part (b) of Figure 1 in our paper. Essentially speaking, Part (a) corresponds to $W^{s,p}$ with larger $s$, while Part (b) corresponds to $W^{s,p}$ with smaller $s$.\n\nFurthermore, we will explain how Sobolev Embedding Theorem above can help us find the space that $f^q$ lies in. Let $r = (\\frac{1}{p}-\\frac{s}{d})^{-1}$. Then we have $f \\in L^{r}(\\Omega)$, which implies that $f^q$ is a $L^{\\frac{r}{q}}$ function. When proving that our estimator is minimax optimal, our main strategy is to embed the influence function $f^{q-1}$ and the error function $f-\\hat{f}$ into \"dual\" $L^p$ spaces. Hence, it's crucial to figure out how to use some $\\hat{f}$ to approximate the original $f$, which is also the main technique contribution we made in the paper.\n\n3. Miscellaneous questions \n\nThe phrase \"Truncated Monte Carlo\" is defined in Subsection 3.2 of our paper. Regarding the introduction of Section 4, we will follow your advice to revise it in the new version.\n\nWe will take all suggestions into account while revising our manuscript. Finally, we would like to express our gratitude to the reviewer for your time and dedication, which helps us improve the quality of our manuscript.\n\n[1] Shi, J., Zhou, Y., Hwang, J., Titsias, M. and Mackey, L., 2022. Gradient estimation with discrete Stein operators. Advances in neural information processing systems, 35, pp.25829-25841.\n\n[2] Angelopoulos, A.N., Bates, S., Fannjiang, C., Jordan, M.I. and Zrnic, T., 2023. Prediction-powered inference. arXiv preprint arXiv:2301.09633.\n\n[3] Liu, H. and Yang, Y., 2020. Regression-adjusted average treatment effect estimates in stratified randomized experiments. Biometrika, 107(4), pp.935-948.\n\n[4] Holzmüller, D. and Bach, F., 2023. Convergence rates for non-log-concave sampling and log-partition estimation. arXiv preprint arXiv:2303.03237.\n\n[5] Belomestny, D., Goldman, A., Naumov, A. and Samsonov, S., 2023. Theoretical guarantees for neural control variates in MCMC. arXiv preprint arXiv:2304.01111."}]

Author rebuttal:

Meta-review: This paper considers estimating moments of a Sobolev function with control variates adjusted by regression. It shows that the resulting estimator satisfies an information theoretic lower-bound, and thus is in a sense unimprovable. Reviewers agreed the paper makes a useful contribution. However, there were also many comments from reviewers on presentation and "story" that might make the paper more accessible and increase impact. The authors are encouraged to take advantage of these.
Decision: Accept (poster)