[everypage]FERMILAB-PUB-23-471-CSAID current page.north east -2in -0.5in

Variance Reduction via Simultaneous Importance Sampling and Control Variates
Techniques Using Vegas

Prasanth Shyamsundar1*, Jacob Scott2, Stephen Mrenna3, Konstantin T. Matchev4,
Kyoungchul Kong2

1 Fermilab Quantum Institute, Fermi National Accelerator Laboratory,

Batavia, IL 60510, USA

2 Department of Physics and Astronomy, University of Kansas,

Lawrence, KS 66045, USA

3 Computational Science and Artificial Intelligence Division, Fermi National
Accelerator Laboratory, Batavia, IL 60510, USA

4 Institute for Fundamental Theory, Physics Department, University of Florida,

Gainesville, FL 32611, USA

*<EMAIL_ADDRESS>
## Abstract

Monte Carlo (MC) integration is an important calculational technique in the
physical sciences. Practical considerations require that the calculations are
performed as accurately as possible for a given set of computational
resources. To improve the accuracy of MC integration, a number of useful
variance reduction algorithms have been developed, including importance
sampling and control variates. In this work, we demonstrate how these two
methods can be applied simultaneously, thus combining their benefits. We
provide a python wrapper, named CoVVVR, which implements our approach in the
Vegas program. The improvements are quantified with several benchmark examples
from the literature.

###### Contents

  1. 1 Introduction
  2. 2 Control Variates and Vegas
    1. 2.1 Naive Monte Carlo Integration
    2. 2.2 Importance Sampling
    3. 2.3 Control Variates
    4. 2.4 Combining Importance Sampling and Control Variates
    5. 2.5 Repurposing Vegas Intermediate Outputs as Control Variates
  3. 3 Results
    1. 3.1 Test Cases
      1. 3.1.1 $d$-dimensional Gaussians
      2. 3.1.2 $d$-dimensional Camel functions
      3. 3.1.3 Entangled circles
      4. 3.1.4 Annulus with cuts
      5. 3.1.5 One-loop scalar box integral
      6. 3.1.6 Polynomial functions
    2. 3.2 Numerical Comparisons
  4. 4 CoVVVR: Control Variates & Vegas Variance Reduction
    1. 5 Summary and Outlook
    2. A Multiple Control Variates

##  1 Introduction

In many fields of science, including particle physics, one has to compute the
value of an integral

$F\equiv\int_{\Omega}d{\mathbf{x}}\,f({\mathbf{x}})\,,$ (1)

over some domain $\Omega$ with volume

$V(\Omega)\equiv\int_{\Omega}d{\mathbf{x}}\,,$ (2)

where ${\mathbf{x}}\in\mathbb{R}^{d}$ is a $d$-dimensional vector of
independent variables. The function $f({\mathbf{x}})$ can be evaluated for a
given ${\mathbf{x}}$, but can be arbitrarily complicated. In high energy
physics, integrations like (1) are ubiquitous, arising when computing total
cross-sections (or differential cross-sections with respect to low-dimensional
event variables), particle lifetimes, convolutions with transfer functions
describing detector effects, etc.

Although widespread, this problem is fundamentally challenging — with the
exception of a few trivial cases (typically found in the textbooks), the
integral cannot be performed analytically, and one has to resort to numerical
methods for its evaluation. Monte Carlo (MC) methods are particularly suited
for high-dimensional integrals, since their accuracy scales as
$\sqrt{\frac{\text{Var}\mathopen{}\mathclose{{}\left[f}\right]}{N_{\text{trials}}}}$,
where $\text{Var}\mathopen{}\mathclose{{}\left[f}\right]$ is the variance of
$f({\mathbf{x}})$ over the domain $\Omega$ and $N_{\text{trials}}$ is the
number of trials [1]. One obvious way to improve the accuracy is to increase
$N_{\text{trials}}$, but this approach soon hits a wall due to resource
limitations. Therefore, much attention has been placed on designing variance-
reducing techniques. Some of the classic variance reduction methods include
importance sampling (IS), stratified sampling, antithetic variables, control
variates (CV), etc. [2]. More recent techniques use quasirandom or low-
discrepancy pointsets in a class of methods known as Quasi Monte Carlo [3, 4],
apply multigrid ideas in Multilevel Monte Carlo estimation [5], or leverage
machine learning (ML) [6]. A parallel research thrust has been the synthesis
of ideas, whereby one tries to apply two such techniques, for example, by
using several control variates [7], combining antithetic variates and control
variates [8, 9], or combining control variates and adaptive importance
sampling [10, 11].

In MC methods, random variates are used to sample the function of interest
over a number of trials. In importance sampling, the probability distribution
of the random variates is adjusted based on the values of the trials to more
closely mirror the behavior of the function. In general, such a change of
variables must be performed on a case-by-case basis and is difficult, if not
impossible, to accomplish for high-dimensional integrands. The Vegas code [12]
and algorithm [13], instead, introduces a parametrization for this change of
variables. The domain of integration is divided into subregions of equal
importance so that the variance of the MC estimate is reduced compared to
other sampling choices. A related method, additionally implemented in Vegas\+
[14], is that of stratified sampling, whereby one uses stratification to
homogenize the resulting weights within each subregion.

The control variate method, instead, aims to reduce the variance by modifying
the integrand. This is accomplished by adding and subtracting a function (the
control variate) that is highly-correlated with the integrand and has a known
value for its integral. The integral of the control variate converges to the
known value, while the variance of the control variate partially cancels the
variance in the original integrand. The challenge is to find a control variate
that is indeed highly correlated (or anti-correlated) with the integrand $f$.
Usually, this can be achieved for only special, low-dimensional cases. Control
variates have been used in various domain-inspired applications in [15, 16,
17, 18, 19, 20, 21].

The thrust of this paper is to develop a general strategy to combine the well-
established importance sampling technique with control variates. We note that
the Vegas method already computes a function that is highly correlated with
the integrand in a general way. We illustrate how this function (or its
evolution as it converges to an optimal value) can serve as a control variate
to improve the accuracy for a given computational budget, or equivalently, to
reduce the computation time for a given target accuracy. The paper is
accompanied with a python code, CoVVVR, short for  Control  variates and
Vegas  variance  reduction. The source code is publicly available at
https://github.com/crumpstrr33/covvvr, together with an introductory tutorial
for its usage.

This paper is organized as follows. Section 2 lays the theoretical groundwork
of our approach. Sections 2.1, 2.2 and 2.3 briefly review the basic ideas of
the Monte Carlo integration, the importance sampling and the control variates
methods, respectively. Section 2.4 introduces the joint application of
importance sampling and control variates, while Section 2.5 explains how one
can leverage the successive Vegas approximations as control variates. Sec. 3
presents the results from our numerical experiments quantifying the achieved
precision improvement on some known benchmark examples. Section 4 contains a
description and usage examples of the CoVVVR code. Section 5 is reserved for a
summary and conclusions.

##  2 Control Variates and Vegas

In this section, we describe the main idea of our method using a one-
dimensional integration example

$I=\int_{a}^{b}\differential{x}f(x)$ (3)

of a real function $f(x)$ from $a$ to $b$, where for definiteness we take
$a<b$. In special cases, $a$ can be $-\infty$ and/or $b$ can be $+\infty$.

###  2.1 Naive Monte Carlo Integration

In the naive Monte Carlo approach, $x$ is sampled uniformly over the domain
$\Omega$, i.e., $x\sim U(\Omega)$, obtaining $N$ independent and uniformly
distributed samples $x_{1},x_{2},\ldots,x_{N}$. An estimate $\widehat{I}_{N}$
of the integral (3) is obtained in terms of the expectation value $E_{x\sim
U(\Omega)}[f]$ of the function $f(x)$ over the domain $\Omega$:

$\displaystyle\widehat{I}_{N}=V(\Omega)~{}\mathrm{E}_{x\sim
U(\Omega)}\big{[}f\big{]}=V(\Omega)\times\frac{1}{N}\sum_{i=1}^{N}f(x_{i})\,,$
(4)

where $x$ is sampled uniformly in the domain $\Omega$, i.e., $x_{i}\sim
U(\Omega)$. In the case of a one-dimensional integral as in (3),
$V(\Omega)=b-a$.

The uncertainty on $\widehat{I}_{N}$ can be estimated in terms of the sample
variance $\text{Var}\mathopen{}\mathclose{{}\left[\widehat{I}_{N}}\right]$

$\delta\widehat{I}_{N}\approx\sqrt{\text{Var}\mathopen{}\mathclose{{}\left[\widehat{I}_{N}}\right]}=V(\Omega)\sqrt{\frac{\text{Var}\mathopen{}\mathclose{{}\left[f}\right]}{N}}\,,$
(5)

which decreases as $N^{-1/2}$. However, increasing $N$ arbitrarily is
infeasible, as it runs into resource limitations. This motivates methods which
attempt to reduce the variance
$\text{Var}\mathopen{}\mathclose{{}\left[f}\right]$ of the integrand function
$f(x)$. Two such methods are discussed in the next two subsections.

###  2.2 Importance Sampling

In importance sampling, we choose a sampling function for the random variates
that resembles $f(x)$ as closely as possible, and whose integral over the
range $(a,b)$ is known. By rescaling with the value of this integral, the new
sampling function can be turned into a unit normalized probability
distribution function (PDF), which we shall denote with $p(x)$:

$\int_{a}^{b}\differential{x}p(x)=1.$ (6)

The integral of interest (3) can be rewritten as

$\displaystyle
I=\int_{a}^{b}\differential{x}f(x)=\int_{a}^{b}\differential{x}\,p(x)\,\frac{f(x)}{p(x)}=\mathrm{E}_{x\sim
p(x)}\mathopen{}\mathclose{{}\left[\frac{f}{p}}\right],$ (7)

which can in turn be estimated with Monte Carlo in $N$ trials as

$\displaystyle\widehat{I}_{N}=\frac{1}{N}\sum_{i=1}^{N}\frac{f(x_{i})}{p(x_{i})}\,,\qquad\text{where
}x_{i}\sim p(x).$ (8)

In analogy to (5) the error on the estimate (8) is given by

$\displaystyle\delta\widehat{I}_{N}\approx\sqrt{\text{Var}\mathopen{}\mathclose{{}\left[\hat{I}_{N}}\right]}$
$\displaystyle=\frac{1}{\sqrt{N}}\mathopen{}\mathclose{{}\left[-I^{2}+\int_{a}^{b}\differential{x}\,p(x)\,\frac{f^{2}(x)}{p^{2}(x)}}\right]^{1/2}.$
(9)

This is minimized when

$\displaystyle p(x)$
$\displaystyle\propto\mathopen{}\mathclose{{}\left|\big{.}f(x)}\right|.$ (10)

In other words, whenever the two functions $f(x)$ and $p(x)$ have similar
shapes, the variance is reduced and the precision of the estimate is improved.
In fact, if we could choose $p(x)$ to be exactly proportional to $f(x)$
everywhere, then $\widehat{I}_{N}$ is exactly equal to $I$ for any $N$. Thus,
importance sampling is most beneficial when the function $p(x)$ mimics $f(x)$
as closely as possible.

###  2.3 Control Variates

In contrast to the importance sampling method, which involves a rescaling of
the integrand, the control variates method adds to $f(x)$ a term The modified
integrand can be taken as

$f_{c}(x)\equiv f(x)+c\bigg{(}g(x)-\frac{G}{b-a}\bigg{)}\,,$ (11)

with

$G\equiv\int_{a}^{b}\differential{x}g(x)\,$ (12)

and where $c$ is a parameter which at this point we are free to choose. It is
easy to see that the modification (11) does not change the value of the
integral, i.e.

$\displaystyle\int_{a}^{b}\differential{x}f(x)=\int_{a}^{b}\differential{x}f_{c}(x)\,,\qquad\forall
c\in\mathbb{R}\,.$ (13)
$\displaystyle\widehat{I}_{N}=V(\Omega)\,\frac{1}{N}\sum_{i=1}^{N}f_{c}(x_{i})\,,\qquad\text{where
}x_{i}\sim U(\Omega).$ (14)

We can now leverage the freedom to choose the value of the parameter $c$ in
order to minimize the variance. Requiring

$\displaystyle\partialderivative{\text{Var}\mathopen{}\mathclose{{}\left[f_{c}}\right]}{c}$
$\displaystyle=0$ (15)

gives the optimal value of $c$ as

$\displaystyle
c^{\ast}=-\frac{\text{Cov}\mathopen{}\mathclose{{}\left(f,g}\right)}{\text{Var}\mathopen{}\mathclose{{}\left[g}\right]}.$
(16)

The resulting variance is

$\displaystyle\text{Var}\mathopen{}\mathclose{{}\left[f_{c^{\ast}}}\right]=\text{Var}\mathopen{}\mathclose{{}\left[f}\right]-\frac{\text{Cov}\mathopen{}\mathclose{{}\left(f,g}\right)^{2}}{\text{Var}\mathopen{}\mathclose{{}\left[g}\right]^{2}}\text{Var}\mathopen{}\mathclose{{}\left[g}\right]=\bigg{[}1-\rho^{2}(f,g)\bigg{]}\text{Var}\mathopen{}\mathclose{{}\left[f}\right],$
(17)

where

$\rho(f,g)\equiv\frac{\text{Cov}\mathopen{}\mathclose{{}\left(f,g}\right)}{\sqrt{\text{Var}\mathopen{}\mathclose{{}\left[f}\right]\text{Var}\mathopen{}\mathclose{{}\left[g}\right]}}$
(18)

is the familiar Pearson correlation coefficient. Note that if
$\absolutevalue{\rho(f,g)}>0$, the variance is necessarily reduced.
Furthermore, the higher the correlation between $g(x)$ and $f(x)$, the larger
the benefit. Therefore, just like in the method of importance sampling, we
desire a function $g(x)$ that i) is highly correlated with $f(x)$ and, ii) has
a known expectation value $\mathrm{E}[g]=G/(b-a)$.

The method can be easily generalized to the case of multiple control variates.
Appendix A contains the derivation for finding the optimal values $c_{i}^{*}$
of the respective coefficients $c_{i}$ in that case.

###  2.4 Combining Importance Sampling and Control Variates

The two methods discussed in the previous two subsections 2.2 and 2.3
(importance sampling and control variates) can be combined together as
follows. Given a known PDF $p(x)$ and a control variate function $g(x)$, we
modify the integrand as

$\displaystyle
I=\int_{a}^{b}\differential{x}\,p(x)\mathopen{}\mathclose{{}\left[\frac{f(x)}{p(x)}+c\mathopen{}\mathclose{{}\left(\frac{g(x)}{p(x)}-E_{p}\mathopen{}\mathclose{{}\left[\frac{g}{p}}\right]}\right)}\right].$
(19)

The corresponding Monte Carlo estimate is

$\displaystyle\widehat{I}_{N}=\frac{1}{N}\sum_{i=1}^{N}\mathopen{}\mathclose{{}\left[\frac{f(x_{i})}{p(x_{i})}+c\mathopen{}\mathclose{{}\left(\frac{g(x_{i})}{p(x_{i})}-E_{p}\mathopen{}\mathclose{{}\left[\frac{g}{p}}\right]}\right)}\right]\,,\qquad\text{where
}x_{i}\sim p(x).$ (20)

We would still need a $p(x)$ that is approximately proportional to $f(x)$ and
a $g(x)$ whose integral is known, so that $f/p$ is correlated with $g/p$. In
this case, the optimal value $c^{\ast}$ is given by

$\displaystyle
c^{\ast}=-\frac{\text{Cov}\mathopen{}\mathclose{{}\left(f/p,g/p}\right)}{\text{Var}\mathopen{}\mathclose{{}\left[g/p}\right]}.$
(21)

###  2.5 Repurposing Vegas Intermediate Outputs as Control Variates

Vegas is an adaptive and iterative Monte Carlo where, at each iteration $i$, a
unit-normalized probability distribution $p_{i}(x)$ is updated to serve as a
probability distribution for importance sampling as described in Section 2.2
[13]. Since Vegas stores and reports these sampling distributions for each
iteration, they can be usefully repurposed as control variates. We can then
apply the result in (20), by choosing the IS function as $p(x)=p_{n}(x)$ for
the final iteration $n$ and the control variate function $g(x)=p_{i}(x)$ from
some previous iteration $i<n$. Note that since $p_{i}(x)$ is unit-normalized
the last term in (20) becomes

$\displaystyle
E_{p}\mathopen{}\mathclose{{}\left[\frac{g}{p}}\right]=E_{p_{n}}\mathopen{}\mathclose{{}\left[\frac{p_{i}}{p_{n}}}\right]=\int_{a}^{b}\differential{x}\,p_{n}(x)\,\frac{p_{i}(x)}{p_{n}(x)}=\int_{a}^{b}\differential{x}\,p_{i}(x)=1\,.$
(22)

In this way, Vegas can be used to simultaneously provide both the sampling
distribution $p(x)$ and the control variate $g(x)$. In principle, different
prescriptions for selecting the previous iteration $i$ to be used in (22) can
be designed. Choosing the optimal among those prescriptions can be done on a
case by case basis, as discussed below in Section 3.

##  3 Results

###  3.1 Test Cases

For comparison with other studies, we compute values for the benchmark
functions presented in Ref. [6]. The domain of integration for each variable
is over the range of 0 to 1, and the corresponding values are given in the
second column of Table III in Ref. [6]. We confirmed those values as shown in
Table 1 and, in what follows, focus on the accuracy of the MC estimate.

For clarity, the definition of the functions used as benchmarks is reproduced
below.

####  3.1.1 $d$-dimensional Gaussians

The $d$-dimensional Gaussians are defined as

$\displaystyle
f_{1}(\mathbf{x})=\frac{1}{(\sigma\sqrt{\pi})^{d}}\exp\mathopen{}\mathclose{{}\left(-\frac{1}{\sigma^{2}}\sum_{i=1}^{d}\mathopen{}\mathclose{{}\left(x_{i}-\mu}\right)^{2}}\right),$
(23)

with mean $\mu=0.5$ and standard deviation $\sigma=0.2$. They serve as a good
starting point and display the effectiveness of control variates for separable
functions. We consider $d=2$, $d=4$, $d=8$ and $d=16$.

####  3.1.2 $d$-dimensional Camel functions

The $d$-dimensional Camel functions are defined by

$\displaystyle
f_{2}(\mathbf{x})=\frac{1}{2(\sigma\sqrt{\pi})^{d}}\mathopen{}\mathclose{{}\left[\exp\mathopen{}\mathclose{{}\left(-\frac{1}{\sigma^{2}}\sum_{i=1}^{d}\mathopen{}\mathclose{{}\left(x_{i}-\mu_{1}}\right)^{2}}\right)+\exp\mathopen{}\mathclose{{}\left(-\frac{1}{\sigma^{2}}\sum_{i=1}^{d}\mathopen{}\mathclose{{}\left(x_{i}-\mu_{2}}\right)^{2}}\right)}\right]$
(24)

in terms of three parameters, $\mu_{1}=1/3$, $\mu_{2}=2/3$, and $\sigma=0.2$.
Unlike the Gaussians case, the integration variables are not separable. We
consider $d=2$, $d=4$, $d=8$ and $d=16$.

####  3.1.3 Entangled circles

This function is given by

$\displaystyle f_{3}(x_{1},x_{2})$
$\displaystyle=x_{2}^{a}\exp[-w\Big{|}(x_{2}-p_{2})^{2}+(x_{1}-p_{1})^{2}-r^{2}\Big{|}\bigg{]}$
$\displaystyle+(1-x_{2})^{a}\exp[-w\Big{|}\big{(}x_{2}-(1-p_{2})\big{)}^{2}+\big{(}x_{1}-(1-p_{1})\big{)}^{2}-r^{2}\Big{|}\bigg{]}$
(25)

It is largely concentrated on two overlapping circles of radius $r$ with
centers at $(p_{1},p_{2})$ and $(1-p_{1},1-p_{2})$, respectively. For the
numerical experiments, we use $p_{1}=0.4$, $p_{2}=0.6$, $r=0.25$, $w=1/0.004$,
and $a=3$.

####  3.1.4 Annulus with cuts

The fourth function is an annulus defined by cuts at $r_{\text{min}}$ and
$r_{\text{max}}$:

$\displaystyle f_{4}(x_{1},x_{2})=\begin{cases}1&\mbox{if
}r_{\text{min}}<\sqrt{x_{1}^{2}+x_{2}^{2}}<r_{\text{max}},\\\
0&\mbox{else}.\end{cases}$ (26)

The cut parameters are chosen to be $r_{\text{min}}=0.2$ and
$r_{\text{max}}=0.45$.

####  3.1.5 One-loop scalar box integral

The fifth function represents a one-loop scalar box integral encountered in
the calculation of $gg\to gh$, an important contribution to Higgs production
from gluon fusion. The integral is

$\displaystyle I_{5}$ $\displaystyle=$ $\displaystyle
S_{\text{Box}}(s_{12},s_{23},s_{1},s_{2},s_{3},s_{4},m_{t}^{2},m_{t}^{2},m_{t}^{2},m_{t}^{2})$
(27) $\displaystyle+$ $\displaystyle
S_{\text{Box}}(s_{23},s_{12},s_{2},s_{3},s_{4},s_{1},m_{t}^{2},m_{t}^{2},m_{t}^{2},m_{t}^{2})$
$\displaystyle+$ $\displaystyle
S_{\text{Box}}(s_{12},s_{23},s_{3},s_{4},s_{1},s_{2},m_{t}^{2},m_{t}^{2},m_{t}^{2},m_{t}^{2})$
$\displaystyle+$ $\displaystyle
S_{\text{Box}}(s_{23},s_{12},s_{4},s_{1},s_{2},s_{3},m_{t}^{2},m_{t}^{2},m_{t}^{2},m_{t}^{2}),$

where

$\displaystyle
S_{\text{Box}}(s_{12},s_{23},s_{1},s_{2},s_{3},s_{4},m_{1}^{2},m_{2}^{2},m_{3}^{4},m_{4}^{2})=\int_{0}^{1}\frac{\differential{x_{1}}\differential{x_{2}}\differential{x_{3}}}{\widetilde{\mathcal{F}}_{\text{Box}}^{2}}$
(28)

and

$\displaystyle\widetilde{\mathcal{F}}_{\text{Box}}=$ $\displaystyle-
s_{12}x_{2}-s_{23}x_{1}x_{3}-s_{1}x_{1}-s_{2}x_{1}x_{2}-s_{3}x_{2}x_{3}-s_{4}x_{3}$
(29)
$\displaystyle+(1+x_{1}+x_{2}+x_{3})(x_{1}m_{1}^{2}+x_{2}m_{2}^{2}+x_{3}m_{3}^{2}+m_{4}^{2}).$

The integrand function $f_{5}(x_{1},x_{2},x_{3})$ in this case is a sum of
four terms of the type $1/{\widetilde{\mathcal{F}}_{\text{Box}}}^{2}$. We test
with the same parameters as in [6]: $s_{12}=130^{2}$, $s_{23}=-s_{12}$,
$s_{1}=s_{2}=s_{3}=0$, $s_{4}=125^{2}$, and $m_{i}=173.9$ for $i=1-4$.

####  3.1.6 Polynomial functions

The final set of integrand functions are quadratic polynomials of $d$
variables

$\displaystyle f_{6}(\mathbf{x})=\sum_{i=1}^{d}x_{i}(1-x_{i}).$ (30)

We consider $d=18$, $d=54$ and $d=96$.

###  3.2 Numerical Comparisons

| Mean | Normalized RMS Error
---|---|---
Function | Dim | True Value | Vegas | CVInt | Vegas | CVInt
Gaussian | 2 | 0.999186 | 0.999190 | 0.999188 | 1.0049e-04 | 9.1006e-05
4 | 0.998373 | 0.998379 | 0.998379 | 0.000145 | 0.000132
8 | 0.996749 | 0.996750 | 0.996866 | 0.000202 | 0.000231
16 | 0.993509 | 0.993507 | 0.993507 | 0.000293 | 0.000265
Camel | 2 | 0.981660 | 0.981618 | 0.981617 | 0.001286 | 0.001283
4 | 0.963657 | 0.963559 | 0.963560 | 0.003180 | 0.003173
8 | 0.928635 | 0.929091 | 0.929223 | 0.010283 | 0.010295
16 | 0.862363 | 0.784645 | 0.784860 | 0.546183 | 0.544413
Entangled Circles | 2 | 0.013680 | 0.013687 | 0.013687 | 0.004349 | 0.004348
Annulus with Cuts | 2 | 0.127627 | 0.127644 | 0.127642 | 0.003773 | 0.003789
Scalar Top Loop | 3 | 1.9374e-10 | 1.9370e-10 | 1.9370e-10 | 0.000328 | 0.000322
Polynomial | 18 | 3.000000 | 2.999998 | 2.999999 | 2.6498e-05 | 2.1906e-05
54 | 9.000000 | 8.999992 | 8.999995 | 4.1208e-05 | 2.9833e-05
96 | 16.000000 | 15.999929 | 15.999954 | 7.2611e-05 | 5.1897e-05

Table 1:  Results for the test case integrals described in Section 3.1 with
$n=50$ iterations, $N_{\text{events}}=5,000$ events per iteration, and
averaged over $N_{\text{runs}}=1,000$ runs. The results displayed in the CVInt
column were obtained with the one control variate that gives maximum variance
reduction. The normalized RMS error is given by (31).  Figure 1: Correlation
between the function $f/p$ and the control variate $g/p$. This is shown for
the 16d Gaussian (left) and the 96d polynomial (right) run for $n=100$
iterations with, at most, $N_{\text{events}}=$ 10000 evaluations per
iteration. The iteration used for the control variate (CV) was chosen
automatically using the parameter cv_iters="auto1" as described in Section 4.
We see that there is a correlation between the functions shown both
qualitatively by the plot and quantitatively by a correlation coefficient of
$\rho\approx 0.42$ (left) and $\rho\approx 0.75$ (right). Additionally, the
expectation is approximately 1, in agreement with (22).

Results from integrating the six types of test functions from Section 3.1 are
shown in Table 1. The first two columns list the name of the function and the
dimension of the integration. The third column shows the expected true answer,
while the next two columns give the MC estimates from Vegas and from CoVVVR,
respectively. For the latter, we use the one control variate which gives the
maximum variance reduction, as explained below, and choose the optimal value
of $c^{\ast}$ according to (21). In the last two columns we show the
normalized RMS error defined as

$\frac{1}{I}\
\sqrt{\frac{1}{N_{\text{runs}}}\sum_{i=1}^{N_{\text{runs}}}(\widehat{I}_{i}-I)^{2}}.$
(31)

We see that, as expected, the accuracy using a CV is comparable or improved.

As implied by eq. (19), the variance reduction results from the presence of a
correlation between the function ratios $f(x)/p(x)$ and $g(x)/p(x)$. This
correlation is illustrated in Figure 1 for the case of the 16-dimensional
Gaussian (left) and the 96-dimensional polynomial (right) function. The
iteration used for the CV was chosen automatically. The correlation is readily
visible by eye, and the correlation coefficient is $\rho\approx 0.42$ (left)
and $\rho\approx 0.75$ (right).

The effect of adding more than one CV is illustrated in Table 2 and Figure 2.
We show results with one (1 CV), two (2 CV) or all 49 intermediate
approximations (All CVs) from Vegas as control variates. The results are
quoted in terms of the variance reduction in percentage (VRP) and the
corresponding computational cost (in seconds, as well as relative to the Vegas
benchmark time). In each case, we select the CV or CVs that reduce the
variance the most. Table 2 confirms that, by construction, the use of CVs
always improves the accuracy of the estimate. The size of the VRP effect
depends on the type of function at hand, and can vary from $\sim 1\%$ to as
much as $50\%$ for one CV and $60\%$ for two CVs. The associated computational
cost is an increase of about $1.5-2.5$ times for one CV and $2-4$ times for 2
CVs. Note that an increase in accuracy could also be achieved by increasing
the number of events used in a standard Vegas calculation. We discuss the
benefits of the CV method later.

| Vegas | 1 CV | 2 CVs | All CVs
---|---|---|---|---
Function | Dim | Time (s) | VRP | Time (s) | VRP | Time (s) | VRP | Time (s)
Gaussian | 2 | 0.07 | 17.02% | 0.09 | (1.4) | 31.40% | 0.14 | (2.1) | 47.15% | 1.76 | (26.2)
4 | 0.10 | 15.86% | 0.14 | (1.5) | 26.55% | 0.22 | (2.3) | 39.78% | 2.91 | (30.1)
8 | 0.15 | 17.28% | 0.25 | (1.7) | 23.11% | 0.38 | (2.6) | 33.95% | 5.74 | (39.4)
16 | 0.23 | 13.95% | 0.49 | (2.1) | 17.22% | 0.79 | (3.4) | 23.87% | 13.34 | (56.9)
Camel | 2 | 0.07 | 0.38% | 0.09 | (1.4) | 0.81% | 0.14 | (2.0) | 1.24% | 1.75 | (25.7)
4 | 0.10 | 0.12% | 0.14 | (1.5) | 0.33% | 0.22 | (2.2) | 0.66% | 2.94 | (30.0)
8 | 0.15 | 0.20% | 0.25 | (1.6) | 0.25% | 0.39 | (2.5) | 0.37% | 5.86 | (37.9)
16 | 0.26 | 0.74% | 0.52 | (2.0) | 3.17% | 0.81 | (3.2) | 7.72% | 13.70 | (53.4)
Entangled Circles | 2 | 0.08 | 0.31% | 0.10 | (1.2) | 0.36% | 0.15 | (1.8) | 0.67% | 1.81 | (22.4)
Annulus with Cuts | 2 | 0.05 | 1.25% | 0.08 | (1.6) | 2.94% | 0.12 | (2.3) | 16.93% | 1.77 | (33.1)
Scalar Top Loop | 3 | 0.12 | 7.31% | 0.14 | (1.2) | 49.33% | 0.21 | (1.8) | 57.91% | 2.30 | (19.9)
Polynomial | 18 | 0.26 | 29.50% | 0.56 | (2.2) | 29.74% | 0.89 | (3.5) | 51.36% | 15.49 | (59.9)
54 | 0.70 | 42.65% | 1.65 | (2.4) | 58.97% | 2.63 | (3.8) | 71.77% | 47.94 | (68.8)
96 | 1.33 | 49.63% | 3.04 | (2.3) | 63.77% | 4.79 | (3.6) | 78.18% | 86.55 | (65.3)

Table 2: Results for the variance reduction in percent (VRP) after using one,
two, or all 49 intermediate approximations from Vegas as control variates. In
each case, we ran $n=50$ iterations and at most $N_{\text{events}}=5000$
events per iteration averaged over $N_{\text{runs}}=100$ independent runs. The
time column shows the corresponding computational cost (in seconds and
relative to the Vegas benchmark time).

Figure 2: The variance reduction due to 1 control variate (along the diagonal)
and 2 control variates for all combinations. We show results for $n=50$
iterations with at most $N_{\text{events}}=$5,000 evaluations per iteration
averaged over $N_{\text{runs}}=$100 runs for a 16d Gaussian (left panel) and
the scalar top loop (right panel).

The dependence of the variance reduction on the choice of CVs is illustrated
in Figure 2. The heat map in each panel depicts the variance reduction due to
1 CV (along the diagonal) and 2 CVs for all combinations. We show results
averaged over 10 runs for a 16d Gaussian (left panel) and the scalar top loop
(right panel). In each run, we perform 50 iterations with at most 25,000
evaluations per iteration. Along each axis, each panel contains a plot of the
normalized variance (red line) and the correlation coefficient between the
target function $f(x)$ and the respective CV (purple line).

Figure 2 shows that the optimal iteration for choosing a CV can vary
significantly — in the case of a 16d Gaussian, it is around 30, while for the
scalar top loop integral, it is at the beginning. When we have the freedom to
choose two CVs, some interesting patterns appear as shown in the heat maps,
and the optimal choice for the Gaussian are the iterations around 25 and 40,
respectively.

| Standard Deviation |
---|---|---
Function | Dim | Events | Total Iter | CV Iter | Vegas | CVInt | VRP
Gaussian | 2 | 150,342 | 32 | 8 | 1.256e-04 | 1.246e-04 | 1.58%
4 | 703,711 | 150 | 37 | 8.288e-05 | 7.976e-05 | 7.37%
8 | 1,435,484 | 311 | 77 | 8.268e-05 | 7.980e-05 | 6.85%
16 | 2,775,411 | 614 | 153 | 8.513e-05 | 7.980e-05 | 12.13%
Camel | 2 | 342,915 | 73 | 18 | 1.053e-03 | 1.053e-03 | 0.00%
4 | 19,379,835 | 4,108 | 1,027 | 3.453e-04 | 3.451e-04 | 0.15%
8 | 47,220,394 | 10,000 | 2,500 | 7.122e-04 | 7.116e-04 | 0.16%
16 | 1,162,234 | 246 | 61 | 7.860e-02 | 7.860e-02 | 0.00%
Entangled Circles | 2 | 48,399,992 | 10,000 | 2,500 | 4.357e-06 | 4.352e-06 | 0.23%
Annulus with Cuts | 2 | 1,862,585 | 389 | 97 | 4.308e-04 | 4.292e-04 | 0.74%
Scalar-top-loop | 3 | 110,748 | 24 | 6 | 6.436e-14 | 6.433e-14 | 0.12%
Polynomial | 18 | 1,807,196 | 400 | 100 | 2.836e-05 | 2.223e-05 | 38.56%
54 | 4,450,434 | 985 | 246 | 7.299e-05 | 5.346e-05 | 46.37%
96 | 13,784,476 | 3,051 | 762 | 1.561e-04 | 1.049e-04 | 54.79%

Table 3: Comparison between native Vegas and the addition of a single control
variate (CV). Number of events and iterations are chosen based on reaching the
precision required in Ref. [6], namely, relative uncertainty of $10^{-4}$ for
the first 11 cases and $10^{-5}$ for the last three. There are a maximum of
$N_{\text{events}}=$5,000 events per iteration. The iteration chosen for the
CV is set at a quarter of the total number of iterations. VRP is the variance
reduction in percentage; it is the reduction in variance between Vegas and
CVInt by percentage.

Since the optimal choice of the iteration is a priori unknown, we show in
Table 3 the corresponding results when the iteration is decided and fixed at
the very beginning. In this case, we pick the CV from the iteration which is
1/4 of the way from the beginning, i.e., $i=\lfloor n/4\rfloor$. The number of
events $N_{\text{events}}$ and total number of iterations $n$ are chosen based
on reaching the precision required in Ref. [6], namely, relative uncertainty
of $10^{-4}$ for the first 11 cases and $10^{-5}$ for the last three. We see
that even when the CV is fixed rather arbitrarily, the variance reduction is
still significant, and can be $\sim 50\%$, as in the case of the polynomial
functions.

##  4 CoVVVR: Control Variates & Vegas Variance Reduction

The paper is accompanied with a python code, CoVVVR, short for  Control
variates & Vegas Variance Reduction, that can be used to reproduce our
results. The source code is publicly available at
https://github.com/crumpstrr33/covvvr. In this section, we provide an
introduction tutorial.

Installation: CoVVVR can be installed via `pip`:

⬇

$ python -m pip install covvvr

\end{result}

Usage: The workflow involves creating a class that inherits the ``Function—
subclass and passing that to the ``CVIntegrator—. The ``Function— class
contains the function to be integrated but also other information such as its
name, the true value of the integration (if available) and parameters of the
function. This class can be a built-in function, such as those used in this
paper, or a custom-made function. The ``CVIntegrator— class does the
integration and stores the results like mean and variance.

Using a Built-In Function:

⬇

from covvvr import CVIntegrator

from covvvr.functions import NGauss

\par# Create 16-dimensional Gaussian

ng = NGauss(16)

\par# Print out all parameters of the class instance

print(ng, ’\n’)

\par# Create integrator classm & use the 20th iteration as control variate

cvi = CVIntegrator(ng, evals=1000, tot_iters=50, cv_iters=20)

\par# Run the integration

cvi.integrate()

\par# Print info

cvi.compare(rounding=5)

This integrates a 16-dimensional Gaussian with 50 iterations and 5000
evaluations per iteration in Vegas. It uses the 20th iteration adaptation from
Vegas as the control variate. The output is

⬇

NGauss(dimension=16, name=16D Gaussian, true_value=0.9935086032227194,

mu=0.5, sigma=0.2)

\par| No CVs | With CVs

———+————-+————-

Mean | 0.99528 | 0.99579

Variance | 4.17808e-06 | 3.54189e-06

St Dev | 0.00204 | 0.00188

VRP | | 15.22696\end{result}

\\\

\par\par\noindent\textbf{Adding User-Defined Function}: The \verb|make_func|
function allows for the definition of a user-defined functions via the
{\ttFunction} subclass. As an example, consider the 2-dimensional function
$f(x_1, x_2) = a x_1^2 + b x_2^2$. It can be defined as\\\

\par\par\begin{minted}{python}

from covvvr import CVIntegrator, make_func

\par# Create function, note that it is vetorized using Numpy slicing

def f(self, x):

return self.a * x[:, 0]**2 + self.b * x[:, 1]

\par# Creating class with name ’WeightedPoly’ and assigning values to the

# parameters in the function

wpoly = make_func(

cname=’WeightedPoly’,

dimension=2,

function=f,

name=’Weighted Polynomial’,

a=0.3,

b=0.6

)

\par# Print out parameters of class (note ‘true_value‘ isn’t shown but it

# can be added if one wants to keep track of that

print(wpoly, ’\n’)

\par# Create integrator class and use multiple control variates

cvi = CVIntegrator(

function=wpoly,

evals=1000,

tot_iters=50,

cv_iters=[10, 15, 20, 25, 30, 35]

)

\par# Run the integration

cvi.integrate()

\par# Print info

cvi.compare(rounding=5)

which outputs

⬇

WeightedPoly(dimension=2, name=Weighted Polynomial, a=0.3, b=0.6)

\par| No CVs | With CVs

———+————-+————-

Mean | 0.39984 | 0.40004

Variance | 5.83781e-08 | 4.90241e-08

St Dev | 2.41616e-04 | 2.21414e-04

VRP | | 16.02309\end{result}

\par\noindentNote that vectorization of the integrand greatly increases the
speed of the computation.

\parTo not have to deal with classes, one can use the \verb|classic_integrate|
function that does the steps above and returns the \verb|CVIntegrator| object.
To run the previous code block, one would use

\par\begin{minted}{python}

cvi = classic_integrate(

function=f,

evals=1000,

tot_iters=50,

bounds=[(0, 1), (0, 1)],

cv_iters=20,

cname=”WeightPoly”,

name=”Weighted Polynomial”,

a=0.3,

b=0.6

)

cvi.compare()

which outputs the same results as before.

⬇

| No CVs | With CVs

———+———–+———–

Mean | 0.400 | 0.400

Variance | 5.826e-08 | 5.025e-08

St Dev | 2.414e-04 | 2.242e-04

VRP | | 13.746\end{result}

\par\noindentNote that in the \verb|classic_integrate| function, \verb|bounds|
is not optional as it is needed in order to define the dimension of the
integral. (In contrast, the \verb|bounds| argument is optional for the

\verb|CVIntegrator| class, since \verb|CVIntegrator| determines the dimension
of the integral from the \verb|Function| class being passed, and then sets the
limits of the integration from 0 to 1 by default.)\\\

\par\par\noindent\textbf{Specifying Control Variate Iterations}: There are
multiple valid arguments that can be passed to \verb|cv_iters| for both the
\verb|CVIntegrator| class and \verb|classic_integrate| which we list here.

\par\begin{itemize}

\itemFor using one control variate, pass an integer representing the iteration
to use.

\itemFor multiple control variates, pass a list of integers representing the
iterations to use.

\itemThe string ‘all’ will use every iteration.

\itemThe string \verb|all\itemThe previous result can be shifted by using
\verb|all\itemThe string \verb|auto1| will estimate the best single CV to use
by sampling each possibility using \verb|auto1_neval| samples specified by the
user.

\end{itemize}

\par\noindent\textbf{Manual Use of Function Class}: To access the function
call manually, use \verb|function| or \verb|f|. Using the second example, one
can run

\par\begin{minted}{python}

wpoly.f([1.2, 1.5]) # returns 1.3319999999999999

wpoly.f([0.2, 1], [0.8, 0.8], [1, 2]) # returns array([0.612, 0.672, 1.5])

This wraps around the private, vectorized ``_function—/``_f— used by Vegas.

##  5 Summary and Outlook

The MC method is widely used in many disciplines, and the issue of variance
reduction in MC estimates has existed as long as the method. A number of
techniques have been developed to reduce the variance of such MC estimates
relative to brute force sampling, including importance sampling and control
variates. In this paper, we combined these two methods by leveraging the
ability of the importance sampling method implemented in Vegas to
automatically provide viable control variate candidates. We demonstrate and
quantify the resulting variance reduction using several benchmark functions
considered previously in a similar context [6]. Our new approach was able to
reduce the variance by ${\mathcal{O}}$(1-50)% using 1 CV and
${\mathcal{O}}$(1-60)% using 2 CVs.

The reduced variance comes at the cost of some computational time. In our
experiments, the main source of additional computational overhead in our
approach is the computation of the control variates for all the datapoints,
using the intermediate Vegas distributions. The inherent theoretical
advantages of our scheme should always be weighed against the possibility to
improve the precision by simply increasing the number of data points. The
recommended use case of our technique is in situations where the latter
approach is comparatively less effective, which will be the case if computing
the integrand $f$ is sufficiently computationally expensive. Potential
examples of such situations in particle physics include inference involving
the slow simulation of the experimental apparatus, or the presence of
additional convolutions in the definition of the integrand $f(\mathbf{x})$
itself, e.g., due to transfer functions.

There are several possible avenues for refining the technique proposed in this
paper. The Vegas algorithm has a few parameters which control how the sampling
distributions are iteratively adapted. The effect of these parameters on the
performance of the Vegas-based control variates can be studied in order to
construct good controls. Furthermore, the Vegas grid adaptation algorithm
could potentially be modified for the specific purpose of providing good
control variates. There are other possible avenues for reducing the variance
of a MC estimate which leverage deep learning [22, 23] or specific domain
knowledge about the integrand [24]. These ideas are left for future work.

## Acknowledgements

The idea for this work arose during the Summer 2022 workshop “Interplay
between Fundamental Physics and Machine Learning” at the Aspen Center for
Physics, which is supported by National Science Foundation grant PHY-1607611.
The authors would like to thank the Aspen Center for Physics for hospitality
during the summer of 2022.

##### Funding information

This work is supported in parts by US DOE DE-SC0024407 and DOE DE-SC0022148.
JS is supported by University of Kansas Research Excellence Initiative Award
and US DOE AI-HEP grant. SM and PS are partially supported by the U.S.
Department of Energy, Office of Science, Office of High Energy Physics
QuantISED program under the grants “HEP Machine Learning and Optimization Go
Quantum”, Award Number 0000240323, and “DOE QuantiSED Consortium QCCFP-
QMLQCF”, Award Number DE-SC0019219. This manuscript has been authored by Fermi
Research Alliance, LLC under Contract No. DEAC02-07CH11359 with the U.S.
Department of Energy, Office of Science, Office of High Energy Physics.

## Code and data availability

The code and data that support the findings of this study are openly available
at the following URL: https://github.com/crumpstrr33/covvvr.

##  Appendix A Multiple Control Variates

For $N_{\text{CV}}>1$ control variates $g_{i}(\mathbf{x})$, the integrand is
modified in analogy to (11)

$\displaystyle f_{c}({\mathbf{x}})$
$\displaystyle=f({\mathbf{x}})+\sum_{i=1}^{N_{\text{CV}}}c_{i}\mathopen{}\mathclose{{}\left(g_{i}({\mathbf{x}})-E[g_{i}]}\right)$
(32)

and its variance is:

$\displaystyle\text{Var}\mathopen{}\mathclose{{}\left[f_{c}}\right]$
$\displaystyle=\text{Var}\mathopen{}\mathclose{{}\left[f({\mathbf{x}})+\sum_{i=1}^{N_{\text{CV}}}c_{i}(g_{i}({\mathbf{x}})-E[g_{i}])}\right]$
(33)
$\displaystyle=\text{Var}\mathopen{}\mathclose{{}\left[f}\right]+2\text{Cov}\mathopen{}\mathclose{{}\left(f,\sum_{i=1}^{N_{\text{CV}}}c_{i}(g_{i}({\mathbf{x}})-E[g_{i}])}\right)+\text{Var}\mathopen{}\mathclose{{}\left[\sum_{i=1}^{N_{\text{CV}}}c_{i}(g_{i}({\mathbf{x}})-E[g_{i}])}\right]$
$\displaystyle=\text{Var}\mathopen{}\mathclose{{}\left[f}\right]+2\text{Cov}\mathopen{}\mathclose{{}\left(f,\sum_{i=1}^{N_{\text{CV}}}c_{i}g_{i}}\right)+\text{Var}\mathopen{}\mathclose{{}\left[\sum_{i=1}^{N_{\text{CV}}}c_{i}g_{i}}\right]$
$\displaystyle=\text{Var}\mathopen{}\mathclose{{}\left[f}\right]+2\sum_{i=1}^{N_{\text{CV}}}c_{i}\text{Cov}\mathopen{}\mathclose{{}\left(f,g_{i}}\right)+\sum_{ij}^{N_{\text{CV}}}c_{i}c_{j}\text{Cov}\mathopen{}\mathclose{{}\left(g_{i},g_{j}}\right)$

where
$\text{Cov}\mathopen{}\mathclose{{}\left(g_{i},g_{i}}\right)=\text{Var}\mathopen{}\mathclose{{}\left[g_{i}}\right]$.
Taking derivatives with respect to the coefficients $c_{j}$ gives

$\displaystyle\partialderivative{\text{Var}\mathopen{}\mathclose{{}\left[f_{c}}\right]}{c_{j}}=2\text{Cov}\mathopen{}\mathclose{{}\left(f,g_{j}}\right)+2\sum_{i=1}^{N_{\text{CV}}}c_{i}\text{Cov}\mathopen{}\mathclose{{}\left(g_{i},g_{j}}\right).$
(34)

Setting that equal to zero implies

$\displaystyle\sum_{i=1}^{N_{\text{CV}}}c_{i}\text{Cov}\mathopen{}\mathclose{{}\left(g_{i},g_{j}}\right)=-\text{Cov}\mathopen{}\mathclose{{}\left(f,g_{j}}\right).$
(35)

If we let $\displaystyle
A_{j}=-\text{Cov}\mathopen{}\mathclose{{}\left(f,g_{j}}\right)$ and
$\displaystyle
B_{ij}=\text{Cov}\mathopen{}\mathclose{{}\left(g_{i},g_{j}}\right)$ then

$\displaystyle\sum_{i=1}^{N_{\text{CV}}}B_{ij}c_{i}=A_{j}\qquad\text{or}\qquad\mathbf{c}=\mathbf{B}^{-1}\mathbf{A}.$
(36)

## References

  * [1] J. Liu,  _Monte Carlo Strategies in Scientific Computing_ ,  Springer Series in Statistics. Springer New York,  ISBN 9780387952307 (2001).
  * [2] L. Lyons,  _Statistics for Nuclear and Particle Physicists_ ,  Cambridge University Press,  ISBN 9780521379342 (1989).
  * [3] X. Wang,  _Variance reduction techniques and quasi-monte carlo methods_ ,  Journal of Computational and Applied Mathematics 132(2), 309 (2001),  https://doi.org/10.1016/S0377-0427(00)00331-9.
  * [4] D. Kroese, T. Taimre and Z. Botev,  _Handbook of Monte Carlo Methods_ ,  Wiley Series in Probability and Statistics. Wiley,  ISBN 9780470177938 (2011).
  * [5] M. B. Giles,  _Multilevel monte carlo path simulation_ ,  Operations Research 56(3), 607 (2008),  10.1287/opre.1070.0496,  https://doi.org/10.1287/opre.1070.0496.
  * [6] C. Gao, J. Isaacson and C. Krause,  _i-flow: High-dimensional Integration and Sampling with Normalizing Flows_ ,  Mach. Learn. Sci. Tech. 1(4), 045023 (2020),  10.1088/2632-2153/abab62,  2001.05486.
  * [7] R. Y. Rubinstein and R. Marcus,  _Efficiency of multivariate control variates in monte carlo simulation_ ,  Operations Research 33(3), 661 (1985),  10.1287/opre.33.3.661,  https://doi.org/10.1287/opre.33.3.661.
  * [8] C. Kwon and J. D. Tew,  _Strategies for combining antithetic variates and control variates in designed simulation experiments_ ,  Management Science 40(8), 1021 (1994),  10.1287/mnsc.40.8.1021,  https://doi.org/10.1287/mnsc.40.8.1021.
  * [9] W.-N. Yang and W.-W. Liou,  _Combining antithetic variates and control variates in simulation experiments_ ,  ACM Trans. Model. Comput. Simul. 6(4), 243–260 (1996),  10.1145/240896.240899.
  * [10] R. Kawai,  _Adaptive importance sampling and control variates_ ,  Journal of Mathematical Analysis and Applications 483(1), 123608 (2020),  https://doi.org/10.1016/j.jmaa.2019.123608.
  * [11] R. Leluc, F. Portier, J. Segers and A. Zhuman,  _A quadrature rule combining control variates and adaptive importance sampling_ (2022), 2205.11890.
  * [12] P. Lepage,  _gplepage/vegas: vegas version 5.4.1_ ,  10.5281/zenodo.7933043 (2023).
  * [13] G. P. Lepage,  _A New Algorithm for Adaptive Multidimensional Integration_ ,  J. Comput. Phys. 27, 192 (1978),  10.1016/0021-9991(78)90004-9.
  * [14] G. P. Lepage,  _Adaptive multidimensional integration: VEGAS enhanced_ ,  J. Comput. Phys. 439, 110386 (2021),  10.1016/j.jcp.2021.110386,  2009.05112.
  * [15] B. Hadzhiyska _et al._ ,  _Mitigating the noise of DESI mocks using analytic control variates_ (2023),  2308.12343.
  * [16] T. Bhattacharya, S. Lawrence and J.-S. Yoo,  _Control variates for lattice field theory_ (2023),  2307.14950.
  * [17] J. DeRose, S.-F. Chen, N. Kokron and M. White,  _Precision redshift-space galaxy power spectra using Zel’dovich control variates_ ,  JCAP 02, 008 (2023),  10.1088/1475-7516/2023/02/008,  2210.14239.
  * [18] N. Chartier and B. D. Wandelt,  _Bayesian control variates for optimal covariance estimation with pairs of simulations and surrogates_ ,  Monthly Notices of the Royal Astronomical Society 515(1), 1296 (2022),  10.1093/mnras/stac1837,  2204.03070.
  * [19] N. Kokron, S.-F. Chen, M. White, J. DeRose and M. Maus,  _Accurate predictions from small boxes: variance suppression via the zel’dovich approximation_ ,  Journal of Cosmology and Astroparticle Physics 2022(09), 059 (2022),  10.1088/1475-7516/2022/09/059,  2205.15327.
  * [20] N. Chartier, B. Wandelt, Y. Akrami and F. Villaescusa-Navarro,  _CARPool: fast, accurate computation of large-scale structure statistics by pairing costly and cheap cosmological simulations_ ,  Mon. Not. Roy. Astron. Soc. 503(2), 1897 (2021),  10.1093/mnras/stab430,  2009.08970.
  * [21] N. Chartier and B. D. Wandelt,  _CARPool covariance: fast, unbiased covariance estimation for large-scale structure observables_ ,  Mon. Not. Roy. Astron. Soc. 509(2), 2220 (2021),  10.1093/mnras/stab3097,  2106.11718.
  * [22] R. Wan, M. Zhong, H. Xiong and Z. Zhu,  _Neural control variates for variance reduction_ (2019), 1806.00159.
  * [23] T. Müller, F. Rousselle, J. Novák and A. Keller,  _Neural control variates_ ,  CoRR abs/2006.01524 (2020),  2006.01524.
  * [24] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli and M. Zaro,  _The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations_ ,  JHEP 07, 079 (2014),  10.1007/JHEP07(2014)079,  1405.0301.