1 Introduction
Machine learning models based on deep neural networks have revolutionized information processing over the last decade. Such models can recognize objects in images krizhevsky2012imagenet ; ResNet ; DenseNet , perform language translation NMT ; GoogleNMT , transcribe spoken language Transcription , and even speak written text Wavenet
at approaching human level. The truly revolutionary aspect of this progress is the generality of deep neural networks: a broad diversity of network architectures can be created from basic building blocks that allow for efficient calculation of gradients via back propagation, and thus efficient optimization through stochastic gradient descent
BackProp. These methods are arbitrarily expressive and can model extremely high dimensional data.
The architecture of a neural network should be designed to process information efficiently, from the input data all the way through to the network’s final output. Indeed, it empirically seems to be the case that networks that process information evenly layerbylayer perform very well. One example of this empirical result is that deep convolutional networks for image processing seem to perform sequentially more abstract operations as a function of depth krizhevsky2012imagenet . Similarly, recurrent networks perform well on time series data, as their recurrent layers naturally describe stepbystep evolution in time mikolov2010recurrent .
The power and generality of deep neural networks has been leveraged across the sciences, and in particular in particle physics. The simplest architecture explored has been the fullyconnected network, which has successfully been applied in a wide variety of contexts, such as in identifying and splitting clusters from multiple particles in the pixel detector Aad:2014yva , in tagging Aad:2015ydr , and in identification Chatrchyan:2012zz . In these basic applications, the neural network optimizes its use of some finite number of relevant physical observables for the task at hand.^{1}^{1}1 For recent work on constructing a basis for neural network inputs, see Datta:2017rhs ; Datta:2017lxt ; Luo:2017ncs , and see Komiske:2017aww for a linear approach that does not require neural network methods.
One drawback of such an approach is that the neural network is limited by the observables it is given. In fact, for these applications, other multivariate methods such as boosted decision trees often have comparable performance using the same inputs, but train faster and can be less sensitive to noise
Gallicchio:2010dq ; ATLPHYSPUB2017004 .As an alternative to feeding a neural network a set of motivated observables, one can feed it raw information. By doing so, one allows the network to take advantage of useful features that physicists have yet to discover. One way of preprocessing the raw data in a fairly unbiased way is through the use of jet images, which contain as pixel intensities the energy deposited by jet constituents in calorimeter cells Cogan:2014oua . Jet images invite the use of techniques from image recognition to discriminate jets of different origins. In Cogan:2014oua
, the pixel intensities in the twodimensional jet image were combined into a vector, and a Fisher linear discriminant was then used to find a plane in the highdimensional space that maximally separates two different jet classes. Treating a 2dimensional jet image as an unstructured collection of pixel intensities, however, ignores the spatial locality of the problem, i.e. that neighboring pixels should have related intensities. Convolutional neural networks (CNNs), which boast reduced complexity by leveraging this spatially local structure, have since been adopted instead, and they generally outperform fullyconnected networks due to their efficient feature detection. In the first applications of CNNs to jet images, on boosted
detection deOliveira:2015xxd and quark/gluon discrimination Komiske:2016rsd , it was indeed found that simple CNNs could generally outperform previous techniques. Since then, a number of studies have aimed to optimize various discrimination tasks using CNNs Komiske:2017ubm ; Kasieczka:2017nvn ; Bhimji:2017qvb ; ATLPHYSPUB2017017 ; Macaluso:2018tck ; Chien:2018dfn .While the twodimensional detector image acts as a natural representation of a jet, especially from an experimental standpoint, the 4momenta of individual jet constituents provide a more fundamental representation for the input to a neural network. One complication in transitioning from the jet image to its list of momenta is that, while the image is a fixedsize representation, the list of momenta will have different sizes for different jets. To avoid this problem, one could truncate the list of momenta in the jet to a fixed size, and zeropad jets smaller than this size
Pearkes:2017hku. Alternatively, there are network architectures, namely recursive (RecNNs) and recurrent neural networks (RNNs), that handle variable length inputs naturally. With such methods, one also has the freedom to choose the order in which constituent momenta are fed into the network. In
Louppe:2017ipp , a RecNN was used to build a fixedsize representation of the jet, and the authors explored various ways of ordering the momenta as input to the network: by jet clustering algorithms, by transverse momentum, and randomly. The resulting representation of the jet was then fed to a fullyconnected neural network for boosted tagging. RecNNs and RNNs have also been used in similar ways for quark/gluon discrimination Cheng:2017rdo , top tagging Egan:2017ojy , and jet charge Fraser:2018ieu . See also Guest:2016iqz ; ATLPHYSPUB2017003 for jet flavor classification using tracks.To date, the majority of applications of machine learning to particle physics employ supervised machine learning techniques. Supervised learning is the optimization of a model to map input to output based on labeled inputoutput pairs in the training data. These training examples are typically simulated by Monte Carlo generators, in which case the labels come from the underlying physical processes being generated. Most of the classification studies mentioned above employ this style of supervised learning, and similar techniques have also been utilized for regression tasks such as pileup subtraction
Komiske:2017ubm . Alternatively, training data can be organized in mixed samples, each containing different proportions of the different underlying processes. In this case, labels correspond to the mixed samples, and learning is referred to as weakly supervised. While full and weak supervision are very similar as computational techniques, the distinction is exceptionally important in particle physics, where the underlying physical processes are unobservable in real collider data. Early studies of weakly supervised learning in particle physics show very promising results: performance comparable to fully supervised methods was found both with lowdimensional inputs Metodiev:2017vrx ; Cohen:2017exh (a few physical observables) and with very highdimensional inputs Komiske:2018oaa (jet images).With supervised learning, there is a notion of absolute accuracy: since every training example is labeled with the desired output, the network predicts this output either correctly or incorrectly. This is in contrast to unsupervised learning, where the machine learns underlying structure that is unlabeled in the training data. Without outputlabeled training examples, there is no notion of absolute accuracy. Several recent studies have employed unsupervised learning techniques in particle physics. In Metodiev:2018ftz , borrowing concepts from topic modelling in text documents, the authors extract observable distributions of underlying quark and gluon jets from two mixed samples. In deOliveira:2017pjk ; Paganini:2017hrr ; Paganini:2017dwg , generative adversarial networks (GANs) are used to efficiently generate realistic jet images and calorimeter showers.
In this work, we explore another approach to unsupervised machine learning in particle physics, in which a deep neural network learns to compute the relative differential cross section of each data point under consideration, or equivalently, the probability distribution generating the data. The power of having access to the probability distribution underlying the data should not be underestimated. For example, likelihood ratios would provide optimal discriminants Neyman289 , and sampling from the probability distribution would provide completely datadriven simulations.
In this paper, we introduce a framework named Junipr: “Jets from UNsupervised Interpretable PRobabilistic models”. We also present a basic implementation of this framework using a deep neural network. This network directly computes the general probability distribution underlying particle collider data using unsupervised learning.
The task of learning the probability distribution underlying collider data comes with challenges due to the complexity of the data. Some past studies have aimed to process collider information efficiently by using neural network architectures inspired by physics techniques already in use Louppe:2017ipp ; Cheng:2017rdo ; Egan:2017ojy ; Fraser:2018ieu ; Guest:2016iqz ; Butter:2017cot . In this paper, we take this idea one step further. We scaffold the neural network architecture around a leadingorder description of the physics underlying the data, from first input all the way to final output. Specifically, we base the Junipr framework on algorithmic jet clustering trees. The tree structure is used, both in processing input information, and in decomposing the network’s output. In particular, Junipr’s output is organized into meaningful probabilities attached to individual nodes in a jet’s clustering tree. In addition to reducing the complexity and increasing the efficiency of the corresponding neural network, this approach also forces the machine to speak a language familiar to physicists, thus enabling its users to interpret the underlying physics it has learned. Indeed, one common downside associated with machine learning techniques in physics is that, though they provide powerful methods to accomplish the tasks learned in training, they do little to clarify the underlying physics that underpins their success. Our approach minimizes this downside.
Let us elaborate on the treebased architecture used for Junipr’s implementation. In particle physics, events at colliders are dominated by the production of collimated collections of particles known as jets. The origin of jets and many of their properties can be understood through the fundamental theory of strong interactions, quantum chromodynamics (QCD). One insight from QCD is that jets have an inherently fractal structure, inherited from the approximate scale invariance of the fundamental theory. The fractal structure is made precise through the notion of factorization, which states that the dynamics in QCD stratify according to soft, collinear, and hard physics Coleman:1965xm ; Collins:1985ue ; Collins:1988ig ; feige2014hard ; feige2013shell , with each sector being separately scale invariant. To capture this structure efficiently in Junipr, we use a kind of factorized architecture, with a dense network to describe local branchings (wellsuited for collinear factorization), and a global RNN superstructure general enough to encode soft coherence and any factorizationviolating effects.
One might naively expect this setup to require knowledge of the sequence of splittings that created the jet. Although there is a sequence of splittings in partonshower simulations, the splittings are only a semiclassical approximation used to model the intensely complex and essentially incalculable distribution of final state particles. Real data is not labelled with any such sequence. In fact, there are many possible sequences which could produce the same event, and the cross section for the event is given by the square of the quantum mechanical sum of all such amplitudes, including effects of virtual particles. A proxy for this fictitious splitting history is a clustering history that can be constructed in a deterministic way using a jetclustering algorithm, such as the algorithm Catani:1993hr ; Ellis:1993tq or the Cambridge/Aachen (C/A) algorithm Dokshitzer:1997in ; Wobisch:1998wt . There is no correct algorithm: each is just a different way to process the momenta in an event. Indeed, there seems to be useful information in the multiple different ways that the same event can be clustered Ellis:2012sn ; Kahawala:2013sba ; Mackey:2015hwa . Any of these algorithms, or any algorithm at all that encodes the momenta of an event into a binary tree, can be used to scaffold a neural network in the Junipr approach.
For practical purposes, Junipr is implemented with respect to a fixed jet clustering algorithm. Without a fixed algorithm, the probability of the finalstate particles constructed through branchings would require marginalization over all possible clustering histories — an extremely onerous computational task. In principle, fixing the algorithm used to implement Junipr should be inconsequential for its output, namely the probability distribution over finalstate momenta, as these momenta are independent of clustering algorithm. To reiterate, the Junipr approach does not require the chosen clustering algorithm to agree with the underlying datageneration process; this is demonstrated in Secs. 5.2 and 5.3 below. On the other hand, the sequence of probabilities assigned to each branching in a clustering tree certainly depends on the algorithm used to define the tree. For example, the same final probability could be reached with one clustering algorithm through the sequence , or with another algorithm through . The key idea is that, if an algorithm is chosen which does correspond to a semiclassical parton shower, the resulting sequence of probabilities may be understandable. This provides avenues for users to interpret what physics the machine learns, and we expect that dissecting Junipr will be useful in such cases. We will demonstrate this throughout the paper.
It is worth emphasizing one fundamental aspect of our approach for clarity. The Junipr framework yields a probabilistic model, not a generative model. The probabilistic model allows us to directly compute the probability density of an individual jet, as defined by its set of constituent particle momenta. To be precise, this is the probability density for those particular momenta to arise in an event, conditioned on the event selection criteria used to select the training data. As a complementary example of this, shower deconstruction Soper:2011cr ; Soper:2014rya provides a theorydriven approach to probabilistic modeling in particle physics, in which probabilities are calculated using QCD rather than a neural network. In contrast, a generative model would output an example jet, taking random noise as input to seed the generation process. Given a distribution of input seeds, the jets output from a generative model should follow the same distribution as the training data. While this means that the probability distribution underlying the data is internally encoded in a generative model, this underlying distribution is hidden from the user. Examples of generative models in particle physics include Monte Carlo event generators and, more recently, GANs used to generate jet images and detector simulations deOliveira:2017pjk ; Paganini:2017hrr ; Paganini:2017dwg .
The direct access to the probability distribution that is enabled by a probabilistic model comes with several advantages. If two different probabilistic models are trained on two different samples of jets, they can be used to compute likelihood ratios that distinguish between the two samples. Likelihood ratios provide theoretically optimal discriminants Neyman289 , which is indeed a major motivation for Junipr’s probabilistic approach. One can also sample from a probabilistic model in order to generate events, though generative models are bettersuited for this application deOliveira:2017pjk ; Paganini:2017hrr ; Paganini:2017dwg . In addition, one can use a probabilistic model to reweight events generated by an imperfect simulator, so that the reweighted events properly agree with data.
In this paper, as a proofofconcept, we use simulated data to train a basic implementation of the Junipr
framework described above. We have not yet attempted to optimize all of this implementation’s hyperparameters; however, we do find that a very simple architecture with no fine tuning is adequate. This is confirmed by its impressive discrimination power and its effective predictivity for a broad class of observables, but more rigorous testing is needed to determine whether this approach can provide stateoftheart results on the most pressing physics problems.
The general probabilistic model, its motivation, and a specific neural network implementation of it are discussed in Sec. 2. A comprehensive discussion of training the model, including the data used and potential subtleties in extending the model are covered in Sec. 3. Results on discrimination, generation, and reweighting are presented in Sec. 4. We provide robustness tests and some conceptually interesting results related to factorization in Sec. 5, including the counterintuitive anti shower generator. There are many ways to generalize our approach, as well as many applications that we do not fully explore in this work. We leave a discussion of some of these possible extensions to Sec. 6, where we conclude.
2 Unsupervised Learning in Jet Physics
To establish the framework clearly and generally, Sec. 2.1 begins by describing Junipr as a general probabilistic model, independent of the specific parametric form taken by the various functions it involves. From this perspective, such a probabilistic model could be implemented in many different ways. Sec. 2.2 then describes the particular neural network implementation of Junipr used in this paper, which has a simple but QCDcustomized architecture and minimal hyperparameter tuning.
2.1 General Probabilistic Model
Consider a set of finalstate 4momenta that we hereafter refer to as “the jet”. Junipr computes the probability density of this set of momenta arising in an event, assuming the event selection criteria used to select the training data. This probability distribution is normalized so that, abstractly,
(1) 
where the integral extends over the physical region of phase space. (In practice, in implementing Junipr we discretized the phase space into cells and assigned a measure of unity to each discrete cell. This results in being a discrete cellsizedependent probability distribution, but this choice is conceptually unimportant here.) A highlevel schematic of Junipr is shown in Fig. 1, which emphasizes that the model does not attempt to learn the quantummechanical evolution that created the jet, but only meaningfully predicts the likelihood of its finalstate momenta.
An unstructured model of the above form would ignore the fact that we know jet evolution is welldescribed by a semiclassical sequence of splittings, due to factorization theorems Coleman:1965xm ; Collins:1985ue ; Collins:1988ig ; feige2014hard ; feige2013shell . A model that ignores factorization would be much more opaque to interpretation, and have many more parameters than needed due to its unnecessary neutrality. Thus, we propose a model that describes a given configuration of finalstate momenta using sequential splittings. Such a sequence is defined by a jet clustering algorithm, which assigns a clustering tree to any set of finalstate momenta, so that a sequential decomposition of the probability distribution can be performed without loss of generality. We imagine fixing a specific algorithm to define the trees, so that there is no need to marginalize over all possible trees in computing a probability, a computation that would be intractable. While a deterministic clustering algorithm cannot directly describe the underlying quantummechanical parton evolution, that is not the goal for this model. With the algorithm set, the model as shown in Fig. 1 becomes that shown in Fig. 2.
We will now formalize this discussion into explicit equations. For the rest of this section we assume that the clustering tree is determined by a fixed jet algorithm (e.g. any of the generalized algorithms Cacciari:2008gp ; Cacciari:2011ma ). The particular algorithm chosen is theoretically inconsequential to the model, as the same probability distribution over final states will be learned for any choice. Practically speaking, however, certain algorithms may have advantages over others. We will discuss the choice of clustering algorithm further in Secs. 5.2 and 5.3.
The application of a clustering algorithm on the jet constituents defines a sequence of “intermediate states” . Here the superscript labels the intermediate state after the branching in the tree (where counting starts at 1) and the subscript enumerates momenta in that state. To be explicit,

the “initial state” consists of a single momentum: ;

at subsequent steps is gotten from by a single momentumconserving branching;

after the final branching, the state is the physical jet: .
In this notation, the probability of the jet (as shown in Fig. 2) can be written as
(2)  
Eq. (2) allows for a natural, sequential description of the jet. However, it obscures the factorization of QCD which predicts an approximately selfsimilar splitting evolution. Thus we decompose the model further, so that each in Eq. (2) is described by a branching function that only indirectly receives information about the rest of the jet. The latter is achieved via an unobserved representation vector of the global state of the jet at step . To be explicit, let denote the branching of a mother into daughters that achieves the transition from to in the clustering tree. Then we can write
(3) 
where is the mother’s discrete index in the intermediate state. We thus have a sequential model that at each step predicts

: probability over binary values for whether or not the tree ends;

: probability over indexing candidate mother momenta;

: probability over possible branchings.
Note that we have left the conditioning on implicit in and , since we will never need to use these functions when . In the product of Eq. (2.1), each subsequent factor is thus conditioned on the outcomes of previous factors, so that breaking up in this way is without loss of generality. In particular, no assumption has been made about the underlying physical processes that generate the data.
With these choices, we force the hidden representation
to encode all global information about the tree, since it must predict whether the tree ends, which momentum branches next, and the branching pattern. In fact, providing with the momenta that directly participate in the branching means that only needs to encode global information. We show that the global structure stored in is crucial for the model to predict the correct branching patterns in Sec. 5.1.2.2 Neural Network Implementation
For a neural network based implementation of the model defined by Eqs. (2) and (2.1), we use an RNN with hidden state augmented by dense neural networks for each of the three probability distributions in Eq. (2.1). The recurrent structure of this implementation is shown in Fig. 3, which emphasizes how the RNN’s hidden representation keeps track of the global state of the jet, by sequentially reading in the momenta that branched most recently.
The fact that learns and remembers the full jet, despite only being shown the two new momenta at step , is ensured by the tasks for which is responsible. These are shown in the detailed network diagram of Fig. 4. There one can see that is the only input into the components of the model that predict when the tree ends and which momentum is next to branch. The domains of the three probability functions in Eq. (2.1) are shown in Fig. 4 as well: is defined over the binary set corresponding to “end” or “not”; is multinomial over the set of candidate mothers; and is defined on the space of possible branchings, which is (a subset of) by momentum conservation. At each step, the model outputs the full probability distributions, which in mathematical notation are , , and .
Fig. 3 and Fig. 4 show how Junipr provides a probability distribution at each step given the momenta emerging from the preceding branching. For clarity, Fig. 5 separately shows how Junipr is used to evaluate the full probability density over finalstate momenta in a jet. At each step , the point in representing whether the tree ends, the point in representing which mother momentum branches, and the point in representing its daughters are plugged into the probability distributions to obtain the probabilities that should be assigned to the jet under consideration. The product of these three probabilities, taken over all steps, leads to .
Let us now go into detail about the neural network architecture used. We use basic RNN cells DeepLearningBook with activation,
(4) 
and found that a hidden representation vector of generic size 100 was sufficient for our needs. We found GRU GRUs and LSTM LSTMs
cells to be unnecessarily complex and highcapacity for the tasks carried out in this paper. This is in contrast to language modelling, for which basic RNN cells are underpowered. To see why this might heuristically be expected, note that a sentence containing 20 words is much more complex than a jet containing 20 momenta, because the words in the sentence are ordered, whereas the momenta in the jet are not. This introduces an additional factor of
to the complexity of language modelling. It is thus reasonable to expect that jet physics will not require all the highpowered tools designed for natural language processing.
For we use a fullyconnected network with
as input, a single hidden layer of size 100 with ReLU activation, and a sigmoid output layer. We use the same setup for
, the only difference being that the output layer is a softmax over the candidate mother momenta, ordered by energy. These choices are generic and not highly tuned. We found that Junipr works well for a very general set of architectures and sizes, so we stick with this simple setup.For the branching function we must describe the probability distribution over all possible configurations of daughter momenta consistent with the mother momentum . For this system, we use coordinates centered around the mother, where is the energy fraction of the softer daughter, () is the opening angle of the softer (harder) daughter, and specifies the plane in which the branching occurs. See Fig. 6 for a visualization of these coordinates.
There are two separate approaches one could take to model the branching function . Firstly, the variables could be treated as discrete, with outputting a softmax probability over discrete cells representing different values. Secondly, one could treat as a continuous variable and use an “energy model” of the form where is a normalizing partition function. In this work we predominantly adopt the former approach, as it is much faster, and most distributions are insensitive to the discretization of . However, we do train an energy model to show that models with continuous are possible, which we discuss in Sec. 3.4.
In the discrete case, we bin the possible values of into a 4dimensional grid with 10 bins per dimension, so that the entire grid has cells. For a given value of , we place a 1 in the bin corresponding to that value, and we place 0’s everywhere else. This 1hot encoding of the possible values of allows us to use a softmax function at the top layer of the neural network describing (see Fig. 4). Furthermore, we use a dense network with a single hidden layer of size 100 and ReLU activation for , just as we did for and . The hidden units in this network receive as input, as well as the mother momentum .
Thus we have a neural network implementation of Eqs. (2) and (2.1), with a representation of the evolving global jet state stored in , and with fullyconnected networks describing , , and . As defined above, the model has a single parameter matrix, mapping the branching function’s 100 dimensional hidden layer to its dimensional output layer, and has parameters elsewhere. One might refer to this implementation as Junipr, as one can imagine many alternative implementations within the Junipr framework that may prove useful in future applications. We will continue to use the term Junipr for brevity, to refer both to the framework and to the basic implementation described here.
3 Training and Validation
We now describe how to train the model outlined in Sec. 2.2. We begin by discussing the training data used, followed by our general approach to training and validation. Finally we discuss an alternative model choice that allows higher resolution on the particle momenta.
3.1 Training Data
To enable proofofconcept demonstrations of Junipr’s various applications, we train the implementation described in Sec. 2.2 using jets simulated in Pythia v8.226 Sjostrand:2006za ; Sjostrand:2014zea and clustered using FastJet v3.2.2 Cacciari:2011ma . We simulated 600k hemisphere jets in Pythia using the process at a centerofmass energy of 1 TeV, with hemispheres defined in FastJet using the exclusive algorithm Catani:1993hr ; Ellis:1993tq , and with an energy window of 450–550 GeV imposed on the jets. To create the deterministic trees that Junipr requires, we reclustered the jets using the C/A clustering algorithm Dokshitzer:1997in ; Wobisch:1998wt , with GeV and . The nonzero values of and make the input to Junipr formally infraredandcollinear safe, but this is by no means necessary. Furthermore, our approach is formally independent of the reclustering algorithm chosen. We demonstrate this by showing results using an absurd reclustering algorithm inspired by a 2D printer in Sec. 5.2, as well as for anti Cacciari:2008gp reclustering in Sec. 5.3.
Thus we have 600k quark jets with GeV and . We use 500k of these jets for training, with 10k set aside as a test set to monitor overfitting, and we use the remaining validation set of 100k jets to make the plots in this paper.
In the applications of Sec. 4, we also make use of several other data sets produced according to the above specifications, with small but important changes. We list these modifications here for completeness. In one case, quark jets from were required to lie in a very tight mass window of 90.7–91.7 GeV. A sample of boosted jets from events was also produced with the same mass cut. And finally, another sample of quark jets was produced, as detailed above, but with the value of in the final state shower changed from Pythia’s default value of 0.1365 to 0.11.
Before being fed to Junipr, jets in these data sets must be clustered, so that each jet becomes a tree of branchings ending in the finalstate momenta of the jet:
(5) 
where the momenta in one column are equal to those of the next column except for a single branching. At each step , only the momenta associated with this branching are fed into Junipr, as detailed in Sec. 2. With this setup, Junipr requires minimal parameters; it learns to update as the tree evolves by focusing only on the stepbystep changes to the jet. Note also that jets of arbitrary length can be considered.
Note that in implementing Junipr, we do not directly evaluate the branching function on the momenta but instead use the parameterization shown in Fig. 6. In fact, we use a nonlinear transformation of this parameterization:
(6) 
This invertible transformation simply maps the range of each coordinate onto , which reduces the amount of global parametric shift required in optimization. Similarly, we perform a transformation on the components of before feeding them into the update rule for in Eq. (4); we do the same for , the input to the branching function . This is a technical point that is not conceptually important.
3.2 Approach to Training
To train Junipr, we maximize the log likelihood over the full set of training data:
(7) 
For a particular jet with finalstate momenta we use Eqs. (2) and (2.1) to compute
(8)  
where is the index of the mother momentum at step in the training example and are its daughters. Maximizing the log likelihood in this way allows the model to learn each step in parallel, providing computational efficiency and stability.
For all models presented in this paper, we use basic stochastic gradient descent with the following learning rate and batch size schedule, where training proceeds from left to right:
Schedule  5 epochs 
5 epochs  5 epochs  5 epochs  5 epochs  5 epochs 

learning rate  
batch size  10  10  10  100  100  100 
We follow such a schedule to slowly increase the resolution and decrease the stochasticity of gradient descent throughout training. Decreasing the learning rate reduces the step size, thereby allowing finer details of the cost surface to be resolved. Increasing the batch size reduces the stochasticity by improving the sample estimates of the true gradients.
3.3 Validation of Model Components
Junipr is constructed as a probabilistic model for jet physics by expanding as a product over steps in the jet’s clustering tree, as shown in Eq. (2). Each step involves three components: the probability that the tree will end, the probability that a given momentum will be the next mother to branch, and the probability over the daughter momenta of the branching, as shown in Eq. (2.1). We now validate each of Junipr’s components using our validation set of 100k previously unseen Pythia jets. In this section, we present histograms of actual outcomes in the Pythia validation set (i.e. frequency distributions) as well as Junipr’s probabilistic output when evaluated on the jets in this data set (i.e. marginalized probability distributions) to check for agreement.
In Fig. 7 we show the probability that the tree should end, as a function of both intermediate state length and maximum particletojetaxis angle. In both cases we see excellent agreement with the validation data, demonstrating a good model fit with low underfitting and no overfitting. Note that Fig. 7 (left) is in onetoone correspondence with the jet constituent multiplicity, and that the shape of Fig. 7 (right) is a direct consequence of C/A clustering with . Indeed, if an opening angle near already exists in an angularordered tree, then there are likely no remaining branchings in the clustering tree.
In Fig. 8 we show the probability that a given candidate will be the next mother to branch in the clustering tree, as a function of both the candidate’s index (which is sorted to be decreasing in energy) and the candidate’s angle from the jet axis. The first of these results is shown in particular for the step in the clustering trees. We observe again that the model fits the validation data well. Note from Fig. 8 (left) that the highest energy branches of the clustering tree are most likely to undergo subsequent branchings, in line with the expectation at leading logarithmic accuracy. Fig. 8 (right) shows consistent predictions, since the highest energy branches also lie at the narrowest angles to the jet axis.
In Fig. 9 we show the branching function , the component of the model that predicts how a mother momentum should split into a pair of daughter momenta. We show the branching function results for and (i.e. with marginalized over the variables not shown) at the first step in the jet evolution , as well as at a later step . (See Fig. 6 for definitions of and and Eq. (6) for their ranges in the data.) This shows the dependency of the branching function on the evolving jet representation , which we will discuss in detail in Sec. 5.1. We see that for these direct predictions, Junipr fits the validation data almost perfectly. Note that in Fig. 9 (top) soft wideangle emissions are the norm at the earliest steps, as expected with the C/A clustering algorithm. In Fig. 9 (bottom) one can see that later in the clustering trees, harder morecollinear branchings are commonplace. It bears repeating that these trends are highly dependent on the chosen clustering algorithm and have no precise connection to the underlying physical processes generating the data.
3.4 Increasing the Branching Function Resolution
In this section, we discuss increasing the resolution of the branching function
(9) 
including the case where is an energy model over continuous . (The coordinates were defined in Fig. 6.) This technical section can easily be skipped without loss of the logical flow of the paper.
We begin by briefly discussing increasing the resolution of the branching function over discrete , the case described in Sec. 2.2. The first thing to note is that with a softmax over 4dimensional , the size of the matrix multiplication required in a dense network is quartic in the number of bins used for each dimension. We generically use 10 bins for each of resulting in an output size of . (In fact we use 10 linearly spaced bins in the transformed coordinates of Eq. (6), and this can be seen on the logarithmic axes of Fig. 9, but this detail is not conceptually important.) Given this quartic scaling, simply increasing the number of discrete cells quickly becomes prohibitively computationally expensive. Potential solutions to this problem include: (i) using a hierarchical softmax HierarchicalSoftmax1 ; HierarchicalSoftmax2
, and (ii) simply interpolating between the discrete bins of the model.
In a hierarchical softmax, a lowresolution probability is predicted first, say with cells, then another celled distribution is predicted inside the chosen lowresolution cell. In principle, this gives resolution at only twice the computational time required for resolution. We briefly implemented the hierarchical softmax, and preliminary tests found it to work efficiently, but perhaps with a decrease in training stability. We chose not to pursue the hierarchical softmax further in this work, primarily because we have not seen the need for resolution much higher than discrete cells.
Due to its ease of use, we do employ linear interpolation between the discrete bins in our baseline model with resolution . This comes at no extra training cost, and removes most of the effects of discretization on the observable distributions generated by sampling from Junipr; see Sec. 4.2.
We now turn to the continuous version of Junipr in which the branching function is given by an undirected energy model:
(10) 
To model , we again use a fullyconnected network with hidden layer of size 100, as used everywhere else, except here the output layer is left to be linear. We perform the integral over using importance sampling:
(11) 
where is the set of ’s sampled from the importance distribution .
Unlike the discrete version of Junipr, where training is relatively straightforward, the continuous version requires a nonstandard technique in training the branching function . This is because, although Eq. (11) provides an unbiased approximation to ,
(12) 
this leads to a biased estimate of the log likelihood, since
(13) 
by Jensen’s inequality. Thus, every gradient step taken is systematically different from the true gradient, and this bias derails training, especially near convergence when the true gradient becomes small.
To overcome this problem, we start by computing the sample variance on our estimate
, which is(14) 
Then the percenterror in our biased estimate of the gradient is approximately
(15) 
This error propagates into the log likelihood, causing the bias in Eq. (13). To mitigate this, we adopt a policy of monitoring during training, and whenever increases above some value (a hyperparameter that we set to 2%) we double the sample size used to compute . This slows down training considerably, but it effectively reduces the bias in our gradient estimates. Note that while generic importance sampling typically fails in higher dimensions, our branching function lives in only 4 dimensions, so this approach is robust using any reasonable importance distribution
. Indeed, we found that a uniform distribution over the transformed coordinates of Eq. (
6) is a fine choice for .In Fig. 10 we show results for Junipr trained with the continuous branching function as described above. In this case, we can use arbitrarily highresolution binning, as Junipr has learned a fully continuous probability density. Fig. 10 can be roughly compared to Fig. 9, where we were required to use 10 bins for each dimension of .
To close this section, we note that in most cases, we expect the discretized branching function with 10 bins per dimension of to be sufficient, especially if one performs a linear interpolation on the output cells. This simple case is certainly faster to train and does not require the technique described here to avoid biased gradient estimates.
4 Applications and Results
With Junipr trained and validated, we turn to some of the most interesting results it enables. Given a jet, Junipr can compute the probability density associated with the momenta inside the jet, conditioned on the criteria used to select the training data. To visualize this, we show a C/Aclustered Pythia jet in Fig. 11 with the Juniprcomputed probability associated with each branching written near that node in the tree. Note that these are small discretized probabilities due to the discretized implementation of Junipr’s branching function described in Sec. 2. This is shown primarily to conceptualize the model, which is constructed to be quite interpretable as it is broken down to compute the probability of each step in the clustering history of a jet.
A direct and powerful application of the Junipr framework, enabled by having access to separate probabilistic models of different data sources, is in discrimination based on likelihood ratios. We discuss discrimination in Sec. 4.1, along with a highly intuitive way of visualizing it. In contrast, an instinctive but indirect use of Junipr as a probabilistic model is in sampling new jets from it. We discuss the observable distributions generated through sampling in Sec. 4.2. However, sampling from a probabilistic model is often inefficient (e.g. slower than Pythia) compared to evaluating probabilities of jets directly. In Sec. 4.3 we discuss reweighting samples from one simulator to match those of another distribution. In principle, this could be used to tweak Pythia samples to match observed collider data simply by reweighting.
4.1 Likelihood Ratio Discrimination
We expect that one of the most exciting applications of Junipr will be in discriminating the underlying physics that could have created a jet.^{2}^{2}2We thank Kyle Cranmer for an early discussion on this topic. For example, suppose we had two sets of jets, one set corresponding to decays of a boosted boson, the other set simply highenergy quarks. We could then train one copy of Junipr on just the boosted sample, giving the probability distribution , and another copy of Junipr on just the quark jets, giving . Finally, for any new jet we could determine whether the jet was initiated by a boosted or by a highenergy quark by looking at the likelihood ratio:
(16) 
where the threshold is set according to the location on the ROC (receiver operating characteristic) curve desired for the discrimination task at hand. In contrast to approaches that try to compute likelihood ratios like this using QCD
Soper:2011cr ; Soper:2014rya , the Junipr approach can learn the separate probability distributions directly from samples of training data.Discrimination based on the likelihood ratio theoretically provides the most statistically powerful discriminant between two hypotheses Neyman289 . Moreover, our setup takes into account all the momenta that define a specific type of jet. Note also that for the task of pairwise discrimination between jet types, this unsupervised approach requires training probabilistic models, whereas a supervised learning approach would require training classifiers. Thus, we expect likelihoodratio discrimination using Junipr to provide a powerful tool.
We note further that we do not even require pure samples of the two underlying processes between which we would like to discriminate Metodiev:2017vrx . Thus, it would be feasible to discriminate based solely on real collider data. In our /quark example above, we would simply train one copy of Junipr on a sample of predominantly boosted jets, and train another copy on predominantly quark jets, and the likelihood ratio of those two models would still be theoretically optimal for /quark discrimination.
In order to get a first look at the potential of likelihoodratio discrimination using Junipr, we continue with the /quark example discussed above. We use Pythia to simulate and events at a centerofmass energy of 1 TeV. We impose a very tight mass window, 90.7 – 91.7 GeV, on the jets in each data set, so that no discrimination power can be gleaned from the jet mass. More details on the generation of the data sets were given in Sec. 3.1. We admit that a more compelling example of discrimination power would be for quark and gluon jets at hadron colliders, but we leave a proper treatment of that important case to future work. The toy scenario studied here serves both to prove that the probabilities output by Junipr are meaningful, and that likelihood ratio discrimination using unsupervised probabilistic models is a promising application of the Junipr framework.
In Fig. 12 we show the /quark separation power achieved by Junipr, both in terms of full likelihood ratio distributions for validation sets of and quark jets, as well as the resulting ROC curve. For comparison, in Fig. 12 we also show the ROC curve achieved using a 2D likelihood ratio discriminant based on 2subjettiness Thaler:2010tr and multiplicity. Junipr’s likelihoodratio discrimination is clearly superior to that based on combining the most natural observables: 2subjettiness, multiplicity, (and keep in mind the tight mass cut). Of course, these observables do not provide stateoftheart discrimination power even in this toy scenario, but we include the comparison in this proofofconcept to provide a sense of scale on the plot.
By design, Junipr naturally processes the information in jets via a recurrent mechanism that tracks the evolution of their clustering trees, and this allows users to peer inside at this structure and access the probabilities at each branching. In particular, we can consider the likelihood ratio at each step in the clustering trees to understand which branchings give rise to the greatest discrimination power. We show this in Fig. 13, where it is clear that Junipr can extract useful discriminatory information at most branchings.
Indeed, visualizing jets as in Fig. 13 can provide a number of insights. Unsurprisingly, we see for the quark jet (on the top) that the likelihood ratio of the first branching is rather extreme, at , since it is unlike the energybalanced first branching associated with boosted jets. However, we also see that almost all subsequent branchings are also unlike those expected in boosted jets, and they combine to provide comparable discrimination power to the first branching alone. Many effects probably contribute to this separation power at later branchings, including that quark jets often gain their mass throughout their evolution instead of solely at the first branching, and that the quark jet is colorconnected to other objects in the global event. Such effects have proven to be useful for discrimination in other contexts Chien:2017xrb .
Similarly, considering the boosted jet on the bottom of Fig. 13 shows that significant discrimination power comes not only from the first branching, but also from subsequent splittings, as the boosted jet evolves as a colorsinglet pair. Note the presence of the predictive secondary emissions sent from one quarksubjet toward the other. This is reminiscent of the pull observable, which has proven useful for discrimination in other contexts Gallicchio:2010sw . More generally, the importance of the energy distribution, opening angles, multiplicity, and branching pattern in highperformance discrimination can be understood from such pictures.
We are very excited by the prospect of visualizing Junipr’s discrimination power on jets, based on the likelihood ratio it assigns at each branching in their clustering trees, as in Fig. 13. Such visualizations could provide intuition that leads to the development of new, humaninterpretable, perhaps calculable observables for discrimination in important contexts.
We would like to make one side note about discrimination, before moving on to the next application of Junipr. The statement that likelihoodratio discrimination is optimal of course only applies in the limit of perfect models. Since this limit is never fully realized, one may worry that discrimination with Junipr may in fact be suboptimal. Since the two probabilistic models we use for discrimination are each trained individually to replicate a certain type of jet, they are not conditioned to focus on the differences between the two jet types, which may be very subtle in the case of a difficult discrimination task. In the realistic case of slightly imperfect models, it may be advantageous for discrimination purposes to instead train the two models to focus on the differences. To be specific, one could train the two models on the two data sets simultaneously, with the goal being to maximize the likelihood ratio on one data set and minimize it on the other. Following this method in the particular example of /quark discrimination used above, one would train the and models on data sets and to maximize the following quantity:
(17) 
Compare this to the approach we have taken above, namley training and to separately maximize the log likelihood of Eq. (7) on their corresponding sets of training data. This alternative training method would correspond to optimizing Junipr for the application of discrimination, leaving intact our ability to visualize discrimination power in clustering trees, but sacrificing the probabilistic interpretation of the model’s output. We have not tested training with Eq. (17), and thus cannot attest to its practicality, but we suspect an approach along these lines may be useful in certain contexts.
4.2 Generation from JUNIPR
We now turn to a more familiar approach to jet physics, but a somewhat less appropriate usage of Junipr models: sampling new jets from the learned probability distribution to generate traditional observable distributions. We include this application here, not only to demonstrate this capability, but also to further validate the distribution learned by Junipr during unsupervised training.
Sampling from Junipr is relatively efficient; one simply samples from the low dimensional distributions at each step and feeds those samples forward as input to subsequent steps. In this way, one generates a full jet in many steps, as detailed in Fig. 14.
We used the baseline implementation of Junipr trained on quark jets, as described in Sec. 3, to generate 100k jets in this way. The resulting jet mass and constituent multiplicity distributions are plotted in Fig. 15 where both distributions sampled from Junipr match those created from our validation set of 100k Pythia jets withheld from training. Reasonable agreement can also be seen in the 2D distributions of Fig. 16.
However, there are two reasons why we do not consider Junipr to be built for generation. (These drawbacks could be avoided with a generative model; see deOliveira:2017pjk ; Paganini:2017hrr ; Paganini:2017dwg .) The first is simply that sampling from probability distributions is generally difficult. As we just showed, it turns out that Junipr is relatively easy to sample from, due to its sequential structure and the fact that distributions are lowdimensional at each step. Despite this, sampling jets from Junipr is still much slower than generation with, for example, Pythia.
The second reason is more fundamental. With a sequential model structured as Junipr is, probability distributions at late steps in generation are highly sensitive to the draws made at earlier steps. Very small defects in the probability distributions at early steps cause feedback in the model that amplifies those errors. Furthermore, as a partially generated jet becomes more misrepresentative of the training data, the resulting probability distributions used at later steps are less trained, which can result in a runaway effect. All of this is to say that, for the purpose of generating jets, Junipr’s accuracy at early steps is disproportionately important. This is in tension with the training method undertaken in Sec. 3.2, namely the maximization of the loglikelihood, which prioritizes all branchings equally. Thus, we should expect that some observable distributions generated by sampling jets from Junipr might agree worse with the validation set of Pythia data than otherwise expected. We mention in passing that this second drawback could be mitigated by reweighting jets after generation, as detailed in Sec. 4.3 below.
In fact, we have found empirically that the Nsubjettiness ratio observables computed by sampling from Junipr do not match the heldout Pythia data perfectly. This can be seen in Fig. 17 with the 2subjettiness distribution, where the difference between the two distributions is more significant.
We consider this disagreement to be both expected and nondiminishing of Junipr’s potential. Indeed, in the next section we will show how to overcome this issue, by generating samples consistent with Junipr’s learned probabilistic model, without ever sampling from it. In particular, the disagreement in Fig. 17 will be rectified in Fig. 18.
4.3 Reweighting Monte Carlo Events
Another application of the Junipr framework is to reweight events. For example, suppose we trained Junipr on data from the Large Hadron Collider (LHC) to yield a probabilistic model . Then one could generate a sample of new events using a relatively accurate Monte Carlo simulator, train another instance of Junipr on that sample to yield , and finally reweight the simulated events by evaluated on an eventbyevent basis. This process yields a sample of events that is theoretically equivalent to the LHC data used in training . The advantage of such an approach is that Junipr can correct the simulated events on different levels, for example using the data reclustered in subjets as we have done in this paper. However, the full simulated event has the complete hadron distributions and can thereby be interfaced with a detector simulation. This is in many ways a simpler approach than trying to improve the simulation directly through the dark art of MonteCarlo tuning.
This reweighting is identical to importance sampling from a proposal distribution given by the simulated data distribution . For example, suppose one wanted to measure the distribution of an observable at the LHC, which is given by
(18) 
where the last approximation is associated with collecting a finite amount of LHC data in order to measure the distribution. (The reader can substitute discretized delta functions appropriate for histogramming if averse to the singular notation used in these equations.) Instead of using real data, if say a public version of were available, then anyone could calculate this observable distribution using only simulated data sampled from as follows:
(19) 
In this way, one could efficiently obtain samples of arbitrary size from by reweighting samples generated by an efficient simulator. The only limitation to this process is that the simulated data must be similar to the actual target data, so that they have overlapping regions of support (formal requirement) and the weights are not too far from unity (efficiency requirement).
As with the likelihoodratio discrimination in Sec. 4.1, here we will show results in a toy scenario as a proofofprinciple. Ideally a model trained on LHC data, with all related complications, would be used to reweight Monte Carlo jets to make the simulated data indiscernible from LHC data; we leave a proper study of this to future work.
Instead, here we use two samples of jets generated using two different versions of Pythia. We reweight jets from one of the samples and demonstrate their agreement with the other sample. In particular, we use our baseline Junipr model trained on Pythia
generated quark jets as our “true distribution”. For the moment, we will refer to this model as
, since its training data was generated using Pythia’s default value of in the final state shower. As our “simulated distribution” we will use , which was trained on quark jets generated with coupling parameter changed to in Pythia’s finalstate shower. (See Sec. 3.1 for a more indepth description of the training data used.) Our goal is to show that reweighting jets from the “simulated distribution” according to the likelihood ratio leads to observables in agreement with the “true distribution”.In Fig. 18 we demonstrate that this is indeed the case. We check this for both the 2subjettiness and 3subjettiness ratio observables, as well as the jet shape observable. On the left side of Fig. 18, one can see that in all cases, the distribution is clearly different from the distribution. On the right side of Fig. 18, one finds that the two distributions come into relatively good agreement once the jets are reweighted by . This also provides further confirmation that Junipr learns subtle correlations between constituent momenta inside jets.
Note that it was the 2subjettiness ratio observable that Junipr struggled to predict well through direct sampling (see Fig. 17), whereas when reweighting another set of samples, Junipr matches the data well on this observable (see topright of Fig. 18). This corroborates the discussion in Sec. 4.2 concerning the difficulties in sampling directly from Junipr.
Before closing this section, let us reiterate one point mentioned above. For the procedure of reweighting events to be practical, the weights used should not be radically different from unity, meaning that the two distributions generating the two samples should not be too different. If this condition is not satisfied, then away from the limit of infinite statistics, a few events with very large weights could vastly overpower the rest of the events, leading to a choppy reweighted distribution with large statistical uncertainties. To avoid this problem in the toy scenario explored in this section, we found it necessary to discard roughly of the jets in the
sample which were outliers with
. These outliers were uncorrelated with the observables shown, and we believe they resulted from imperfections in the trained model. It is clear that much more needs to be understood about the application of reweighting, but this would perhaps be more effectively done in the context of a specific task of interest involving LHC data.5 Factorization and JUNIPR
In the previous section, we showed some preliminary but very exciting results for likelihoodratio discrimination and for the generation of observables by reweighting simulated jets. Both of these applications require access to an unsupervised probabilistic model. Next we discuss some of the more subtle internal workings of Junipr, which are intimately related to the underlying physics of factorization.
In particular, we show that the hidden representation indeed stores important global information about intermediate states of jets in Sec. 5.1. We then discuss the clusteringalgorithm independence of Junipr by considering two distinct clustering algorithms: a “printer” algorithm in Sec. 5.2, where momenta are processed lefttoright and toptobottom as if by an inkjet printer; and the anti algorithm in Sec. 5.3, which allows us to present another counterintuitive result, the anti shower generator.
5.1 The Encoding of Global Information
We have constructed Junipr so that all global information about the jet is contained in the RNN’s hidden state . Only the branching function receives the local branching information in addition to . This forces to contain all the information needed to predict when the shower should end, , to predict which momentum should branch next, , and to inform the branching function of the relevant global structure. As the primary feature vector for all three of these distinct tasks, must learn an effective representation of the jet at evolution step .
To explicitly show that stores important global information about the intermediate jet state at step , we train a new model on our baseline quark jet data (see Sec. 3.1) with the difference that we remove as an input to the branching function . We expect that such a “local” branching model will not evolve correctly as the global jet structure evolves, since all global information is being withheld. This is indeed what we find, as can be seen in Fig. 19. On the left side of that figure, the evolution of the distribution (defined in Fig. 6) from to is shown using 100k Pythia jets from our heldout set of validation data. There we see the gradual decrease in angle as expected for C/A trees. On the right side of Fig. 19, the evolution of the branching function is shown for the “local” branching model, and the disagreement between this damaged model and Pythia is clear. Note that this prediction of incorrect distributions at intermediate branchings in the C/A tree will inevitably lead to an incorrect probability distribution over finalstate momenta.
While we do not show the corresponding results from our baseline (global) model in Fig. 19 to avoid clutter, the agreement with Pythia is essentially perfect, as one would expect from the similar check performed in Fig. 9. This confirms the success of the jet representation in supplying the branching function with important information about the global structure.
5.2 Clustering Algorithm Independence
Another subtle aspect of Junipr is its theoretical clustering algorithm independence. In principle, the model as described in Sec. 2.1 is indeed independent of the chosen algorithm, which is fixed simply to avoid a sum over all possible trees consistent with the finalstate momenta. That is, for each clustering procedure chosen by the user, a different model is learned, but one that describes the same probability distribution over finalstate momenta, at least formally.
However, it is not guaranteed that a given neuralnetwork implementation of Junipr will work well for every clustering algorithm. We have chosen an architecture that stores the global jet physics in the RNN’s hidden state and the local branching physics in the branching function . This architecture is motivated by the factorizing structure of QCD, and thus Junipr will most easily learn jet trees that are most similar to QCD — our primary reason for predominantly using the C/A algorithm. Consequently, though the model described in Sec. 2.1 is formally independent of clustering algorithm, the particular implementation adopted in Sec. 2.2 may weakly depend on the chosen algorithm by virtue of the ease with which it can learn the data.
To put this to the test, we have introduced a jet clustering algorithm that is nothing like QCD, but more like a 2D printer.^{3}^{3}3We thank Eric Metodiev for this suggestion. The “printer” clustering algorithm scans the 2D jet image (i.e. the cross sectional image perpendicular to the jet axis) from righttoleft and bottomtotop, clustering particles as it encounters them. Run in reverse (i.e. as a shower) particles are emitted from the jet core from lefttoright and toptobottom; this is how a jet image would be printed by an inkjet printer with a single printing tip. In Fig. 20 we show a single Pythia jet clustered using the printer algorithm. As can be seen in the jet image on the right side of Fig. 20, momenta are indeed emitted toptobottom. On the left side of Fig. 20, we see that any collinear branching structure is completely absent from the clustering tree; instead, particles are steadily emitted upandtotheleft.
Though Junipr’s neural network architecture is not optimized for the informational structure of the printer algorithm, it is still able to learn the structure, by relying much more heavily on the the jet representation . We demonstrate this by training Junipr on our data set of Pythiagenerated quark jets (see Sec. 3.1) clustered with the printer algorithm, thus yielding the probabilistic model . Indeed, in Fig. 21 one can see a jet sampled from , which correctly follows the printer structure.
As expected, however, the distributions sampled from are not quite as good as our C/A results. On the left side of Fig. 22 we show the 2dimensional distribution over jet mass and constituent multiplicity generated using 100k jets sampled directly from . Comparing to the distribution generated by Pythia (see the left side of Fig. 16) this distribution matches well. However, for the 2subjettiness ratio observable on the right side of Fig. 22 we get a somewhat worse match to the Pythia validation data; compare this to the results of the C/A model in Fig. 17. Of course, we discussed in Sec. 4.2 why we do not expect direct sampling from Junipr to be perfectly reliable (and we discussed a way around this in Sec. 4.3), but it is still clear that such distributions are comparably worse when using the printer clustering algorithm, instead of the more natural C/A algorithm.
5.3 Anti Shower Generator
Reassured by the results of the previous section, we next consider Junipr trained on Pythia jets reclustered with anti Cacciari:2008gp . Like the printer algorithm, anti does not approximate the natural collinear structure of QCD. Unlike the printer algorithm, however, anti is a very commonly used tool. For the latter reason we explore anti jets here.
Perhaps the most interesting result associated with an anti version of Junipr is that it provides access to an anti shower generator. Generating an anti shower is counterintuitive, because the anti algorithm generally clusters soft emissions onebyone with the hard jet core. Thus, a generator must remember where previous emissions landed in order to send subsequent emissions nearby. This is required to reproduce the correct collinear structure in the distribution of finalstate of momenta. Said in another way, since the collinear factorization of QCD is not built into the anti clustering algorithm, a local (or factorized) anti generator could not produce emissions with the correct collinear distribution. Thus, we should expect that, in an anti version of Junipr, higher demands will be placed on the jet representation to monitor all the radiation in the jet. This is certainly possible, but not the task for which our neural network architecture is optimized.
To see to what extent an anti implementation of Junipr relies on the global information stored in , we trained two models on Pythiagenerated quark jets clustered with anti (see Sec. 3.1 for more details on the training data used). One model, , has the baseline architecture outlined in Sec. 2. The other, , is a local branching model like the one used in Sec. 5.1, in which the global representation is withheld as input to the branching function.
In Fig. 23 (bottom) we show a jet sampled from . In this case, though the tree itself does not properly guide the collinear structure of emissions, one can see that the emission directions are highly correlated with one another, demonstrating the success of the jet representation in tracking the global branching pattern. In Fig. 23 (top) we show for comparison a jet sampled from , in which the branching function does not receive . In the latter case, all correlation between the emission directions is lost. This shows that the global representation is crucial for a successful anti branching model.
In Fig. 24 we show the 2dimensional distribution over jet mass and constituent multiplicity, as well as the 2subjettiness distribution, generated with . One can see that the former distribution is consistent with the distribution generated by Pythia in Fig. 16. Mild disagreement between ’s 2subjettiness distribution and Pythia’s can be seen on the right side of Fig. 24. This is on par with the agreement obtained by sampling from the C/A model in Fig. 17.
In Sec. 5.1 we saw that the RNN’s hidden state manages the global information in Junipr’s neural network architecture. This is an efficient and natural way to characterize QCDlike jets, and therefore also C/A clustering trees. Though Junipr is formally independent of jet algorithm (i.e. in the infinitecapacity and perfecttraining limit), we might expect Junipr’s performance to degrade somewhat when paired with clustering algorithms that require significantly more information to be stored in . This was explored in Secs. 5.2 and 5.3 using two separate nonQCDlike clustering algorithms, namely the “printer” and anti algorithms. Despite these clustering algorithms being unnatural choices for Junipr, we were able to demonstrate conceptually interesting and novel results, such as the anti shower generator. This further demonstrates that Junipr can continue to function well, even when the clustering algorithm chosen for implementation bears little resemblance to the underlying physical processes that generate the data.
6 Conclusions and Outlook
In this paper, we have introduced Junipr as a framework for unsupervised machine learning in particle physics. The framework calls for a neural network architecture designed to efficiently describe the leadingorder physics of splittings, alongside a representation of the global jet physics. This requires the momenta in a jet to be clustered into a binary tree. The choice of clustering algorithm is not essential to Junipr’s performance, but choosing an algorithm that has some correspondence with an underlying physical model, such as the angularordered parton shower in quantum chromodynamics, gives improved performance and allows for intrerpretability of the network. At Junipr’s core is a recurrent neural network with three interconnected components. It moves along the jet’s clustering tree, evaluating the likelihood of each branching. More generally, Junipr is a function that acts on a set of 4momenta in an event to compute their relative differential cross section, i.e. the probability density for this event to occur, given the event selection criteria used to select the training sample. One of the appealing features of Junipr is its interpretability: it provides a desconstruction of the probability density into contributions from each point in the clustering history.
There are many promising applications of Junipr, and we have only been able to touch on a few proofofconcept tests in this introductory work. One exciting use case is discrimination. In contrast to supervised models which directly learn to discriminate between two samples, Junipr learns the features of the samples separately. It then discriminates by comparing the likelihood of a given event with respect to alternative models of the underlying physics. The resulting likelihood ratio provides theoretically optimal statistical power. As an example, we showed that Junipr can discriminate between boosted bosons and quark jets (in a very tight mass window around ) in events when trained on the two samples separately. With Junipr, it is not only possible to perform powerful discrimination using unsupervised learning, but the discrimination power can be visualized over the entire clustering tree of each jet, as in Fig. 13. This opens new avenues for physicists to gain intuition about the physics underlying highperformance discrimination. Such studies might even inspire the construction of new calculable observables.
Another exciting potential application of Junipr is the reweighting of Monte Carlo events, in order to improve agreement with real collider data. A proofofconcept of this idea was given in Fig. 18, where jets generated with one Pythia tune were reweighted to match jets generated with another. The reason this application is important is that current Monte Carlo event generators do an excellent job of simulating events on average, but are limited by the models and parameters within them. It may be easier to correct for systematic bias in event generation by a small reweighting factor appropriate for a particular data sample, rather than by trying to isolate and improve faulty components of the model. In this context, Junipr can be thought of as providing small but highly granular tweaks to simulations in order to improve agreement with data.
The Junipr framework was used here to compute the likelihood that a given set of particle momenta will arise inside a jet. One can also imagine more general models that act on all the momenta in an entire event, including particle identification tags, or even lowlevel detector data. A particularly interesting direction would be to consider applying Junipr to heavy ion collisions, in which the medium where the jets are produced and evolve is not yet well understood. In this case, comparing the probabilities in data to those of simulation could give insights into how to improve the simulations, or more optimistically, to improve our understanding of the underlying physics.
Acknowledgements.
We benefited from interesting discussions with D. Barber, E. Bernton, A. Botev, Y.T. Chien, K. Cranmer, R. Elayavalli, M. Freytsis, B. Gaujac, R. Habib, P. Komiske, E. Metodiev, B. Nachman, and J. Thaler. AA, CF, and MDS are supported in part by the Department of Energy under contract DESC0013607. Support for AA and CF was provided in part by the Harvard Data Science Initiative.References
 (1) A. Krizhevsky, I. Sutskever and G. E. Hinton, Imagenet classification with deep convolutional neural networks, in Advances in neural information processing systems, pp. 1097–1105, 2012.
 (2) K. He, X. Zhang, S. Ren and J. Sun, Deep residual learning for image recognition, pp. 770–778, 2016. 1512.03385.
 (3) G. Huang, Z. Liu and K. Q. Weinberger, Densely connected convolutional networks, 2017. 1608.06993.
 (4) D. Bahdanau, K. Cho and Y. Bengio, Neural machine translation by jointly learning to align and translate, 2014. 1409.0473.
 (5) Y. Wu, M. Schuster, Z. Chen, Q. V. Le, M. Norouzi, W. Macherey et al., Google’s neural machine translation system: Bridging the gap between human and machine translation, 1609.08144.
 (6) A. Graves and N. Jaitly, Towards endtoend speech recognition with recurrent neural networks, 2014.
 (7) A. Van Den Oord, S. Dieleman, H. Zen, K. Simonyan, O. Vinyals, A. Graves et al., Wavenet: A generative model for raw audio, 2016. 1609.03499.
 (8) D. E. Rumelhart, G. E. Hinton and R. J. Williams, Learning representations by backpropagating errors, nature 323 (1986) 533.
 (9) T. Mikolov, M. Karafiát, L. Burget, J. Černockỳ and S. Khudanpur, Recurrent neural network based language model, in Eleventh Annual Conference of the International Speech Communication Association, 2010.
 (10) ATLAS collaboration, G. Aad et al., A neural network clustering algorithm for the ATLAS silicon pixel detector, JINST 9 (2014) P09009, [1406.7690].
 (11) ATLAS collaboration, G. Aad et al., Performance of Jet Identification in the ATLAS Experiment, JINST 11 (2016) P04008, [1512.01094].
 (12) CMS collaboration, S. Chatrchyan et al., Performance of taulepton reconstruction and identification in CMS, JINST 7 (2012) P01001, [1109.6034].
 (13) K. Datta and A. Larkoski, How Much Information is in a Jet?, JHEP 06 (2017) 073, [1704.08249].
 (14) K. Datta and A. J. Larkoski, Novel Jet Observables from Machine Learning, JHEP 03 (2018) 086, [1710.01305].
 (15) H. Luo, M.x. Luo, K. Wang, T. Xu and G. Zhu, Quark jet versus gluon jet: deep neural networks with highlevel features, 1712.03634.
 (16) P. T. Komiske, E. M. Metodiev and J. Thaler, Energy flow polynomials: A complete linear basis for jet substructure, 1712.07124.
 (17) J. Gallicchio, J. Huth, M. Kagan, M. D. Schwartz, K. Black and B. Tweedie, Multivariate discrimination and the Higgs + W/Z search, JHEP 04 (2011) 069, [1010.3698].
 (18) ATLAS Collaboration collaboration, Identification of HadronicallyDecaying W Bosons and Top Quarks Using HighLevel Features as Input to Boosted Decision Trees and Deep Neural Networks in ATLAS at = 13 TeV, Tech. Rep. ATLPHYSPUB2017004, CERN, Geneva, Apr, 2017.

(19)
J. Cogan, M. Kagan, E. Strauss and A. Schwarztman,
JetImages: Computer Vision Inspired Techniques for Jet Tagging
, JHEP 02 (2015) 118, [1407.5675]. 
(20)
L. de Oliveira, M. Kagan, L. Mackey, B. Nachman and A. Schwartzman,
Jetimages — deep learning edition
, JHEP 07 (2016) 069, [1511.05190].  (21) P. T. Komiske, E. M. Metodiev and M. D. Schwartz, Deep learning in color: towards automated quark/gluon jet discrimination, JHEP 01 (2017) 110, [1612.01551].
 (22) P. T. Komiske, E. M. Metodiev, B. Nachman and M. D. Schwartz, Pileup Mitigation with Machine Learning (PUMML), JHEP 12 (2017) 051, [1707.08600].
 (23) G. Kasieczka, T. Plehn, M. Russell and T. Schell, Deeplearning Top Taggers or The End of QCD?, JHEP 05 (2017) 006, [1701.08784].
 (24) W. Bhimji, S. A. Farrell, T. Kurth, M. Paganini, Prabhat and E. Racah, Deep Neural Networks for Physics Analysis on lowlevel wholedetector data at the LHC, in 18th International Workshop on Advanced Computing and Analysis Techniques in Physics Research (ACAT 2017) Seattle, WA, USA, August 2125, 2017, 2017. 1711.03573.
 (25) ATLAS Collaboration collaboration, Quark versus Gluon Jet Tagging Using Jet Images with the ATLAS Detector, Tech. Rep. ATLPHYSPUB2017017, CERN, Geneva, Jul, 2017.
 (26) S. Macaluso and D. Shih, Pulling Out All the Tops with Computer Vision and Deep Learning, 1803.00107.
 (27) Y.T. Chien and R. Kunnawalkam Elayavalli, Probing heavy ion collisions using quark and gluon jet substructure, 1803.03589.
 (28) J. Pearkes, W. Fedorko, A. Lister and C. Gay, Jet Constituents for Deep Neural Network Based Top Quark Tagging, 1704.02124.
 (29) G. Louppe, K. Cho, C. Becot and K. Cranmer, QCDAware Recursive Neural Networks for Jet Physics, 1702.00748.
 (30) T. Cheng, Recursive Neural Networks in Quark/Gluon Tagging, 1711.02633.
 (31) S. Egan, W. Fedorko, A. Lister, J. Pearkes and C. Gay, Long ShortTerm Memory (LSTM) networks with jet constituents for boosted top tagging at the LHC, 1711.09059.
 (32) K. Fraser and M. D. Schwartz, Jet Charge and Machine Learning, 1803.08066.
 (33) D. Guest, J. Collado, P. Baldi, S.C. Hsu, G. Urban and D. Whiteson, Jet Flavor Classification in HighEnergy Physics with Deep Neural Networks, Phys. Rev. D94 (2016) 112002, [1607.08633].
 (34) ATLAS Collaboration collaboration, Identification of Jets Containing Hadrons with Recurrent Neural Networks at the ATLAS Experiment, Tech. Rep. ATLPHYSPUB2017003, CERN, Geneva, Mar, 2017.
 (35) E. M. Metodiev, B. Nachman and J. Thaler, Classification without labels: Learning from mixed samples in high energy physics, JHEP 10 (2017) 174, [1708.02949].
 (36) T. Cohen, M. Freytsis and B. Ostdiek, (Machine) Learning to Do More with Less, JHEP 02 (2018) 034, [1706.09451].
 (37) P. T. Komiske, E. M. Metodiev, B. Nachman and M. D. Schwartz, Learning to Classify from Impure Samples, 1801.10158.
 (38) E. M. Metodiev and J. Thaler, On the Topic of Jets, 1802.00008.
 (39) L. de Oliveira, M. Paganini and B. Nachman, Learning Particle Physics by Example: LocationAware Generative Adversarial Networks for Physics Synthesis, Comput. Softw. Big Sci. 1 (2017) 4, [1701.05927].
 (40) M. Paganini, L. de Oliveira and B. Nachman, Accelerating Science with Generative Adversarial Networks: An Application to 3D Particle Showers in Multilayer Calorimeters, Phys. Rev. Lett. 120 (2018) 042003, [1705.02355].
 (41) M. Paganini, L. de Oliveira and B. Nachman, CaloGAN : Simulating 3D high energy particle showers in multilayer electromagnetic calorimeters with generative adversarial networks, Phys. Rev. D97 (2018) 014021, [1712.10321].
 (42) J. Neyman and E. S. Pearson, Ix. on the problem of the most efficient tests of statistical hypotheses, Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 231 (1933) 289–337, [http://rsta.royalsocietypublishing.org/content/231/694706/289.full.pdf].
 (43) A. Butter, G. Kasieczka, T. Plehn and M. Russell, Deeplearned Top Tagging with a Lorentz Layer, 1707.08966.
 (44) S. Coleman and R. Norton, Singularities in the physical region, Nuovo Cim. 38 (1965) 438–442.
 (45) J. C. Collins, D. E. Soper and G. F. Sterman, Factorization for Short Distance Hadron  Hadron Scattering, Nucl.Phys. B261 (1985) 104.
 (46) J. C. Collins, D. E. Soper and G. F. Sterman, Soft Gluons and Factorization, Nucl.Phys. B308 (1988) 833.
 (47) I. Feige and M. D. Schwartz, Hardsoftcollinear factorization to all orders, Physical Review D 90 (2014) 105020.
 (48) I. Feige and M. D. Schwartz, An onshell approach to factorization, Physical Review D 88 (2013) 065021.
 (49) S. Catani, Y. L. Dokshitzer, M. H. Seymour and B. R. Webber, Longitudinally invariant clustering algorithms for hadron hadron collisions, Nucl. Phys. B406 (1993) 187–224.
 (50) S. D. Ellis and D. E. Soper, Successive combination jet algorithm for hadron collisions, Phys. Rev. D48 (1993) 3160–3166, [hepph/9305266].
 (51) Y. L. Dokshitzer, G. D. Leder, S. Moretti and B. R. Webber, Better jet clustering algorithms, JHEP 08 (1997) 001, [hepph/9707323].
 (52) M. Wobisch and T. Wengler, Hadronization corrections to jet crosssections in deep inelastic scattering, in Monte Carlo generators for HERA physics. Proceedings, Workshop, Hamburg, Germany, 19981999, pp. 270–279, 1998. hepph/9907280.
 (53) S. D. Ellis, A. Hornig, T. S. Roy, D. Krohn and M. D. Schwartz, Qjets: A NonDeterministic Approach to TreeBased Jet Substructure, Phys. Rev. Lett. 108 (2012) 182003, [1201.1914].
 (54) D. Kahawala, D. Krohn and M. D. Schwartz, Jet Sampling: Improving Event Reconstruction through Multiple Interpretations, JHEP 06 (2013) 006, [1304.2394].
 (55) L. Mackey, B. Nachman, A. Schwartzman and C. Stansbury, Fuzzy Jets, JHEP 06 (2016) 010, [1509.02216].
 (56) D. E. Soper and M. Spannowsky, Finding physics signals with shower deconstruction, Phys. Rev. D84 (2011) 074002, [1102.3480].
 (57) D. E. Soper and M. Spannowsky, Finding physics signals with event deconstruction, Phys. Rev. D89 (2014) 094005, [1402.1189].
 (58) M. Cacciari, G. P. Salam and G. Soyez, The Antik(t) jet clustering algorithm, JHEP 04 (2008) 063, [0802.1189].
 (59) M. Cacciari, G. P. Salam and G. Soyez, FastJet User Manual, Eur. Phys. J. C72 (2012) 1896, [1111.6097].
 (60) I. Goodfellow, Y. Bengio and A. Courville, Deep Learning. MIT Press, 2016.
 (61) K. Cho, B. van Merrienboer, C. Gulcehre, D. Bahdanau, F. Bougares, H. Schwenk et al., Learning phrase representations using rnn encoder–decoder for statistical machine translation, in Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing (EMNLP), pp. 1724–1734, Association for Computational Linguistics, 2014. DOI.
 (62) S. Hochreiter and J. Schmidhuber, Long shortterm memory, Neural Comput. 9 (Nov., 1997) 1735–1780.
 (63) T. Sjostrand, S. Mrenna and P. Z. Skands, PYTHIA 6.4 Physics and Manual, JHEP 05 (2006) 026, [hepph/0603175].
 (64) T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten et al., An Introduction to PYTHIA 8.2, Comput. Phys. Commun. 191 (2015) 159–177, [1410.3012].
 (65) Theano Development Team collaboration, R. AlRfou, G. Alain, A. Almahairi, C. Angermueller, D. Bahdanau, N. Ballas et al., Theano: A Python framework for fast computation of mathematical expressions, arXiv eprints abs/1605.02688 (May, 2016) .
 (66) F. Morin and Y. Bengio, Hierarchical probabilistic neural network language model, in AISTATS’05, pp. 246–252, 2005.
 (67) A. Mnih and G. E. Hinton, A scalable hierarchical distributed language model, in Advances in Neural Information Processing Systems 21 (D. Koller, D. Schuurmans, Y. Bengio and L. Bottou, eds.), pp. 1081–1088. Curran Associates, Inc., 2009.
 (68) J. Thaler and K. Van Tilburg, Identifying Boosted Objects with Nsubjettiness, JHEP 03 (2011) 015, [1011.2268].
 (69) Y.T. Chien, A. Emerman, S.C. Hsu, S. Meehan and Z. Montague, Telescoping jet substructure, 1711.11041.
 (70) J. Gallicchio and M. D. Schwartz, Seeing in Color: Jet Superstructure, Phys. Rev. Lett. 105 (2010) 022001, [1001.5027].