QCD-aware recursive neural networks for jet physics

Recent progress in applying machine learning for jet physics has been built upon an analogy between calorimeters and images. In this work, we present a novel class of recursive neural networks built instead upon an analogy between QCD and natural languages. In the analogy, four-momenta are like words and the clustering history of sequential recombination jet algorithms is like the parsing of a sentence. Our approach works directly with the four-momenta of a variable-length set of particles, and the jet-based tree structure varies on an event-by-event basis. Our experiments highlight the flexibility of our method for building task-specific jet embeddings and show that recursive architectures are significantly more accurate and data efficient than previous image-based networks. We extend the analogy from individual jets (sentences) to full events (paragraphs), and show for the first time an event-level classifier operating on all the stable particles produced in an LHC event.


Introduction
By far the most common structures seen in collisions at the Large Hadron Collider (LHC) are collimated sprays of energetic hadrons referred to as 'jets'. These jets are produced from the fragmentation and hadronization of quarks and gluons as described by quantum chromodynamics (QCD). Several goals for the LHC are centered around the treatment of jets, and there has been an enormous amount of effort from both the theoretical and experimental communities to develop techniques that are able to cope with the experimental realities while maintaining precise theoretical properties. In particular, the communities have converged on sequential recombination jet algorithms, methods to study jet substructure, and grooming techniques to provide robustness to pileup.
One compelling physics challenge is to search for highly boosted standard model particles decaying hadronically. For instance, if a hadronically decaying W boson is highly boosted, then its decay products will merge into a single fat jet with a characteristic substructure. Unfortunately, there is a large background from jets produced by more mundane QCD processes. For this reason, several jet 'taggers' and variables sensitive to jet substructure have been proposed. Initially, this work was dominated by techniques inspired by our intuition and knowledge of QCD; however, more recently there has been a wave of approaches that eschew this expert knowledge in favor of machine learning techniques. In this paper, we present a hybrid approach that leverages the structure of sequential recombination jet algorithms and deep neural networks.
Recent progress in applying machine learning techniques for jet physics has been built upon an analogy between calorimeters and images [1][2][3][4][5][6][7][8]. These methods take a variablelength set of 4-momenta and project them into a fixed grid of η − φ towers or 'pixels' to produce a 'jet image'. The original jet classification problem, hence, reduces to an image classification problem, lending itself to deep convolutional networks and other machine learning algorithms. Despite their promising results, these models suffer from the fact that they have many free parameters and that they require large amounts of data for training. More importantly, the projection of jets into images also loses information, which impacts classification performance. The most obvious way to address this issue is to use a recurrent neural network to process a sequence of 4-momenta as they are. However, it is not clear how to order this sequence. While p T ordering is common in many contexts [5], it does not capture important angular information critical for understanding the subtle structure of jets.
In this work, we propose instead a solution for jet classification based on an analogy between QCD and natural languages, as inspired by several works from natural language processing [9][10][11][12][13][14]. Much like a sentence is composed of words following a syntactic structure organized as a parse tree, a jet is also composed of 4-momenta following a structure dictated by QCD and organized via the clustering history of a sequential recombination jet algorithm. More specifically, our approach uses 'recursive' networks where the topology of the network is given by the clustering history of a sequential recombination jet algorithm, which varies on an event-by-event basis. This event-by-event adaptive structure can be contrasted with the 'recurrent' networks that operate purely on sequences (see e.g., [15]). The network is therefore given the 4-momenta without any loss of information, in a way that also captures substructures, as motivated by physical theory.
It is convenient to think of the recursive neural network as providing a 'jet embedding', which maps a set of 4-momenta into R q . This embedding has fixed length and can be fed into a subsequent network used for classification or regression. Thus the procedure can be used for jet tagging or estimating parameters that characterize the jet, such as the masses of resonances buried inside the jet. Importantly, the embedding and the subsequent network can be trained jointly so that the embedding is optimized for the task at hand.
Extending the natural language analogy paragraphs of text are sequence of sentences, just as event are sequence of jets. In particular, we propose to embed the full particle content of an event by feeding a sequence of jet-embeddings into a recurrent network. As before, this event-level embedding can be fed into a subsequent network used for classification or regression. To our knowledge, this represents the first machine learning model operating on all the detectable particles in an event.
The remainder of the paper is structured as follows. In Sec. 2, we formalize the classification tasks at the jet-level and event-level. We describe the proposed recursive network architectures in Sec. 3 and detail the data samples and preprocessing used in our experiments in Sec. 4. Our results are summarized and discussed first in Sec. 5 for experiments on a jet-level classification problem, and then in Sec. 6 for experiments on an event-level classification problem. In Sec. 7, we relate our work to close contributions from deep learning, natural language processing, and jet physics. Finally, we gather our conclusions and directions for further works in Sec. 8.

Problem statement
We describe a collision event e ∈ E as being composed of a varying number of particles, indexed by i, and where each particle is represented by its 4-momentum vector v i ∈ R 4 , such that e = {v i |i = 1, . . . , N }.
The 4-momenta in each event can be clustered into jets with a sequential recombination jet algorithm that recursively combines (by simply adding their 4-momenta) the pair i, i that minimize while d α ii is less than min(p 2α ti , p 2α ti ) [16,17]. These sequential recombination algorithms have three hyper-parameters: R, p t,min , α, and jets with p t < p t,min are discarded. At that point, the jet algorithm has clustered e into M jets, each of which can be represented by a binary tree t j ∈ T indexed by j = 1, . . . , M with N j leaves (corresponding to a subset of the v i ). In the following, we will consider the specific cases where α = 1, 0, −1, which respectively correspond to the k t , Cambridge-Aachen and anti-k t algorithms.
In addition to jet algorithms, we consider a 'random' baseline that corresponds to recombining particles at random to form random binary trees t j , along with 'asc-p T ' and 'desc-p T ' baselines, which correspond to degenerate binary trees formed from the sequences of particles sorted respectively in ascending and descending order of p T .
For jet-level classification or regression, each jet t j ∈ T in the training data comes with labels or regression values y j ∈ Y jet . In this framework, our goal is to build a predictive model f jet : T → Y jet minimizing some loss function L jet . Similarly, for event-level classification or regression, we assume that each collision event e l ∈ E in the training data comes with labels or regression values y l ∈ Y event , and our goal is to build a predictive model f event : E → Y event minimizing some loss function L event .

Individual jets
Let us first consider the case of an individual jet whose particles are topologically structured as a binary tree t j , e.g., based on a sequential recombination jet clustering algorithm or a simple sequential sorting in p T . Let k = 1, . . . , 2N j − 1 indexes the node of the binary tree t j , and let the left and right children of node k be denoted by k L and k R respectively. Let also k L always be the hardest child of k. By construction, we suppose that leaves k map to particles i(k) while internal nodes correspond to recombinations. Using these notations, we recursively define the embedding h jet k ∈ R q of node k in t j as Classifier Jet embedding Figure 1. QCD-motivated recursive jet embedding for classification. For each individual jet, the embedding h jet 1 (t j ) is computed recursively from the root node down to the outer nodes of the binary tree t j . The resulting embedding is chained to a subsequent classifier, as illustrated in the top part of the figure. The topology of the network in the bottom part is distinct for each jet and is determined by a sequential recombination jet algorithm (e.g., k t clustering).
and b u ∈ R q form together the shared parameters to be learned, q is the size of the embedding, σ is the ReLU activation function [18], and g is a function extracting the kinematic features p, η, θ, φ, E, and p T from the 4-momentum o k . When applying Eqn. 3.1 recursively from the root node k = 1 down to the outer nodes of the binary tree t j , the resulting embedding, denoted h jet 1 (t j ), effectively summarizes the information contained in the particles forming the jet into a single vector. In particular, this recursive neural network (RNN) embeds a binary tree of varying shape and size into a vector of fixed size. As a result, the embedding h jet 1 (t j ) can now be chained to a subsequent classifier or regressor to solve our target supervised learning problem, as illustrated in Figure 1. All parameters (i.e., of the recursive jet embedding and of the classifier) are learned jointly using backpropagation through structure [9] to minimize the loss L jet , hence tailoring the embedding to the specific requirements of the task. Further implementation details, including an efficient batched computation over distinct binary trees, can be found in Appendix C. ...
... f event (e) Classifier Figure 2. QCD-motivated event embedding for classification. The embedding of an event is computed by feeding the sequence of pairs (v(t j ), h jet 1 (t j )) over the jets it is made of, where v(t j ) is the unprocessed 4-momentum of the jet t j and h jet 1 (t j ) is its embedding. The resulting event-level embedding h event M (e) is chained to a subsequent classifier, as illustrated in the right part of the figure.
In addition to the recursive activation of Eqn. 3.1, we also consider and study its extended version equipped with reset and update gates (see details in Appendix A). This gated architecture allows the network to preferentially pass information along the left-child, right-child, or their combination.
While we have not performed experiments, we point out that there is an analogous style of architectures based on jet algorithms with 2 → 3 recombinations [17,19,20].

Full events
We now embed entire events e of variable size by feeding the embeddings of their individual jets to an event-level sequence-based recurrent neural network.
As an illustrative example, we consider here a gated recurrent unit [21] (GRU) operating on the p T ordered sequence of pairs (v(t j ), h jet 1 (t j )), for j = 1, . . . , M , where v(t j ) is the unprocessed 4-momentum of the jet t j and h jet 1 (t j ) is its embedding. The final output h event M (e) (see Appendix B for details) of the GRU is chained to a subsequent classifier to solve an event-level classification task. Again, all parameters (i.e., of the inner jet embedding function, of the GRU, and of the classifier) are learned jointly using backpropagation through structure [9] to minimize the loss L event . Figure 2 provides a schematic of the full classification model. In summary, combining two levels of recurrence provides a QCDmotivated event-level embedding that effectively operates at the hadron-level for all the particles in the event.
In addition and for the purpose of comparison, we also consider the simpler baselines where i) only the 4-momenta v(t j ) of the jets are given as input to the GRU, without augmentation with their embeddings, and ii) the 4-momenta v i of the constituents of the event are all directly given as input to the GRU, without grouping them into jets or providing the jet embeddings.

Data, Preprocessing and Experimental Setup
In order to focus attention on the impact of the network architectures and the projection of input 4-momenta into images, we consider the same boosted W tagging example as used in Refs. [1,2,4,6]. The signal (y = 1) corresponds to a hadronically decaying W boson with 200 < p T < 500 GeV, while the background (y = 0) corresponds to a QCD jet with the same range of p T .
We are grateful to the authors of Ref. [6] for sharing the data used in their studies. We obtained both the full-event records from their PYTHIA benchmark samples, including both the particle-level data and the towers from the DELPHES detector simulation. In addition, we obtained the fully processed jet images of 25×25 pixels, which include the initial R = 1 anti-k t jet clustering and subsequent trimming, translation, pixelisation, rotation, reflection, cropping, and normalization preprocessing stages detailed in Ref. [2,6].
Our training data was collected by sampling from the original data a total of 100,000 signal and background jets with equal prior. The testing data was assembled similarly by sampling 100,000 signal and background jets, without overlap with the training data. For direct comparison with Ref. [6], performance is evaluated at test time within the restricted window of 250 < p T < 300 and 50 ≤ m ≤ 110, where the signal and background jets are re-weighted to produce flat p T distributions. Results are reported in terms of the area under the ROC curve (ROC AUC) and of background rejection (i.e., 1/FPR) at 50% signal efficiency (R =50% ). Average scores reported include uncertainty estimates that come from training 30 models with distinct initial random seeds. About 2% of the models had technical problems during training (e.g., due to numerical errors), so we applied a simple algorithm to ensure robustness: we discarded models whose R =50% was outside of 3 standard deviations of the mean, where the mean and standard deviation were estimated excluding the five best and worst performing models.
For our jet-level experiments we consider as input to the classifiers the 4-momenta v i from both the particle-level data and the DELPHES towers. We also compare the performance with and without the projection of those 4-momenta into images. While the image data already included the full pre-processing steps, when considering particle-level and tower inputs we performed the initial R = 1 anti-k t jet clustering to identify the constituents of the highest p T jet t 1 of each event, and then performed the subsequent translation, rotation, and reflection pre-processing steps (omitting cropping and normalization). When processing the image data, we inverted the normalization that enforced the sum of the squares of the pixel intensities be equal to one. 1 For our event-level experiments we were not able to use the data from Ref. [6] because the signal sample corresponded to pp → W (→ J)Z(→ νν) and the background to pp → jj. Thus the signal was characterized by one high-p T jet and large missing energy from Z(→ νν) which is trivially separated from the dijet background. For this reason, we generated our own PYTHIA and DELPHES samples of pp → W → W (→ J)Z(→ J) and QCD background such that both the signal and background have two high-p T jets. We use m W = 700 GeV and restrictp t of the 2 → 2 scattering process to 300 <p t < 350 GeV. Our focus is to demonstrate the scalability of our method to all the particles or towers in an event, and not to provide a precise statement about physics reach for this signal process. In this case each event e was clustered by the same anti-k t algorithm with R = 1, and then the constituents of each jet were treated as in Sec. 3.1 (i.e., reclustered using k t or a sequential ordering in p T to provide the network topology for a non-gated embedding). Additionally, the constituents of each jet were pre-processed with translation, rotation, and reflection as in the individual jet case. Training was carried out on a dataset of 100,000 signal and background events with equal prior. Performance was evaluated on an independent test set of 100,000 other events, as measured by the ROC AUC and R =80% of the model predictions. Again, average scores are given with uncertainty estimates that come from training 30 models with distinct initial random seeds.
In both jet-level and event-level experiments, the dimension of the embeddings q was set to 40. Training was conducted using Adam [22] as an optimizer for 25 epochs, with a batch size of 64 and a learning rate of 0.0005 decayed by a factor of 0.9 after every epoch. These parameters were found to perform best on average, as determined through an optimization of the hyper-parameters. Performance was monitored during training on a validation set of 5000 samples to allow for early stopping and prevent from overfitting.

Performance studies
We carried out performance studies where we varied the following factors: the projection of the 4-momenta into an image, the source of those 4-momenta, the topology of the RNN, and the presence or absence of gating.
Impact of image projection The first factor we studied was whether or not to project the 4-momenta into an image as in Refs. [2,6]. The architectures used in previous studies required a fixed input (image) representation, and cannot be applied to the variable length set of input 4-momenta. Conversely, we can apply the RNN architecture to the discretized image 4-momenta. Table 1 shows that the RNN architecture based on a k t topology performs almost as well as the MaxOut architecture in Ref. [6] when applied to the image pre-processed 4-momenta coming from DELPHES towers. Importantly the RNN architecture is much more data efficient. While the MaxOut architecture in Ref. [6] has 975,693 parameters and was trained with 6M examples, the non-gated RNN architecture has 8,481 parameters and was trained with 100,000 examples only. Next, we compare the RNN classifier based on a k t topology on tower 4-momenta with and without image preprocessing. Table 1 and Fig. 3 show significant gains in not using jet images, improving ROC AUC from 0.8321 to 0.8807 (resp., R =50% from 12.7 to 24.1) in the case of k t topologies. In addition, this result outperforms the MaxOut architecture operating on images by a significant margin. This suggests that the projection into an image loses information and impacts classification performance. We suspect the loss of information to be due to some of the construction steps of jet images (i.e., pixelisation, rotation, zooming, cropping and normalization). In particular, all are applied at the imagelevel instead of being performed directly on the 4-momenta, which might induce artefacts due to the lower resolution, particle superposition and aliasing. By contrast, the RNN is able to work directly with the 4-momenta of a variable-length set of particles, without any loss of information. For completeness, we also compare to the performance of a classifier based purely on the single n-subjettiness feature τ 21 := τ2 /τ1 and a classifier based on two features (the trimmed mass and τ 21 ) [23]. In agreement with previous results based on deep learning [2,6], we see that our RNN classifier clearly outperforms this variable.
Measurements of the 4-momenta The second factor we varied was the source of the 4-momenta. The towers scenario, corresponds to the case where the 4-momenta come from the calorimeter simulation in DELPHES. While the calorimeter simulation is simplistic, the granularity of the towers is quite large (10 • in φ) and it does not take into account that tracking detectors can provide very accurate momenta measurements for charged particles that can be combined with calorimetry as in the particle flow approach. Thus, we also consider the particles scenario, which corresponds to an idealized case where the 4-momenta come from perfectly measured stable hadrons from PYTHIA. Table 1 and Fig. 3 show that further gains could be made with more accurate measurements of the 4-momenta, improving e.g. ROC AUC from 0.8807 to 0.9185 (resp., R =50% from 24.1 to 68.3) in the case of k t topologies. We also considered a case where the 4-momentum came from the DELPHES particle flow simulation and the data associated with each particle was augmented with a particle-flow identifier distinguishing ± charged hadrons, photons, and neutral hadrons. This is similar in motivation to Ref. [7], but we did not observe any significant gains in classification performance with respect to the towers scenario.
Topology of the binary trees The third factor we studied was the topology of the binary tree t j described in Sections 2 and 3.1 that dictates the recursive structure of the RNN. We considered binary trees based on the anti-k t , Cambridge-Aachen (C/A), and k t sequential recombination jet algorithms, along with random, asc-p T and desc-p T binary trees. Table 1 and Fig. 4 show the performance of the RNN classifier based on these various topologies. Interestingly, the topology is significant.
For instance, k t and C/A significantly outperform the anti-k t topology on both tower and particle inputs. This is consistent with intuition from previous jet substructure studies where jets are typically reclustered with the k t algorithm. The fact that the topology is important is further supported by the poor performance of the random binary tree topology. We expected however that a simple sequence (represented as a degenerate binary tree) based on ascending and descending p T ordering would not perform particularly well, particularly since the topology does not use any angular information. Surprisingly, the simple descending p T ordering slightly outperforms the RNNs based on k t and C/A topologies. The descending p T network has the highest p T 4-momenta near the root of the tree, which we expect to be the most important. We suspect this is the reason that the descending . Jet classification performance for various input representations of the RNN classifier, using k t topologies for the embedding. The plot shows that there is significant improvement from removing the image processing step and that significant gains can be made with more accurate measurements of the 4-momenta.
p T outperforms the ascending p T ordering on particles, but this is not supported by the performance on towers. A similar observation was already made in the context of natural languages [24][25][26], where tree-based models have at best only slightly outperformed simpler sequence-based networks. While recursive networks appear as a principled choice, it is conjectured that recurrent networks may in fact be able to discover and implicitly use recursive compositional structure by themselves, without supervision.
Gating The last factor that we varied was whether or not to incorporate gating in the RNN. Adding gating increases the number of parameters to 48,761, but this is still about 20 times smaller than the number of parameters in the MaxOut architectures used in previous jet image studies. Table 1 shows the performance of the various RNN topologies with gating. While results improve significantly with gating, most notably in terms of R =50% , the trends in terms of topologies remain unchanged.
Other variants Finally, we also considered a number of other variants. For example, we jointly trained a classifier with the concatenated embeddings obtained over k t and antik t topologies, but saw no significant performance gain. We also tested the performance of recursive activations transferred across topologies. For instance, we used the recursive activation learned with a k t topology when applied to an anti-k t topology and observed a significant loss in performance. We also considered particle and tower level inputs with an additional trimming preprocessing step, which was used for the jet image studies, but we  . This plot shows that topology is significant, as supported by the fact that results for k t , C/A and desc-p T topologies improve over results for anti-k t , asc-p T and random binary trees. Best results are achieved for C/A and desc-p T topologies, depending on the metric considered.
saw a significant loss in performance. While the trimming degraded classification performance, we did not evaluate the robustness to pileup that motivates trimming and other jet grooming procedures.

Infrared and Collinear Safety Studies
In proposing variables to characterize substructure, physicists have been equally concerned with classification performance and the ability to ensure various theoretical properties of those variables. In particular, initial work on jet algorithms focused on the Infrared-Collinear (IRC) safe conditions: • Infrared safety. The model is robust to augmenting e with additional particles {v N +1 , . . . , v N +K } with small transverse momentum.
• Collinear safety. The model is robust to a collinear splitting of a particle, which is represented by replacing a particle v j ∈ e with two particles v j 1 and v j 2 , such that The sequential recombination algorithms lead to an IRC-safe definition of jets, in the sense that given the event e, the number of jets M and their 4-momenta v(t j ) are IRC-safe.
An early motivation of this work is that basing the RNN topology on the sequential recombination algorithms would provide an avenue to machine learning classifiers with some theoretical guarantee of IRC safety. If one only wants to ensure robustness to only one soft particle or one collinear split, this could be satisfied by simply running a single iteration of the jet algorithm as a pre-processing step. However, it is difficult to ensure a more general notion of IRC safety on the embedding due to the non-linearities in the network. Nevertheless, we can explicitly test the robustness of the embedding or the subsequent classifier to the addition of soft particles or collinear splits to the input 4-momenta. Table 2 shows the results of a non-gated RNN trained on the nominal particle-level input when applied to testing data with additional soft particles or collinear splits. The collinear splits were uniform in the momentum fraction and maintained the small invariant mass of the hadrons. We considered one or ten collinear splits on both random particles and the highest p T particles. We see that while the 30 models trained with a descending p T topology very slightly outperform the k t topology for almost scenarios, their performance in terms of R =50% decreases relatively more rapidly when collinear splits are applied (see e.g., the collinear10-max scenarios where the performance of k t decreases by 4%, while the performance of p T decreases by 10%). This suggests a higher robustness towards collinear splits for recursive networks based on k t topologies.
We also point out that the training of these networks is based solely on the classification loss for the nominal sample. If we are truly concerned with the IRC-safety considerations, then it is natural to augment the training of the classifiers to be robust to these variations. A number of modified training procedures exist, including e.g., the adversarial training procedure described in Ref. [27].

Experiments with event-level classification
As in the previous section, we carried out a number of performance studies. However, our goal is mainly to demonstrate the relevance and scalability of the QCD-motivated approach we propose, rather than making a statement about the physics reach of the signal process. Results are discussed considering the idealized particles scenario, where the 4-momenta come from perfectly measured stable hadrons from PYTHIA. Experiments for the towers scenario (omitted here) reveal similar qualitative conclusions, though performance was slightly worse for all models, as expected.
Number of jets The first factor we varied was the maximum number of jets in the sequence of embeddings given as input to the GRU. While the event-level embedding can be computed over all the jets it is constituted by, QCD suggests that the 2 highest p T jets hold most of the information to separate signal from background events, with only marginal discriminating information left in the subsequent jets. As Table 3 and Fig. 5 show, there is indeed significant improvement in going from the hardest jet to the 2 hardest jets, while there is no to little gain in considering more jets. Let us also emphasize that the event-level models have only 18,681 parameters, and were trained on 100,000 training examples.
Topology of the binary trees The second factor we studied was the architecture of the networks used for the inner embedding of the jets, for which we compare k t against descend- Table 2. Performance of pre-trained RNN classifiers (without gating) applied to nominal and modified particle inputs. The collinear1 (collinear10) scenarios correspond to applying collinear splits to one (ten) random particles within the jet. The collinear1-max (collinear10-max) scenarios correspond to applying collinear splits to the highest p T (ten highest p T ) particles in the jet. The soft scenario corresponds to adding 200 particles with p T = 10 −5 GeV uniformly in 0 < φ < 2π and −5 < η < 5.

Scenario
Architecture ing p T topologies. As in the previous section, best results are achieved with descending p T topologies, though the difference is only marginal.
Other variants Finally, we also compare with baselines. With respect to an event-level embedding computed only from the 4-momenta v(t j ) (for j = 1, . . . , M ) of the jets, we find that augmenting the input to the GRU with jet-level embeddings yields significant improvement, e.g. improving ROC AUC from 0.9606 to 0.9875 (resp. R =80% from 21.1 to 174.5) when considering the 2 hardest jets case. This suggests that jet substructures are important to separate signal from background events, and correctly learned when nesting embeddings. Similarly, we observe that directly feeding the GRU with the 4-momenta v i , for i = 1, . . . , N , of the constituents of the event performs significantly worse. While performance remains decent (e.g., with a ROC AUC of 0.8925 when feeding the 50 4momenta with largest p T ), this suggests that the recurrent network fails to leverage some of the relevant information, which is otherwise easier to identify and learn when inputs to the GRU come directly grouped as jets, themselves structured as trees. In contrast to our previous results for jet-level experiments, this last comparison underlines the fact that integrating domain knowledge by structuring the network topology is in some cases crucial for performance. Overall, this study shows that event embeddings inspired from QCD and produced by nested recurrence, over the jets and over their constituents, is a promising avenue for building effective machine learning models. To our knowledge, this is the first classifier Table 3. Summary of event classification performance. Best results are achieved through nested recurrence over the jets and over their constituents, as motivated by QCD.

Related work
Neural networks in particle physics have a long history. They have been used in the past for many tasks, including early work on quark-gluon discrimination [28,29], particle identification [30], Higgs tagging [31] or track identification [32]. In most of these, neural networks appear as shallow multi-layer perceptrons where input features were designed by experts to incorporate domain knowledge. More recently, the success of deep convolutional networks has triggered a new body of work in jet physics, shifting the paradigm from engineering input features to learning them automatically from raw data, e.g., as in these works treating jets as images [1][2][3][4][5][6][7][8]. Our work builds instead upon an analogy between QCD and natural languages, hence complementing the set of algorithms for jet physics with techniques initially developed for natural language processing [9][10][11][12][13][14]. In addition, our approach does not delegate the full modeling task to the machine. It allows to incorporate domain knowledge in terms of the network architecture, specifically by structuring the recursion stack for the embedding directly from QCD-inspired jet algorithms (see Sec. 3) Between the time that this work appeared on the arXiv preprint server and submitted for publication there has been a flurry of activity connecting deep learning techniques and jet physics (for reviews see Refs. [33][34][35]). In particular the method described here was also used for quark/gluon tagging in Ref. [36] and a variant of this method was used to reconstruct a jet's charge [37]. The authors of Ref. [38] used the tree structure defined by the jet clustering history to define a substructure ordering scheme for use with a sequential recurrent neural network. Going the opposite direction, graph neural networks and message passing neural networks have also been applied to the the same jet-level classification problem and data described in this work [39]. There has also been a spate of recent work on using QCD-inspired variables, enforcing physical constraints into neural networks, and ensuring infrared safety of neural network based approaches to jet physics [40][41][42][43][44]. Exploring a complementary direction, several authors have developed ways to train machine learning techniques using real data to avoid sensitivity to systematic effects in the simulation [45][46][47][48]. These recent training techniques are agnostic to the network architecture and can be paired with the RNN architectures described here. Finally, deep learning techniques are now being studied as generative models for jets, where the tree-based model mimics the parton shower and can be trained on real data [49]. Learning generative modes for both signal and background classes of jets can be used to define a classifier in which each branch of the tree can be interpreted as a contribution to a likelihood ratio discriminant [49].

Conclusions
Building upon an analogy between QCD and natural languages, we have presented in this work a novel class of recursive neural networks for jet physics that are derived from sequential recombination jet algorithms. Our experiments have revealed that preprocessing steps applied to jet images during their construction (specifically the pixelisation) loses information, which impacts classification performance. By contrast, our recursive network is able to work directly with the four-momenta of a variable-length set of particles, without the loss of information due to discretization into pixels. Our experiments indicate that this results in significant gains in terms of accuracy and data efficiency with respect to previous image-based networks. Finally, we also showed for the first time a hierarchical, event-level classification model operating on all the hadrons of an event. Notably, our results showed that incorporating domain knowledge derived from jet algorithms and encapsulated in terms of the network architecture led to improved classification performance.
While we initially expected recursive networks operating on jet recombination trees to outperform simpler p T -ordered architectures, our results still clearly indicate that the topology has an effect on the final performance of the classifier. However, our initial studies indicate that architectures based on jet trees are more robust to infrared radiation and collinear splittings than the simpler p T -ordered architectures, which may outweigh what at face value appears to be a small loss in performance. Accordingly, it would be natural to include robustness to pileup, infrared radiation, and collinear splittings directly in the training procedure [27]. Moreover, it is compelling to think of generalizations in which the optimization would include the topology used for the embedding as learnable component instead of considering it fixed a priori. An immediate challenge of this approach is that a discontinuous change in the topology (e.g., from varying α or R) makes the loss non-differentiable and rules out standard back propagation optimization algorithms. Nevertheless, solutions for learning composition orders have recently been proposed in NLP, using either explicit supervision [50] or reinforcement learning [51]; both of which could certainly be adapted to jet embeddings. Another promising generalization is to use a graph-convolutional network that operates on a graph where the vertices correspond to particle 4-momenta v i and the edge weights are given by d α ii or a similar QCD-motivated quantity [39,[52][53][54][55][56][57][58]. In conclusion, we feel confident that there is great potential in hybrid techniques like this that incorporate physics knowledge and leverage the power of machine learning.
[58] T. N. Kipf  A Gated recursive jet embedding The recursive activation proposed in Sec. 3.1 suffers from two critical issues. First, it assumes that left-child, right-child and local node information h jet k L , h jet k R , u k are all equally relevant for computing the new activation, while only some of this information may be needed and selected. Second, it forces information to pass through several levels of nonlinearities and does not allow to propagate unchanged from leaves to root. Addressing these issues and generalizing from [12][13][14], we recursively define a recursive activation equipped with reset and update gates as follows: and b u ∈ R q form together the shared parameters to be learned, σ is the ReLU activation function and denotes the element-wise multiplication. Intuitively, the reset gates r L , r R and r N control how to actively select and then merge the left-child embedding h jet k L , the right-child embedding h jet k R and the local node information u k to form a new candidate activationh jet k . The final embedding h jet k can then be regarded as a choice among the candidate activation, the left-child embedding, the right-child embedding and the local node information, as controlled by the update gates z H , z L , z R and z N . Finally, let us note that the proposed gated recursive embedding is a generalization of Section 3.1, in the sense that the later corresponds to the case where update gates are set to z H = 1, z L = 0, z R = 0 and z N = 0 and reset gates to r L = 1, r R = 1 and r N = 1 for all nodes k.

B Gated recurrent event embedding
In this section, we formally define the gated recurrent event embedding introduced in Sec. 3.2. Our event embedding function is a GRU [21] operating on the p T ordered sequence of pairs (v(t j ), h jet 1 (t j )), for j = 1, . . . , M , where v(t j ) is the unprocessed 4-momentum (φ, η, p T , m) of the jet t j and h jet 1 (t j ) is its embedding. Its final output h event j=M is recursively defined as follows: where W hx ∈ R r×4+q , W hh ∈ R r×r , b h ∈ R r , W rx ∈ R r×4+q , W rh ∈ R r×r , b r ∈ R r , W zx ∈ R r×4+q , W zh R r×r and b z ∈ R r are the parameters of the embedding function, r is the size of the embedding, σ is the ReLU activation function, and h event 0 = 0. In the experiments of Sec. 6, only the 1, 2 or 5 hardest jets are considered in the sequence j = 1, . . . , M , as ordered by ascending values of p T .

C Implementation details
While tree-structured networks appear to be a principled choice in natural language processing, they often have been overlooked in favor of sequence-based networks on the account of their technical incompatibility with batch computation [50]. Because tree-structured networks use a different topology for each example, batching is indeed often impossible in standard implementations, which prevents them from being trained efficiently on large datasets. For this reason, the recursive jet embedding we introduced in Sec. 3.1 would undergo the same technical issues if not implemented with caution.
In our implementation, we achieve batch computation by noticing that activations from a same level in a recursive binary tree can be performed all at once, provided all necessary computations from the deeper levels have already been performed. This principle extends to the synchronized computation across multiple trees, which enables the batch computation of our jet embeddings across many events. More specifically, the computation of jet embeddings is preceded by a traversal of the recursion trees for all jets in the batch, and whose purpose is to unroll and regroup computations by their level of recursion. Embeddings are then reconstructed level-wise in batch, in a bottom-up fashion, starting from the deepest level of recursion across all trees.
Finally, learning is carried out through gradients obtained by the automatic differentiation of the full model chain on a batch of events (i.e., the recursive computation of jet embeddings, the sequence-based recurrence to form the event embeddings, and the forward pass through the classifier). The implementation is written in native Python code and makes use of Autograd [59] for the easy derivation of the gradients over dynamic structures. Code is available at 2 under BSD license for further technical details.