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

TL;DR

This paper explores deep generative models for physics event simulation, demonstrating that a density buffer improves the accuracy and speed of event generation, with applications in density estimation, anomaly detection, and importance sampling.

Contribution

It introduces a density information buffer for VAEs that significantly enhances event generation accuracy and efficiency in physics simulations.

Findings

Standard GANs and VAEs struggle to learn distributions precisely.

Buffering density information enables accurate, fast event sampling.

Applications include density estimation, anomaly detection, and phase space integration.

Abstract

We present a study for the generation of events from a physical process with deep generative models. The simulation of physical processes requires not only the production of physical events, but also to ensure these events occur with the correct frequencies. We investigate the feasibility of learning the event generation and the frequency of occurrence with Generative Adversarial Networks (GANs) and Variational Autoencoders (VAEs) to produce events like Monte Carlo generators. We study three processes: a simple two-body decay, the processes $e^+e^-\to Z \to l^+l^-$ and $p p \to t\bar{t} $ including the decay of the top quarks and a simulation of the detector response. We find that the tested GAN architectures and the standard VAE are not able to learn the distributions precisely. By buffering density information of encoded Monte Carlo events given the encoder of a VAE we are able to construct a prior for the sampling of new events from the decoder that yields distributions that are in very good agreement with real Monte Carlo events and are generated several orders of magnitude faster. Applications of this work include generic density estimation and sampling, targeted event generation via a principal component analysis of encoded ground truth data, anomaly detection and more efficient importance sampling, e.g. for the phase space integration of matrix elements in quantum field theories.

Tables

Table 1. Table 1 : The combinations of dim 𝐳 , B , α = 1 dimension 𝐳 𝐵 𝛼 1 \dim\mathbf{z},B,\alpha=1 and γ 𝛾 \gamma giving the top-15 performance w.r.t. δ 𝛿 \delta with the corresponding δ OF subscript 𝛿 OF \delta_{\text{OF}} .

$(dim𝐳,B,\alpha ,\gamma )$	$\delta$	${\delta }_{\text{OF}}$
$(20,{10}^{-6},1,0.01)$	0.0076	3.69
$(20,{10}^{-7},1,0.01)$	0.0090	3.81
$(20,{10}^{-6},1,0.05)$	0.0090	1.01
$(16,{10}^{-7},1,0.01)$	0.0095	4.29
$(16,{10}^{-6},1,0.01)$	0.0101	3.30
$(16,{10}^{-6},1,0.05)$	0.0122	0.51
$(16,{10}^{-5},1,0.05)$	0.0137	0.46
$(24,{10}^{-5},1,0.05)$	0.0138	0.65
$(16,{10}^{-5},1,0.01)$	0.0148	1.23
$(20,{10}^{-5},1,0.01)$	0.0148	1.18
$(24,{10}^{-5},1,0.01)$	0.0149	1.06
$(28,{10}^{-7},1,0.01)$	0.0155	4.22
$(24,{10}^{-7},1,0.01)$	0.0156	3.63
$(24,{10}^{-6},1,0.01)$	0.0165	3.68
$(28,{10}^{-6},1,0.01)$	0.0176	3.71

Table 2. Table 2 : The combinations of dim 𝐳 , B , α dimension 𝐳 𝐵 𝛼 \dim\mathbf{z},B,\alpha and γ = 0 𝛾 0 \gamma=0 giving the top-15 performance w.r.t. δ 𝛿 \delta with the corresponding δ OF subscript 𝛿 OF \delta_{\text{OF}} .

$(dim𝐳,B,\alpha ,\gamma )$	$\delta$	${\delta }_{\text{OF}}$
$(28,{10}^{-7},1,0)$	0.0066	94.00
$(24,{10}^{-7},1,0)$	0.0074	97.92
$(20,{10}^{-6},1,0)$	0.0075	17.48
$(20,{10}^{-7},1,0)$	0.0084	106.64
$(20,{10}^{-7},5,0)$	0.0088	5.67
$(16,{10}^{-7},5,0)$	0.0093	7.96
$(16,{10}^{-7},1,0)$	0.0094	133.48
$(16,{10}^{-6},1,0)$	0.0102	14.02
$(16,{10}^{-7},10,0)$	0.0102	2.05
$(20,{10}^{-7},10,0)$	0.0112	1.65
$(24,{10}^{-7},5,0)$	0.0144	4.55
$(24,{10}^{-6},1,0)$	0.0156	15.45
$(28,{10}^{-6},1,0)$	0.0162	13.52
$(28,{10}^{-7},5,0)$	0.0166	4.60
$(28,{10}^{-6},5,0)$	0.0189	0.68

Table 3. Table 3 : Performance of the B-VAE with different latent space density estimation techniques.

Model	$\delta$	${\delta }_{\text{OF}}$
KDE	0.2631	7.30
GMM, 50	0.0297	11.54
GMM, 100	0.0266	12.70
GMM, 1000	0.0302	7.64
S-VAE, B=1	1.4452	840.60
S-VAE, B=0.1	1.1438	756.56
S-VAE, B=0.01	0.0689	12.73
S-VAE, B=${10}^{-3}$	0.0438	3.36
S-VAE, B=${10}^{-4}$	0.0992	9.26
S-VAE, B=${10}^{-5}$	0.3179	12.33
S-VAE, B=${10}^{-6}$	0.8753	238.33

Table 4. Table 4 : δ 𝛿 \delta and δ OF subscript 𝛿 OF \delta_{\text{OF}} for several GAN architectures.

GAN model	$\delta$	${\delta }_{\text{OF}}$
DijetGAN	0.3477	19.70
LSGAN	0.3592	3.03
MMD-GAN	0.9454	642.20
WGAN-GP	1.1605	723.79
WGAN	1.0672	840.63

Table 5. Table 5 : KDE and GMM model performance evaluated on figures of merit δ 𝛿 \delta and δ OF subscript 𝛿 OF \delta_{\text{OF}} .

Model	$\delta$	${\delta }_{\text{OF}}$
KDE	0.6038	4.99
GMM, 50	0.1078	16.33
GMM, 100	0.0948	20.43
GMM, 1000	0.0874	12.09
5 % Smearing	0.0093	3.53
10 % Smearing	0.0192	3.89

Equations

$$\theta =\arccos \frac{p_z}{\sqrt{p_x^2+p_y^2+p_z^2}},\varphi =\arctan \frac{p_y}{p_x}.$$

$$E^2-{\mathbf{p}}^2-m^2=0.$$

$$p_{\theta }(\mathbf{x})=\int p_{\theta }(\mathbf{z})p_{\theta }(\mathbf{x}|\mathbf{z})d\mathbf{z}$$

$$\log p_{\theta }({\mathbf{x}}^1,\dots ,{\mathbf{x}}^N)=\sum _{i=1}^N\log p_{\theta }({\mathbf{x}}^i)$$

$$\log p({\mathbf{x}}^i)=\log E_{p(\mathbf{z}|{\mathbf{x}}^i)}\left[\frac{p({\mathbf{x}}^i,\mathbf{z})}{p(\mathbf{z}|{\mathbf{x}}^i)}\right]$$

$$\log p({\mathbf{x}}^i)\ge E_{p(\mathbf{z}|{\mathbf{x}}^i)}\left[\log \frac{p({\mathbf{x}}^i,\mathbf{z})}{p(\mathbf{z}|{\mathbf{x}}^i)}\right].$$

$$\log p_{\theta }({\mathbf{x}}^i)=D_{\text{KL}}(q_{\varphi }(\mathbf{z}|{\mathbf{x}}^i)\Vert p_{\theta }(\mathbf{z}|{\mathbf{x}}^i))+\mathcal{L}(\theta ,\varphi ;{\mathbf{x}}^i).$$

$$\begin{array}{rl}\mathcal{L}(\theta ,\varphi ;\mathbf{x})& =E_{q_{\varphi }(\mathbf{z}|\mathbf{x})}\left[-\log q_{\varphi }(\mathbf{z}|\mathbf{x})+\log p_{\theta }(\mathbf{x},\mathbf{z})\right]\\ & =-D_{\text{KL}}(q_{\varphi }(\mathbf{z}|\mathbf{x})\Vert p_{\theta }(\mathbf{z}))\\ & +E_{q_{\varphi }(\mathbf{z}|\mathbf{x})}\left[\log (p_{\theta }(\mathbf{x}|\mathbf{z})\right].\end{array}$$

$$\log q_{\varphi }(\mathbf{z}|\mathbf{x})=\log \mathcal{N}(\mathbf{z};\mu ,{\sigma }^{2}1)$$

$$\begin{array}{rl}-D_{\text{KL}}(& q_{\varphi }(\mathbf{z}|\mathbf{x})\Vert p_{\theta }(\mathbf{z}))\\ & =\frac{1}{2}\sum _{j=1}^{\dim \mathbf{z}}1+\log ({\sigma }_j^2)-{\mu }_j^2-{\sigma }_j^2.\end{array}$$

$$\mathcal{L}(\theta ,\varphi ;\mathbf{X})\simeq {\mathcal{L}}^M(\theta ,\varphi ;{\mathbf{X}}^M)=\frac{N}{M}\sum _{i=1}^M\tilde{\mathcal{L}}(\theta ,\varphi ;{\mathbf{x}}^i),$$

$$\begin{array}{rl}\tilde{\mathcal{L}}(\theta ,\varphi ;{\mathbf{x}}^i)=& -D_{\text{KL}}(q_{\varphi }(\mathbf{z}|{\mathbf{x}}^i)\Vert p_{\theta }(\mathbf{z}))\\ & +\log (p_{\theta }({\mathbf{x}}^i|g_{\varphi }({ϵ}^i,{\mathbf{x}}^i)))\end{array}.$$

$$L=\frac{1}{M}\sum _{i=1}^M(1-B)\cdot \text{MSE}+B\cdot D_{\text{KL}}.$$

$$\dim \mathbf{z}=\{16,20,24,28\},$$

$$B=\{10^{-7},10^{-6},10^{-5},10^{-4}\},$$

$$\hat{p}(y)=\sum _{i=1}^{N}K\left(\frac{y-x_i}{h}\right)$$

$$K(x;h)\propto \exp \left(-\frac{x^2}{2h^2}\right),$$

$$p_{\varphi ,{\mathbf{X}}_{\text{L}}}(\mathbf{z})=\sum _{i=1}^m{q}_{\varphi }(\mathbf{z}|{\mathbf{x}}^i)p({\mathbf{x}}^i)\text{with}p({\mathbf{x}}^i)=\frac{1}{m}.$$

$$\mathbf{z}^{i}\sim\left\{\begin{array}[]{ll}\mathcal{N}(\mu^{i},\alpha\sigma^{2,i})&\text{if }\sigma<\sigma_{T}\\ \mathcal{N}(\mu^{i},\sigma^{2,i})&\text{else}\end{array}\right.,$$

$${\left(z^i\right)}_{j=1}^{\dim z}\sim {\left(\mathcal{N}\left({\mu }_j^i,{\alpha }_j{\sigma }_j^{2,i}+{\gamma }_j\right)\right)}_{j=1}^{\dim \mathbf{z}},$$

$${\left(z^i\right)}_{j=1}^{\dim \mathbf{z}}\sim {\left(\mathcal{N}\left({\mu }_j^i,{\sigma }_j^{2,i}+\gamma \right)\right)}_{j=1}^{\dim \mathbf{z}}.$$

$$\Delta R=\sqrt{(\Delta \varphi {)}^2+(\Delta \eta {)}^2}$$

$$\begin{array}{rl}\delta & =\frac{1}{N_{\text{hist}}}\sum _i\sum _{j<i}{\chi }_{ij}^2\end{array}$$

$${\chi }^2=\sum _{u=1}^{N_{\text{bins}}}\frac{{\left(p_u-p_u^{\text{MC}}\right)}^2}{p_u+p_u^{\text{MC}}}$$

$${\delta }_{\text{OF}}={\int }_0^{1000}d\sqrt{N_{\text{bins}}}\left|f_{\text{e}}\left(\sqrt{N_{\text{bins}}}\right)-f_{\text{e}}^{\text{MC}}\left(\sqrt{N_{\text{bins}}}\right)\right|.$$

$$\delta (N,A,B,C,D)=\frac{A}{\log (BN+C)}+D$$

$$\gamma$$

$$\alpha$$

$$p_{\theta }({\tilde{\mathbf{x}}}^i|\mathbf{z}({\mathbf{x}}^i))=\delta (\mathbf{x}-{\mathbf{x}}^i).$$

Full text

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 $10$ 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 $\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 $\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(\mathbf{x})$. The main challenge we tackle is to create a model that learns a transformation from a random variable $\mathbf{z}\to\mathbf{x}$ such that the distribution of $\mathbf{x}$ follows $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 $Z$ boson production from $e^{+}e^{-}$ collisions and its decay to two leptons, $e^{+}e^{-}$ and $\mu^{+}\mu^{-}$, with four 4-vectors per data point of which two are always zero.

• •

the 26-dimensional $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 $Z$ bosons reveals that standard variational autoencoders can’t reliably model the process but confirms good agreement for the B-VAE. For $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 $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

• •

provide a brief overview of GANs,

• •

explain VAEs and our method, the B-VAE,

• •

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 $t\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 $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. $Z$-boson production in $e^{+}e^{-}$ collisions and 3. top quark production and decay in proton collisions, i.e. $pp\to t\bar{t}$. Here, we describe the procedures for obtaining the data sets.

II.1 10-dimensional toy model

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

[TABLE]

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 $E_{1}$, $E_{2}$ of the daughter particles, the phase space components $p_{x}$, $p_{y}$, $p_{z}$ for each particle and their masses $m_{1}$ and $m_{2}$. This introduces a degeneracy with the goal of checking whether the generative models learn the relativistic dispersion relations

[TABLE]

II.2 16-dimensional $e^{+}e^{-}\to Z\to l^{+}l^{-}$

We generate $10^{6}$ events of the $e^{+}e^{-}\to Z\to l^{+}l^{-}$ ($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: $\left(E_{e^{-}},p_{x,e^{-}},p_{y,e^{-}},p_{z,e^{-}},E_{e^{+}},p_{x,e^{+}},p_{y,e^{+}},p_{z,e^{+}},\right.$ $\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^{-}$ or $\mu^{+}\mu^{-}$.

II.3 26-dimensional $pp\to t\bar{t}$

We generate $1.2\cdot 10^{6}$ events of $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,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.

III.1 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).

Generative Adversarial Networks

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 $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 $\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 $\mathbf{z}$. With this in mind one may think of event generation as a two-step process: (1) sampling from a parameterized prior $p_{\mathbf{\theta}}(\mathbf{z})$ (2) sampling $\mathbf{x}^{i}$ from the conditional distribution $p_{\mathbf{\theta}}(\mathbf{x}^{i}|\mathbf{z})$, representing the likelihood. For deep neural networks the marginal likelihood

[TABLE]

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_{\mathbf{\theta}}(\mathbf{z}|\mathbf{x})$ is created: a probabilistic encoder $q_{\mathbf{\phi}}(\mathbf{z}|\mathbf{x})$. Given a data point $x^{i}$ it will produce a distribution over the latent code $\mathbf{z}$ from which the data point might have been generated. Similarly, a probabilistic decoder $p_{\mathbf{\theta}}(\mathbf{x}|\mathbf{z})$ is introduced that produces a distribution over possible events $\mathbf{x}^{i}$ given some latent code $\mathbf{z}$. In this approach the encoder and decoder are deep neural networks whose parameters $\mathbf{\phi}$ and $\mathbf{\theta}$ are learned jointly.

The marginal likelihood

[TABLE]

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

[TABLE]

and applying Jensen’s inequality, one finds that

[TABLE]

For this situation, one must substitute $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_{\phi}(\mathbf{z}|\mathbf{x}^{i})$. From here one can derive the variational lower bound $\mathcal{L}(\theta,\phi;\mathbf{x}^{i})$ and find that

[TABLE]

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

[TABLE]

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

Variational Autoencoder

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_{\phi}(\mathbf{z}|\mathbf{x})$ one reparametrizes $\hat{\mathbf{z}}\sim q_{\phi}(\mathbf{z}|\mathbf{x})$ using a differentiable transformation $g_{\phi}(\epsilon,\mathbf{x}),\epsilon\sim\mathcal{N}(0,1)$ with an auxiliary noise variable $\epsilon$. Choosing $\mathbf{z}\sim p(\mathbf{z}|\mathbf{x})=\mathcal{N}(\mu,\sigma^{2})$ with a diagonal covariance structure, such that

[TABLE]

where $\mu$ and $\sigma^{2}$ are outputs of the encoding deep neural network. Reparametrizing $\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:

[TABLE]

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 $\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>100$ it is sufficient to sample $\epsilon$ once which is adopted in our implementation. By calculating the lower bound for a batch of $M$ samples $\mathbf{X}^{M}\subset\mathbf{X}$ they construct the estimator of $\mathcal{L}$:

[TABLE]

where

[TABLE]

We use the gradients $\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 $D_{\text{KL}}$ and the MSE such that the loss function of the VAE $L\propto D_{\text{KL}}+MSE$. In our approach we introduce a multiplicative factor $B$ for $D_{\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 $\beta>1$ such that each dimension is more closely related to features of the output. In contrast we choose $B\ll 1$ to emphasize a good reconstruction. The loss function of the VAE can subsequently be written as

[TABLE]

This however also implies that $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=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\to l^{+}l^{-}$) or six hidden layers ($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\cdot 10^{-6}$ for the toy model and $Z\to l^{+}l^{-}$. The number of latent space dimensions are 9 for the toy model and 10 for $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\to l^{+}l^{-}$ we employ a learning rate scheduling for $7\times 80$ epochs and SWATS Shirish Keskar and Socher (2017), i.e. switching from Adam to SGD during training.

For $pp\to t\bar{t}$ we perform a scan over hyperparameters with

[TABLE]

We perform this scan on a small training data set with $10^{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}$ 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\to t\bar{t}$ the setup is identical except for the number of epochs: we train $4\times 240$ epochs with Adam and then for $4\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(\mathbf{z})=\mathcal{N}(0,1)$ used to sample $p_{\theta}(\mathbf{x}|\mathbf{z})$ isn’t identical to the distribution over the latent code $\mathbf{z}$ resulting from the encoding of true observations $q_{\phi}(\mathbf{z}|\mathbf{X})$. The generated distribution over $\mathbf{x}$ given $p_{\theta}(\mathbf{x}|\mathbf{z})$ therefore doesn’t match the reference when assuming a unit Gaussian over $\mathbf{z}$. We address this issue by estimating the prior $p(\mathbf{z})$ for the probabilistic decoder from data using a strategy similar to the Empirical Bayes method Robbins (1956). We collect observations $\mathbf{Z}=\{\mathbf{z}_{1},\ldots,\mathbf{z}_{m}\}$ by sampling $q_{\phi}(\mathbf{z}|\mathbf{X}_{\text{L}})$ where $\mathbf{X}_{\text{L}}\subset\mathbf{X}$ is a subset of physical events. $\mathbf{Z}$ is then used as the data for another density estimation to create a generative model for $p(\mathbf{z})$: this is what we call the density information buffer. This buffer is used in several ways to construct $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.

Kernel Density Estimation

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

[TABLE]

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

[TABLE]

and have optimised $h$. 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 $h$ using GridSearchCV and a 5-fold cross validation for 20 samples for $h$ distributed uniformly on a log-scale between 0.1 and 10.

Gaussian Mixture Models

Since the VAE also minimizes $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_{\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\}$ with the full covariance matrix and $10^{5}$ encodings $\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_{\eta}(\mathbf{z}|\mathbf{z}^{\prime})$ from latent code observations $\mathbf{Z}$ that is sampled using $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},\ldots,1\}$. We slightly modify the loss function from eq. (13) and remove the $(1-B)$ in front of the MSE term because we want to test higher values of $B$ of up to $B=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\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 $\mathbf{z}^{\prime}$ to then sample $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(\mathbf{z})$ and $q_{\phi}(\mathbf{z})$ is to explicitly construct a prior $p_{\phi,\mathbf{X}_{\text{L}}}(\mathbf{z})$ by aggregating (a subset of) the encodings of the training data:

[TABLE]

Practically this is done by saving all $\mu^{i}$ and $\sigma^{2,i}$ for all $m$ events in $\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 $\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 $B$ 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 $\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 $\alpha=\{1,5,10\}$ and only apply $\alpha$ if $\sigma<\sigma_{T}=0.05$, such that

[TABLE]

with $\mu^{i}$ and $\sigma^{2,i}$ being the Gaussian parameters in latent space corresponding to events $i=1,\ldots,m$ in $\mathbf{X}_{\text{L}}$. It is straightforward to expand this approach to have more freedom in $\alpha$, e.g. by optimizing $(\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 $\gamma_{j}$ for the standard deviation such that

[TABLE]

we have $2\cdot\dim\mathbf{z}$ additional hyperparameters. More generally we can try to learn a vector-valued function $\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 $\gamma=\gamma_{j}=\{0.01,0.05,0.1\}$ to be the minimal standard deviation, such that

[TABLE]

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 $\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\cdot 10^{6}$ Monte Carlo samples, $10^{5}$ of those samples are used as training data for the generative models. We then produce sets of $1.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 $E$, $p_{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

[TABLE]

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 $N_{\text{2D,bins}}=25^{2}$. We then define

[TABLE]

with $i$ and $j$ summing over observables and $N_{\text{hist}}=528$ which represent averages over all histograms of the test statistic

[TABLE]

where $p_{u}$ and $p_{u}^{\text{MC}}$ are the normalized bin contents of bins $u$ 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 $\eta_{j_{1}}$ vs. $\eta_{j_{2}}$, $\eta_{i}\in[-2.5,2.5]$ with a variable number of bins $N_{\text{bins}}$ and measure the fraction of empty bins $f_{\text{e}}\left(\sqrt{N_{\text{bins}}}\right)$. As $N_{\text{bins}}$ is increased, the histogram granularity probes the rate of overfitting in increasingly smaller regions of the phase space. The function $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 $f^{\text{MC}}_{\text{e}}\left(\sqrt{N_{\text{bins}}}\right)$. We therefore compute the figure of merit as

[TABLE]

After searching the best model with respect to $\eta_{j_{1}}$ vs. $\eta_{j_{2}}$, we independently check its behavior in $\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 $\delta_{\text{OF}}$ for $1.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\cdot 10^{5}$ events and then extrapolate as shown in Fig. 5(c). We expect that an ideal generator has $\delta\approx 0.0003$ and $\delta_{\text{OF}}\approx 0.5$ for $1.2\cdot 10^{6}$ events. For $\delta$, this value is obtained by fitting the parameters $A,B,C,D$ to the empirically motivated function

[TABLE]

for values of $\delta$ obtained by evaluating the MC dataset with $N=\{50000,100000,\ldots,600000\}$ and extrapolating to $N=1200000$ using the $curve\_fit$ function in scipy. The expected $\delta_{\text{OF}}$ is assumed to be roughly flat around 0.5. We select our best model by requiring $\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 $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\mathbf{z}$ and $B$ 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 $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.

IV.1 Two-body decay toy model

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 $p_{x}$, $p_{y}$ and $p_{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 $\phi_{j_{1}}\sim U(-\pi,\pi)$. It is questionable whether this technique can be generalized to higher multiplicities.

IV.2 $e^{+}e^{-}\rightarrow Z\rightarrow l^{+}l^{-}$

Fig. 1(b) and 2 show the results for the $Z$ events, where the $Z$ 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 $p_{T}$, $\theta$ and $\phi$ and good agreement for invariant mass $M_{\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 $p_{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 $pp\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 $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 $b$-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 ($p_{T}$) of 20 GeV in the Monte Carlo simulation. We first explore the hyperparameter space of the B-VAE in $\dim\mathbf{z}$, $B$, $\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\mathbf{z},B,\alpha,\gamma)$ combinations evaluated on the figures of merit defined in section III.3. For all possible combinations of $\dim\mathbf{z}$ and $B$ as defined in section III.1 we have separately investigated

[TABLE]

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

It is not surprising that the best performance in $\delta_{\text{1D}}$ is attained by the B-VAE with the highest latent code dimensionality, the lowest $B$ and $\alpha=1,\gamma=0$. The downside however is a very poor performance in $\delta_{\text{OF}}$. Comparing to the values for the $5\%$ Gaussian smearing of events in Table 5, they are very similar in $\delta$ but even worse in $\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 $\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 $\delta_{\text{OF}}$, whereas models with $B>10^{-5}$ perform poorly in $\delta$. For illustrative purposes we proceed to show and discuss details for the model we consider best: $\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,p_{T},\eta$ and $\phi$ for all four jets and the leading lepton,

• •

MET and MET$\phi$,

• •

$\Delta\phi$ between MET and leading lepton,

• •

$\Delta$R between leading and subleading jets and

• •

the invariant mass $M_{\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 $10^{5}$ samples that were used to generate $1.2\cdot 10^{6}$ events which are compared to $1.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 $\Delta$R we see almost perfect agreement.

Generating $10^{7}$ $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 ($\phi_{1}$ vs $\phi_{2}$) that were created using the B-VAE with $\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 $[2,3]\times[2,3]$. The first row displays the training data consisting of $10^{5}$ events. The second and third row of Fig. 4 show $1.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 $10^{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 $\alpha=1$ and $\gamma=0$ and

• •

generator artefacts appear around $(\pm\pi,0)$ and $(0,\pm\pi)$. For $E$, $p_{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 $\eta_{j_{1}}$ vs. $\eta_{j_{2}}$ and $\phi_{j_{1}}$ vs. $\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)$. 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 $B$. This is an expected result since the encoder $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(\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 $\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 $\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 $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 $10^{4}$ samples and then perform a KDE with $h_{\text{opt}}$ on $10^{5}$ samples. Additionally, we create a GMM of those $10^{5}$ samples with $50,100$ and $1000$ components with a maximum of 500 iterations. Subsequently we generate $1.2\cdot 10^{6}$ samples from the KDE and the three GMM models and evaluate them with our figures of merit $\delta$ and $\delta_{\text{OF}}$ as presented in Table 5. Additionally we take $10^{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 $\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 $\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 $\delta_{\text{OF}}$.

When $\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

[TABLE]

Since $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 $\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 $t\bar{t}$ events of a VAE with a 20-dimensional latent space and $B=10^{-6}$ and the reconstructions of the same VAE for $\mathbf{x}_{\text{test}}\sim U(0,1)$, where we clearly see that the VAE does not simply learn the identity function. The parameters $\alpha,B,\gamma$ and $\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\times 8$ grid in polar coordinates $r$ 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.5$ and $\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 $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 $\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_{\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 $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 $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.