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

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 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/ √ 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 (Ω). Given N uniform samples, the integral of f (x) can be approximated by, and the uncertainty is determined by the standard deviation of the mean, where V is the volume encompassed by Ω and 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 highdimensional 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 a uniform distribution, one samples from a distribution g(x) that ideally has the same shape as the integrand f (x). Using the transformation dx = dG(x)/g(x), with G(x) the cumulative distribution function of g(x), one obtains In the ideal case when g(x) → 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 [13,14] 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.
where f : R D → R and f i : R → 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. [35]). The resulting uncertainty can be reduced further by applying stratified sampling, in addition to the VEGAS algorithm, after the binning [36]. Foam [15] 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 [37], 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 [38].

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 FIG. 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 . 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 x: Based on the chain rule, the output x K 's probability distribution, g K , can be inferred given the base probability distribution g 0 from which x is drawn: 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 [28,29], then generalized by [30,31], 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 shows the basic structure of a coupling layer, which is a special design of the bijective mapping c. For each map, the input variable x = {x 1 , .., x D } is partitioned into two subsets, x A and x B which can be determined arbitrarily so long as neither is the empty set. This arbitrary partitioning will be referred to as masking. Without loss of generality, one simple partitioning is given by x A = {x 1 , .., x d } and 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 The NN takes x A as inputs and outputs m( x A ) that represents the parameters of the invertible "Coupling Transform", C, that will be applied to x B . The inverse map is given by which leads to the simple Jacobian Note that Eq. (9) does not require the computation of the gradient of m( x A ), which would scale as O(D 3 ) with D the number of dimensions. In addition, taking ∂C/∂ 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 A that at most 2 log 2 D number of Coupling Layers are required in order to express non-separable structures of the integrand.

A. Number of Coupling Layers
The minimum number of coupling layers required to capture all the correlations between every dimension of the integration variable, n min , depends on the dimensionality of the integral, D [30]. In the cases of D = 2 and D = 3, each dimension is transformed once based on the other dimension(s) and thus n min = 2 and n min = 3, respectively. This way of counting n min could be generalized to higher D. In fact, this is what autoregressive flows are based on [39]. Here we show that the number of coupling layers required to capture all the correlations is bounded below by 3 and above by 2 log 2 D .
Let x be a D-dimensional input variable that gets split into two subsets, (x A ; x B ) = (x 1 , . . . , x d ; x d+1 , . . . , x D ). After applying three coupling layers, the dependence on other dimensions for each input dimension is given by: Let x be a D-dimensional input variable that gets split into two subsets. After applying three coupling layers, all the input dimensions have been transformed based on the information from one another, at least in an indirect way. However, since the transformation introduced by a single coupling layer is not expressive enough for a generic integrand, the actual number of coupling layers that we use is larger. In addition, the normalizing flow shown above transforms each dimension of x A only once, while the ones in x B are transformed twice. To approximate the integrand equally well in all dimensions, we require that all of them are transformed the same number of times, which increases the minimum number of coupling layers for D ≥ 4 to n min = 4. To systematically construct coupling layers that ensure a direct dependence between all the input dimensions on each other, we require that no set of two (or more) variables are always transformed together and that all the variables must be transformed equal number of times. To meet these conditions, the minimum number of coupling layers can be shown to be n max = 2 log 2 D . The proof is as follows. Since the dimension can either be transformed or not, the first step will be to number each dimension from 0 to D − 1 and take the binary representation of that number with exactly as many bits required to represent the number D − 1. The first layer is found by taking the most significant bit of each number, and transforming those with an "1" and passing through those with a "0". The second layer is done the same as the first, but transforming those with a "0" and passing through those with an "1", which ensures that each dimension is transformed equally. The rest of the layers are found by continuing to iterate through the bits. Thus, the number of layers is given by twice as many bits required to represent the number D − 1 in binary, or 2 log 2 D . See Table I for an example.

B. Using i-flow
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.

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 batch , D), where n batch is the number of points to sample per training step. Secondly, the function must return an array with shape (n 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.

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 A.
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 [30]. This provides the user the flexibility to achieve the expressiveness required for their specific problem.

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 [40]. However, the code can use any optimizer implemented within TensorFlow.

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 [40]), 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 [30], 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. If the latter is used, the number of input bins has to be specified. 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.

Putting it all together
The networks are trained by sampling a fixed number of points using the current state of g(x) [41]. 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 five 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: The result of integrating f 1 from zero to one is given by: In the following, we use α = 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: The result of integrating f 2 from zero to one is given by: In the following, we use α = 0.2. The third case is given by This function has two circles with shifted centers, varying thickness and height. Also, the function exhibits nonfactorizing behavior. The integral of f 3 between 0 and 1 can be computed numerically using Mathematica [42], which is 0.0136848 ± (5 · 10 −9 ), with p 1 = 0.4, p 2 = 0.6, r = 0.25, w = 1/0.004 and a = 3 [43].  The fourth case is a ring function with hard cuts: This function demonstrates how i-flow learns hard, non-factorizing cuts. The result of integrating f 4 from zero to one is given by: π 0.45 2 − 0.2 2 = 0.1625π. 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 → 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 [44], the integral of interest is given by 1 The result of integrating f 5 from zero to one can be obtained through the use of LoopTools [45], which gives a numerical result of 1.9369640238 · 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. Further applications to event generation of high-energy particle collisions is discussed in [46] and also in [47]. These papers investigate using normalizing flows to improve upon phase space integration for event simulation at particle colliders.

V. RESULTS
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 set the number of bins to 100 and use stratified sampling as implemented in [36]. The setup of Foam requires a number of points per cell, here we will fix this to be 1000. In the setup of i-flow, we use 2 log 2 D number of coupling layers with the masking discussed in Sec. III A, and the coupling transform C taken to be a Piecewise Rational Quadratic spline (Appendix A C). 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 · (3n bins + 1), where d is the number of dimensions to be transformed. We further set the number of bins (n bins ) in each output dimension to 16. The learning rate was set to 1 · 10 −3 in all cases except for the Ring function, see details further below.  We used in total five million (5M) points to train the integrators. For VEGAS and i-flow, this number is split into 1000 epochs with 5000 points per epoch, for Foam we used 5000 cells with 1000 points per cell. After we completed the training, we used one million (1M) points to evaluate the integral and report the result in Tab. II. Since this way of presenting results highlights the precision of the integrator, we show in Tab. III the relative deviation from the true value in units of standard deviations, i.e. the accuracy of the result. The relative deviation is defined as where I code is the result from VEGAS, Foam, or i-flow, I true is the true value of the integral, and the ∆I terms signify the uncertainty in the integral. Note, that ∆I true is only non-zero for the case of the entangle circles for which it is 5 · 10 −9 . It is clear from these numbers that for the low-dimensional integrals (D ≤ 4), all three integrators achieve reasonable results. However, i-flow is the only one that shows consistent performance up to 16 dimensions. For example, VEGAS fails in the integration of 16D Camel function completely (missing one of the peaks) and Foam has a large uncertainty on the final result. 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 5000 cells. Therefore, when dealing with high-dimensional integrals, Foam is the least efficient integrator.
Not only does i-flow perform better in high-dimensional integrals, it also yields a smaller uncertainty compared to VEGAS and Foam after the entire training of using 5M points, even though the latter converge much faster at the initial stage. This is a common behaviour observed in i-flow despite the type of functions that undergo training, as illustrated in Figure 3. Interestingly, out of all the examples plotted, VEGAS performs the worst for the 4D Camel and the entangled circles, which is precisely due to its inability to deal with non-factorizable functions. As mentioned before, VEGAS assumes that the integrand must factorize which explains why the uncertainty for the camel function, entangled circles, ring with cuts, and scalar-top-loop all have larger uncertainties than either Foam or i-flow. It is interesting to note, that in the case of the ring with cuts, Foam slightly outperforms i-flow in terms of precision, but not in terms of accuracy.
As discussed in the earlier sections, i-flow also allows efficient sampling once it "learns" the integral up to small uncertainty. As an example, Fig. 4 shows a sample distribution after training i-flow with 5M points (1000 epochs with 5000 points per epoch) on the Ring function of eq. (15). The resulting value for the integral is in Tab. II. For training, we used a learning schedule with exponential decay. An initial learning rate of 2 · 10 −3 is halved every 250 epochs. The cut efficiency is 89.6%. Figure 5 shows the weights of 10000 points sampled after training with 10M points on the Entangled Circles of eq. (14). In the ideal case of g → 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.

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.
The setup suggested in [48] would allow the extension of i-flow to discrete distributions, which has also applications in HEP [46,47]. Another way to implement this type of information is utilizing Conditional Normalizing Flows [49].
The implementation of transflow-learning, which was suggested in [50], 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 [51] have the potential to reduce the memory usage substantially, therefore allowing more points to be used at each training step or larger NN architectures. Currently, if too many points per training step are used, the program crashes due to insufficient space in RAM for such an array. These cases can be circumvented by breaking up the large set of points into several smaller mini batches. The gradients are then computed for each mini batch consecutively before they are applied.
The setup presented in [52], utilizing 1 × 1 convolutions, showed an improved performance over the vanilla implementation of the normalizing flows, which possibly also applies to our case. Additionally, this would modify the maximum number of coupling layers required by having more expressive permutations.  (14). g is a flat distribution before training and approximately resembles the shape of f after training.

A. COUPLING LAYER DETAILS
The implementation of the layers available in i-flow are detailed below. The layers are based on the work of [30,31] and are reproduced here for the convenience of the reader.

A. Piecewise Linear
For the piecewise linear coupling layer [30], given K bins of width w, the probability density function (PDF) is defined as: . . .
The cumulative distribution function (CDF) is defined by the integral giving: where b is the bin in which The inverse CDF is given by: The Jacobian for this network is straight forward to calculate, and gives: The piecewise linear layers require fixed bin widths in each layer. For details on why this is required, see Appendix B of [30].

B. Piecewise Quadratic
For the piecewise quadratic coupling layer [30], given K bins with widths W ik , with K + 1 vertex heights given by V ik , the PDF is defined as: Integrating the above equation leads to the CDF: where b is defined as the solution to [53] b−1 is the relative position of x B i in bin b. Inverting the CDF leads to: where b is defined as the solution to and β is the relative position of C i in the bin b, and is given by:
First, we define the bin widths (w (k) = x (k+1) − x (k) ) and the slopes ( ). We next obtain the fractional distance (ξ) between the two knots that the point of interest (x) lies (ξ = x−x (k) w (k) , where k is the bin x lies in). The CDF is given by: where the details of α(ξ) and β(ξ) can be found in [54], but simplifies to: which is noted to be less prone to numerical issues [54]. The inverse can be found by solving a quadratic equation [31]: where the coefficients are given in [31], solving this equation for the solution that gives a monotonically increasing x results in: where the second form is numerically more precise when 4ac is small, and is also valid for a = 0 [31].