A Quadrature Rule combining Control Variates and Adaptive Importance Sampling
Authors: Rémi Leluc
François Portier
Johan Segers
Aigerim Zhuman
Venue: NeurIPS 2022 Accept (2022)
Keywords: Monte Carlo
quadrature rule
control variates
adaptive importance sampling
Bayesian inference

Abstract: Driven by several successful applications such as in stochastic gradient descent or in Bayesian computation, control variates have become a major tool for Monte Carlo integration. However, standard methods do not allow the distribution of the particles to evolve during the algorithm, as is the case in  sequential simulation methods. Within the standard adaptive importance sampling framework, a simple weighted least squares approach is proposed to improve the procedure with control variates. The procedure takes the form of a quadrature rule with adapted quadrature weights to reflect the information brought in by the control variates. The quadrature points and weights do not depend on the integrand, a computational advantage in case of multiple integrands. Moreover, the target density needs to be known only up to a multiplicative constant. Our main result is a non-asymptotic bound on the probabilistic error of the procedure. The bound proves that for improving the estimate's accuracy, the benefits from adaptive importance sampling and control variates can be combined. The good behavior of the method is illustrated empirically on synthetic examples and real-world data for Bayesian linear regression.

A Quadrature Rule combining Control Variates and Adaptive Importance Sampling
Rémi Lelue
François Portier CREST
LTCI, Télécom Paris
Institut Polytechnique de Paris, France
ENSAL France
remi Jeluc@@teleccompprrbee
francois ponaccis.portiir@maall....o
Aigerim Zhuman
Johan Segers LIDAM, ISBA
LIDAM, ISBA
UCLouvain, Belgium aige digerim.m. zhuman@uclouvain. aiggrrin..huunnneecc.ovvv.a be
UCLouvain, Belgium johan. wwraansseerr@uuloaaai.....
Abstract
Driven by several successful applications such in stochastic gradient descent or in Bayesian computation, control variates have become a major tool for Monte Carlo integration. However, standard methods do not allow the distribution of the particles to evolve during the algorithm, as is the case in sequential simulation methods. Within the standard adaptive importance sampling framework, simple weighted least squares approach is proposed to improve the procedure with control variates. The procedure takes the form of a quadrature rule with adapted quadrature weights to reflect the information brought in by the control variates. The quadrature points and weights do not depend on the integrand, computational advantage in case of multiple integrands. Moreover, the target density needs to be known only up a multiplicative constant. Our main result is non-asymptotic bound on the prob- abilistic error of the procedure. The bound proves that for improving the estimate's accuracy, the benefits from adaptive importance sampling and control variates can be combined. The good behavior of the method is illustrated empirically on synthetic examples and real-world data for Bayesian linear regression.
Introduction
In recent years, sequential simulation has emerged a leading approach to compute multidimensional integrals, A key object in sequential simulation is the sequence of distributions, called the policy, from which generate the random variables, called particles, used approximate the integrals of interest. The policy is designed to evolve in the course of the algorithm to mimic the target density, which may itself be known only up to proportionality constant. While the design of algorithms with adaptive policies has been of major interest recently, only few studies have focused on using control variates to reduce the variance. This paper provides new method to incorporate control variates within standard sequential algorithms. The proposed approach significantly improves the accuracy of the initial algorithm, both theoretically and in practice.
The sequential framework. Consider the problem of approximating the integral f gf dA Fau (x) dr, where À is the d-dimensional Lebesgue measure, is a probability density on Rd and the integrand g is real-valued function on gd For instance, one may think of as the posterior density in Bayesian inference. Let (q:) 120 be the policy of the algorithm, i.c., sequence of proba- bility densities which evolves adaptively depending on previous outcomes. The particles (X,)221 are generated sequentially- at iteration i, particle X, is drawn from 9:-1. The integral f of dA is
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

estimated by the normalized sum (Ei=1 wig(Xi))/(Eit ws), where wi = f( (Xx() are the importance weights. The normalization En= Wi allows to deal with situations where the target density f is known only up to a proportionality constant.
Such an algorithm is part of the adaptive importance sampling (AIS) framework. Many different ways have been investigated to update the densities qi adaptively, Early works that inspired such sequential schemes include [13, 22, 29] where the sampling policy is chosen out of a parametric family. The parametric approach has been further extended by the Population Monte Carlo framework [4, 5. 26]. Various asymptotic results have been obtained in [6. 10. 34]. In [7, 9, 23, 391, nonparametric importance sampling based on kernel smoothing is studied. The latter bears resemblance sequential Monte Carlo methods [8, 6]. in which the target distribution f changes in the course of the algorithm.
Control variates. Let h = each k the integral J f dA is known. Without loss of generality, suppose that / h dA = 0. The functions ha are called control variates and can be obtained in different ways. In Bayesian statisties, Stein control variates [28] are constructed by applying the second-order Stein operator to functions satisfying certain regularity conditions [27]. Other control variates might be created by re-weighting a function h' that satisfies f hn dA- = 0 via h = /f The use of control variates is a well studied technique [14, 30J. The benefits can be established theoretically in terms of error bounds [28, 24]. weak convergence [35], the excess risk [2] and even uniform error bounds over large classes of integrands [33]. In practice, the control variates framework has led to efficient procedures in reinforcement learning [18, 25] and optimization [38], to name a few. Importance sampling and control variates in case of a Gaussian target density explored in [20]. The procedure in [21] incorporates control variates and is said to involve adaptive importance sampling, but in fact
be vector of real-valued functions on Rd such that for
the particles are always sampled from the uniform distribution on the unit cube. To the best of our knowledge, the existing control variate methods do not account for sequential changes in the particle distribution as the case in AIS.
AISCV estimate. The proposed approach to use control variates within the sequential AIS framework
relies on the ordinary least squares expression of control variates (see for instance [35]). To take care of the policy changes, some re-weighting must be applied. The AISCV estimate of the integral gf dA is defined as the first coordinate of the solution to the weighted least squares problem
n
(ōn, Bn) NER.EE RM
arg min
Σ We (g(X.) - a- b -- b(x))))
with wi the importance weights from before. The AISCV estimate ån has several interesting properties: (a) whenever is of the form a + h for some a E IR and 3 E R the error is zero, i... ån = a f gf dA; (b) the estimate takes the form of a quadrature rule å, = Σi= (x,,)(X,).. for quadrature weights tn that do not depend on the function g and that can be computed by single weighted least squares procedure; and (c) it can be computed even when f is known only up to a multiplicative constant. Point (a) suggests that when the linear combinations of the functions hę span a rich function class, the integration error is likely to be small. Point (b) implies that multiple integrals can be computed just easily Bayesian computations. In addition, the control variates can be brought into play in a post-hoc scheme, after generation of the particles and importance weights, and this for any AIS algorithm.
h
single one. Point (c) shows that the approach is applicable for
Main result. The main theoretical result of the paper is a probabilistic, non- con-asymptotic bound on ån a. Under appropriate conditions, the bound scales as T T/V where -² is the scale constant in a sub-Gaussian tail condition on the error variable E = g - a - ,3 h for (a, 3) = arg min,, I' bT hT)² dd. Note that has the smallest possible variance one could get using control variates h. As a consequence, when the space of control variates is well suited for approximating g. the AISCV estimate will be highly accurate. Also, our bound depends only the linear function space spanned by the control variates /i,..., hmm not on the particular basis chosen in that space. The results rely on martingale theory, in particular on a concentration inequality for norm-subGaussian martingales in [19]. In the course of the proof, we develop a novel bound on the smallest eigenvalue of certain random matrices, extending inequality from [37] to the martingale case.
Outline. Section 2 introduces the general framework of adaptive importance sampling and control variates. Next, Section 3 presents the AISCV estimate and the associated quadrature rule. Section 4 contains the statements of the theoretical results while Section gathers practical considerations, including the construction of control variates. Numerical experiments are presented in Section 6.
2

---

2 Preliminaries on Monte Carlo integration
The aim of this section is to present the required mathematical framework for Monte Carlo integration and the variance reduction methods of interest, namely adaptive importance sampling and the control variate technique. Recall that g : R - R is an integrand and f probability density on Rd. The aim is to compute Eg = of dA.
daaptive importance sampling, In adaptive importance sampling (AIS), Eyl is estimated by a weighted mean over sample of random particles XI...., X, in R. Since appropriate sampling densities naturally depend on 9 and J. we generally cannot simulate from them. They are then approximated in an adaptive manner by family of tractable densities (g:):20 that often evolve towards a density lopt that optimizes some criterion. While the starting density do is fixed, the density qi for i 2 is determined in function of the particles X1 , X; already sampled; think for instance of a parametric family, where the parameter of qi is a function of XI,..., X,. Given the particles X1, XI,..,,', Xi, random vectors on I& defined on some probability space (92 $0 PP. The distribution of the sequence (Xt)121 is specified by its policy as defined below.
the next particle, Xi+1, is then drawn from q. Formally, let (Xi)121 be sequence of
Definition 1 (Policy). A policy is a random sequence of probability density functions (g: (9:).20 on Rd adapted to the o-field (Fi))zn defined by Fn = = (2) and P; = oX1,.... X;) for i 2 1. The sequence (g: lizo is the policy of (Xi)121 whenever X; has density qi-1 conditionally on F:
The (normalized) adaptive importance sampling estimate of [g is then defined as f(X,)
r(ais) (g) Ei= (Xg(X.) where
(x,) for i = ,,...,
(1)
Ei=1 W/
The sampling weights w; reflect the fact that X; has been sampled from qi rather than from f. The division by 2im w; rather than by n has two benefits: first, the integration is exact for constant integrands, and second, f needs be known only up a proportionality constant, an advantage for Bayesian inference.
Since updating the density qi at each iteration may be computationally expensive, it is customary to hold it fixed over a pre-determined number of iterations. Writing n = n2 + +... + ny in terms of positive integers (n=)==1 called the allocation policy, the AIS estimate then becomes
(nis) (g) Ei=1 E (X.s)
T
((x,s)
where We,i 2e (X,,)
(2)
E W3,5
nl. At stage I, the particles .e,,,..., Xr,,,, are sampled indepen- dently from Qt-1. while all particles sampled up to and including stage i are used to determine the sampling density qe for stage t + 1. It is easy to see that the two formulations of the AIS estimate are equivalent: (1) arises from (2) by setting n, = 1 for all t, while (2) can be obtained from (1) by constructing the policy in such a way that the densities q: do not change within integer intervals of the form (O...., My - }. (MI,... ...., , + n2 1}, and so on While the shorter representation (1) is more convenient for theoretical purposes, formulation (2) is the one used in practice (see Section 6).
for / - 1,..., and i
Interestingly, the AIS estimate (1) may be seen as a weighted least-squares estimate minimizing the loss function + E (g(X,) - - This perspective is key to understand control variates.
Control variates. The control variates method is a variance reduction technique that consists in incorporating a new piece of information- the known values of the integrals of some control functions-in basic Monte Carlo framework. Control variates are simply functions h,, ... Am Lz(f) with known integrals. Without loss of generality, assume that E [hy] = 0 for all j = 1..... m. Let h = (hj...., hm)² denote the -value function with the m control variates as elements. For any coefficient vector 3 E R.m, we have R, (g - ,3 M] = F [g]. Given an independent random sample XI,..., Xn, from f, any B E 1n therefore results in an unbiased estimator of Ey[g] by
(cv) (g. B) = - Σ 33(X,))3 h(x)))
Provided the covariance matrix G == [hh J is invertible, there is unique coefficient vector 3* E IR for which the variance of f(ov) (g) is minimal and it is given by
8 = (EE 1hh [) -1 Ey [hg).
(4)
3

---

This vector being generally unknown, it needs to be estimated from the particles .I..... Casting the problem in ordinary least squares framework leads to the control variate estimate
T((0i (g) I((v) (g. ((xr))
a(cv)
where
(5)
(a/o*) B(ev)) E arg min
Σ (o(X;) - -bb))(x))) x
The estimator (g) is well-defined provided the minimizer (ev) to (5) is unique. This is the case if and only if there does not exist E IR such that 6 h(X,) = for all I....,n.
The asymptotic distribution of 11xx) 3* were used in (3). In particular, the asymptotic variance I((x) (g) is o (g)/n where
(g) as 3o is the same as if the variance- variance-minnmizing vector
o (g) min E, [(g [ 3 hh
SER"
Interestingly, when using only the first out of control variates, where l E {0, 1,..., m). we have o (g) I (g). In terms of asymptotic variance, therefore never harms to add more control variates. Their construction will be addressed in Section 5.1.
3 Combining adaptive importance sampling with control variates
AISCV estimator. Consider the same integration problem Ey[g] = gf dA as Section With the idea of performing variance reduction when calculating integrals with respect to the posterior
density in Bayesian inference, we incorporate control variates into the AIS estimate. Let the particles (X1)121 be generated according vector of control variates, i.e., hj E L2()) and [y = 0 for every j = 1,..., m. Combining (1) and (3), the proposed estimate takes the form
policy (2))20 in Definition 1. Let h (hi,.... .... T be
E wi -
A (aisev) (g.
(6)
Σx W/
where B E R remains to be determined. To do so, the ordinary least-squares problem in (5) is replaced by a weighted one, yielding the novel AISCV estimator
(aiscv))
(g) = 11 ((aisov) (g, Pn) = ān where
n
(ān, Br) arg min (a.b)ERxR"
Σ wi (o(X,) - -a-b-(xx,)))
The estimator is well-defined only the minimizer &r is unique-the minimizer 3n need not be. We will come back to this in the paragraph.
As in (2), the policy may be divided into T stages in order to reduce the number of times the sampler needs to be updated. Stage / = 1,..., T has length nı, with Σ = n.. Within each stage, the sampling density remains constant. In practice, this leads to the AISCV estimate in Algorithm 1.
Quadrature rule. The AIS estimate (1) is a quadrature rule with quadrature points X, and quadrature weights proportional to the sampling weights w,. The AISCV estimate (7) has the same property, but with adapted quadrature weights. Let en = (e,,i) be the vector of residuals resulting from the weighted least-squares regression of the constant vector 1, 1)) E R" on the control variates but without intercept:
Cn,i =1 1- (1..) ((X,))
where
Br (1,) E arg min
Σ (m,(1 - (1-bT)(()))
(8)
DERm
Even though the vector Ba (1,) is not necessarily unique, the weighted least squares fit ((x.)) always According to the next proposition, the quadrature weights are proportional to (W)en.))=

---

Algorithm Adaptive Importance Sampling with Control Variates (AISCV)
Require: integrand J, target density f (up to a proportionality constant), number of stages 7 E No, allocation policy (n4) initial density qu- update rule for the sampling policy
1: for = I...., T do
2: Generate independent random sample Xt,1, Xe,n, from 1:-1
3: Compute the vector of weights (urr)) where W4, = f(2))// (Xe,)
4: Construct the matrix of control variates H2 = (X,(),,)) 1=
i=1,..,,
5: Evaluate the integrand in the particles: (g(X),)))
6: Update the sampler de based on all previous particles (X,,² : 1,... ,,...,,;; ..,...
7: end for
8: Compute (á.r. (r) arg min (RRRR""
{ E wiss (g(Xe,)) - -- -xT(()))))
9: return 7/ (aisev) (g) Örr.
Proposition 1 (AISCV quadrature rule). The minimizer on, in (7) is unique if and only en, # 0 in (8). In that case, the AISCV estimate is
(g) = ân Ei= wien i9(e
T(alsov) n
(9)
E WiCn,i
If en = 0, then there exists 6 R such that 6T (T,,)) for all L,..., n. In that case, the minimizer o, in (7) is not unique and the AISCV estimate is not well-defined. To remedy this, one can for instance reduce the number of control variates. This issue already occurs with the ordinary control variate estimator in (3).
Rather than requiring different weighted least squares problem for every integrand g as in (7), the quadrature rule in (9) only involves single weighted least squares problem (8), whatever Given the quadrature weights, calculating the AISCV estimate for novel integrand only requires the evaluations of that function on the sampled particles, making the whole procedure post-hoe scheme. The steps case the sampling policy is divided into 7" stages are given in Algorithm 2, which gives the same result as Algorithm 1, but with less effort if multiple integrands are into play.
Algorithm 2 Quadrature Rule - AISCV post-hoc scheme
Require: integrand g. T Nº allocation policy (no) weights (wy (w:) with wi (w,, 1 matrices (Ho with He =
-=I,.... particles (X,, : L .,, T; .1)
i=l...,mm
1: Compute 3. (1.,) arg min Σ 2t=1 Ei Whi (1 - 1
2: Compute u, diag(u (w) H, (1n) fort 1....,T
3: Compute = ΣT Ei1 U4
4: Compute weights U/₁i ut,i/8 s for 1 1,..., Tandi- i 1,...,me
5: return 1 ((assvv) (g) = Σ Ena (Xess)
Theoretical properties of the AISCV estimate
Here we point out several theoretical properties of the novel AISCV estimate. A first point is that the integration rule is exact on the linear span of the control variates and the constant function.
Proposition 2 (Exact integration). For integrands of the form + B h for E IR and B3 E R
the AISCV estimate exact: [ (g) a Eg
A second property is that we may apply arbitrary invertible linear transformations to the control variates without changing the AISCV estimate. This can be advantageous computationally, to make the underlying weighted least squares problem more stable numerically. Also, it means that without
5

---

loss of generality, we may assume that the control variates are uncorrelated and have unit variance, which simplifies the theoretical performance analysis.
Proposition 3 (Invariance). If the matrix A E Im is invertible, then the AISCV estimate based on the control variates Ah the same as the one based on h.
Our main result is a non-asymptote bound on the error of the AISCV estimate for J gf dA when J g dA is finite. First, we introduce some assumptions and definitions.
The first condition that is required concerns the policy given by the AIS part of the algorithm. It is supposed that any element from the policy should dominate the function J.
Assumption (Dominated measures), There exists c 2 1 such that, for all x E Rd and for any
i i= 1, ni, we have fxx S c cx))
This assumption represents a safe approach to importance sampling, as the policy will always allow to sample in places where f is positive. A well-known and well-spread [16, 30, 9] technique to achieve such a defensive strategy is to a use mixture density q: (1 n) + 1/00 where 1) (0, 1) and where has sufficiently heavy tails to dominate f. Such a mixture allows to choose the densities fi with some flexibility using in principle any AIS algorithm. Second, the control variates shall be linearly independent and bounded.
Assumption 2 (Control variates). We have supx:(()>>0 [b)(z) co for all j Iy..., m. The matrix G = f hhh T f dd is invertible.
The previous condition allows to define the standardized vector of control variates as 1- By Proposition 3, this change does not affect the AISCV estimate. The orthonormal control variates will play a key role through the following quantity
B sup [|xx))]1
2=
The quadratic form |m<=)] Tx=TTT->(()) is referred as the leverage function in ordinary linear regression as it quantifies the influence of a training point x on the prediction of the observed response. It is invariant with respect to invertible linear transformations of the control variate vector.
Assumption and the fact that the integrand g is square integrable with respect to f allows to define the residual function = 9 J of dλ - h 3 where 3 has been introduced in (4) as minimizer of the residual variance. Since we work in the space [J)()), we assume without loss of generality that g and h vanish outside {x : fzx) > 0} and we put =(x)= = 0 for x E Fd such that f(r) = 0. The residual function should satisfy the following tail condition.
Assumption 3 (Residual tail). There exists T > 0 such that, for all 1 > 0 and all integer i 2 1, we have PP P[[wee (X:)) 1 S (exp((((((22))
The previous assumption concerns both the function e and the policy sequence (gi) Since EX,, F.-I = 0, it is implied by the so-called sub-Gaussian condition [3] that Elexp(Auu;e (X,)) [-,] I exp(-X?
for any À R. In the proof of Theorem 1, As- sumption 3 allows to derive concentration bounds on eisidual-based sums using recent results from [19, 24]. We are now in position to state our main result on the error of the AISCV estimate.
Theorem (Concentration inequality for AISCV estimate). If Assumptions 1, 2 and 3 hold, then, for any d (0, 1) and for all 2 Cic² log (10m/0), we have, with probability at least 1 - 8, that
[ g(x)f(z) dx S CCTV / Rd
[oo (10/8)
+ C3CBB CCc log(10m/5))
r(aisev)(g) n
(g)
where C1. C2. C2 are universal constants specified in the proof.
Remark 1 (Understanding 7). The quantity in Assumption 3 is related to the conditional variance They actually coincide when wea(X.) is Gaussian. For a policy satisfying Assumption 1, Elwfe (X,) _ S co which for certain combinations of integrands and control functions scales as m s/d [35] where the parameter s represents the degree of smoothness of g.
Elwiee (X.)) F.-i,,
Remark 2 (Convergence rates). Consider an asymptotic regime where the number of control variates m lends to infinity with the sumple size The AISCV estimate improves upon the AIS method (m = 0), which has rate 1/ /T, as soon as T + T log(m)/ - 0. To recover the same order of an oracle estimate with rate 7 / nn, one must have B log(m) O(vn) as n CO.

---

5 Practical considerations
This section presents ways to build control variates using either families of polynomials general functions based on Stein's method, with highlight on computations in the Bayesian framework.
5.1 Control variate constructions
Orthogonal polynomials. When the target density f can be decomposed as a product of univariate densities f = P1 @ .... @ Pd, multidimensional control functions may be constructed based on univariate ones. This happens for instance for the uniform distribution over the unit cube [o, 1] or with uncorrelated Gaussian distributions on I Such univariate control variates may be easily constructed using families of polynomials [12]. such as Legendre polynomials for the uniform distribution on [0. 1] and Hermite polynomials for the Gaussian distribution on R. This technique can also be used when f is dominated by another density f having the said product form by transforming zero-mean control variates h" with respect to f via h = h ff.
Let (h,, ................................ hx)) be vector of univariate control functions with respect to density p. i.e., 1yp] 0 for all j 1, 1....,,k k. Let ho = 1 denote the constant function equal to one. For a multi-index == (6,---,, fa) in {0..... ......... kjd a {(0...... 0)), multivariate controls with respect to pdd are built by forming tensor products of the form (1 Fa) = he, (#1) ... hea (xa)) yielding a total number of m = (k: + 1)d - 1 control functions. Alternative approaches yielding smaller control spaces consist of imposing L 0 for all but a small number of coordinates j 1,.... d by the constraint & + .....42 + Ea I Q for some Q 2 1.
Stein control variates. In the general case where one has only access to the evaluations of ʃ. control variates may be constructed using Stein's ss method. The technique relies on the gradient V, log f(x) which can either be directly computed (see the example of Bayesian regression below) or which may be available through automatic differentiation provided in popular API's such as Tensorflow and Py [1, 31J. Let A V V. denote the Laplace operator. By definition, the second-order Stein operator L [36, 15] associated to the density f is defined by:
Vy E & (RR,,,,
(L)) = -n<(()) + V V><(<) v log fxx)
VV,
The transformation guarantees that [ 0 for all op with weak regularity conditions [27]. There- fore, we can build infinitely many control variates h2 = Ly from given functions s. One simple way is to let yp be polynomial with bounded total degree: for a degree vector (a1..... aa) E Nd with a1 + ...... + ad S Q. define (a(x) = =i ...x." Given the dimension d and the total degree Q, there are m (d= 1 such degree vectors, yielding the associated control variates ha hep. For fast computation, note that, writing (x) = (x)
aa) and D2(x) = diag(an 1)/aq,.... .alaa - 1)/=2), we have Vasa (x) Dí (z))aa (x) and Arta (I) I (D)()))a In practice, all combinations of a are stored in matrix A E Nmxd
(r) L. (1 = diag(on/21
5.2 Bayesian inference
Given data Ø and a parameter of interest 8 E θ C Fd posterior integrals take the form Jae (())p(() ( Ø) dé, where p((]9) & (()) is the posterior distribution. proportional to prior π (0) and a likelihood function (())). For instance, when g(0) = 0, the integral above recovers the posterior mean. Stein control variates involve the computation of the gradient of the log- posterior v, log 2())2), which implicitly relies on the score function V log ((2)0). We point out two common examples linear and logistic regression- where these functions are easy to compute.
Bayesian linear regression. Consider a linear regression problem comprised of observations X R with labels y RV.. In the Gaussian fixed design setting, the predictor z; produces the response y/1 = xz. 0 + & where E1,...,EN, ~ (0,, o²) are centered Gaussian noises. The likelihood (X,y(0) is proportional to (a2) -N/2 (xp(-(y- X0) (y- - (2))() yielding score function V log e(X) = X (y X2))(2o²).
yielding the
Bayesian logistic regression. Next, consider the logistic regression problem comprised of observa- tions X E N xd with associated binary labels y E {0, 1}.V. Letting (s)) = 1/(1 + e- denote the sigmoid function, the likelihood function is E (() = II,, (0 );;) (1 a00 =;)))-u, The score function is simply Vo log ((X,y0)) = X (y (X0)).
7

---

Numerical illustration
To compare the finite-sample performance of the AIS and AISCV estimators, we first present in Section 6.1 synthetic data examples involving the integration problem over the unit cube [0, I and then with respect to some Gaussian mixtures as in [4]. The goal is to compute f of dλ for vectors of integrands 9 R + Iº. We consider various dimensions > 1 and several choices for the number of control variates m. Section 6.2 deals with real-world datasets in the context of Bayesian inference. For ease of reproducibility, the code, numerical details and additional results are available in the supplementary material.
Parameters. all simulations, the sampling policy is taken within the family of multivariate Student & distributions of degree v denoted by {4,,E% μ E RR) with Eo Lala(v 2)/v and v > 2, og > O. Similarly to [34], the mean μ4, updated at each stage / -,.... by the generalized method of moments (GMM), leading to fı = (E+ Ei= W === Dn w. The allocation policy is fixed n, 1000 and the number of stages is T E {5; 10; 20; 30; 50}. The different Monte Carlo estimates are compared by their mean squared error (MSE) obtained over 100 independent replications.
6.1 Synthetic examples
Integration on [0, 112. We seek to integrate functions g with respect to the uniform density
f()) = 1 for x E [0, in dimensions d E {4:8}. We rely on Legendre polynomials for the control variates. Consider the integrands (xx) = 1 + sin E2 2; - 1))- 92(x) =
IT (2/x)\/2x, -e- log(2.)* -ebgtxi and ga(x) = IT log(2)21~>, all of which integrate to 1 on [0, 1| None of the integrands is linear combination of the control variates. The policy parameters are and σ 0.1. The control variates of Legendre polynomials where the degree l, equals 0 for all but two coordinates, leading to a total number of m kd + h2 d(d 1)/2 control variates. The maximum degree in each variable k = 6, 240 and m 1056 control variates in dimensions d and d 8 respectively. Figure presents the boxplots of the AIS and AISCV estimates. The error reduction obtained thanks the control variates is huge: the AISCV estimate has mean squared error smaller than the of the AIS estimate by factor at least 10 and up to 100 (see Table in the supplement).
No = (0.5.... 0.5) E Rd. y
built out of tensor products
yielding
1.06
AB
1.06
1.04
1.09
ABCV
AISCV
AISCV
1.07
AISCV
1.04
1.04
1.02
1.05
1.02
1.02
1.03
1.00
1.00
1.01
0.98
0.98
0.98
0.97
0.99
0.96
0.96
0.95
0.96
0.93
0.91
5
10
20
50
5
10
30
50
5
20
30
50
10
20
30 50
Total Sample Size (xle3)
Total Sample Size
Total Sample Size
Sample Size (x1e3)
(a) gn (d 4)
(b) gi (el 8)
(c) 92 4)
(d) gs (d = 8)
(g) and / (aisev) (g) with integrands 91 f2, yo in dimensions (4:8) obtained over 100 replications. The true integral equals
Integration on [0, 1] boxplots of estimates /(als)
Figure 1:
Gaussian mixture f and Stein control variates. In this setting we assume we only have access to the evaluations of the target density f. We consider the classical example introduced in [4] where f is a mixture two Gaussian distributions. The control variates are built using Stein's 's method (Section 5.1) out of polynomials of total degree at most Q {2; 3), leading number of control variates m E {14; 34} in dimension d = 4 and m E (44; 164) in dimension d = 8 respectively. We consider two cases: an isotropic and anisotropic one.
Isotropic case. Let (z) = 0.50 14) (x μ1) where (1, (11...,1)T/2/443 = Ia/d and Φ5; is the multivariate normal density function with zero mean and covariance matrix Σ. The Euclidean distance between the two mixture centers is 1, independently of d. The initial density go is
the multivariate Student distribution with mean (1, 11,11, ....................... .), 0) Vd and variance (5/d)Ia. The

---

initial mean value differs from the null vector to prevent the naive algorithm using the initial density from having good results due to the symmetrical set-up.
Anisotropic case. In this case, the mixture is unbalanced and each Gaussian is anisotropic. The target density is (z) = 0. v (x μ²) + v.25 (z + μt) where = (I....,1) /)T/2vd and V= = dia 11.,,) l).d. The initial density go is the same as for the isotropic case.
Figure 2 presents the evolution of the logarithm of the mean squared error [}(()) T(g)]]3 Once again, the AISCV estimators are the clear winners with a mean squared error smaller by factor up to 1000 for the anisotropic case (see Table in the supplement).
bo (ifo)- 100111
Log of MSE (11gg) 11g0113
ABDV2
NSSVV
Log ((1g])
Log MSE liño)
-10
-8-
A/S0/3
ASOO
-9
-8
-10
-12
-12
AISEV2
AISDV2
-13
-13
MBCC
AISCV3
5
10
20
30
50
5
10
20
30
50
10 20
30
50
10
20
30
50
lotal Sample Sizee Ixlex)
Total Sample size Odles)
Total Sample Size (3te3)
Total Sample Size (xle.))
(a) fx 4)
(b) fx = 8)
(c) (d 4)
(d) fv 8)
Figure 2: Gaussian mixture density: Logarithm of | - (g) for g(x) = x with target isotropic
fy and anisotropic fv in dimensions d E 8} obtained over 100 replications.
6.2 Real-world examples
We place ourselves in the framework of Bayesian linear regression (Section 5.2) with fea- tures X E R and continuous responses y E TRN. The posterior distribution p(0) ()) in- volves a Gaussian prior (0) ~ N Ex) and a likelihood function {(2(0) proportional to (o²)
N/2 (o/2 exp(-(y - TO) (v- - (20)( The noise level is fixed and taken sufficiently large at o = 50 to account for general priors. The posterior distribution is Gaussian too: N (Ps, Eg} with -,,))= Σバ')) and the control variates are built with the Stein operator (Section 5.1) out of monomials with total degree Q {1: leading to the AISCVI and AISCV2 estimators respectively.
H4 -x(a-~XXT y + Σ (aa) and Σ, 2 X +
The integrand is g(0) = - 63
Datasets and parameters. Classical datasets from [11] are considered housing (N 506; d 13; m. E (12; 104)); abalone (N = 4177; d = 8; m E {7:44)) 44)); red wine (N = 1599; d = 11; m E {10; 77)); and white wine (N = 4898; d = 11; m E {10; 77))) ). The initial density is the multivariate Student distribution with y = 10 degrees of freedom, zero mean and covariance matrix Es.
Boopot 11gihh00
0.15
AISCV1
Bonplot Ret- Highing)
0.30
AB
ABEVE
Boophett for (llg) High)/ll
0.15
AFOOV]
AIS
Bropritt for lig) lighttig)
0.15
0.20
MSCVI
0.10
AIS022
0.20
AISC22
AFFOV2
0.10
ASSC2
0.05
0.10
0.05
0.05
0.00
0.00
0.00
-0.05
-0.05
-0.55
320
0.20
-0.10
2013
-0.10
-0.0
201
.0.30
-0.15
-0.15
20
30
50
10
5
10
20
30
50
10
20
30
50
Total Sample Size (dle3)
Total Sample Size (cie3)
Total Sample Size (xte3)
Total Sampte Size Oxle3)
(a) Housing
(b) Abalone
(c) Red Wine
(d) White Wine
Figure 3: Bayesian linear regression: boxplots of (i(s) - I(g))//()) for g(0) = =0))=
Results. Figure 3 presents the boxplots of the relative error ((()))((g))/[())) revealing the benefits of control variates even with polynomials degree Q = 1. When Q = 2, the error of the AISCV2 estimate is virtually zero (see Table 3 in the supplement), in line with Proposition 2. The mean squared error of the AISCVI estimate smaller than that of the AIS estimate by a factor ranging between 2 and 10.

---

7 Discussion
While control variates are a well-known tool for Monte Carlo integration, standard methods do not allow the distribution of particles to evolve throughout the algorithm, as is the case for sequential methods. Within the standard adaptive importance sampling framework, we have developed weighted least-squares procedure to improve numerical integration by incorporating control variates. The underlying adapted weights of this quadrature rule do not depend on the integrand and our non-asymptotic bound highlights the benefits of combining adaptive importance sampling with control variates. Different ways for constructing control variates are proposed. The method is for computing integrals with respect to the posterior density in Bayesian analysis, as the target density only needs to be known up to a multiplicative constant.
A limitation of the combined AISCV approach is that it requires the user make quite some design choices, notably the sampling policy for the AIS part and the control variates for the CV part. These culminate into the factor in Assumption 3, which appears prominently in the error bound in Theorem and which can be interpreted roughly as the standard deviation of we, where w the importance weight well behaved when the policy is well-chosen in relation to the target density and where E is some residual function - well behaved when the control variates are well-chosen with respect the integrand. Further, choosing too many control variates may result in ill-conditioned empirical Gram matrix or in overfitting The least-squares solution could become unstable, requiring some kind of regularization, such the LASSO [24].
Acknowledgements
The authors grateful to the Area chair and three anonymous Reviewers for their valuable comments and interesting suggestions. Aigerim Zhuman gratefully acknowledges a research grant from the National Bank of Belgium.
References
[1] Abadi, M., Barham, P., Chen, J., Chen, Z., Davis, A., Dean, J., Devin, M., Ghemawat, S., Irving, G., Isard, M., et al. (2016). TensorFlow: A system for Large-Scale machine learning. In 12th USENIX symposium on operating systems design and implementation (OSDI 16), pages 265-283.
[2] Belomestny, D., Iosipoi, L., Paris, Q., and Zhivotovskiy, N. (2022). Empirical variance minimiza- tion with applications in variance reduction and optimal control. Bernoulli, 28(2):138--140)
[3] Boucheron, S., Lugosi, G., and Massart, P. (2013). Concentration Inequalities: A Nonasymptotic Theory of Independence. Oxford University Press.
[4] Cappé, O., Douc, R., Guillin, A., Marin, J.-M., and Robert, C. (2008). Adaptive importance sampling in general mixture classes. Statistics and Computing, 18(4):447-459,
[5] Cappé, O., Guillin, A., Marin, J.-M., and Robert, C. (2004). Population Monte Carlo. Journal of Computational and Graphical Statisties, 13(4):907-929,
[6] Chopin, N. (2004). Central limit theorem for sequential Monte Carlo methods and its application to bayesian inference. The Annals of Statistics, 328652233-
[7] Dai, B., He, N.. Dai, H., and Song. (2016). Provable Bayesian inference via particle mirror descent. In Artificial Intelligence and Statistics, pages 985-994, PMLR.
[8] Del Moral, P., Doucet, A., and Jasra, A. (2006). Sequential Monte Carlo samplers. Journal of the Royal Statistical Society. Series B (Statistical Methodology), 68(3):411
[9] Delyon, B. and Portier, F. (2021). Safe adaptive importance sampling: mixture approach. The Annals of Statistics, (912))885-917.
[10] Douc, R. and Moulines, F. (2008). Limit theorems for weighted samples with applications to sequential Monte Carlo methods. The Annals of Statisties, pages 2344-2376.
10

---

[11] Dua, D. and Graff, C. (2019). Uci Machine Learning Repository [http://acchie. ics. uci. edu/ml].
irvine, ca: University of california. School of Information and Computer Science, 25:27.
[12] Gautschi, W. (2004). Orthogonal polynomials: computation and approximation. OUP Oxford.
[13] Geweke, (1989). Bayesian inference in econometric models using Monte Carlo integration. Econometrica: Journal of the Econometric Society, pages 1317-1339.
[14) Glynn, P. W. and Szechtman, R. (2002). Some new perspectives the method of control variates. In Monte Carlo and Quasi-Monte Carlo Methods 2000, pages 27-49. Springer.
[15] Gorham, J. and Mackey, L. (2015). Measuring sample quality with Stein's 's method. Advances in Neural Information Processing Systems, 28.
[16] Hesterberg, T. (1995). Weighted average importance sampling and defensive mixture distribu- tions. Technometrics, 37(2):185-194,
[17] Hoffman, M. D., Gelman, A., et al. (2014). The no-u-turn sampler: adaptively setting path lengths hamiltonian monte carlo. J. Mach. Learn. Res., 15(1):1593-1633.
[18] Jie, T. and Abbeel, P. (2010). On a connection between importance sampling and the likelihood ratio policy gradient. In Advances in Neural Information Processing Systems, volume 23. Curran Associates, Inc.
[19] Jin, C., Netrapalli, P. Ge, R., Kakade, S. M., and Jordan, M. I. (2019). A short note on concen- tration inequalities for random vectors with subgaussian norm. arXiv preprint arXiv: 190.. 03736.
[20] Jourdain, B. (2009). Adaptive variance reduction techniques in finance. Radon Series Comp. Appl. Math, 8:1-18.
[21] Kawai, R. (2020). Adaptive importance sampling and control variates. Journal of Mathematical Analysis and Applications, 483(1):123308)
[22] Kloek, T. and Van Dijk, H. K. (1978). Bayesian estimates of equation system parameters: an application of integration by Monte Carlo. Econometrica: Journal of the Econometric Society, pages 1-19.
[23] Korba, A. and Portier, F. (2022). Adaptive importance sampling meets mirror descent: a biri-variance tradeoff. In International Conference on Artificial Intelligence and Statistics, pages 1 11503-11527. PMLR.
[24] Leluc, R., Portier, F,, and Segers, J. (2021). Control variate selection for Monte Carlo integration. Statistics and Computing, 31.
[25] Liu, H., Feng, Y., Mao, Y., Zhou, D., Peng, J., and Liu, Q. (2017). Action- control variates for policy optimization via Stein's 's identity. arXiv preprint arXiv:
[26] Martino, I.., Elvira, V., Luengo, D., and Corander, J. (2017). Layered adaptive importance sampling. Statistics and Computing, 27(3):599-622 623.
[27] Mira, A., Solgi, R., and Imparato, D. (2013). Zero variance Markov Chain Monte Carlo for Bayesian estimators. Statistics and Computing, 23(5):653-662.
[28] Oates, C. J., Girolami, M., and Chopin, N. (2017). Control functionals for Monte Carlo integration. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 79(3):695- 718.
[29] Oh, M.-S. and Berger, J. O. (1992), Adaptive importance sampling in monte carlo integration. Journal of Statistical Computation and Simulation, 41(3-4):143-168
[30] Owen, A. and Zhou, Y. (2000). Safe and effective importance sampling. J. Amer. Statist. Assoc., 95(449):135-143
[31] Paszke, A., Gross, S., Chintala, S., Chanan, G., Yang, E., DeVito, 7., Lin, Z., Desmaison, A., Antiga, L., and Lerer, A. (2017). Automatic differentiation in pytorch.
11

---

[32] Patil, A., Huard, D., and Fonnesbeck, C. J. (2010). Pyme: Bayesian stochastic modelling in python. Journal of statistical software, 35(4):1.
[33] Plassier, V., Portier, F., and Segers, J. (2020). Risk bounds when learning infinitely many response functions by ordinary linear regression. arXiv preprint arXiv:2006.0922 To appear in Annales de Institttt Henri Poincaré (B) Probabilités et Statistiques.
[34] Portier, F. and Delyon, B. (2018). Asymptotic optimality of adaptive importance sampling. Advances in Neural Information Processing Systems, 31.
[35] Portier, F. and Segers, J. (2019). Monte Carlo integration with growing number of control variates. Journal of Applied Probability, 56(4): 168-1186.
[36] Stein, C. (1972). A bound for the error in the normal approximation to the distribution of sum dependent random variables. In Proceedings of the sixth Berkeley symposium on mathematical statistics and probability, volume 2: Probability theory, volume 6, pages 583-603. University of California Press.
[37] Tropp. (2015). An introduction to matrix concentration inequalities. arXiv:1501001571.
[38] Wang, C., Chen, X., Smola, A., and Xing. E. (2013). Variance reduction for stochastic gradient optimization. Advances in Neural Information Processing Systems, volume 26. Curran Associates, Inc.
[39] Zhang, P. (1996), Nonparametric importance sampling. J. Amer. Statist. Assoc., 91(435): 1245- 1253.
12


Reviews:
[{"review_id": "QHrUBibw5qN", "reviewer": "Reviewer_KHDF", "summary": "The authors present AISCV estimator that combined Adaptive Importance Sampling (AIS) and control variates (CV). The paper presents a theoretical proof of its convergence rates and presents practical considerations that are really helpful. ", "questions": "Would it be worth having a baseline of just CV estimators? It is interesting to see how big the reduction is from CV -> AISCV. Especially, how would it compare on Figure 3. My thinking is the \"virtually zero\" variance is a CV effect rather than an AISCV effect. Is this correct?", "limitations": "The authors did not address this issue. I do not see any reason why the authors would need to address it either, the paper does not touch any sensitive area of society.", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The theoretical side of this paper is extremely well executed. Everything is carefully described and seems technically sound (I didn't carefully check details in the proof of Theorem 1, but seems valid based on what we know from AIS and CV. I did check the outline of the proof). \n\nI would like the authors to reflect on if AISCV is just a combination of AIS and CV, or if they think it is greater than \"the sum of its parts\".  \n\nCould you Figure 1(a) and 1(b) have the same cutoff on the y-axis for easier visual comparison between input dimensions. Same where appropriate in Figure 2.\n\nThis paper has enough quality to be accepted to NeurIPS. The drawdown is, I think, that the short format of NeurIPS is making this paper poorer than it could be in a longer format. The experimental section is likely a bit sub-par to many other papers at this venue. Yet, the authors have more experiments in the appendix that I did not read carefully.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "6IehXvQ3aGG", "reviewer": "Reviewer_tPbu", "summary": "This paper proposes a novel method for approximating intractable integrals that combines adaptive importance sampling and control variates. The authors prove a probabilistic bound on the error of their proposed estimator and demonstrate its superiority to the simple AIS estimator. ", "questions": "Could it ever be the case that poorly (or even adversarially) constructed control variates result in unboundedly worse performance than AIS? Or is Assumption 2 a sufficient safeguard against this kind of behavior? \n\nCan the authors speak to why their AISCV methods in Figure 2 demonstrate a higher variance than AIS in this setting? It seems odd to me that control variants, a method originally proposed as a variance reduction technique, could cause this behavior. ", "limitations": "Related to my questions above, it seems that the use of control variates introduces a whole host of design choices that have the potential to result in poor performance. While the theoretical results speak to what requirements the control variates need to satisfy, in practice it seems like this could result in undesirable behavior if the requirements are not met.", "rating": 7, "confidence": 2, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Overall, I found this paper an enjoyable read although I confess, many of the technical details surrounding the theoretical results escaped me. \n\n1. In terms of originality, the primary proposed idea is not particularly inspired, in that it combines two commonly-used, existing methods. \n2. The level of analysis, both theoretical and empirical, is compelling. \n3. The paper is extremely well written and clear. \n4. The work seems likely to be significant as the task of estimating integrals is commonplace and the proposed method appears sufficiently general so as to be widely applicable. ", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "5NvdPn9BS5D", "reviewer": "Reviewer_KrYe", "summary": "In this paper, the authors combine adaptive importance sampling (AIS) with control variates (CV) leading to a novel Monte Carlo estimator which they call AISCV. Non-asymptotic bounds were developed for the novel estimator and the experiments section demonstrated the superiority of AISCV versus a standard AIS estimate. ", "questions": "Below I will list general questions/comments\n\n1. In line 126, I think there is a typo. The paper says that $\\mathbb{E}_f[g - \\beta^\\top h] = 0$ but it should be $\\mathbb{E}_f[g - \\beta^\\top h] = \\mathbb{E}_f[g]$.", "limitations": "There is no section in the paper describing the limitations of the approach. While presumably the limitations of AISCV are similar to AIS and CV, it would be beneficial to state them in the paper.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "# Strengths\nI can not express how cool the proposed approach is. While standard CV estimators require $\\beta$ to be recomputed for every target function $g$, AISCV sidesteps this completely! Instead---just like AIS---quadrature weights are computed using the particles and their corresponding importance weights, where, crucially, the quadrature weights **are not a function of** $g$! Thus, this process is done once and can be used for any $g$. Moreover, the bounds on the AISCV estimate are very impressive as well!\n\n# Weaknesses\n\nWhile the experiments section clearly demonstrated the superiority of AISCV, I wish a more practical experiment was done. AISCV was motivated by wanting to compute posterior expectations. Thus, it was disappointing that the only example done was Bayesian linear regression where the posterior is available analytically. Instead, it would have been more compelling if the AISCV estimate was compared against Bayesian linear regression, which can be sampled by using HMC for instance. \n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept