Event Generation and Statistical Sampling for Physics with Deep Generative Models and a Density Information Buffer

Sydney Otten, Sascha Caron, Wieske de Swart, Melissa van Beekveld, Luc Hendriks, Caspar van Leeuwen, Damian Podareanu, Roberto Ruiz de Austri, Rob Verheyen

I Introduction

The simulation of physical and other statistical processes is typically performed in two steps: first, one samples (pseudo)random numbers; in the second step an algorithm transforms these random numbers into simulated physical events. Here, physical events are high energy particle collisions. This is known as the Monte Carlo (MC) method. Currently, a fundamental problem with these numerical simulations is their immense need for computational resources. As such, the corresponding scientific progress is restricted due to the speed of and budget for simulation. As an example, the full pipeline of the MC event generation in particle physics experiments including the detector response may take up to 1010 minutes per event Alwall et al. 2011; Sjöstrand et al. 2015; de Favereau et al. 2014; Corcella et al. 2001; Gleisberg et al. 2009; Belyaev et al. 2013; Kilian et al. 2011 and largely depends on non-optimal MC sampling algorithms such as VEGAS Lepage 1980. Accelerating the event generation pipeline with the help of machine learning can provide a significant speed up for signal studies allowing e.g. broader searches for signals of new physics. Another issue is the inability to exactly specify the properties of the events the simulation produces. Data analysis often requires the generation of events which are kinematically similar to events seen in the data. Current event generators typically accommodate this by generating a large number of events and then selecting the interesting ones with a low efficiency. Events that were not selected in that procedure are often discarded. It is of interest to investigate ways in which the generation of such events can be avoided.

Most of the efforts of the machine learning community regarding generative models are typically not directly aimed at learning the correct frequency of occurrence. So far, applications of generative ML approaches in particle physics focused on image generation Paganini et al. 2018; Erdmann et al. 2018; de Oliveira et al. 2017; Erbin and Krippendorf 2018 due to the recent successes in unsupervised machine learning with generative adversarial networks (GANs) Goodfellow et al. 2014; Radford et al. 2015; Brock et al. 2018 to generate realistic images according to human judgement Salimans et al. 2016; Heusel et al. 2017. GANs were applied to the simulation of detector responses to hadronic jets and were able to accurately model aggregated pixel intensities as well as distributions of high level variables that are used for quark/gluon discrimination and merged jets tagging Musella and Pandolfi 2018. The authors start from jet images and use an Image-to-Image translation technique Isola et al. 2016 and condition the generator on the particle level content.

Soon after the initial preprint of the present article, two relevant papers appeared that model the event generation with GANs Di Sipio et al. 2020; Hashemi et al. 2019. There the authors have achieved an approximate agreement between the true and the generated distributions. Since those papers looked at processes involving two objects such that the generator output was 7 or 8 dimensional, it is still an open question which generative models are able to reliably model processes with a larger number of objects. Additionally, in both papers the authors report difficulties with learning the azimuthal density ϕ\phi which we also target in our studies. In Di Sipio et al. 2020 the authors circumvent the trouble of learning ϕ\phi explicitly with their GAN by learning only Δϕ\Delta\phi, manually sampling ϕj1\phi_{j_{1}} from a uniform distribution and processing the data with an additional random rotation of the system. This further reduces the dimensionality of the studied problems.

In this article we outline an alternative approach to the MC simulation of physical and statistical processes with machine learning and provide a comparison between traditional methods and several deep generative models. All of these processes are characterized by some outcome x\mathbf{x}. The data we use to train the generative models is a collection of such outcomes and we consider them as samples drawn from a probability density p(x)p(\mathbf{x}). The main challenge we tackle is to create a model that learns a transformation from a random variable z→x\mathbf{z}\to\mathbf{x} such that the distribution of x\mathbf{x} follows p(x)p(\mathbf{x}) and enables us to quickly generate more samples.

We investigate several GAN architectures with default hyperparameters and Variational Autoencoders (VAEs) Kingma and Welling 2013 and provide more insights that pave the way towards highly efficient modeling of stochastic processes like the event generation at particle accelerators with deep generative models. We present the B-VAE, a setup of the variational autoencoder with a heavily weighted reconstruction loss and a latent code density estimation based on observations of encoded ground truth data. We also perform a first exploration of its hyperparameter space to optimize the generalization properties.

To test our setup, three different types of data with increasing dimensionality and complexity are generated. In a first step we construct generative models for a 10-dimensional two-body decay toy-model and compare several distributions in the real MC and the generated ML model data. We confirm the recent findings that both GANs and VAEs are generally able to generate events from physical processes. Subsequently, we study two more complex processes:

the 16-dimensional ZZ boson production from e+e−e^{+}e^{-} collisions and its decay to two leptons, e+e−e^{+}e^{-} and μ+μ−\mu^{+}\mu^{-}, with four 4-vectors per data point of which two are always zero.

the 26-dimensional ttˉt\bar{t} production from proton collisions, where at least one of the top quarks is required to decay leptonically with a mixture of five or six final state objects.

The study on ZZ bosons reveals that standard variational autoencoders can’t reliably model the process but confirms good agreement for the B-VAE. For ttˉt\bar{t} we find that by using the B-VAE we are able to produce a realistic collection of events that follows the distributions present in the MC event data. We search for the best B-VAE architecture and explore different possibilities of creating a practical prior by trying to learn the latent code density of encoded ground truth data. We also present results for several GAN architectures with the recommended hyperparameters.

We perform a principal component analysis (PCA) Pearson 1901; Shlens 2014 of encoded ground truth data in the latent space of the VAE for ttˉt\bar{t} production and demonstrate an option to steer the generation of events. Finally, we discuss several further applications of this work including anomaly detection and the utilization for the phase space integration of matrix elements.

In short, the structure of the paper is as follows: In section II. we briefly explain how we create the training data. In section III. we present the methodology. We

present several methods to assess the density of the latent code of a VAE and

define figures of merit that are subsequently used to evaluate our generative models.

In section IV. we present the results. We show

the two-body decay toy model and the leptonic Z-decay,

the ttˉ→4j+1or2lt\bar{t}\to 4j+1\rm{or}2l, where we optimize for several hyperparameters, assess different ways of utilizing latent code densities and show how several GAN architectures with default hyperparameters perform and

two sanity checks on ttˉt\bar{t}: a) Gaussian smearing, creating Gaussian Mixture Models and Kernel Density Estimators for events and b) investigating whether the B-VAE learns the identity function.

In section V. we propose several applications of the B-VAE and provide our conclusions in section VI.

II Monte Carlo Data

We study the generation of physical events using three different sets of generated events: 1. a simple toy-model, 2. ZZ-boson production in e+e−e^{+}e^{-} collisions and 3. top quark production and decay in proton collisions, i.e. pp→ttˉpp\to t\bar{t}. Here, we describe the procedures for obtaining the data sets.

For the toy model we assume a stationary particle with mass MM decaying into two particles with masses m1m_{1} and m2m_{2} and calculate their momentum 4-vectors by sampling m1m_{1}, m2m_{2}, θ\theta and ϕ\phi from uniform distributions 10610^{6} times. θ\theta and ϕ\phi are the polar and azimuthal angle of the direction into which particle 1 travels:

These angles and momentum conservation fix the direction of particle 2. The quantities of the model that are used as training data for the generative models are the energies E1E_{1}, E2E_{2} of the daughter particles, the phase space components pxp_{x}, pyp_{y}, pzp_{z} for each particle and their masses m1m_{1} and m2m_{2}. This introduces a degeneracy with the goal of checking whether the generative models learn the relativistic dispersion relations

II.2 16-dimensional e+​e−→Z→l+​l−e^{+}e^{-}\to Z\to l^{+}l^{-}

We generate 10610^{6} events of the e+e−→Z→l+l−e^{+}e^{-}\to Z\to l^{+}l^{-} (l≡e,μl\equiv e,\mu) process at matrix element level with a center-of-mass energy of 91 GeV using MG5_aMC@NLO v6.3.2 Alwall et al. 2011. The four-momenta of the produced leptons are extracted from the events given in LHEF format Alwall et al. 2007, and are directly used as input data for the generative models. The dimensionality of the input and output data is therefore 16: (Ee−,px,e−,py,e−,pz,e−,Ee+,px,e+,py,e+,pz,e+,\left(E_{e^{-}},p_{x,e^{-}},p_{y,e^{-}},p_{z,e^{-}},E_{e^{+}},p_{x,e^{+}},p_{y,e^{+}},p_{z,e^{+}},\right. Eμ−,px,μ−,py,μ−,pz,μ−,Eμ+,px,μ+,py,μ+,pz,μ+)\left.E_{\mu^{-}},p_{x,\mu^{-}},p_{y,\mu^{-}},p_{z,\mu^{-}},E_{\mu^{+}},p_{x,\mu^{+}},p_{y,\mu^{+}},p_{z,\mu^{+}}\right) and will always contain 8 zeros, since the events consist of e+e−e^{+}e^{-} or μ+μ−\mu^{+}\mu^{-}.

II.3 26-dimensional p​p→t​t¯pp\to t\bar{t}

We generate 1.2⋅1061.2\cdot 10^{6} events of pp→ttˉpp\to t\bar{t}, where at least one of the top-quarks is required to decay leptonically. We used MG5_aMC@NLO v6.3.2 Alwall et al. 2011 for the matrix element generation, using the NNPDF PDF set Ball et al. 2017. Madgraph is interfaced to Pythia 8.2 Sjöstrand et al. 2015, which handles showering and hadronization. The matching with the parton shower is done using the MLM merging prescription Mangano et al. 2003. Finally, a quick detector simulation is done with Delphes 3 de Favereau et al. 2014; Cacciari et al. 2012, using the ATLAS detector card. For all final state objects we use (E,pT,η,ϕ)(E,p_{T},\eta,\phi) as training data and also include MET and METϕ\phi. We have 5 or 6 objects in the final state, four jets and one or two leptons, i.e. our generative models have a 26 dimensional input and output, while those with only one lepton contain 4 zeros at the position of the second lepton.

III Methods

This section summarizes the methodology used to investigate deep and traditional generative models to produce a realistic collection of events from a physical process. We present the generative techniques we have applied to the data sets, GANs and VAEs, with a focus on our technique: an explicit probabilistic model, the B-VAE, which is a method combining a density information buffer with a variant of the VAE. Subsequently we discuss several traditional methods to learn the latent code densities and finally, we present figures of merit to assess the performance of the generative models.

In this section we give a brief description of GANs and a thorough description of our B-VAE technique. For the latter, we provide the details of the corresponding architecture as well as hyperparameters and training procedures, and also show how the density information buffer is created and how it is utilized to generate events. The GANs and VAEs are trained on an Nvidia Geforce GTX 970 and a Tesla K40m GPU using tensorflow-gpu 1.14.0 Abadi et al. 2016, Keras 2.2.5 Chollet et al. and cuDNN 7.6.1 Chetlur et al. 2014.

GANs learn to generate samples from a data distribution by searching for the global Nash equilibrium in a two-player game. The two players are neural networks: one that tries to generate samples that convinces the other, a discriminator that tries to distinguish real from fake data. There are many possibilities to realize this, accompanied by large hyperparameter spaces. We try to create event generators with several recent GAN architectures:

the regular Jensen-Shannon GAN Goodfellow et al. 2014 that was only applied to the first toy model dataset,

several more recent GAN architectures that have been applied to the ttˉt\bar{t} dataset: Wasserstein GAN (WGAN) Arjovsky et al. 2017, WGAN with Gradient Penalty (WGAN-GP) Gulrajani et al. 2017, Least Squares GAN (LSGAN) Mao et al. 2016, Maximum Mean Discrepancy GAN (MMDGAN) Li et al. 2017.

We use the recommended hyperparameters from the corresponding papers. Note here that an extensive hyperparameter scan may yield GANs that perform better than those reported in this paper. The performance of the GAN models we present serve as baselines.

Explicit probabilistic models

Consider that our data, N particle physics events X={xi}i=1N\mathbf{X}=\{\mathbf{x}^{i}\}_{i=1}^{N}, are the result of a stochastic process and that this process is not known exactly. It depends on some hidden variables called latent code z\mathbf{z}. With this in mind one may think of event generation as a two-step process: (1) sampling from a parameterized prior pθ(z)p_{\mathbf{\theta}}(\mathbf{z}) (2) sampling xi\mathbf{x}^{i} from the conditional distribution pθ(xi∣z)p_{\mathbf{\theta}}(\mathbf{x}^{i}|\mathbf{z}), representing the likelihood. For deep neural networks the marginal likelihood

is often intractable. Learning the hidden stochastic process that creates physical events from simulated or experimental data requires us to have access to an efficient approximation of the parameters θ\mathbf{\theta}. To solve this issue an approximation to the intractable true posterior pθ(z∣x)p_{\mathbf{\theta}}(\mathbf{z}|\mathbf{x}) is created: a probabilistic encoder qϕ(z∣x)q_{\mathbf{\phi}}(\mathbf{z}|\mathbf{x}). Given a data point xix^{i} it will produce a distribution over the latent code z\mathbf{z} from which the data point might have been generated. Similarly, a probabilistic decoder pθ(x∣z)p_{\mathbf{\theta}}(\mathbf{x}|\mathbf{z}) is introduced that produces a distribution over possible events xi\mathbf{x}^{i} given some latent code z\mathbf{z}. In this approach the encoder and decoder are deep neural networks whose parameters ϕ\mathbf{\phi} and θ\mathbf{\theta} are learned jointly.

can be written as a sum of the likelihood of individual data points. Using that

and applying Jensen’s inequality, one finds that

For this situation, one must substitute p(z∣xi)→qϕ(z∣xi)p(\mathbf{z}|\mathbf{x}^{i})\rightarrow q_{\phi}(\mathbf{z}|\mathbf{x}^{i}) since we don’t know the true posterior but have the approximating encoder qϕ(z∣xi)q_{\phi}(\mathbf{z}|\mathbf{x}^{i}). From here one can derive the variational lower bound L(θ,ϕ;xi)\mathcal{L}(\theta,\phi;\mathbf{x}^{i}) and find that

DKLD_{\text{KL}} measures the distance between the approximate and the true posterior and since DKL≥0D_{\text{KL}}\geq 0, L(θ,ϕ;xi)\mathcal{L}(\theta,\phi;\mathbf{x}^{i}) is called the variational lower bound of the marginal likelihood

We optimize L(θ,ϕ;xi)\mathcal{L}(\theta,\phi;\mathbf{x}^{i}) with respect to its variational and generative parameters ϕ\phi and θ\theta.

Using the Auto-Encoding Variational Bayes Algorithm Kingma and Welling 2013 (AEVB) a practical estimator of the lower bound is maximized: for a fixed qϕ(z∣x)q_{\phi}(\mathbf{z}|\mathbf{x}) one reparametrizes z^∼qϕ(z∣x)\hat{\mathbf{z}}\sim q_{\phi}(\mathbf{z}|\mathbf{x}) using a differentiable transformation gϕ(ϵ,x),ϵ∼N(0,1)g_{\phi}(\epsilon,\mathbf{x}),\epsilon\sim\mathcal{N}(0,1) with an auxiliary noise variable ϵ\epsilon. Choosing z∼p(z∣x)=N(μ,σ2)\mathbf{z}\sim p(\mathbf{z}|\mathbf{x})=\mathcal{N}(\mu,\sigma^{2}) with a diagonal covariance structure, such that

where μ\mu and σ2\sigma^{2} are outputs of the encoding deep neural network. Reparametrizing z=μ+σ⊙ϵ\mathbf{z}=\mu+\sigma\odot\epsilon yields the Variational Autoencoder (VAE) Kingma and Welling 2013. In that case the first term in eq. (8) can be calculated analytically:

The second term in eq. (8) corresponds to the negative reconstruction error that, summed over a batch of samples, is proportional to the mean squared error (MSE) between the input xi\mathbf{x}^{i} and its reconstruction given the probabilistic encoder and decoder. The authors in Kingma and Welling 2013 state that for batch-sizes M>100M>100 it is sufficient to sample ϵ\epsilon once which is adopted in our implementation. By calculating the lower bound for a batch of MM samples XM⊂X\mathbf{X}^{M}\subset\mathbf{X} they construct the estimator of L\mathcal{L}:

We use the gradients ∇θ,ϕLM(θ,ϕ;XM,ϵ)\nabla_{\theta,\phi}\mathcal{L}^{M}(\theta,\phi;\mathbf{X}^{M},\epsilon) for the SWATS optimization procedure Shirish Keskar and Socher 2017, beginning the training with the Adam optimizer Kingma and Ba 2014 and switching to stochastic gradient descent. Practically, the maximization of the lower bound is turned into the minimization of the positive DKLD_{\text{KL}} and the MSE such that the loss function of the VAE L∝DKL+MSEL\propto D_{\text{KL}}+MSE. In our approach we introduce a multiplicative factor BB for DKLD_{\text{KL}} to tune the relative importance of both terms. The authors in Burgess et al. 2018 introduce a similar factor β\beta, but their goal is to disentangle the latent code by choosing β>1\beta>1 such that each dimension is more closely related to features of the output. In contrast we choose B≪1B\ll 1 to emphasize a good reconstruction. The loss function of the VAE can subsequently be written as

This however also implies that DKL(qϕ(z∣x)∥pθ(z))D_{\text{KL}}(q_{\phi}(\mathbf{z}|\mathbf{x})\|p_{\theta}(\mathbf{z})) is less important, i.e. there is a much smaller penalty when the latent code distribution deviates from a standard Gaussian. This incentivizes narrower Gaussians because the mean squared error for a single event reconstruction grows as the sampling of the Gaussians in latent space occur further from the mean. Note that for B=0B=0 one obtains the same loss function as for a standard autoencoder Rumelhart et al. 1986: the reconstruction error between input and output. Although the standard deviations will be small, the VAE will still maintain its explicit probabilistic character because it contains probabilistic nodes whose outputs are taken to be the mean and logarithmic variance, while the standard autoencoder doesn’t.

The encoders and the decoders of our VAEs have the same architectures consisting of four (toy model and Z→l+l−Z\to l^{+}l^{-}) or six hidden layers (pp→ttˉpp\to t\bar{t}) with 128 neurons each and shortcut connections between every other layer He et al. 2015; Huang et al. 2016. We choose B=3⋅10−6B=3\cdot 10^{-6} for the toy model and Z→l+l−Z\to l^{+}l^{-}. The number of latent space dimensions are 9 for the toy model and 10 for Z→l+l−Z\to l^{+}l^{-}. For the toy model we use a simple training procedure using the Adam optimizer with default values for 100 epochs. For Z→l+l−Z\to l^{+}l^{-} we employ a learning rate scheduling for 7×807\times 80 epochs and SWATS Shirish Keskar and Socher 2017, i.e. switching from Adam to SGD during training.

For pp→ttˉpp\to t\bar{t} we perform a scan over hyperparameters with

We perform this scan on a small training data set with 10510^{5} samples. We use a batch-size of 1024 and the exponential linear unit (ELU) Clevert et al. 2015 as the activation function of hidden layers. The output layer of the decoder is a hyperbolic tangent such that we need to pre- and postprocess the input and output of the VAE. We do this by dividing each dimension of the input by the maximum of absolute values found in the training data. We apply this pre- and post-processing in all cases. We initialize the hidden layers following a normal distribution with mean 0 and a variance of (1.55/128)0.5(1.55/128)^{0.5} such that the variance of the initial weights is approximately equal to the variance after applying the activation function on the weights Clevert et al. 2015. For pp→ttˉpp\to t\bar{t} the setup is identical except for the number of epochs: we train 4×2404\times 240 epochs with Adam and then for 4×1204\times 120 epochs with SGD. Due to the increasing complexity of the data sets we perform more thorough training procedures.

III.2 Latent Code Density Estimation

In the case of VAEs the prior p(z)=N(0,1)p(\mathbf{z})=\mathcal{N}(0,1) used to sample pθ(x∣z)p_{\theta}(\mathbf{x}|\mathbf{z}) isn’t identical to the distribution over the latent code z\mathbf{z} resulting from the encoding of true observations qϕ(z∣X)q_{\phi}(\mathbf{z}|\mathbf{X}). The generated distribution over x\mathbf{x} given pθ(x∣z)p_{\theta}(\mathbf{x}|\mathbf{z}) therefore doesn’t match the reference when assuming a unit Gaussian over z\mathbf{z}. We address this issue by estimating the prior p(z)p(\mathbf{z}) for the probabilistic decoder from data using a strategy similar to the Empirical Bayes method Robbins 1956. We collect observations Z={z1,…,zm}\mathbf{Z}=\{\mathbf{z}_{1},\ldots,\mathbf{z}_{m}\} by sampling qϕ(z∣XL)q_{\phi}(\mathbf{z}|\mathbf{X}_{\text{L}}) where XL⊂X\mathbf{X}_{\text{L}}\subset\mathbf{X} is a subset of physical events. Z\mathbf{Z} is then used as the data for another density estimation to create a generative model for p(z)p(\mathbf{z}): this is what we call the density information buffer. This buffer is used in several ways to construct p(z)p(\mathbf{z}): we apply Kernel Density Estimation Parzen 1962, Gaussian Mixture Models using the expectation maximization algorithm College and Dellaert 2002, train a staged VAE Dai and Wipf 2019 and directly use the density information buffer. Note that the Kernel Density Estimation and the Gaussian Mixture Models are also used in another context, namely in the attempt to construct such a traditional generative model that is optimized on physical events instead of the latent code as suggested here.

Given N samples from an unknown density pp the kernel density estimator (KDE) p^\hat{p} for a point yy is constructed via

where the bandwidth hh is a smoothing parameter that controls the trade-off between bias and variance. Our experiments make use of N={104,105}N=\{10^{4},10^{5}\}, a Gaussian kernel

and have optimised hh. We use the KDE implementation of scikit-learn Pedregosa et al. 2011 that offers a simple way to use the KDE as a generative model and optimize the bandwidth hh using GridSearchCV and a 5-fold cross validation for 20 samples for hh distributed uniformly on a log-scale between 0.1 and 10.

Gaussian Mixture Models

Since the VAE also minimizes DKL(qϕ(z∣x)∥N(0,1))D_{\text{KL}}(q_{\phi}(\mathbf{z}|\mathbf{x})\|\mathcal{N}(0,1)) it’s incentivized that even with low values of β\beta the latent code density qϕ(z∣XL)q_{\phi}(\mathbf{z}|\mathbf{X}_{\text{L}}) is similar to a Gaussian. It therefore appears promising that a probabilistic model that assumes a finite set of Gaussians with unknown parameters can model the latent code density very well. We use the Gaussian Mixture Model as implemented in scikit-learn, choosing the number of components to be {50,100,1000}\{50,100,1000\} with the full covariance matrix and 10510^{5} encodings zi∈Z\mathbf{z}_{i}\in\mathbf{Z}.

Two-Stage VAE

The idea of the two-stage VAE (S-VAE) is to create another probabilistic decoder pη(z∣z′)p_{\eta}(\mathbf{z}|\mathbf{z}^{\prime}) from latent code observations Z\mathbf{Z} that is sampled using p(z′)=N(0,1)p(\mathbf{z}^{\prime})=\mathcal{N}(0,1) Dai and Wipf 2019. We use a lower neural capacity for this VAE with three hidden layers with 64 neurons each without shortcut connections for each neural network, and use B={10−6,10−5,…,1}B=\{10^{-6},10^{-5},\ldots,1\}. We slightly modify the loss function from eq. (13) and remove the (1−B)(1-B) in front of the MSE term because we want to test higher values of BB of up to B=1B=1 and don’t want to completely neglect the MSE. Every other hyperparameter including the training procedure is identical to those in the VAE for pp→ttˉpp\to t\bar{t}. It is straightforward to expand this even further and also apply KDE, or create a density information buffer from the latent codes z′\mathbf{z}^{\prime} to then sample pη(z∣z′)p_{\eta}(\mathbf{z}|\mathbf{z}^{\prime}) with the data-driven prior.

Density Information Buffer

Another way to take care of the mismatch between p(z)p(\mathbf{z}) and qϕ(z)q_{\phi}(\mathbf{z}) is to explicitly construct a prior pϕ,XL(z)p_{\phi,\mathbf{X}_{\text{L}}}(\mathbf{z}) by aggregating (a subset of) the encodings of the training data:

Practically this is done by saving all μi\mu^{i} and σ2,i\sigma^{2,i} for all mm events in XL\mathbf{X}_{\text{L}} to a file, constituting the buffer. The advantage of this procedure is that the correlations are explicitly conserved by construction for the density information buffer while the KDE, GMM and the staged VAE may only learn an approximation of the correlations in z∼qϕ(z∣XL)\mathbf{z}\sim q_{\phi}(\mathbf{z}|\mathbf{X}_{\text{L}}). A disadvantage of this approach is that the resulting density is biased towards the training data, in the sense that the aggregated prior is conditioned on true observations of the latent code for the training data and has a very low variance when BB in eq. (13) is small. One can interpret this as overfitting to the data with respect to the learned density. To counter this effect, we introduce a smudge factor α\alpha such that we sample zi∼N(μi,ασ2,i) ∀ xi∈XL\mathbf{z}^{i}\sim\mathcal{N}(\mu^{i},\alpha\sigma^{2,i})\,\forall\,\mathbf{x}^{i}\in\mathbf{X}_{\text{L}}. In our experiments we investigate α={1,5,10}\alpha=\{1,5,10\} and only apply α\alpha if σ<σT=0.05\sigma<\sigma_{T}=0.05, such that

with μi\mu^{i} and σ2,i\sigma^{2,i} being the Gaussian parameters in latent space corresponding to events i=1,…,mi=1,\ldots,m in XL\mathbf{X}_{\text{L}}. It is straightforward to expand this approach to have more freedom in α\alpha, e.g. by optimizing (αj)j=1dim⁡z(\alpha_{j})_{j=1}^{\dim\mathbf{z}}, a smudge factor for each latent code dimension. One can include more hyperparameters that can be optimized with respect to figures of merit. By introducing a learnable offset γj\gamma_{j} for the standard deviation such that

we have 2⋅dim⁡z2\cdot\dim\mathbf{z} additional hyperparameters. More generally we can try to learn a vector-valued function γ(ρ(z))\mathbf{\gamma}(\rho(\mathbf{z})) that determines the offset depending on the local point density in latent space. While all of these approaches may allow a generative model to be optimized, it introduces a trade-off by requiring an additional optimization step that increases in complexity with increasing degrees of freedom. In our experiments we only require γ=γj={0.01,0.05,0.1}\gamma=\gamma_{j}=\{0.01,0.05,0.1\} to be the minimal standard deviation, such that

III.3 Figures of Merit

Having discussed a number of candidate generative models, we now define a method of ranking these models based on their ability to reproduce the densities encoded in the training data. While work that was done so far predominantly relies on χ2\chi^{2} between observable distributions and pair-wise correlations Di Sipio et al. 2020; Hashemi et al. 2019, we aim to capture the generative performance more generally. Starting from a total of 1.2⋅1061.2\cdot 10^{6} Monte Carlo samples, 10510^{5} of those samples are used as training data for the generative models. We then produce sets of 1.2⋅1061.2\cdot 10^{6} events with every model and compare to the Monte Carlo data.

The comparison is carried out by first defining a number of commonly used phenomenological observables. They are the MET, METϕ\phi and EE, pTp_{T}, η\eta and ϕ\phi of all particles, the two-, three- and four jet, four jet plus one and two lepton mass, the angular distance between the leading and subleading jet

and the azimuthal distance between the leading lepton and the MET. An object is considered to be leading if it has the highest energy in its object class. All possible 2D histograms of these 33 observables are then set up for all models, and are compared with those of the Monte Carlo data. We create 2D histograms with N2D,bins=252N_{\text{2D,bins}}=25^{2}. We then define

with ii and jj summing over observables and Nhist=528N_{\text{hist}}=528 which represent averages over all histograms of the test statistic

where pup_{u} and puMCp_{u}^{\text{MC}} are the normalized bin contents of bins uu Porter 2008. We have verified that the Kullback-Leibler divergence and the Wasserstein distance lead to identical conclusions. This figure of merit is set up to measure correlated performance in bivariate distributions.

A second figure of merit is included with the goal of measuring the rate of overfitting on the training data. If any amount of overfitting occurs, the generative model is expected to not properly populate regions of phase space that are uncovered by the training data. However, given a large enough set of training data, these regions may be small and hard to identify. To test for such a phenomenon, we produce 2D histograms in ηj1\eta_{j_{1}} vs. ηj2\eta_{j_{2}}, ηi∈[−2.5,2.5]\eta_{i}\in[-2.5,2.5] with a variable number of bins NbinsN_{\text{bins}} and measure the fraction of empty bins fe(Nbins)f_{\text{e}}\left(\sqrt{N_{\text{bins}}}\right). As NbinsN_{\text{bins}} is increased, the histogram granularity probes the rate of overfitting in increasingly smaller regions of the phase space. The function fe(Nbins)f_{\text{e}}\left(\sqrt{N_{\text{bins}}}\right) also depends on the underlying probability distribution, and is thus only meaningful in comparison to the Monte Carlo equivalent feMC(Nbins)f^{\text{MC}}_{\text{e}}\left(\sqrt{N_{\text{bins}}}\right). We therefore compute the figure of merit as

After searching the best model with respect to ηj1\eta_{j_{1}} vs. ηj2\eta_{j_{2}}, we independently check its behavior in δOF(ϕj1,ϕj2)\delta_{\text{OF}}(\phi_{j_{1}},\phi_{j_{2}}). The performance, indicated by both figures of merit, is better the lower the value is. We derive our expectations for the values of δ\delta and δOF\delta_{\text{OF}} for 1.2⋅1061.2\cdot 10^{6} events by evaluating these figures of merit on MC vs. MC data. Since we only have 1.2 million events in total, we can only evaluate the figures of merit for up to 6⋅1056\cdot 10^{5} events and then extrapolate as shown in Fig. 5(c). We expect that an ideal generator has δ≈0.0003\delta\approx 0.0003 and δOF≈0.5\delta_{\text{OF}}\approx 0.5 for 1.2⋅1061.2\cdot 10^{6} events. For δ\delta, this value is obtained by fitting the parameters A,B,C,DA,B,C,D to the empirically motivated function

for values of δ\delta obtained by evaluating the MC dataset with N={50000,100000,…,600000}N=\{50000,100000,\ldots,600000\} and extrapolating to N=1200000N=1200000 using the curve_fitcurve\_fit function in scipy. The expected δOF\delta_{\text{OF}} is assumed to be roughly flat around 0.5. We select our best model by requiring δOF⪅1\delta_{\text{OF}}\lessapprox 1 while minimizing δ\delta.

IV Results

We study the behavior of generative models on the three different data sets described in II. Most of the conducted studies focus on the ttˉt\bar{t} dataset beginning in IV.3 and we only present short, preliminary studies on the two-body and the leptonic Z decay in IV.1 and IV.2. Our study finds that by using the B-VAE, we are able to capture the underlying distribution such that we can generate a collection of events that is in very good agreement with the distributions found in MC event data with 12 times more events than in the training data. Our study finds that many GAN architectures with default parameters and the standard VAE do not perform well. The best GAN results in this study are achieved by the DijetGAN Di Sipio et al. 2020 with the implementation as delivered by the authors and the LSGAN Mao et al. 2016 with the recommended hyperparameters. The failure of the standard VAE is accounted to the fact that the distributions of encoded physical events in latent space is not a standard normal distribution. We find that the density information buffer can circumvent this issue. To this end we perform a brief parameter scan beyond dim⁡z\dim\mathbf{z} and BB for the smudge factors α\alpha and offsets γ\gamma. The performance of the optimized B-VAE is presented in Figures 3 to 7. Additionally, we investigate whether improvements to the density information buffer can be achieved by performing a Kernel Density Estimation, creating a Gaussian Mixture Model or learning the latent code density with another VAE. Finally, we perform sanity checks with the ttˉt\bar{t} dataset in IV.4, obtain benchmark performances from traditional methods and test whether our proposed method is trivial, i.e. whether it is only learning the identity function.

The comparison of the generative model performances for the toy model in Fig. 1(a) indicates that the B-VAE with an adjusted prior, given in Eq. 16, is the best investigated ML technique that is able to reliably model the pxp_{x}, pyp_{y} and pzp_{z} distributions when compared to regular GANs and VAEs with a standard normal prior, although these models still give good approximations. We find that all models learn the relativistic dispersion relation which underlines the findings in Wu and Tegmark 2018; Iten et al. 2018. It is noteworthy that for this data set, we only try regular GANs with small capacities and find that they can already model the distributions reasonably well. We confirm the findings in Hashemi et al. 2019; Di Sipio et al. 2020 that it is problematic for GANs to learn the uniform distribution in ϕ\phi. While it is one of the few deviations that occur in Hashemi et al. 2019, Di Sipio et al. 2020 circumvents the issue with ϕ\phi by only learning Δϕ\Delta\phi between the two jets and manually sampling ϕj1∼U(−π,π)\phi_{j_{1}}\sim U(-\pi,\pi). It is questionable whether this technique can be generalized to higher multiplicities.

IV.2 e+​e−→Z→l+​l−e^{+}e^{-}\rightarrow Z\rightarrow l^{+}l^{-}

Fig. 1(b) and 2 show the results for the ZZ events, where the ZZ boson decays leptonically. Here we find that the B-VAE is able to accurately generate events that respect the probability distribution of the physical events. We find very good agreement between the B-VAE and physical events for distributions of pTp_{T}, θ\theta and ϕ\phi and good agreement for invariant mass MinvM_{\text{inv}} of the lepton pair around 91 GeV. While the standard VAE fails for the momentum conservation beside having a peak around 0, the B-VAE is much closer to the true distribution. When displaying ϕ\phi, θ\theta and the transverse momentum pTp_{T} of lepton 1 against lepton 2 (Fig. 2), we find good agreement for the B-VAE, while the standard VAE results in a smeared out distribution. In addition, it can be seen that the events generated by the standard VAE are not always produced back to back but are heavily smeared. We conclude that if we do not use density information buffering, a standard VAE is not able to accurately generate events that follow the Monte Carlo distributions. In particular, events with four leptons are sometimes generated if no buffering is used.

IV.3 p​p→t​t¯→4​jets+1​or​ 2​leptonspp\rightarrow t\bar{t}\rightarrow 4\,\rm{jets+}1\,\rm{or}\,2\,\rm{leptons}

Here we present and discuss the results for the more complicated ttˉt\bar{t} production with a subsequent semi-leptonic decay. We train the generative models on events that have four jets and up to two leptons in the final state such that their input and output dimension is 26. For simplicity we do not discriminate between bb-jets and light-flavored jets, nor between different kinds of leptons. A jet is defined as a clustered object that has a minimum transverse momentum (pTp_{T}) of 20 GeV in the Monte Carlo simulation. We first explore the hyperparameter space of the B-VAE in dim⁡z\dim\mathbf{z}, BB, α\alpha, γ\gamma and recommend a best practice for the creation of a generative model for physical events. Subsequently we investigate various methods to learn the latent code density of encoded ground truth data. Finally we try to create a generative model for physical events with several GAN architectures.

Tables 1 and 2 show the top-15 performances of (dim⁡z,B,α,γ)(\dim\mathbf{z},B,\alpha,\gamma) combinations evaluated on the figures of merit defined in section III.3. For all possible combinations of dim⁡z\dim\mathbf{z} and BB as defined in section III.1 we have separately investigated

For the γ\gamma-study we fixed α=1\alpha=1 and for the α\alpha-study we fixed γ=0\gamma=0. Tables 1 and 2 show the ranking in δ1D\delta_{\text{1D}} for the studies on γ\gamma and α\alpha respectively.

It is not surprising that the best performance in δ1D\delta_{\text{1D}} is attained by the B-VAE with the highest latent code dimensionality, the lowest BB and α=1,γ=0\alpha=1,\gamma=0. The downside however is a very poor performance in δOF\delta_{\text{OF}}. Comparing to the values for the 5%5\% Gaussian smearing of events in Table 5, they are very similar in δ\delta but even worse in δOF\delta_{\text{OF}} and thus, this model provides no advantage over simple smearing without using machine learning techniques: it essentially learns to reproduce the training data. We observe similar patterns for the ranking in δ\delta: the models that perform best only provide a small advantage. Other models do provide a bigger advantage but there is a trade-off between performance in δ\delta and δOF\delta_{\text{OF}} that can in principle be weighted arbitrarily. By introducing the factor α\alpha we smear the B-VAE events in latent space. Models with neither smearing nor an offset perform poorly in δOF\delta_{\text{OF}}, whereas models with B>10−5B>10^{-5} perform poorly in δ\delta. For illustrative purposes we proceed to show and discuss details for the model we consider best: dim⁡z=20,B=10−6,α=1,γ=0.05\dim\mathbf{z}=20,B=10^{-6},\alpha=1,\gamma=0.05. Fig. 3 shows the comparison between B-VAE events and ground truth data in 29 one-dimensional histograms for this model:

E,pT,ηE,p_{T},\eta and ϕ\phi for all four jets and the leading lepton,

Δϕ\Delta\phi between MET and leading lepton,

Δ\DeltaR between leading and subleading jets and

the invariant mass MinvM_{\text{inv}} for 2, 3 and 4 jets and 4 jets + 1 and 2 leptons.

Note that the training data and the density information buffer consist of the same 10510^{5} samples that were used to generate 1.2⋅1061.2\cdot 10^{6} events which are compared to 1.2⋅1061.2\cdot 10^{6} ground truth samples. We observe that the ground truth and generated distributions generally are in good agreement. For the invariant masses we again observe deviations in the tail of the distribution. For MET, METϕ\phi, Δϕ\Delta\phi and Δ\DeltaR we see almost perfect agreement.

Generating 10710^{7} ttˉt\bar{t} events with the VAE has taken 177.5 seconds on an Intel i7-4790K and is therefore several orders of magnitude faster than the traditional MC methods.

Fig. 4 shows eight histograms of ϕ\phi of the leading jet vs. ϕ\phi of the next to leading jet (ϕ1\phi_{1} vs ϕ2\phi_{2}) that were created using the B-VAE with dim⁡z=20,B=10−6,α=1,γ=0.05\dim\mathbf{z}=20,B=10^{-6},\alpha=1,\gamma=0.05. The left column shows the histogram for the full range [−π,π]×[−π,π][-\pi,\pi]\times[-\pi,\pi] whereas the right column shows the same histogram zoomed in on ×\times. The first row displays the training data consisting of 10510^{5} events. The second and third row of Fig. 4 show 1.2⋅1061.2\cdot 10^{6} ground truth and B-VAE events respectively allowing for a comparison of how well the B-VAE generalizes considering it was trained on only 10510^{5} events. The amount of empty bins (holes) present for the ground truth and B-VAE events is very similar. Also the general features of the generated distribution are in very good agreement with the ground truth. However, one can spot two shortcomings:

the presented model smears the detector granularity that is visible in ϕ\phi due to the γ\gamma parameter which would be learned for α=1\alpha=1 and γ=0\gamma=0 and

generator artefacts appear around (±π,0)(\pm\pi,0) and (0,±π)(0,\pm\pi). For EE, pTp_{T} and η\eta we observe larger deviations in the tails of the distributions while for ϕ\phi we only observe slightly more events produced around ±π\pm\pi.

The first effect is most likely due to the γ\gamma parameter and the second effect was already expected from the deviations in the one-dimensional azimuthal distributions around ±π\pm\pi.

Fig. 5 shows how the fraction of empty bins evolves with respect to the number of bins in 2D histograms of ηj1\eta_{j_{1}} vs. ηj2\eta_{j_{2}} and ϕj1\phi_{j_{1}} vs. ϕj2\phi_{j_{2}} for several models including the ground truth. One can see that our chosen model, whose performance was presented in Figures 3 and 4, also accurately follows the fraction of empty bins of the Monte Carlo data.

As discussed in section III.2 we compare four different methods for constructing a prior for the generative model. We compare a KDE, three GMMs and several S-VAEs to the explicit latent code density of encoded ground truth data. To demonstrate this we choose the same B-VAE model as in the preceding paragraph: (20,10−6,1,0.05)(20,10^{-6},1,0.05). Figure 7 shows histograms of all 20 latent code dimensions coming from the different approaches. We observe that all dimensions are generally modelled well by all approaches, except for the S-VAEs with extreme values of BB. This is an expected result since the encoder qϕ(z∣x)q_{\phi}(\mathbf{z}|\mathbf{x}) transforms the input into multivariate Gaussians for which a density estimation is much easier than for such non-Gaussian densities present in physical events. Table 3 shows the performance of the different approaches.

It is remarkable that the KDE and GMM models of the prior p(z)p(\mathbf{z}) provide such good performance in δ\delta, especially the GMM with 1000 components. A drawback for all of the models that try to learn the latent code density is that the resulting performance in δ\delta and δOF\delta_{\text{OF}} is very poor when compared to the explicit use of the density information buffer.

We compare several state of the art GAN architectures in Table 4.

Table 4 shows the evaluation of the GAN models on our figures of merit. However, we find that no GAN architecture we tried is able to provide a satisfactory performance with respect to δ\delta and that all of the tried architectures perform worse than traditional methods such as KDE and GMM except for the LSGAN. The best GAN we find is the LSGAN that, in contrast to all GANs we try otherwise, outperforms all traditional and several B-VAE models with respect to δOF\delta_{\text{OF}}. Fig. 6 shows the loss curves for the GAN architectures that are also shown in Table 4 and the B-VAE.

Considering the GAN literature, the results found are not surprising; the authors in Di Sipio et al. 2020; Hashemi et al. 2019 report difficulties when trying to learn ϕ\phi. Several other papers report that it is very difficult or technically unfeasible to learn densities with GANs Mescheder et al. 2017; Arora and Zhang 2017; Fisher et al. 2018; Abbasnejad et al. 2019. Some of these papers even show that the regular GAN and the WGAN can even fail to learn a combination of 2D Gaussians and that they are not suited to evaluate densities by design Abbasnejad et al. 2019.

Note that all the GAN models we have tried here were trained using the hyperparameters that were recommended in the corresponding papers. However, each of these models is accompanied by large hyperparameter spaces that impact the performance of the generator. The poor performance we find for most GAN models therefore does not imply that GANs are ruled out as potential generative models.

IV.4 Sanity Checks

We perform two sanity checks: (1) we show that two traditional density learning algorithms, Kernel Density Estimation and Gaussian Mixture Models, do not work well when applied directly on the events. (2) we check whether the VAE learns the identity function. Both checks are performed on the ttˉt\bar{t} data.

We perform a KDE with an initial grid-search as described in III.2 to find the optimal bandwidth on a reduced data set with 10410^{4} samples and then perform a KDE with hopth_{\text{opt}} on 10510^{5} samples. Additionally, we create a GMM of those 10510^{5} samples with 50,10050,100 and 10001000 components with a maximum of 500 iterations. Subsequently we generate 1.2⋅1061.2\cdot 10^{6} samples from the KDE and the three GMM models and evaluate them with our figures of merit δ\delta and δOF\delta_{\text{OF}} as presented in Table 5. Additionally we take 10510^{5} events and smear them by sampling from a Gaussian around these events. To this end, we pre-process them in the same way as above and multiply every dimension of every event with N(1,σ2={0.05,0.1})\mathcal{N}\left(1,\sigma^{2}=\{0.05,0.1\}\right) and sample 12 times per event. Table 5 generally shows a poor performance of all models, especially for δOF\delta_{\text{OF}}. Only the smearing shows good performance for δ\delta This procedure however does not respect the correlations in the data and therefore also performs poorly for δOF\delta_{\text{OF}}.

When dim⁡z\dim\mathbf{z} is greater than or equal to the number of dimensions of the training data, it becomes questionable whether a VAE is merely learning the identity function, i.e. whether

Since qϕ(z∣xi)q_{\phi}(\mathbf{z}|\mathbf{x}^{i}) always had non-zero variance, no delta functions occur practically. However, one can notice a bias in some variables when feeding random uniform noise xtest∼U(0,1)\mathbf{x}_{\text{test}}\sim U(0,1) into the VAE. This is no surprise since the encoder and decoder are constructed to learn a function that can reconstruct the input. In Figure 8 we show the reconstructions for the 26-dimensional ttˉt\bar{t} events of a VAE with a 20-dimensional latent space and B=10−6B=10^{-6} and the reconstructions of the same VAE for xtest∼U(0,1)\mathbf{x}_{\text{test}}\sim U(0,1), where we clearly see that the VAE does not simply learn the identity function. The parameters α,B,γ\alpha,B,\gamma and dim⁡z\dim\mathbf{z} allow one to tune how the B-VAE generalizes.

V Applications

We have found that the B-VAE as a deep generative model can be a good generator of collision data. In this section we discuss several further applications of this work such as anomaly detection and improved MC integration. We demonstrate the option of how one can utilize the B-VAE to steer the event generation.

To steer the event generation we need to find out which regions in latent space correspond to which events generated by the decoder, i.e. we want to find a mapping from relevant latent space volumes to phase space volumes. To this end, we perform a principal component analysis of the latent space representation of physical events. The PCA is an orthogonal transformation of the data that defines new axes such that the first component accounts for most of the variance in the dataset. We look at the first two principal components, sample a grid in these components and apply the inverse PCA transformation to get back to a latent space representation. We choose 64 points in latent space that were picked after finding that physical events in PCA space are distributed on an approximately circular area. Because of that finding we created an equidistant 8×88\times 8 grid in polar coordinates rr and ϕ\phi. The grid in PCA space is then transformed back to a latent space representation and used as input for the decoder to generate events that are being displayed in Fig. 9. The 64 chosen points on a polar grid correspond to the events in Fig. 9. This is effectively a two-dimensional PCA map of latent space. Observing the event displays reveals that we are in fact able to capture where we find events with what number of jets and leptons, what order of MET and what kind of orientations. In case one wants to produce events that e.g. look like event 62, one can do this by sampling around r=3.5r=3.5 and ϕ=225°\phi=225\degree in PCA space, then transform these events back to a latent space representation and to use that as input for the decoder. This will offer the possibility to narrow down the characteristics of the events even further and many iterations of this procedure will finally allow the generation of events with arbitrarily precise characteristics. Alternatively, one could create a classifier that defines boundaries of a latent space volume and corresponds to the desired phase space volume.

Having found that the B-VAE can be used to sample highly complex probability distributions, one possible application may be to provide a very efficient method for the phase space integration of multi-leg matrix elements. Recent work has shown that machine learning approaches to Monte Carlo integration of multidimensional probability distributions Bendavid 2017 and phase space integration of matrix elements Klimek and Perelstein 2018 may be able to obtain much better rejection efficiency than the current widely used methods Lepage 1980. We point out that event weights can be obtained from the B-VAE in similar fashion to the above papers.

The reconstruction of noise and test events in Fig. 8 clearly shows that ttˉt\bar{t} events beyond the training data are a) embedded well in latent space and b) reconstructed very well when compared to the reconstruction of noise. This suggests that one can use the (relative) reconstruction loss histograms or the (relative) reconstruction losses to detect anomalies, i.e. departures from the training data in terms of single events or their frequency of occurrence. The obvious use case of this is to train a B-VAE on a mixture of standard model events to detect anomalies in experimental data that correspond to new physics similarly to Nachman and Shih 2020. The B-VAE makes it possible to increase the ability to reconstruct the training and test data compared to a normal VAE, so it may be a better anomaly detector.

VI Discussion

We have provided more evidence for the capability of deep generative models to learn physical processes. To compare the performance of all of the investigated models, we have introduced two figures of merit, δ\delta and δOF\delta_{\text{OF}}. In particular, we describe and optimize a method for this task: the B-VAE. Several GAN architectures with recommended hyperparameters and the VAE with a standard normal prior fail to correctly produce events with the right frequency of occurrence. By creating a density information buffer with encoded ground truth data we presented a way to generate events whose probabilistic characteristics are in very good agreement with those found in the ground truth data. We identified the relevant hyperparameters of the B-VAE that allow for the optimization of its generalization properties and performed a first exploration of that hyperparameter space. We find that the dimensionality of the latent space should be smaller than, but close to, the input dimension. We find that it is necessary to heavily weight the reconstruction loss to create an accurate generative model and to tune the underestimated variance of the latent code. We have tested several traditional density estimation methods to learn the latent code density of encoded ground truth data, and concluded that the explicit use of the density information buffer with the parameters α\alpha and γ\gamma performs better. In a final step, we have investigated several GAN architectures with default hyperparameters but failed to create a model that successfully generates physical events with the right densities. Improvements could be made by performing a stricter model selection and to sweep through the full hyperparameter space beyond the hyperparameter recommendations given in the corresponding GAN papers. More generally, the GAN training procedure may be improved because the simultaneous gradient ascent that is currently used to find local Nash equilibria of the two-player game has issues that may be overcome by other objectives like the consensus optimization Mescheder et al. 2017 or by approaches such as the generative adversarial density estimator Abbasnejad et al. 2019.

By performing a principal component analysis of the latent space representations of MC events and a subsequent exploration of the corresponding PCA space we introduced an option to steer the event generation. In section IV we demonstrate that the statistics generalize to some degree. In future work it will be necessary to identify to which degree the implicit interpolation of the density pθ(x)p_{\theta}(\mathbf{x}) generalizes beyond the observed ground truth - and to maximize it. Another missing piece to complete the puzzle, is to find which generative models can describe processes that contain both events with very high and low multiplicities with up to twenty or more final state objects. Independent of what the outcome will be, potential improvements to all presented techniques can be made by incorporating auxiliary features as in Musella and Pandolfi 2018; Hashemi et al. 2019. Furthermore, improvements can be made by adding regression terms to the loss function that penalize deviations from the desired distributions in the generated data as in Hashemi et al. 2019 and by utilizing classification terms that force the number of objects and the corresponding object type to be correct. Another promising class of methods to create generative models are flow based models Rezende and Mohamed 2015; Kingma et al. 2016; Germain et al. 2015 and a thorough comparison of all available methods would be useful.

All in all, the results of this investigation indicate usefulness of the hereby proposed method not only for particle physics but for all branches of science that involve computationally expensive Monte Carlo simulations, that have the interest to create a generative model from experimental data, or that have the need to sample from high-dimensional and complex distributions.

VII Acknowledgements

This work was partly funded by and carried out in the SURF Open Innovation Lab project ”Machine learning enhanced high performance computing applications and computations” and was partly performed using the Dutch national e-infrastructure. S. C. and S. O. thank the support by the Netherlands eScience Center under the project iDark: The intelligent Dark Matter Survey. M. v. B. and R. V. acknowledge support by the Foundation for Fundamental Research of Matter (FOM), program 156, ”Higgs as Probe and Portal”. R. RdA, thanks the support from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 674896, the “SOM Sabor y origen de la Materia” MEC projects and the Spanish MINECO Centro de Excelencia Severo Ochoa del IFIC program under grant SEV-2014-0398.

VIII Author contributions

S. O. contributed to the idea, wrote most of the paper, invented the B-VAE and performed the majority of trainings and the analysis, as well as the data creation for all figures except figure 6 and the creation of figures 4, 5 and 7. S. C. contributed to the idea, invented the B-VAE and discussed every section and, also intermediate, results in great detail. W. d. S. contributed to the data generation of the toy model and the initial GAN for the toy model. M. v. B. contributed the ttˉt\bar{t} data and figures 1, 2 and 9, as well as the discussion and editing of the initial preprint. L. H. contributed to the idea, discussion and editing with his machine learning expertise, figure 6, as well as the training of several GANs. C. v. L. and D. P. contributed the HPC processing of the data and to discussions. R. R. d. A. contributed the Z-decay data and to discussions of the manuscript. R. V. contributed discussions and implementations of the figures of merit as well as discussions and editing of the manuscript, in particular section III.3 and figures 3 and 8.

IX Competing Interests

The authors declare no competing interests.

X Data Availability Statement

The ttˉt\bar{t} dataset that was used to obtain the results in sections IV.4 and IV.3 is available under the doi:10.5281/zenodo.3560661, https://zenodo.org/record/3560661#.XeaiVehKiUk. All other data that were used to train or that were generated with one of the trained generative models are available from the corresponding author upon request.

XI Code Availability Statement

The custom code that was created during the work that led to the main results of this article is published in a public GitHub repository: https://github.com/SydneyOtten/DeepEvents.

References