FERMILAB-PUB-20-010-T

i-flow: High-dimensional Integration and Sampling with Normalizing Flows

Christina Gao1, Joshua Isaacson1, and Claudius Krause1

1 Theoretical Physics Department, Fermi National Accelerator Laboratory,
Batavia, IL, 60510, USA

In many fields of science, high-dimensional integration is required. Numerical
methods have been developed to evaluate these complex integrals. We introduce
the code i-flow, a python package that performs high-dimensional numerical
integration utilizing normalizing flows. Normalizing flows are machine-
learned, bijective mappings between two distributions. i-flow can also be used
to sample random points according to complicated distributions in high
dimensions. We compare i-flow to other algorithms for high-dimensional
numerical integration and show that i-flow outperforms them for high
dimensional correlated integrals. The i-flow code is publicly available on
gitlab at https://gitlab.com/i-flow/i-flow.

##  I Introduction

Simulation based on first principles is an important practice, because it is
the only way that a theoretical model can be checked against experiments or
real-world data. In high-energy physics (HEP) experiments, a thorough
understanding of the properties of known physics forms the basis of any
searches that look for new effects. This can only be achieved by an accurate
simulation, which in many cases boils down to performing an integral and
sampling from it. Often high-dimensional phase space integrals with non-
trivial correlations between dimensions are required in important theory
calculations. Monte-Carlo (MC) methods still remain as the most important
techniques for solving high-dimensional problems across many fields, including
for instance: biology Hobolth2008 ; MC.applications , chemistry
1955JChPh..23..356R , astronomy MacGillivray1982 , medical physics
2006PMB….51R.287R , finance jackel2002monte  and image rendering
kalos2008monte . In high-energy physics, all analyses at the Large Hadron
Collider (LHC) rely strongly on multipurpose Monte Carlo event generators
Webber:1986mc ; Buckley:2011ms  for signal or background prediction. However,
the extraordinary performance of the experiments requires an amount of
simulated data that soon cannot be delivered with current algorithms and
computational resources Buckley:2019wov ; ATLASCS .

A main endeavour in the field of MC methods is to improve the error estimate.
In particular, stratified sampling — dividing the integration domain in sub-
domains, and importance sampling — sampling from non-uniform distributions
James:1968gu  are two ways of reducing the variance. Currently, the most
widely used numerical algorithm that exploits importance sampling is the VEGAS
algorithm Lepage:1977sw ; Lepage:1980dq . But VEGAS assumes the
factorizability of the integrand, which can be a bad approximation if the
variables have complex correlations amongst one another. Foam Jadach:2002kn
is a popular alternative that tries to address this issue. It uses an adaptive
strategy to attempt to model correlations, but requires exponentially large
samples in high dimensions.

Lately, the burgeoning field of machine learning (ML) has brought new
techniques into the game. For the following discussion, we restrict ourselves
to focus on progress made in the field of high-energy physics, see
Bourilkov:2019yoi  for a recent review. However, these techniques are also
widely applied in other areas of research. Concerning event generation,
Bendavid:2017zhk  used boosted decision trees and generative adversarial
networks (GANs) to improve MC integration. Reference Klimek:2018mza  proposed
a novel idea that uses a dense neural network (DNN) to learn the phase space
directly and shows promising results. In principle, once the neural
network(NN)-based algorithm for MC integration is trained, one can invert the
network and use it for sampling. However, the inversion of the NN requires
evaluating its Jacobian, which incurs a computational cost that scales as
$\mathcal{O}(D^{3})$ for $D$-dimensional integrals 111An $N$ particle final
state phase space is a $D\approx 4N-3$ dimensional integral, when including
recursive multichannel selection in the integral.. Therefore, it is extremely
inefficient to use a standard NN-based algorithm for sampling.

In addition to generating events from scratch, it is possible to generate
additional events from a set of precomputed events. References Otten:2019hhl ;
Hashemi:2019fkn ; DiSipio:2019imz ; Butter:2019cae ; Carrazza:2019cnt ;
SHiP:2019gcl ; Butter:2019eyo  used GANs and Variational Autoencoders (VAEs)
to achieve this goal. While their work is promising, they have a few
downsides. The major advantage of this approach is the drastic speed
improvement over standard techniques. They report improvements in generation
of a factor around 1000. However, this approach requires a significant number
of events already generated which may be cost prohibitive for interesting,
high-multiplicity problems. Furthermore, these approaches can only generate
events similar to those already generated. Therefore, this would not improve
the corners of distributions Matchev:2020tbw  and can even result in incorrect
total cross-sections. Yet another approach to speed up event generation is to
use NN as interpolator and learn the Matrix Element Bishara:2019iwh .

Our goal is to explore NN architectures that allow both efficient MC
integration and sampling. A ML algorithm based on _normalizing flows_ (NF)
provides such a candidate. The idea was first proposed by _non-linear
independent components estimation_ (NICE) DBLP:journals/corr/DinhKB14 ;
Dinh2016DensityEU , and generalized in rezende2015variational ;
DBLP:journals/corr/abs-1808-03856 ; 2019arXiv190604032D , for example. They
introduced _coupling layers_ (CL) allowing the inclusion of NNs in the
construction of a bijective mapping between the target and initial
distributions such that the $\mathcal{O}(D^{3})$ evaluation of the Jacobian
can be reduced to an analytic expression. This expression can now be evaluated
in $\mathcal{O}(D)$ time. These techniques have also been combined with Markov
Chain Monte Carlo methods, showing promising results 2017arXiv170607561S ;
2017arXiv171109268L ; 2019arXiv190303704H .

Our contribution is a complete, openly available implementation of normalizing
flows into TensorFlow tensorflow2015-whitepaper , to be used for any high-
dimensional integration problem at hand. Our code includes the original
proposal of DBLP:journals/corr/abs-1808-03856  and the additions of
2019arXiv190604032D . We further include various different loss functions,
based on the class of $f$-divergences 2014ISPL…21…10N . The paper is organized
in the following way. The basic principles of MC integration and importance
sampling are reviewed in Section II. In Section III, we review the concept of
normalizing flows and work done on CL-based flow by
DBLP:journals/corr/DinhKB14 ; Dinh2016DensityEU ;
DBLP:journals/corr/abs-1808-03856 ; 2019arXiv190604032D . We investigate the
minimum number of CLs required to capture the correlations between every other
input dimension. Section IV sets up the stage for a comparison between our
code, VEGAS, and Foam on various trial functions, of which we give results in
Section V. This comparison is based on several criteria, allowing a potential
user to judge whether it might be worth trying out. Section VI contains our
conclusion and outlook.

##  II Monte Carlo Integrators

While techniques exist for accurate one-dimensional integration, such as
double exponential integration takahasi1974double , using them for high
dimensional integrals requires repeated evaluation of one dimensional
integrals. This leads to an exponential growth in computation time as a
function of the number of dimensions. This is often referred to as the _curse
of dimensionality_. In other words, when the dimensionality of the integration
domain increases, the points become more and more sparse and no statistically
significant statement can be made without increasing the number of points
exponentially. This can be seen in the ratio of the volume of a
$D$-dimensional hypersphere to the $D$-dimensional hypercube, which vanishes
as $D$ goes to infinity. However, Monte-Carlo techniques are statistical in
nature and thus always converge as $1/\sqrt{N}$ for any number of dimensions.

Therefore, MC integration is the most important technique in solving high-
dimensional integrals numerically. The naïve MC approach samples uniformly on
the integration domain ($\Omega$). Given $N$ uniform samples, the integral of
$f\left(x\right)$ can be approximated by,

$I\approx\frac{V}{N}\sum_{i=1}^{N}f(x_{i})\equiv V\,\langle f\rangle_{x}\,,$
(1)

and the uncertainty is determined by the standard deviation of the mean,

$\sigma_{I}=\sqrt{\text{Var}}\approx
V\,\sqrt{\frac{\langle{f}^{2}\rangle_{x}-\langle f\rangle_{x}^{2}}{N-1}}\;,$
(2)

where $V$ is the volume encompassed by $\Omega$ and $\langle\ \rangle_{x}$
indicates that the average is taken with respect to a uniform distribution in
$x$. While this works for simple or low-dimensional problems, it soon becomes
inefficient for high-dimensional problems. This is what our work is concerned
with. In particular, we are going to focus on improving current methods for MC
integration that are based on importance sampling.

In importance sampling, instead of sampling from an uniform distribution, one
samples from a distribution $g(x)$ that ideally has the same shape as the
integrand $f(x)$. Using the transformation ${\rm d}x={\rm d}G(x)/g(x)$, with
$G(x)$ the cumulative distribution function of $g(x)$, one obtains

$I=\int_{\Omega}\frac{f(x)}{g(x)}\,\text{d}G(x)=V\,\langle
f/g\rangle_{G}\;,\qquad\sigma_{I}=V\,\sqrt{\frac{\langle(f/g)^{2}\rangle_{G}-\langle
f/g\rangle_{G}^{2}}{N-1}}\;.$ (3)

In the ideal case when $g(x)\to f(x)/I$, Eq. (3) would be estimated with
vanishing uncertainty. However, this requires already knowing the analytic
solution to the integral! The goal is thus to find a distribution $g(x)$ that
resembles the shape of $f(x)$ most closely, while being integrable and
invertible such as to allow for fast sampling. We review the current MC
integrators that are widely used, especially in the field of high-energy
physics.

VEGAS Lepage:1977sw ; Lepage:1980dq  approximates all 1-dimensional
projections of the integrand using a histogram and an adaptive algorithm. This
algorithm adjusts the bin widths such that the area of the bins are roughly
equal. To sample a random point from VEGAS can be done in two steps. First,
select a bin randomly for each dimension. Second, sample a point from each bin
according to a uniform distribution. However, this algorithm is limited
because it assumes that the integrand factorizes, i.e.

$f\left(\vec{x}\right)=f_{1}\left(x_{1}\right)\cdots f_{D}\left(x_{D}\right),$
(4)

where $f\colon\mathbb{R}^{D}\mapsto\mathbb{R}$ and
$f_{i}\colon\mathbb{R}\mapsto\mathbb{R}$. High-dimensional integrals with non-
trivial correlations between integration variables, that are often needed for
LHC data analyses, cannot be integrated efficiently with the VEGAS algorithm
(c.f. Hoeche:2019rti ). The resulting uncertainty can be reduced further by
applying stratified sampling, in addition to the VEGAS algorithm, after the
binning VEGAS .

Foam Jadach:2002kn  uses a cellular approximation of the integrand and is
therefore able to learn correlations between the variables. In the first phase
of the algorithm, the so-called exploration phase, the cell grid is built by
subsequent binary splits of existing cells. Since the first cell consists of
the full integration domain, all regions of the integration space are explored
by construction. The second phase of the algorithm uses this grid to generate
points either for importance sampling or as an event generator. In this work
we use the implementation of Hoeche:foam , which implemented an additional
reweighting of the cells at the end of the optimization.

However, both Foam and VEGAS are based on histograms, whose edge effects would
be detrimental to numerical analyses that demand high precision. As we will
explain below, our code i-flow uses a spline approximation which does not
suffer from these effects. These edge effects are an important source of
uncertainty for high-precision physics Dulat:2017prg .

##  III Importance Sampling with Normalizing Flows

As we detailed in the previous section, importance sampling requires finding
an approximation $g(x)$ that can easily be integrated and subsequently
inverted, so that we can use it for sampling. Mathematically, this corresponds
to a coordinate transformation with an inverse Jacobian determinant that is
given by $g(x)$. General ML algorithms incorporate NNs in learning the
transformation, which inevitably involve evaluating the Jacobian of the NNs.
This results in inefficient sampling. Coupling Layer-based Normailzing Flow
algorithms precisely circumvent this problem. To begin, let us review the
concept of a normalizing flow (NF).

Let $c_{k}$, with $k=1,...,K$, be a series of bijective mappings on the random
variable $\vec{x}$:

$\vec{x}_{K}=\;c_{K}(c_{K-1}(\cdots c_{2}(c_{1}(\vec{x})))).$ (5)

Based on the chain rule, the output $\vec{x}_{K}$’s probability distribution,
$g_{K}$, can be inferred given the base probability distribution $g_{0}$ from
which $\vec{x}$ is drawn:

${g}_{K}(\vec{x}_{K})=\;{g}_{0}(\vec{x}_{0})\prod_{k=1}^{K}\,\left|\frac{\partial
c_{k}(\vec{x}_{k-1})}{\partial\vec{x}_{k-1}}\right|^{-1}\;,\quad\text{where}\quad\begin{cases}\vec{x}_{0}=\vec{x}\\\
\vec{x}_{k}=c_{k}(\vec{x}_{k-1})\end{cases}\;.$ (6)

One sees that the target and base distributions are related by the inverse
Jacobian determinant of the transformation. For practical uses, the Jacobian
determinant must be easy to compute, restricting the allowed functional forms
of $c_{k}$. However, with the help of _coupling layers_ , first proposed by
DBLP:journals/corr/DinhKB14 ; Dinh2016DensityEU , then generalized by
DBLP:journals/corr/abs-1808-03856 ; 2019arXiv190604032D , one can incorporate
NNs into the construction of $c_{k}$, thus greatly enhancing the level of
complexity and expressiveness of NF without introducing any intractable
Jacobian computations.

Figure 1: Structure of a Coupling Layer. $m$ is the output of a neural network
and defines the Coupling Transform, $C$, that will be applied to $x_{B}$. See
Eqs. (7), (8), (9) for the mathematical description of a Coupling Layer.

Figure 1 shows the basic structure of a coupling layer, which is a special
design of the bijective mapping $c$. For each map, the input variable
$\vec{x}=\\{x_{1},..,x_{D}\\}$ is partitioned into two subsets, $\vec{x}_{A}$
and $\vec{x}_{B}$ which can be determined arbitrarily so long as neither is
the empty set. This arbitrary partitioning will be referred to as a _masking_.
Without loss of generality, one simple partitioning is given by
$\vec{x}_{A}=\\{x_{1},..,x_{d}\\}$ and $\vec{x}_{B}=\\{x_{d+1},..,x_{D}\\}$.
Different maskings can be achieved via permutations of the simple example
above. Under the bijective map, $C$, the resulting variable transforms as

$\begin{split}x_{A}^{\prime}&=x_{A},\quad\quad\quad\quad\quad\quad
A\in[1,d]\;,\\\ x_{B}^{\prime}&=C(x_{B};m(\vec{x}_{A})),\quad
B\in[d+1,D]\;.\end{split}$ (7)

The NN takes $x_{A}$ as inputs and outputs $m(\vec{x}_{A})$ that represents
the parameters of the invertible “Coupling Transform”, $C$, that will be
applied to $x_{B}$. We detail various choices for $C$, like piecewise linear,
piecewise quadratic, or piecewise rational quadratic spline functions in
Appendix A. The inverse map is given by

$\begin{split}x_{A}&=x_{A}^{\prime}\;,\\\
x_{B}&=C^{-1}(x_{B}^{\prime};m(\vec{x}_{A}^{\prime}))=C^{-1}(x_{B}^{\prime};m(\vec{x}_{A}))\;,\end{split}$
(8)

which leads to the simple Jacobian

$\left|\frac{\partial
c(\vec{x})}{\partial\vec{x}}\right|^{-1}=\left|\left(\begin{array}[]{cc}\vec{1}&0\\\
\frac{\partial C}{\partial m}\frac{\partial
m}{\partial\vec{x}_{A}}&\frac{\partial
C}{\partial\vec{x}_{B}}\end{array}\right)\right|^{-1}=\left|\frac{\partial
C(\vec{x}_{B};m(\vec{x}_{A}))}{\partial\vec{x}_{B}}\right|^{-1}\;.$ (9)

Note that Eq. (9) does not require the computation of the gradient of
$m(\vec{x}_{A})$, which would scale as $\mathcal{O}(D^{3})$ with $D$ the
number of dimensions. In addition, taking $\partial C/\partial\vec{x}_{B}$ to
be diagonal further reduces the computation complexity of the determinant to
be linear with respect to the dimensionality of the problem. Linear scaling
makes this approach tractable even for high dimensional problems. In summary,
the NN learns the parameters of a transformation and not the transformation
itself, thus the Jacobian can be calculated analytically.

To construct a complete Normalizing Flow, one simply compounds a series of
Coupling Layers with the freedom of choosing any of the outputs of the
previous layer to be transformed in the subsequent layer. We show in Sec III.1
that $2\lceil\log_{2}\mathrm{D}\rceil$ number of Coupling Layers are required
in order to express all non-separable structures of the integrand.

###  III.1 Number of Coupling Layers

The minimum number of coupling layers required to capture all possible
correlations between every dimension of the integration variable,
$n_{\text{min}}$, depends on the dimensionality of the integral, $D$
DBLP:journals/corr/abs-1808-03856 . In the cases of $D=2$ and $D=3$, each
dimension is transformed once based on the other dimension(s) and thus
$n_{\text{min}}=2$ and $n_{\text{min}}=3$, respectively. This way of counting
$n_{\text{min}}$ could be generalized to higher $D$. In fact, this is what
_autoregressive flows_ are based on 2017arXiv170507057P . Here we show that
the number of coupling layers required to capture all the correlations is
$2\lceil\log_{2}D\rceil$ for $D>5$, and $D$ layers for $D\leq 5$. This can be
considered the minimum number of layers required in order to capture all
correlations, adding an additional layer will not add any new information, and
similar effects should be achieved with increasing the depth of the network
associated with each layer. On the other hand, this can be considered the
maximum number of layers needed to capture all the correlations. If a function
has fewer correlations, then all the correlations can be captured with less
than $2\lceil\log_{2}D\rceil$.

######  Theorem 1.

Given a set of correlated random variables $x$, if a transformation exists
that takes the variables $x$ to $z$, such that the correlation between the
variables $z$ is zero, then a composition of normalizing flows can create such
a transformation. Given a set of infinitely wide NNs that are universal
function approximators, and requiring that all variables are transformed equal
number of times, it is possible to represent all the correlations between
variables in a normalizing flow using $2\lceil\log_{2}D\rceil$ layers for
$D>5$. When $D\leq 5$ it is possible to represent all correlations with $D$
layers.

###### Proof.

Given the random variables $x_{1},\ldots x_{D}$, with means
$\mu_{1},\ldots,\mu_{D}$ and joint probability distribution
$f\left(x_{1},\ldots,x_{D}\right)$, the correlation between all the variables
is given by:

$\langle\left(x_{1}-\mu_{1}\right)\ldots\left(x_{D}-\mu_{D}\right)\rangle=\int_{0}^{1}dx_{1}\ldots
dx_{D}\left(x_{1}-\mu_{1})\ldots(x_{D}-\mu_{D}\right)f\left(x_{1},\ldots,x_{D}\right).$
(10)

Using two layers of a normalizing flow network, which can be seen as a
universal function approximator, defines a transformation $T:x\mapsto y$, with
the bounds of integration being mapped such that
$y_{i}\left(T\left(x=0\right)\right)=0$ and
$y_{i}\left(T\left(x=1\right)\right)=1\ \forall i\in[1,D]$, with the sets
$\\{y_{a}\\}=\\{y_{i}|\,i\equiv 1\pmod{2},\,i\in[1,D]\\}$ and
$\\{y_{b}\\}=\\{y_{i}|\,i\equiv 0\pmod{2},\,i\in[1,D]\\}$, such that
$f\left(x_{1},\ldots x_{D}\right)\mapsto
g\left(\\{y_{a}\\}\right)h\left(\\{y_{b}\\}\right)$, and with the Jacobian
$J(y,x)$. The transformation also maps the means: $\mu\mapsto\mu^{y}$. This
decomposition is possible following from the arguments of
(papamakarios2019normalizing, , Sec. 2.2). Applying the transformation to Eq.
10 gives:

$\langle\left(x_{1}-\mu_{1}\right)\ldots\left(x_{D}-\mu_{D}\right)\rangle=\int_{0}^{1}dy_{1}\ldots
dy_{D}J(y,x)\left(y_{1}-\mu^{y}_{1})\ldots(y_{D}-\mu^{y}_{D}\right)g\left(\\{y_{a}\\}\right)h\left(\\{y_{b}\\}\right).$
(11)

If we now consider the correlation between the variables $y$, we obtain:

$\displaystyle\langle\left(y_{1}-\mu^{y}_{1}\right)\ldots\left(y_{D}-\mu^{y}_{D}\right)\rangle=$
$\displaystyle\int_{0}^{1}dy_{1}\ldots
dy_{D}\left(y_{1}-\mu^{y}_{1}\left)\ldots\right(y_{D}-\mu^{y}_{D}\right)g\left(\\{y_{a}\\}\right)h\left(\\{y_{b}\\}\right)$
$\displaystyle=$
$\displaystyle\int_{0}^{1}\prod_{\\{y_{a}\\}}\left(dy_{a}\left(y_{a}-\mu^{y}_{a}\right)\right)g\left(\\{y_{a}\\}\right)\times\int_{0}^{1}\prod_{\\{y_{b}\\}}\left(dy_{b}\left(y_{b}-\mu^{y}_{b}\right)\right)h\left(\\{y_{b}\\}\right)$
$\displaystyle=$
$\displaystyle\bigg{\langle}\prod_{\\{y_{a}\\}}\bigg{(}y_{a}-\mu^{y}_{a}\bigg{)}\bigg{\rangle}\bigg{\langle}\prod_{\\{y_{b}\\}}\bigg{(}y_{b}-\mu^{y}_{b}\bigg{)}\bigg{\rangle}.$
(12)

The result of the transformation shows that the variables $\\{y_{a}\\}$ are
now not correlated with $\\{y_{b}\\}$. We can construct a subsequent
transformation $T^{\prime}:y\mapsto z$, with
$z_{i}\left(T^{\prime}\left(y=0\right)\right)=0$ and
$z_{i}\left(T^{\prime}\left(y=1\right)\right)=1\ \forall i\in[1,D]$, and the
sets $\\{z_{a}\\}=\\{z_{i}|\,i\equiv 1\pmod{4},\,i\in[1,D]\\}$,
$\\{z_{b}\\}=\\{z_{i}|\,i\equiv 2\pmod{4},\,i\in[1,D]\\}$,
$\\{z_{c}\\}=\\{z_{i}|\,i\equiv 3\pmod{4},\,i\in[1,D]\\}$, and
$\\{z_{d}\\}=\\{z_{i}|\,i\equiv 0\pmod{4},\,i\in[1,D]\\}$, such that
$g\left(\\{y_{a}\\}\right)h\left(\\{y_{b}\\}\right)\mapsto
g^{\prime}\left(\\{z_{a},z_{b}\\}\right)h^{\prime}\left(\\{z_{c},z_{d}\\}\right)$
with the constraint that such a transformation does not introduce new
correlations between the variables that have already been decorrelated. In
other words, the composition of $T$ and $T^{\prime}$ can be defined as a
transformation $T^{\prime\prime}:x\mapsto z$, with
$z_{i}(T^{\prime\prime}(x)=0)=0$ and $z_{i}(T^{\prime\prime}(x)=1)=1\ \forall
i\in[1,D]$, such that:

$f\left(x_{1},\ldots,x_{N}\right)\mapsto
g_{1}\left(\\{z_{a}\\}\right)g_{2}\left(\\{z_{b}\\}\right)g_{3}\left(\\{z_{c}\\}\right)g_{4}\left(\\{z_{d}\\}\right),$

and the means are mapped from $\mu$ to $\mu^{z}$. Thus, the correlation
between the variables $z$ is given by:

$\displaystyle\langle\left(z_{1}-\mu^{z}_{1}\right)\ldots\left(z_{D}-\mu^{z}_{D}\right)\rangle=$
$\displaystyle\int_{0}^{1}dz_{1}\ldots
dz_{D}\left(z_{1}-\mu^{z}_{1})\ldots(z_{D}-\mu^{z}_{D}\right)g_{1}\left(\\{z_{a}\\}\right)g_{2}\left(\\{z_{b}\\}\right)g_{3}\left(\\{z_{c}\\}\right)g_{4}\left(\\{z_{d}\\}\right)$
$\displaystyle=$
$\displaystyle\left\langle\prod_{\\{z_{a}\\}}\left(z_{a}-\mu^{z}_{a}\right)\right\rangle\left\langle\prod_{\\{z_{b}\\}}\left(z_{b}-\mu^{z}_{b}\right)\right\rangle\left\langle\prod_{\\{z_{c}\\}}\left(z_{c}-\mu^{z}_{c}\right)\right\rangle\left\langle\prod_{\\{z_{d}\\}}\left(z_{d}-\mu^{z}_{d}\right)\right\rangle.$
(13)

The above transformations can be iterated until all the variables are
decorrelated. A method of determining the mapping for each step can be
obtained by the following procedure:

  1. 1.

Reindex the dimension numbers from $[1,D]$ to $[0,D-1]$

  2. 2.

Convert all dimensions to their binary representation, using the minimum
number of bits required to represent the number $D-1$

  3. 3.

Consider the least significant bit for each dimension, and define the
transformation as $T:x\mapsto y$, with $f(\\{x_{0}\\},\\{x_{1}\\})\mapsto
g(\\{x_{0}\\})h(\\{x_{1}\\})$, where $\\{x_{0}\\}(\\{x_{1}\\})$ is the set of
variables with a 0 (1) for the least significant bit

  4. 4.

Repeat the third step taking the next least significant bit, until the most
significant bit is reached

See Table 1 for an example of the steps above. In that example, transformation
1 would groups the first 8 dimensions in $g(\\{x_{0}\\})$ and the last 4 in
$h(\\{x_{1}\\})$ etc.

Dimension |  0 |  1 |  2 |  3 |  4 |  5 |  6 |  7 |  8 |  9 |  10 |  11
---|---|---|---|---|---|---|---|---|---|---|---|---
Transformation 1 |  0 |  1 |  0 |  1 |  0 |  1 |  0 |  1 |  0 |  1 |  0 |  1
Transformation 2 |  0 |  0 |  1 |  1 |  0 |  0 |  1 |  1 |  0 |  0 |  1 |  1
Transformation 3 |  0 |  0 |  0 |  0 |  1 |  1 |  1 |  1 |  0 |  0 |  0 |  0
Transformation 4 |  0 |  0 |  0 |  0 |  0 |  0 |  0 |  0 |  1 |  1 |  1 |  1
Table 1: Finding the unique masking to capture all correlations in a $D=12$
space, using the procedure detailed above.

The number of steps for this procedure can easily be seen to be
$\lceil\log_{2}\left(D\right)\rceil$. However, since we need two layers per
transformation to ensure that all variables are equally transformed by the
network leads to a requirement of $2\lceil\log_{2}\left(D\right)\rceil$.

In the situation of $D\leq 5$, the $i^{\text{th}}$ coupling layer can be
defined to take the $i^{\text{th}}$ variable and transform the variable, such
that it is not correlated with any other variable. This leads to a requirement
of $D$ layers. ∎

The requirement on the above theorem is that we require that a transformation
exists in order to perform the above mapping. However, even if a
transformation does not exist for the integrand itself, with the use of
importance sampling, it is only necessary to find a function $g$ which is as
close to the integrand as possible. In i-flow splines are used to create the
function $g$, and according to the Stone-Weierstrass Theorem PINKUS20001 ;
10.2307/1989788 ; 10.2307/3029337 , it is possible to represent $g$ such that
it is $\varepsilon$ close to $f$. A corollary of the Stone-Weierstrass Theorem
for $\mathbb{R}^{n}$ can be expressed as:

######  Corollary 1.

Given a function $f:\mathbb{R}^{n}\mapsto\mathbb{R}$, $\varepsilon>0$, and
$C(\mathbb{R}^{n},\mathbb{R})$: the space of all real-valued continuous
functions in $\mathbb{R}^{n}$, there exists a polynomial spline $g\in
C(\mathbb{R}^{n},\mathbb{R})$ such that:

$|f(x)-g(x)|<\varepsilon,$ (14)

for all $x\in\mathbb{R}^{n}$.

###### Proof.

The space $R^{n}$ is a subset of the spaces proved in the Stone-Weierstrass
Theorem, and thus the proof of the corollary follows directly from the Stone-
Weierstrass Theorem. ∎

Furthermore, this can be extended to a sum of piecewise polynomials, such that
any continuous and bounded function $f$ can be represented by an infinite
series of polynomials (see Theorem C from PINKUS20001 ). In i-flow, we will
consider the case of discontinuous functions, but these can be approximated as
a continuous function with a slope of $1/\varepsilon$ in the region of
discontinuity. This will lead to some difference between $f$ and $g$, but
since the goal is to find a function $g$ as close to $f$ as possible, then
this is acceptable and should still allow for high precision importance
sampling.

###  III.2 Using i-flow

Figure 2: Illustration of one step in the training of i-flow. Users need to
provide a normalizing flow network, a function $f$ to integrate, and a loss
function. $\tilde{I}$ stands for the Monte-Carlo estimate of the integral
using the sample of points $\vec{x}_{i}$, and $g(\vec{x}_{i})$ is the
probability of a given point occurring in the i-flow sampling.

The i-flow package requires three pieces of information from the user: the
function to be integrated, the normalizing flow network, and the method of
optimizing the network. Figure 2 shows schematically how one step in the
training of i-flow works. The code is publicly available on gitlab at
https://gitlab.com/i-flow/i-flow. Running the script iflow_test.py will
produce the results presented in section V.

####  III.2.1 Integrand

The function to be integrated has very few requirements on how it is
implemented in the code. Firstly, the function must accept an array of points
with shape $(n_{\text{batch}},D)$, where $n_{\text{batch}}$ is the number of
points to sample per training step. Secondly, the function must return an
array with shape $(n_{\text{batch}})$ to be used to estimate the integral.
Finally, the number of dimensions in the integral is required to be at least
2. However, one dimension can be treated as a dummy dimension integrated from
0 to 1, which will not change any result.

####  III.2.2 Normalizing flow network

A normalizing flow network consists of a series of coupling layers compounded
together. To construct each coupling layer, one needs to specify the choice of
coupling transform $C$ (cf App. A), the number of coupling layers, the masking
for each level, and the neural network $m(x_{A})$ that constitutes the
coupling transform. We provide the ability to automatically generate the
masking and number of layers according to Sec. III.1.

The neural networks $m(x_{A})$ must be provided by the user. However, we
provide examples for a dense network and the U-shape network of
DBLP:journals/corr/abs-1808-03856 . This provides the user the flexibility to
achieve the expressiveness required for their specific problem.

####  III.2.3 Optimizing the network

To uniquely define the optimization algorithm of the network, two pieces of
information are required. Firstly, the loss function to be minimized is
required. We supply a large set of loss functions from the set of
$f$-divergences, which can be found in App. B. By default, the i-flow code
uses the exponential loss function. Secondly, an optimizer needs to be
supplied. In the examples we used the ADAM optimizer kingma2014adam . However,
the code can use any optimizer implemented within TensorFlow.

####  III.2.4 Hyperparameters

The setup we presented here has several hyperparameters that can be adjusted
for better performance. However, i-flow has the flexibility for the user to
implement additional features in each section beyond what is discussed below.
This would come with additional hyperparameters as well.

The first group concerns the architecture of the NNs $m(x_{A})$. Once the
general type of network (dense or U-shape) is set, the number of layers and
nodes per layer have to be specified. In the case of the U-shape network, the
user can specify the number of nodes in the first layer and the number of
“downward” steps.

The second group of hyperparameters concerns the optimization process. Apart
from setting an optimizer (e.g. ADAM kingma2014adam ), a learning schedule
(e.g. constant or exponentially decaying), an initial learning rate, and a
loss function have to be specified. Some of these options come with their own,
additional set of hyperparameters. The number of points per training epoch and
the number of epochs have to be set as well.

The third group of hyperparameters concerns the setup of i-flow directly. As
was discussed in DBLP:journals/corr/abs-1808-03856 , there are two ways to
pass $x_{A}$ into $m(x_{A})$: either directly or with one-blob encoding.
i-flow supports both of these options. One-blob encoding
DBLP:journals/corr/abs-1808-03856  is a generalization of one-hot encoding.
The input $x_{A}$ is passed through a Gaussian kernel and several adjacent
bins are activated. If one-blob encoding is used, the number of input bins has
to be specified, the width of the Gaussian is set to the inverse of the number
of bins. Further, the type of coupling function $C(x_{B},m(x_{A}))$, the
number of output bins, the number of CLs and the maskings have to be set.

####  III.2.5 Putting it all together

The networks are trained by sampling a fixed number of points using the
current state of $g(x)$ 222Since we initialize the last layer of each network
with vanishing bias and weights, in the first sampling $g(x)$ is constant.. We
use one of the statistical divergences as a measure for how much the
distribution $g(x)$ resembles the shape of the integrand $f(x)$, and an
optimizer to minimize it. Because we can generate an infinite set of random
numbers and evaluate the target function for each of the points, this approach
corresponds to supervised learning with an infinite dataset. Drawing a new set
of points at every training epoch automatically also ensures that the networks
cannot overfit.

##  IV Integrator Comparison

To illustrate the performance of i-flow and compare it to VEGAS and Foam, we
present a set of six test functions, each highlighting a different aspect of
high-dimensional integration and sampling. These functions demonstrate how
each algorithm handles the cases of a purely separable function, functions
with correlations, and functions with non-factorizing hard cuts. In most
cases, an analytic solution to the integral is known.

The first test function is an $n$-dimensional Gaussian, serving as a sanity
check:

$f_{1}(\vec{x})=(\alpha\sqrt{\pi})^{-n}\exp{\\{-\frac{\sum_{i}(x_{i}-0.5)^{2}}{\alpha^{2}}\\}}\;.$
(15)

The result of integrating $f_{1}$ from zero to one is given by:

$\int_{0}^{1}d^{n}\vec{x}\,f_{1}(\vec{x})=\text{erf}\left(\frac{1}{2\alpha}\right)^{n}\;.$
(16)

In the following, we use $\alpha=0.2$.

The second test function is an $n$-dimensional Camel function, which would
show how i-flow learns correlations that VEGAS (without stratified sampling)
would not learn:

$f_{2}(\vec{x})=\frac{1}{2}(\alpha\sqrt{\pi})^{-n}\left(\exp{\\{-\frac{\sum_{i}(x_{i}-\tfrac{1}{3})^{2}}{\alpha^{2}}\\}}+\exp{\\{-\frac{\sum_{i}(x_{i}-\tfrac{2}{3})^{2}}{\alpha^{2}}\\}}\right)\;.$
(17)

The result of integrating $f_{2}$ from zero to one is given by:

$\int_{0}^{1}d^{n}\vec{x}\,f_{2}(\vec{x})=\left(\frac{1}{2}\left(\text{erf}\left(\frac{1}{3\alpha}\right)+\text{erf}\left(\frac{2}{3\alpha}\right)\right)\right)^{n}\;.$
(18)

In the following, we use $\alpha=0.2$.

The third case is given by

$\displaystyle\begin{aligned}
f_{3}(x_{1},x_{2})&=x_{2}^{a}\exp{\\{-w|(x_{2}-p_{2})^{2}+(x_{1}-p_{1})^{2}-r^{2}|\\}}\\\
&+(1-x_{2})^{a}\exp{\\{-w|(x_{2}-1+p_{2})^{2}+(x_{1}-1+p_{1})^{2}-r^{2}|\\}}\;.\end{aligned}$
(19)

This function has two circles with shifted centers, varying thickness and
height. Also, the function exhibits non-factorizing behavior. The integral of
$f_{3}$ between 0 and 1 can be computed numerically using Mathematica
Mathematica , which is $0.0136848\pm(5\cdot 10^{-9})$, with
$p_{1}=0.4,p_{2}=0.6,r=0.25,w=1/0.004$ and $a=3$ 333There is no known analytic
solution to this given function.

The fourth case is an annulus function with hard cuts:

$f_{4}(x_{1},x_{2})=\left\\{\begin{array}[]{cc}1&0.2<\sqrt{x_{1}^{2}+x_{2}^{2}}<0.45\\\
0&\text{else}\end{array}\right.\;.$ (20)

This function demonstrates how i-flow learns hard, non-factorizing cuts. The
result of integrating $f_{4}$ from zero to one is given by:
$\pi\left(0.45^{2}-0.2^{2}\right)=0.1625\pi$.

The fifth case is motivated by high energy physics, and is a one-loop scalar
box integral representative of an integral required for the calculation of
$gg\to gh$ in the Standard Model. This calculation is an important
contribution for the total production cross-section of the Higgs boson. As
explained in App. C, after Feynman parametrisation and sector decomposition
Binoth:2002xh , the integral of interest is given by

$\begin{split}f_{5}=&S_{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})\\\
&+S_{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})\\\
&+S_{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})\\\
&+S_{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})\\\
S_{Box}&(s_{12},s_{23},s_{1},s_{2},s_{3},s_{4},m_{1}^{2},m_{2}^{2},m_{3}^{2},m_{4}^{2})=\int_{0}^{1}dt_{1}dt_{2}dt_{3}\frac{1}{\tilde{\mathcal{F}}_{Box}^{2}}\\\
\tilde{\mathcal{F}}_{Box}&=(-s_{12})t_{2}+(-s_{23})t_{1}t_{3}+(-s_{1})t_{1}+(-s_{2})t_{1}t_{2}+(-s_{3})t_{2}t_{3}\\\
&\quad\quad+(-s_{4})t_{3}+(1+t_{1}+t_{2}+t_{3})(t_{1}m_{1}^{2}+t_{2}m_{2}^{2}+t_{3}m_{3}^{2}+m_{4}^{2})\end{split}$
(21)

The result of integrating $f_{5}$ from zero to one can be obtained through the
use of LoopTools HAHN1999153 , which gives a numerical result of
$1.9369640238\cdot 10^{-10}$ for the inputs
$s_{12}=130^{2},s_{23}=-130^{2},s_{1}=0,s_{2}=0,s_{3}=0,s_{4}=125^{2},m_{t}=175$.

As a sixth test function, we consider the polynomial

$f_{6}(x_{1},\dots,x_{n})=\sum_{i=1}^{n}-x_{i}^{2}+x_{i}.$ (22)

The result of integrating $f_{6}$ from zero to one is given by:

$\int_{0}^{1}\,dx_{1}\dots dx_{n}\,f_{6}(x_{1},\dots,x_{n})=\frac{n}{6}$ (23)

This function can easily be integrated in a high number of dimensions and,
unlike the Gaussian or Camel functions, has support in almost all of the
integration domain. It therefore does not suffer that much from the curse of
dimensionality.

Further applications to event generation of high-energy particle collisions is
discussed in Gao:2020zvv  and also in Bothmann:2020ywa . These papers
investigate using normalizing flows to improve upon phase space integration
for event simulation at particle colliders. The integral dimension for
processes with $n_{f}$ particles in the final state is $D=4n_{f}-3$. In
Gao:2020zvv , we studied processes with $n_{f}\leq 6$.

##  V Results

| Dim | VEGAS | Foam | i-flow
---|---|---|---|---
Gaussian | 2 | $\mathbf{164,436}$ | $6,259,812$ | $2,310,000$
4 | $\mathbf{631,874}$ | $24,094,679$ | $2,285,000$
8 | $\mathbf{1,299,718}$ | $>50,000,000\;\dagger$ | $3,095,000$
16 | $\mathbf{2,772,216}$ | $>50,000,000\;\dagger$ | $7,230,000$
Camel | 2 | $\mathbf{421,475}$ | $5,619,646$ | $2,225,000$
4 | $24,139,889$ | $21,821,075$ | $\mathbf{8,220,000}$
8 |  $>50,000,000\;\dagger$ |  $>50,000,000$ |  $\mathbf{19,460,000}$
16 | $993,294\;\dagger$ | $>50,000,000\;\dagger$ |  $32,145,000\;\dagger$fn:camel16
Entangled circles | 2 | $43,367,192$ | $\mathbf{17,499,823}$ | $23,105,000$
Annulus w. cuts | 2 | $4,981,080\;\dagger$ | $\mathbf{11,219,498}$ | $17,435,000$
Scalar-top-loop | 3 | $\mathbf{152,957}$ | $5,290,142$ | $685,000$
Polynomial | 18 | $42,756,678$ | $>50,000,000$ | $\mathbf{585,000}$
54 | $>50,000,000$ | $>21,505,000\;*$ | $\mathbf{685,000}$
96 | $>50,000,000\;\dagger$ | $>10,235,000\;*$ | $\mathbf{1,145,000}$
Table 2: Number of functional calls to reach a total relative uncertainty of $10^{-4}$ (for the first 11 cases) or $10^{-5}$ (for the last 3 cases). The total relative uncertainty is defined as the inverse-variance weighted combination of the uncertainties of each optimization iteration divided by the true integral value. The integrator with the fewest functional calls, which also is within 5 standard deviations of the true result, is highlighted in boldface. We set an upper cut-off of $5\cdot 10^{7}$ calls. A $\dagger$ indicates that the algorithm did not converge to the true integral value within 5 standard deviations (see tab. 3), a $*$ indicates cases where the algorithm ran out of memory before the cut-off was reached. $\alpha=0.2$ for Gaussian and Camel functions. | Dim | VEGAS | (pull) | Foam | (pull) | i-flow | (pull) | true value
---|---|---|---|---|---|---|---|---
Gaussian | 2 | $0.99925(10)$ | $0.7$ | $0.99925(10)$ | $0.6$ | $\mathbf{0.99919(10)}$ | $\mathbf{0.1}$ | $0.999186$
4 | $0.99861(10)$ | $2.4$ | $\mathbf{0.99835(10)}$ | $\mathbf{-0.2}$ | $0.99841(10)$ | $0.4$ | $0.998373$
8 | $0.99694(10)$ | $1.9$ | $0.99439(37)\;\dagger$ | $-6.4$ | $\mathbf{0.99684(10)}$ | $\mathbf{0.9}$ | $0.996749$
16 | $0.99357(10)$ | $0.6$ | $0.54986(235)\;\dagger$ | $-188$ | $\mathbf{0.99354(10)}$ | $\mathbf{0.4}$ | $0.993509$
Camel | 2 | $0.98175(10)$ | $0.9$ | $0.98163(10)$ | $-0.3$ | $\mathbf{0.98165(10)}$ | $\mathbf{-0.1}$ | $0.98166$
4 | $0.96345(10)$ | $-2.2$ | $0.96361(10)$ | $-0.5$ | $\mathbf{0.96365(10)}$ | $\mathbf{-0.02}$ | $0.963657$
8 | $0.92495(28)\;\dagger$ | $-13$ | $0.92798(19)\;\dagger$ | $-3.5$ | $\mathbf{0.92843(9)}$ | $\mathbf{-2.2}$ | $0.928635$
16 | $0.43137(9)$ | $-5001$ | $0.76921(129)\;\dagger$ | $-72$ |  $\mathbf{0.85940(9)}$ fn:camel16  | $\mathbf{-34}$ | $0.862363$
Entangled circles | 2 | $0.0136798(14)$ | $-3.6$ | $\mathbf{0.0136838(14)}$ | $\mathbf{-0.7}$ | $0.0136829(14)$ | $-1.4$ | $0.0136848$
Annulus w. cuts | 2 | $0.509813(51)$ | $-14$ | $0.510559(51)$ | $1.0$ | $\mathbf{0.510511(51)}$ | $\mathbf{0.1}$ | $0.510508$
Scalar-top-loop | 3 |  $1.93711(19)\cdot 10^{-10}$ | $0.7$ |  $\mathbf{1.93708(19)\cdot 10^{-10}}$ | $\mathbf{0.6}$ |  $1.93677(19)\cdot 10^{-10}$ | $-1.0$ |  $1.936964\cdot 10^{-10}$
Polynomial | 18 | $2.99989(3)$ | $-3.6$ | $2.99986(12)\;\dagger$ | $-1.1$ | $\mathbf{2.99997(3)}$ | $\mathbf{-1.1}$ | $3$
54 | $8.99972(19)\;\dagger$ | $-1.5$ | $9.00013(32)\;*$ | $0.4$ | $\mathbf{9.00001(9)}$ | $\mathbf{0.2}$ | $9$
96 | $0.15547(52)\;\dagger$ | $-30683$ | $16.0004(3)\;*$ | $1.7$ | $\mathbf{15.9998(2)}$ | $\mathbf{-1.2}$ | $16$
Table 3: Integral estimate and uncertainty of the runs of Table 2 together with their relative deviations (“pull”), defined in Eq. (25). A $\dagger$ indicates that the algorithm reached a cut-off of $5\cdot 10^{7}$ function calls before the target uncertainty was reached, a $*$ indicates cases where the algorithm ran out of memory before the cut-off was reached. The result with the smallest relative deviation is boldfaced. $\alpha=0.2$ for Gaussian and Camel functions. | Dim | VEGAS | Foam | i-flow
---|---|---|---|---
Gaussian | 2 | $\mathbf{7\cdot 10^{-4}}$ | $3\cdot 10^{-3}$ | $2\cdot 10^{-3}\;*$
4 | $\mathbf{1.5\cdot 10^{-3}}$ | $3\cdot 10^{-3}$ | $\mathbf{1.5\cdot 10^{-3}\;*}$
8 | $2.5\cdot 10^{-3}$ | $3\cdot 10^{-2}$ | $\mathbf{1.5\cdot 10^{-3}\;*}$
16 |  $3.5\cdot 10^{-3}$ |  $2\cdot 10^{-2}$ |  $\mathbf{2.5\cdot 10^{-3}\;*}$
Camel | 2 | $\mathbf{2\cdot 10^{-3}}$ | $\mathbf{2\cdot 10^{-3}}$ | $\mathbf{2\cdot 10^{-3}\;*}$
4 | $8\cdot 10^{-3}$ | $1\cdot 10^{-2}$ | $\mathbf{4\cdot 10^{-3}}$
8 | $4\cdot 10^{-2}$ | $1.6\cdot 10^{-2}$ | $\mathbf{5\cdot 10^{-3}}$
16 | $\dagger$ | $1.5\cdot 10^{-1}$ | $\mathbf{5\cdot 10^{-3}}$
Entangled circles | 2 | $1\cdot 10^{-2}$ | $\mathbf{4\cdot 10^{-3}}$ | $5\cdot 10^{-3}\;*$
Annulus w. cuts | 2 | $\mathbf{3\cdot 10^{-3}}$ | $4\cdot 10^{-3}\;*$ | $5\cdot 10^{-3}$
Scalar-top-loop | 3 | $7\cdot 10^{-4}$ | $\mathbf{5\cdot 10^{-4}}$ | $\mathbf{5\cdot 10^{-4}\;*}$
Polynomial | 18 | $1.5\cdot 10^{-3}$ | $1.5\cdot 10^{-3}\;*$ | $\mathbf{8\cdot 10^{-5}\;*}$
54 | $3\cdot 10^{-3}$ | $9\cdot 10^{-4}\;*$ | $\mathbf{8\cdot 10^{-5}\;*}$
96 | $\dagger$ | $8\cdot 10^{-4}\;*$ | $\mathbf{1\cdot 10^{-4}\;*}$
Table 4: Relative uncertainty on the integral estimate of the last iteration
of the runs of Table 2, based on a sample of 5000 points. The integrator that
adapted best to the integrand is boldfaced. A $*$ indicates when the value was
still decreasing and had not yet converged, a $\dagger$ is in place where the
algorithm did not converge to the true integrand.

In this section we show the performance of i-flow and compare it to VEGAS and
Foam based on the test functions we introduced in Sec. IV. For the VEGAS
algorithm, we use the default parameters as implemented in VEGAS . This
includes the use of stratified sampling and a maximum of 1000 bins per axis.
We further set the number of points per iteration to 5000. However, the
implementation in VEGAS  uses this number as a maximum, so we monitor the
actual number of function calls separately. The setup of Foam requires a
number of points per cell, which we fix to 5000. In the setup of i-flow, we
use $2\lceil\log_{2}D\rceil$ number of coupling layers with the masking
discussed in Sec. III.1, and the coupling transform $C$ taken to be a
Piecewise Rational Quadratic spline (Appendix A.3). The neural network in each
CL is taken to be a DNN of 5 layers with 32 nodes in each of the first four
layers. The number of nodes in the last layer depends on the coupling
transform $C$ and the dimensionality of the integrand. For the case of
Piecewise Rational Quadratic splines, the number of nodes is given by
$d\cdot(3n_{\text{bins}}+1)$, where $d$ is the number of dimensions to be
transformed. We further set the number of bins ($n_{\text{bins}}$) in each
output dimension to 16. The learning rate was set to $1\cdot 10^{-3}$ in all
cases. We use the exponential divergence, see eq. (44), as loss function.

To compare the integrators, we set a relative uncertainty on the integral
estimate as target. We then optimize the algorithms until the standard
deviation of the inverse-variance weighted combination of the estimates of
each optimization iteration (epoch) reaches this target. The inverse-variance
weighted combination is defined as:

$\mu=\frac{\sum_{i}\mu_{i}/\sigma_{i}^{2}}{\sum_{i}1/\sigma_{i}^{2}},\quad\sigma^{2}=\frac{1}{\sum_{i}1/\sigma_{i}^{2}},$
(24)

where $\mu_{i}(\sigma_{i})$ is the mean (standard deviation) of the
$i^{\text{th}}$ epoch and $\mu(\sigma)$ is the combination. The relative
uncertainty is defined as the uncertainty of the given estimate, normalized to
the true value of the integral 444Note that Foam directly gives the
uncertainty including all sampled points up to the given iteration and no
combination is needed.. Given this setup, there are three metrics that we use
to compare the integrators: 1) the number of function calls needed to reach
the target uncertainty; 2) how close the estimated integral value is to the
true value; 3) the uncertainty of the estimates in the last iterations. Each
of those highlights a different aspect of the integrator and we detail them
below. The results are shown in Tables 2–4. We chose a relative target
uncertainty of $10^{-4}$ for the non-polynomial test functions and $10^{-5}$
for the polynomials. For Gaussian and Camel functions, we use $\alpha=0.2$. In
addition, we set a cut-off at $5\cdot 10^{7}$ function calls.

Number of function calls. This number shows how often the integrand was
evaluated by the algorithm until the target uncertainty was reached. Having
fewer function calls is especially important when the function is numerically
expensive to evaluate and the computational overhead of the integration
algorithm becomes subleading. The results are shown in Tab. 2. We highlight
the entry with the fewest calls in boldface. In addition, we mark entries in
which the final integral estimate differs by more than 5 standard deviations
from the true result with a $\dagger$ and entries in which too much memory was
required by a $*$.

Integral estimate and uncertainty. This obviously shows how well the
integrator estimated the value of the integral. We show our results in Tab. 3
and compare them to the true, known results. We highlight in boldface the
entry with the smallest relative deviation (“pull”), defined as

$\frac{I_{\text{code}}-I_{\text{true}}}{\sqrt{(\Delta
I_{\text{code}}^{2}+\Delta I_{\text{true}}^{2})}}\;.$ (25)

Here, $I_{\text{code}}$ is the result from VEGAS, Foam, or i-flow,
$I_{\text{true}}$ is the true value of the integral, and the $\Delta I$ terms
signify the uncertainty in the integral. Note, that $\Delta I_{\text{true}}$
is only non-zero for the case of the entangle circles for which it is $5\cdot
10^{-9}$. Cases in which the cut-off for function calls was reached (see Tab.
2) are marked with a $\dagger$, cases that ran into memory problems are marked
with a $*$.

Relative uncertainty on the integral estimate in the last iterations. The
uncertainty on the integral estimate after adaptation is a measure for how
well the algorithm adapted to the integrand. Once the algorithm is fully
adapted, the uncertainty of a single integral estimate will be constant and
the combination of all iterations will follow the $1/\sqrt{N}$ scaling law for
MC estimates based on $N$ points. A better adapted algorithm introduces a
smaller coefficient for that scaling and therefore require fewer function
calls to reach a smaller uncertainty. We show our results in Tab. 4. Cases in
which VEGAS failed to converge to the right integral value are marked with a
$\dagger$, a $*$ shows entries that still showed a downward trend at the end
of the optimization, indicating that the algorithm was still adapting to the
integrand. We highlight the integrator with the smallest uncertainty in
boldface.

For the Gaussians, VEGAS always has the fewest calls. This is expected, since
the integrand factorizes. However, the number grows rapidly for increasing
integrand dimension, whereas the number for i-flow grows slower. Foam is not
able to reach the target uncertainty for $D=8,16$ before the cut-off of
$5\cdot 10^{7}$ function calls. i-flow has adapted best to all Gaussians of
$D>2$, as can be seen in Tab. 4. This means that if a sufficiently small
target uncertainty is required, i-flow would potentially need fewer function
calls to reach it. The fact that the optimization of i-flow was not complete
when the target uncertainty was reached can also be seen in Fig. 3(a), where
the accumulated uncertainty of i-flow (red line) was falling quicker than
$1/\sqrt{N}$ (dashed gray lines). In almost all of the cases, the integral
estimate of i-flow was closest to the true integral value.

For the Camel functions, VEGAS only has the fewest calls for $D=2$, in higher
dimensions i-flow needs fewer calls. Note that in 16 dimensions, VEGAS
completely misses one of the two peaks, yielding an estimate that is off by a
factor of two. Since in this case the integrand is like a Gaussian, VEGAS
converges quicker than the other algorithms. The integral estimate of i-flow
seems also off fn:camel16 , but this is due to the fact that it needs roughly
200 epochs to “see” the structure of the integrand and all of those iterations
contribute to the final number in Tab. 2. Again, Foam needs too many points
for $D>4$ to reach the target uncertainty, the integral estimates of i-flow
are closest to the true value, and i-flow has adapted best to the integrand.
The latter can be seen by the small relative uncertainties in Tab. 4 and the
scaling in the 4 dimensional case shown in Fig. 3(b).

We discuss the entangled circles and the annulus after the polynomials. For
the scalar top loop, VEGAS needs the fewest function calls, but all 3
integrators estimate the true value within one standard deviation. It is,
however, interesting to see that both VEGAS and Foam seem to be fully adapted,
whereas the uncertainties of the estimates of i-flow were improving much
faster than the $1/\sqrt{N}$ expectation, see Fig. 3(d).

The polynomials show the strength of i-flow. It has no problems adapting to
the high-dimensional integrand, as can be seen in Tab. 4. Therefore, i-flow
needs comparatively few function calls to reach the target uncertainty. Since
the polynomial does not factorize, VEGAS does not adapt well, or in the case
of $D=96$ not at all. The difference between the adaptation of the algorithms
is also visible in Fig. 3(c). There, however, we see an interesting pattern in
the accumulated uncertainty of Foam that we want to comment on. First, since
Foam estimates the integral and uncertainty of a given iteration always on the
points of all previous iterations, the uncertainty can grow for a growing
number of points if the central value shifts. Second, due to the symmetry of
the polynomial integrand, we see a periodic pattern that we can understand as
follows. We start with an uncertainty based on the first 5000 points. Adding
more points at this initial stage lets the algorithm “see” more structure of
the integrand and the uncertainty grows. A large cell with a large spread of
functional values within it is then further split consecutively into many
smaller cells. That reduces the spread of functional values per cell and
therefore the uncertainty of the integral estimate. Once the uncertainty drops
below a certain value, Foam stops splitting these (smaller) cells and returns
to one of the “bigger” cells it did not split in the beginning and starts
splitting it. This initially increases the uncertainty again, because the
spread of functional values in the large cell is larger than it was in the
smaller cells. The result is the oscillating pattern we see in Fig. 3(c). Note
that the minima of this pattern follow the $1/\sqrt{N}$ scaling.

The entangled circles are best integrated by Foam, as it is only 2
dimensional, yet non-factorizable. i-flow is slightly worse, but not by much.
VEGAS, however, does not perform well. Similar statements can be made about
the annulus function with hard cuts. VEGAS does the worst because of the non-
factorizing structure of the integrand and Foam does well because it is only a
2-dimensional problem. As discussed in the earlier sections, i-flow also
allows efficient sampling once it “learns” the integral up to small
uncertainty, we therefore use these test functions to illustrate the sampling
performance of i-flow. As an example, Fig. 4 shows a sample distribution after
training i-flow with $5\cdot 10^{6}$ points (1000 epochs with 5000 points per
epoch) on the annulus function of Eq. (20). For training, we used a learning
schedule with exponential decay. An initial learning rate of $2\cdot 10^{-3}$
is halved every 250 epochs. The cut efficiency, defined as the fraction of the
generated points that pass the hard cut, is $89.6\%$. Figure 5 shows the
weights of 10000 points sampled after training with 10M points on the
Entangled Circles of Eq. (19). In the ideal case of $g\to f/I$, we expect the
weight distribution to approach a delta function. In Fig. 5, we see that the
trained results are much more like a delta function than the flat prior,
showing significant improvement in the ability to draw samples from this
function.

It is clear from these considerations that for the low-dimensional integrals
($D\leq 4$), all three integrators achieve reasonable results. If the target
uncertainty is not very small, VEGAS or Foam provide the best integrator,
depending on the integrand at hand. If, however, a very small target
uncertainty is needed, i-flow is the better option as it adapts really well to
the shape of the integrand. It is only the fact that i-flow adapts slower than
VEGAS that makes i-flow lose in the beginning, as illustrated in Figure 3. For
higher-dimensional integrands ($D\geq 4$) i-flow requires fewer function calls
because it adapts better to the integrand. For example, VEGAS fails in the
integration of 16-dimensional Camel function completely (missing one of the
peaks) and Foam has a large uncertainty on the final result, even though it
has much more function calls. Foam also performs poorly in the case of the
Gaussian in 16 dimensions. In both of these cases, Foam approximately requires
$b^{D}$ number of cells to map out all the features of a function, where $b$
is the average number of bins in each dimension. If $b$ is taken to be 2, for
16 dimensions, the number of cells required is at least $2^{16}$, which is far
greater than our set cut-off of 10000 cells. Therefore, when dealing with
high-dimensional integrals, Foam is the least efficient integrator.

To quantify the computational overhead of i-flow in comparison to VEGAS, we
trained both for 100 iterations with 5000 points per iteration on the
polynomial function. It took VEGAS consistently 2 seconds for 2, 4, 8, 16, and
32 dimensions, and it took i-flow 14.7, 37.2, 80.1, 176.4, and 359.2 seconds
respectively on a Laptop with Intel(R) Core(TM) i7-7700HQ CPU @ 2.80GHz. This
increase is due to needing more Coupling Layers and therefore increasing the
number of trainable parameters in higher dimensions. Working out the time for
the 32 dimensions, we find that if the function evaluation takes much longer
than 720 $\mu\text{s}$ then the overhead starts to become unimportant.
Additionally, if the difference in function evaluations to reach a target
precision are taken into account, the time for function evaluation is even
smaller in order for the additional overhead of i-flow to become
insignificant.

To summarize, i-flow provides the best integrator for integrals in 4 or more
dimensions, especially if a high precision is needed and/or the integrand is
numerically expensive and slow to evaluate.

(a) 4-dimensional Gaussian

(b) 4-dimensional Camel

(c) 54-dimensional Polynomial

(d) Scalar top loop

Figure 3: Accumulated (total) relative uncertainty of the integral, defined as
the inverse-variance weighted combination of the uncertainties (normalized to
the integral value) per iteration, as a function of number of points used in
training for four different test functions. The dashed lines indicate a
$1/\sqrt{N}$ scaling.

Figure 4: A set of 7500 points sampled after training i-flow with 5M points on
the Ring function. 6720 are inside (blue), 780 outside (red).

Figure 5: Weights of 10000 points, sampled after training i-flow on Entangled
circles (19). $g$ is a flat distribution before training and approximately
resembles the shape of $f$ after training.

##  VI Conclusion and Outlook

As shown in the previous section, i-flow tends to do better than both VEGAS
and Foam for all the test cases provided. However, i-flow comes with a few
downsides. Since i-flow has to learn all the correlations of the function, it
takes significantly longer to achieve optimal performance compared to the
other integrators. This can be seen in Fig. 3. This obviously translates to
longer training times. Additionally, the memory footprint required for i-flow
is much larger due to requiring storage for quicker parameter updates within
the NNs. Both of these can be overcome with future improvements.

There are several directions in which we plan to improve the presented setup
in the future. So far, we only used simple NN architectures in our coupling
layers. Using convolutional NNs instead might improve the convergence of the
normalizing flow for complicated integrands, as these networks have the
ability to learn complicated shapes in images with fewer parameters than dense
networks.

The setup suggested in 2019arXiv190510347T  would allow the extension of
i-flow to discrete distributions, which also has applications in HEP
Gao:2020zvv ; Bothmann:2020ywa . Another way to implement this type of
information is utilizing Conditional Normalizing Flows winkler2019learning .

The implementation of transflow-learning, which was suggested in
gambardella2019transflow , would allow the use of a trained normalizing flow
on different, but similar problems without retraining the network. Such
problems arise in HEP when new-physics effects from high energy scales modify
scattering properties at low energies slightly and are described in an
effective field theory framework. Another application for transflow-learning
would be to train one network for a given dimensionality and adapt the network
for another problem with the same dimensionality.

Using techniques like gradient checkpointing 2016arXiv160406174C  have the
potential to reduce the memory usage substantially, therefore allowing more
points to be used at each training step or larger NN architectures.

The setup presented in 2018arXiv180703039K , which introduces invertible
$1\times 1$ convolutions, showed an improved performance over the vanilla
implementation of the normalizing flows, which possibly also applies to our
case. These $1\times 1$ convolutions are generalizations of permutation
operators acting on the inputs. Additionally, this would modify the maximum
number of coupling layers required by having more expressive permutations.

## Acknowledgements

We thank Joao M. Goncalves Caldeira, Felix Kling, Luisa Lucie-Smith, Tilman
Plehn, Holger Schulz, Nhan Tran, Paddy Fox, William Jay, David Shih, and the
participants of the Aspen workshop “The Energy Frontier Beyond the LHC Run 2”
for their comments and discussions. We further thank Stefan Höche for helpful
discussions, comments, and for his Foam implementation Hoeche:foam .
This manuscript has been authored by Fermi Research Alliance, LLC under
Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of
Science, Office of High Energy Physics. This work was performed in part at
Aspen Center for Physics, which is supported by National Science Foundation
grant PHY-1607611. C.K. acknowledges the support of the Alexander von Humboldt
Foundation.

##  Appendix A Coupling Layer details

The implementation of the layers available in i-flow are detailed below. The
layers are based on the work of DBLP:journals/corr/abs-1808-03856 ;
2019arXiv190604032D  and are reproduced here for the convenience of the
reader.

###  A.1 Piecewise Linear

For the piecewise linear coupling layer DBLP:journals/corr/abs-1808-03856 ,
given $K$ bins of width $w$, the probability density function (PDF) is defined
as:

$q_{i}(t)=\begin{cases}Q_{i1}/w&t<w\\\ Q_{i2}/w&w\leq t<2w\\\ \quad\vdots&\\\
Q_{iK}/w&1-w\leq t<1\\\ \end{cases}\;.$ (26)

The cumulative distribution function (CDF) is defined by the integral giving:

$C(x_{i}^{B};Q)=\alpha Q_{ib}+\sum_{k=1}^{b-1}Q_{ik},$ (27)

where $b$ is the bin in which $x_{i}^{B}$ occurs ($(b-1)w\leq x_{i}^{B}<bw$),
and $\alpha=\frac{x_{i}^{B}-(b-1)w}{w}$. Alternatively, we can define $b$ as
the maximal $b$ for which $\left(C_{i}-\sum_{k=1}^{b-1}Q_{ik}\right)>0$. The
inverse CDF is given by:

$x_{i}^{B}(C_{i};Q)=\frac{w\left(C_{i}-\sum_{k=1}^{b-1}Q_{ik}\right)}{Q_{ib}}+(b-1)w.$
(28)

The Jacobian for this network is straightforward to calculate, and gives:

$\left|\frac{\partial C}{\partial x_{B}}\right|=\prod_{i}Q_{ib}/w\;.$ (29)

The piecewise linear layers require fixed bin widths in each layer. For
details on why this is required, see Appendix B of
DBLP:journals/corr/abs-1808-03856 .

###  A.2 Piecewise Quadratic

For the piecewise quadratic coupling layer DBLP:journals/corr/abs-1808-03856 ,
given $K$ bins with widths $W_{ik}$, with $K+1$ vertex heights given by
$V_{ik}$, the PDF is defined as:

$q_{i}(t)=\begin{cases}\frac{V_{i2}-V_{i1}}{W_{i1}}t+V_{i1}&t<W_{i1}\\\
\frac{V_{i3}-V_{i2}}{W_{i2}}\left(t-W_{i1}\right)+V_{i2}&W_{i1}\leq
t<W_{i1}+W_{i2}\\\ \quad\vdots&\\\
\frac{V_{i(K+1)}-V_{iK}}{W_{iK}}\left(t-\sum_{k=1}^{K-1}W_{ik}\right)+V_{iK}&\sum_{k=1}^{K-1}W_{ik}\leq
t<1\\\ \end{cases}$ (30)

Integrating the above equation leads to the CDF:

$C(x_{i}^{B};W,V)=\frac{\alpha^{2}}{2}\left(V_{ib+1}-V_{ib}\right)W_{ib}+V_{ib}W_{ib}\alpha+\sum_{k=1}^{b-1}\frac{V_{ik+1}+V_{ik}}{2}W_{ik},$
(31)

where $b$ is defined as the solution to 555Note that this definition means
$b\in[1,K]$. $\sum_{k=1}^{b-1}W_{ik}\leq x_{i}^{B}<\sum_{k=1}^{b}W_{ik}$, and
$\alpha=\frac{x_{i}^{B}-\sum_{k=1}^{b-1}W_{ik}}{W_{ib}}$ is the relative
position of $x_{i}^{B}$ in bin $b$. Inverting the CDF leads to:

$x_{i}^{B}(C_{i};W,V)=W_{ib}\left(\frac{-V_{ib}}{V_{ib+1}-V_{ib}}+\sqrt{\frac{V_{ib}^{2}}{\left(V_{ib+1}-V_{ib}\right)^{2}}+2\beta}\right)+\sum_{k=1}^{b-1}W_{ik},$
(32)

where $b$ is defined as the solution to

$\sum_{k=1}^{b-1}\frac{V_{ik}+V_{ik+1}}{2}W_{ik}\leq
C_{i}<\sum_{k=1}^{b}\frac{V_{ik}+V_{ik+1}}{2}W_{ik},$ (33)

and $\beta$ is the relative position of $C_{i}$ in the bin $b$, and is given
by:

$\beta=\frac{C_{i}-\sum_{k=1}^{b-1}\frac{V_{ik}+V_{ik+1}}{2}W_{ik}}{\left(V_{ib+1}-V_{ib}\right)W_{ib}}.$
(34)

###  A.3 Piecewise Rational Quadratic

For the piecewise rational quadratic coupling layer 2019arXiv190604032D ,
given $K+1$ knot points
$\left\\{\left(x^{(k)},y^{(k)}\right)\right\\}_{k=0}^{K}$ that are
monotonically increasing, with $(x^{(0)},y^{(0)})=(0,0)$ and
$(x^{(K)},y^{(K)})=(1,1)$, and $K+1$ non-negative derivatives
$\left\\{d^{(k)}\right\\}_{k=0}^{K}$, the CDF can be calculated using the
algorithm from 10.1093/imanum/2.2.123 , which is roughly reproduced below.

First, we define the bin widths ($w^{(k)}=x^{(k+1)}-x^{(k)}$) and the slopes
($s^{(k)}=\frac{y^{(k+1)}-y^{(k)}}{w^{(k)}}$). We next obtain the fractional
distance ($\xi$) between the two knots that the point of interest ($x$) lies
($\xi=\frac{x-x^{(k)}}{w^{(k)}}$, where $k$ is the bin $x$ lies in). The CDF
is given by:

$g(x)=\frac{\alpha(\xi)}{\beta(\xi)},$ (35)

where the details of $\alpha(\xi)$ and $\beta(\xi)$ can be found in
10.1093/imanum/2.2.123 , but simplifies to:

$g(x)=y^{(k)}+\frac{\left(y^{(k+1)}-y^{(k)}\right)\left[s^{(k)}\xi^{2}+d^{(k)}\xi\left(1-\xi\right)\right]}{s^{(k)}+\left[d^{(k+1)}+d^{(k)}-2s^{(k)}\right]\xi\left(1-\xi\right)},$
(36)

which is noted to be less prone to numerical issues 10.1093/imanum/2.2.123 .
The inverse can be found by solving a quadratic equation 2019arXiv190604032D :

$q(x)=\alpha(\xi)-y\beta(\xi)=ax^{2}+bx+c=0,$ (37)

where the coefficients are given in 2019arXiv190604032D , solving this
equation for the solution that gives a monotonically increasing $x$ results
in:

$x=\frac{-b+\sqrt{b^{2}-4ac}}{2a}=\frac{2c}{-b-\sqrt{b^{2}-4ac}},$ (38)

where the second form is numerically more precise when $4ac$ is small, and is
also valid for $a=0$ 2019arXiv190604032D .

##  Appendix B Loss functions

We implemented several different divergences that can be used as loss
functions. They differ in $p\leftrightarrow q$ symmetry, relative weight
between small and large deviations, treatment of $p=0$ case (also in the
derivative), and numerical complexity. All of them are from the class of
$f$-divergences 2014ISPL…21…10N .
Pearson $\chi^{2}$ divergence:

$D_{\chi^{2}}=\int\frac{(p(x)-q(x))^{2}}{q(x)}\;dx$ (39)

Kullback–Leibler divergence:

$D_{KL}=\int p(x)\log{\left(\frac{p(x)}{q(x)}\right)}\;dx$ (40)

squared Hellinger distance:

$D_{H^{2}}=\int 2\left(\sqrt{p(x)}-\sqrt{q(x)}\right)^{2}\;dx$ (41)

Jeffreys divergence:

$D_{J}=\int(p(x)-q(x))(\log{p(x)}-\log{q(x)})\;dx$ (42)

Chernoff’s $\alpha$-divergence:

$D_{C\alpha}=\frac{4}{1-\alpha^{2}}\left(1-\int
p(x)^{\frac{1-\alpha}{2}}q(x)^{\frac{1+\alpha}{2}}\;dx\right)$ (43)

exponential divergence:

$D_{e}=\int p(x)\log{\left(\frac{p(x)}{q(x)}\right)}^{2}\;dx$ (44)

$(\alpha,\beta)$-product divergence:

$D_{\alpha\beta}=\frac{2}{(1-\alpha)(1-\beta)}\int\left(1-\left(\frac{q(x)}{p(x)}\right)^{\frac{1-\alpha}{2}}\right)\left(1-\left(\frac{q(x)}{p(x)}\right)^{\frac{1-\beta}{2}}\right)p(x)\;dx$
(45)

Jensen–Shannon divergence:

$D_{JS}=\frac{1}{2}\int
p(x)\log{\left(\frac{2p(x)}{p(x)+q(x)}\right)}+q(x)\log{\left(\frac{2q(x)}{p(x)+q(x)}\right)}\;dx$
(46)

##  Appendix C Sector decomposition of scalar loop integrals

Following Binoth:2002xh , we give the integral representations of triangle and
box functions in 4 dimensions using the Feynman parametrisation. To begin
with, the triangle integral with external particles of energy
$\sqrt{s_{1}},\sqrt{s_{2}},\sqrt{s_{3}}$ and internal propagators of masses
$m_{1},m_{2},m_{3}$ is given by

$\begin{split}I_{3}(s_{1},s_{2},s_{3},m_{1}^{2},m_{2}^{2},m_{3}^{2})=&\int\frac{d^{4}k}{i\pi^{2}}\frac{1}{[(k-r_{1})^{2}-m_{1}^{2}][(k-r_{2})^{2}-m_{2}^{2}][k^{2}-m_{3}^{2}]}\\\
=&-\int_{0}^{\infty}d^{3}x\,\delta(1-x_{123})\frac{x_{123}^{-1}}{\mathcal{F}_{Tri}}\\\
\mathcal{F}_{Tri}=(-s_{1})x_{1}x_{3}+(-s_{2})&x_{1}x_{2}+(-s_{3})x_{2}x_{3}+x_{123}(x_{1}m_{1}^{2}+x_{2}m_{2}^{2}+x_{3}m_{3}^{2})-i\epsilon\\\
x_{123}=x_{1}+x_{2}+x_{3}\quad\quad\quad&\end{split}$ (47)

The 3-dimensional integral is further split into 3 sectors by the
decomposition:

$1=\Theta(x_{1}>x_{2},x_{3})+\Theta(x_{2}>x_{1},x_{3})+\Theta(x_{3}>x_{1},x_{2})$
(48)

For example, when $x_{3}>x_{1},x_{2}$, after the variable transformation
$t_{i}=x_{i}/x_{3}(i=1,2)$, the integral simplifies to

$\begin{split}S_{Tri}&(s_{1},s_{2},s_{3},m_{1}^{2},m_{2}^{2},m_{3}^{2})=\int_{0}^{1}dt_{1}dt_{2}\frac{(1+t_{1}+t_{2})^{-1}}{\tilde{\mathcal{F}}_{Tri}(s_{1},s_{2},s_{3},m_{1}^{2},m_{2}^{2},m_{3}^{2},t_{1},t_{2})}\\\
&\tilde{\mathcal{F}}_{Tri}(s_{1},s_{2},s_{3},m_{1}^{2},m_{2}^{2},m_{3}^{2},t_{1},t_{2})=(-s_{1})t_{1}+(-s_{2})t_{1}t_{2}+(-s_{3})t_{2}\\\
&\qquad\qquad\qquad\qquad\qquad\qquad\qquad+(1+t_{1}+t_{2})(t_{1}m_{1}^{2}+t_{2}m_{2}^{2}+m_{3}^{2})\end{split}$
(49)

Therefore,

$\begin{split}I_{3}=S_{Tri}(s_{1},s_{2},s_{3},m_{1}^{2},m_{2}^{2},m_{3}^{2})+S_{Tri}(1\,2\,3\,\to
2\,3\,1\,)+S_{Tri}(1\,2\,3\,\to 3\,1\,2\,)\end{split}$ (50)

One can perform the same trick to treat the box integral with 4 external
fields and 4 propagators. After sector decomposition, one gets

$\begin{split}I_{4}=&S_{Box}(s_{12},s_{23},s_{1},s_{2},s_{3},s_{4},m_{1}^{2},m_{2}^{2},m_{3}^{2},m_{4}^{2})+\textbf{permutations}\\\
S_{Box}&(s_{12},s_{23},s_{1},s_{2},s_{3},s_{4},m_{1}^{2},m_{2}^{2},m_{3}^{2},m_{4}^{2})=\int_{0}^{1}dt_{1}dt_{2}dt_{3}\frac{1}{\tilde{\mathcal{F}}_{Box}^{2}}\\\
\tilde{\mathcal{F}}&{}_{Box}=(-s_{12})t_{2}+(-s_{23})t_{1}t_{3}+(-s_{1})t_{1}+(-s_{2})t_{1}t_{2}+(-s_{3})t_{2}t_{3}\\\
&\quad\quad+(-s_{4})t_{3}+(1+t_{1}+t_{2}+t_{3})(t_{1}m_{1}^{2}+t_{2}m_{2}^{2}+t_{3}m_{3}^{2}+m_{4}^{2})\end{split}$
(51)

## References

  * (1) Asger Hobolth, Marcy K Uyenoyama, and Carsten Wiuf.  Importance sampling for the infinite sites model.  Statistical Applications in Genetics and Molecular Biology, 7(1), January 2008.  URL: https://doi.org/10.2202/1544-6115.1400, doi:10.2202/1544-6115.1400.
  * (2) Charles J. Mode et al.  Applications of Monte Carlo Methods in Biology, Medicine and Other Fields of Science.  IntechOpen, 2011.  doi:10.5772/634.
  * (3) Marshall N. Rosenbluth and Arianna W. Rosenbluth.  Monte Carlo Calculation of the Average Extension of Molecular Chains.  The Journal of Chemical Physics, 23(2):356–359, Feb 1955.  doi:10.1063/1.1741967.
  * (4) H. T. MacGillivray and R. J. Dodd.  Monte-carlo simulations of galaxy systems.  Astrophysics and Space Science, 86(2):419–435, Sep 1982.  URL: https://doi.org/10.1007/BF00683346, doi:10.1007/BF00683346.
  * (5) D. W. O. Rogers.  REVIEW: Fifty years of Monte Carlo simulations for medical physics.  Physics in Medicine and Biology, 51(13):R287–R301, Jul 2006.  doi:10.1088/0031-9155/51/13/R17.
  * (6) Peter Jäckel.  Monte Carlo methods in finance, volume 71.  Wiley, 2002.
  * (7) László Szirmay-Kalos.  Monte Carlo methods in global illumination : photo-realistic rendering with randomization.  VDM Verl. Dr. Müller, Saarbrücken, 2008.
  * (8) B.R. Webber.  Monte Carlo Simulation of Hard Hadronic Processes.  Ann. Rev. Nucl. Part. Sci., 36:253–286, 1986.
  * (9) Andy Buckley et al.  General-purpose event generators for LHC physics.  Phys. Rept., 504:145–233, 2011.  arXiv:1101.2599, doi:http://dx.doi.org/10.1016/j.physrep.2011.03.005.
  * (10) Andy Buckley.  Computational challenges for MC event generation.  In 19th International Workshop on Advanced Computing and Analysis Techniques in Physics Research: Empowering the revolution: Bringing Machine Learning to High Performance Computing (ACAT 2019) Saas-Fee, Switzerland, March 11-15, 2019, 2019.  arXiv:1908.00167.
  * (11) The ATLAS collaboration.  Computing and Software - Public Results.  https://twiki.cern.ch/twiki/bin/view/AtlasPublic/ComputingandSoftwarePublicResults.
  * (12) F. James.  Monte-Carlo phase space.  1968\.  CERN-68-15.
  * (13) G. Peter Lepage.  A New Algorithm for Adaptive Multidimensional Integration.  J. Comput. Phys., 27:192, 1978.  doi:10.1016/0021-9991(78)90004-9.
  * (14) G. Peter Lepage.  VEGAS: AN ADAPTIVE MULTIDIMENSIONAL INTEGRATION PROGRAM.  1980\.
  * (15) S. Jadach.  Foam: A General purpose cellular Monte Carlo event generator.  Comput. Phys. Commun., 152:55–100, 2003.  arXiv:physics/0203033, doi:10.1016/S0010-4655(02)00755-5.
  * (16) Dimitri Bourilkov.  Machine and Deep Learning Applications in Particle Physics.  2019\.  arXiv:1912.08245.
  * (17) Joshua Bendavid.  Efficient Monte Carlo Integration Using Boosted Decision Trees and Generative Deep Neural Networks.  2017\.  arXiv:1707.00028.
  * (18) Matthew D. Klimek and Maxim Perelstein.  Neural Network-Based Approach to Phase Space Integration.  2018\.  arXiv:1810.11509.
  * (19) An $N$ particle final state phase space is a $D\approx 4N-3$ dimensional integral, when including recursive multichannel selection in the integral.
  * (20) Sydney Otten, Sascha Caron, Wieske de Swart, Melissa van Beekveld, Luc Hendriks, Caspar van Leeuwen, Damian Podareanu, Roberto Ruiz de Austri, and Rob Verheyen.  Event Generation and Statistical Sampling for Physics with Deep Generative Models and a Density Information Buffer.  2019\.  arXiv:1901.00875.
  * (21) Bobak Hashemi, Nick Amin, Kaustuv Datta, Dominick Olivito, and Maurizio Pierini.  LHC analysis-specific datasets with Generative Adversarial Networks.  2019\.  arXiv:1901.05282.
  * (22) Riccardo Di Sipio, Michele Faucci Giannelli, Sana Ketabchi Haghighat, and Serena Palazzo.  DijetGAN: A Generative-Adversarial Network Approach for the Simulation of QCD Dijet Events at the LHC.  JHEP, 08:110, 2020.  arXiv:1903.02433, doi:10.1007/JHEP08(2019)110.
  * (23) Anja Butter, Tilman Plehn, and Ramon Winterhalder.  How to GAN LHC Events.  SciPost Phys., 7:075, 2019.  arXiv:1907.03764, doi:10.21468/SciPostPhys.7.6.075.
  * (24) Stefano Carrazza and Frédéric A. Dreyer.  Lund jet images from generative and cycle-consistent adversarial networks.  2019\.  arXiv:1909.01359.
  * (25) C. Ahdida et al.  Fast simulation of muons produced at the SHiP experiment using Generative Adversarial Networks.  2019\.  arXiv:1909.04451.
  * (26) Anja Butter, Tilman Plehn, and Ramon Winterhalder.  How to GAN Event Subtraction.  2019\.  arXiv:1912.08824.
  * (27) Konstantin T. Matchev and Prasanth Shyamsundar.  Uncertainties associated with GAN-generated datasets in high energy physics.  2020\.  arXiv:2002.06307.
  * (28) Fady Bishara and Marc Montull.  (Machine) Learning amplitudes for faster event generation.  2019\.  arXiv:1912.11055.
  * (29) Laurent Dinh, David Krueger, and Yoshua Bengio.  NICE: non-linear independent components estimation.  In 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Workshop Track Proceedings, 2015\.  URL: http://arxiv.org/abs/1410.8516.
  * (30) Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio.  Density estimation using real nvp.  ArXiv, abs/1605.08803, 2016.
  * (31) Danilo Jimenez Rezende and Shakir Mohamed.  Variational inference with normalizing flows, 2015.  arXiv:1505.05770.
  * (32) Thomas Müller, Brian Mcwilliams, Fabrice Rousselle, Markus Gross, and Jan Novák.  Neural importance sampling.  ACM Trans. Graph., 38(5), October 2019.  URL: https://doi.org/10.1145/3341156, doi:10.1145/3341156.
  * (33) Conor Durkan, Artur Bekasov, Iain Murray, and George Papamakarios.  Neural Spline Flows.  arXiv e-prints, page arXiv:1906.04032, Jun 2019.  arXiv:1906.04032.
  * (34) Jiaming Song, Shengjia Zhao, and Stefano Ermon.  A-NICE-MC: Adversarial Training for MCMC.  arXiv e-prints, Jun 2017.  arXiv:1706.07561.
  * (35) Daniel Levy, Matthew D. Hoffman, and Jascha Sohl-Dickstein.  Generalizing Hamiltonian Monte Carlo with Neural Networks.  arXiv e-prints, Nov 2017.  arXiv:1711.09268.
  * (36) Matthew Hoffman, Pavel Sountsov, Joshua V. Dillon, Ian Langmore, Dustin Tran, and Srinivas Vasudevan.  NeuTra-lizing Bad Geometry in Hamiltonian Monte Carlo Using Neural Transport.  arXiv e-prints, Mar 2019.  arXiv:1903.03704.
  * (37) Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dandelion Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015\.  Software available from tensorflow.org.  URL: https://www.tensorflow.org/.
  * (38) Frank Nielsen and Richard Nock.  On the Chi Square and Higher-Order Chi Distances for Approximating f -Divergences.  IEEE Signal Processing Letters, 21(1):10–13, Jan 2014.  arXiv:1309.3029, doi:10.1109/LSP.2013.2288355.
  * (39) Hidetosi Takahasi and Masatake Mori.  Double exponential formulas for numerical integration.  Publications of the Research Institute for Mathematical Sciences, 9(3):721–741, 1974.
  * (40) Stefan Höche, Stefan Prestel, and Holger Schulz.  Simulation of vector boson plus many jet final states at the high luminosity LHC.  2019\.  arXiv:1905.05120.
  * (41) G. Peter Lepage.  Vegas documentation release 3.4.  URL: https://vegas.readthedocs.io/en/latest/index.html.
  * (42) S. Hoeche.  private implementation of foam integrator.  URL: {https://gitlab.com/shoeche/foam}.
  * (43) Falko Dulat, Bernhard Mistlberger, and Andrea Pelloni.  Differential Higgs production at N3LO beyond threshold.  JHEP, 01:145, 2018.  arXiv:1710.03016, doi:10.1007/JHEP01(2018)145.
  * (44) George Papamakarios, Theo Pavlakou, and Iain Murray.  Masked Autoregressive Flow for Density Estimation.  arXiv e-prints, page arXiv:1705.07057, May 2017.  arXiv:1705.07057.
  * (45) George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshminarayanan.  Normalizing flows for probabilistic modeling and inference, 2019.  arXiv:1912.02762.
  * (46) Allan Pinkus.  Weierstrass and approximation theory.  Journal of Approximation Theory, 107(1):1 – 66, 2000.  URL: http://www.sciencedirect.com/science/article/pii/S0021904500935081, doi:https://doi.org/10.1006/jath.2000.3508.
  * (47) M. H. Stone.  Applications of the theory of boolean rings to general topology.  Transactions of the American Mathematical Society, 41(3):375–481, 1937.  URL: http://www.jstor.org/stable/1989788.
  * (48) Marshall H. Stone.  The generalized weierstrass approximation theorem.  Mathematics Magazine, 21(5):237–254, 1948.  URL: http://www.jstor.org/stable/3029337.
  * (49) Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization, 2014.  arXiv:1412.6980.
  * (50) Since we initialize the last layer of each network with vanishing bias and weights, in the first sampling $g(x)$ is constant.
  * (51) Wolfram Research, Inc.  Mathematica, Version 12.0.  Champaign, IL, 2019.  URL: https://www.wolfram.com/mathematica.
  * (52) There is no known analytic solution to this given function.
  * (53) T. Binoth, G. Heinrich, and N. Kauer.  A Numerical evaluation of the scalar hexagon integral in the physical region.  Nucl. Phys., B654:277–300, 2003.  arXiv:hep-ph/0210023, doi:10.1016/S0550-3213(03)00052-X.
  * (54) T. Hahn and M. Pérez-Victoria.  Automated one-loop calculations in four and d dimensions.  Computer Physics Communications, 118(2):153 – 165, 1999.  URL: http://www.sciencedirect.com/science/article/pii/S0010465598001738, doi:https://doi.org/10.1016/S0010-4655(98)00173-8.
  * (55) Christina Gao, Stefan Höche, Joshua Isaacson, Claudius Krause, and Holger Schulz.  Event Generation with Normalizing Flows.  2020\.  arXiv:2001.10028.
  * (56) Enrico Bothmann, Timo Janßen, Max Knobbe, Tobias Schmale, and Steffen Schumann.  Exploring phase space with Neural Importance Sampling.  2020\.  arXiv:2001.05478.
  * (57) The estimate of i-flow only deviates from the true value because the estimates from all iterations are combined and the first $\sim 200$ epochs only “see” one of the two peaks. Combining 15 of the last epochs yields $0.86377(136)$, which is closer to the true value.
  * (58) Note that Foam directly gives the uncertainty including all sampled points up to the given iteration and no combination is needed.
  * (59) Dustin Tran, Keyon Vafa, Kumar Krishna Agrawal, Laurent Dinh, and Ben Poole.  Discrete Flows: Invertible Generative Models of Discrete Data.  arXiv e-prints, page arXiv:1905.10347, May 2019.  arXiv:1905.10347.
  * (60) Christina Winkler, Daniel Worrall, Emiel Hoogeboom, and Max Welling.  Learning likelihoods with conditional normalizing flows.  2019\.  arXiv:1912.00042.
  * (61) Andrew Gambardella, Atılım Güneş Baydin, and Philip H. S. Torr.  Transflow learning: Repurposing flow models without retraining.  2019\.  arXiv:1911.13270.
  * (62) Tianqi Chen, Bing Xu, Chiyuan Zhang, and Carlos Guestrin.  Training Deep Nets with Sublinear Memory Cost.  arXiv e-prints, page arXiv:1604.06174, Apr 2016.  URL: https://github.com/cybertronai/gradient-checkpointing, arXiv:1604.06174.
  * (63) Diederik P. Kingma and Prafulla Dhariwal.  Glow: Generative Flow with Invertible 1x1 Convolutions.  arXiv e-prints, page arXiv:1807.03039, Jul 2018.  arXiv:1807.03039.
  * (64) Note that this definition means $b\in[1,K]$.
  * (65) J. A. GREGORY and R. DELBOURGO.  Piecewise Rational Quadratic Interpolation to Monotonic Data.  IMA Journal of Numerical Analysis, 2(2):123–130, 04 1982.  URL: https://doi.org/10.1093/imanum/2.2.123, arXiv:http://oup.prod.sis.lan/imajna/article-pdf/2/2/123/2267745/2-2-123.pdf, doi:10.1093/imanum/2.2.123.