Binary JUNIPR: an interpretable probabilistic model for discrimination

TL;DR

Binary JUNIPR enhances jet classification in particle physics by refining probabilistic models for discrimination, achieving state-of-the-art results while providing interpretable insights into jet differences.

Contribution

It introduces Binary JUNIPR, a refined probabilistic model optimized for classification tasks like quark/gluon discrimination and top-tagging, with improved performance and interpretability.

Findings

Achieves state-of-the-art quark/gluon discrimination.

Provides physical insights into jet classification.

Differentiates between gluon jets from different simulations.

Abstract

JUNIPR is an approach to unsupervised learning in particle physics that scaffolds a probabilistic model for jets around their representation as binary trees. Separate JUNIPR models can be learned for different event or jet types, then compared and explored for physical insight. The relative probabilities can also be used for discrimination. In this paper, we show how the training of the separate models can be refined in the context of classification to optimize discrimination power. We refer to this refined approach as Binary JUNIPR. Binary JUNIPR achieves state-of-the-art performance for quark/gluon discrimination and top-tagging. The trained models can then be analyzed to provide physical insight into how the classification is achieved. As examples, we explore differences between quark and gluon jets and between gluon jets generated with two different simulations.

Equations

$$P_{\mathcal{J}}(\text{jet})$$

$$\displaystyle\times P^{(n)}\big{(}\text{end}\big{|}k_{1}^{(n)}\!\ldots\,k_{n}^{(n)}\big{)}$$

$$P^{(n)}\big{(}\text{end}\big{|}k_{1}^{(n)}\ldots k_{n}^{(n)}\big{)}=P_{\text{end}}\big{(}\text{true}\big{|}h^{(n)}\big{)}$$

$$\displaystyle P^{(t)}\big{(}k^{(t+1)}_{1}\!\ldots\big{|}k^{(t)}_{1}\!\ldots\big{)}$$

$$\displaystyle\times P_{\text{mother}}\big{(}m^{(t)}\big{|}h^{(t)}\big{)}$$

$$\displaystyle\times P_{\text{branch}}\big{(}k_{d_{1}}^{(t+1)}k_{d_{2}}^{(t+1)}\big{|}k_{m}^{(t)}\,h^{(t)}\big{)}$$

$$\displaystyle\text{log likelihood}=\sum_{\text{jets}}\;\log P_{\mathcal{J}}\big{(}\text{jet}\big{)}$$

$$P(a|\text{jet})=\frac{P(\text{jet}|a)P(a)}{P(\text{jet})}$$

$$P(a|\text{jet})=\frac{P(\text{jet}|a)P(a)}{P(\text{jet}|a)P(a)+P(\text{jet}|b)P(b)}$$

$$\mathcal{L}$$

$$+\sum _{\text{b-jets}}\log \frac{P_{\mathcal{J}}(\text{jet}|b)f_b}{P_{\mathcal{J}}(\text{jet}|a)f_a+P_{\mathcal{J}}(\text{jet}|b)f_b}$$

Full text

Modern machine learning has already made impressive contributions to particle physics. Convolutional deOliveira:2015xxd ; Baldi:2016fql ; Komiske:2016rsd ; Komiske:2017ubm ; Kasieczka:2017nvn ; Macaluso:2018tck , recurrent and recursive networks Butter:2017cot ; Louppe:2017ipp ; Cheng:2017rdo ; Egan:2017ojy ; Fraser:2018ieu , autoencoders Collins:2018epr ; Farina:2018fyg ; Blance:2019ibf ; Roy:2019jae , adversarial networks Paganini:2017dwg ; deOliveira:2017pjk ; Paganini:2017hrr and more have been shown effective in applications including quark/gluon jet discrimination, top-tagging and pileup removal. A key question that is beginning to be addressed is: what is the optimal representation of the information in an event? Is it through analogy with images deOliveira:2015xxd ; Baldi:2016fql , natural-language processing Louppe:2017ipp ; Fraser:2018ieu , or set theory Komiske:2018cqr ; Qu:2019gqs ? In many of these approaches, there is a competition between effectiveness in some task (e.g. pileup removal, jet classification) and interpretability of the neural network. An approach to machine learning for particle physics called Junipr JUNIPR builds a separate network for each jet type using a physical representation of the information in the jet: the jet clustering tree. In JUNIPR a method for construction and training of such a network was introduced. In this paper, we show how the Junipr framework can be used in discrimination tasks, achieving state-of-the art classification power while maintaining physical interpretability.

Junipr begins by taking each jet in some sample and clustering it into a binary tree according to some deterministic algorithm. See Fig. 2 below for an example of such a tree. The algorithm can be physically motivated (like the $k_{T}$ Ellis:1993tq or Cambridge/Aachen Dokshitzer:1997in algorithms) but does not have to be. In such a tree, the momenta of each mother branch is the sum of the momenta of her daughters. We denote the momenta of the particles in the jet by $\{p_{1}\ldots p_{n}\}$ and the momenta in the clustering tree by $\{k_{1}^{(t)}\!\ldots\,k_{t}^{(t)}\}$ at branching step $t$. To be concrete, at $t=1$ we have $k_{1}^{(1)}=p_{1}+\cdots+p_{n}$, at $t=n$ we have $\{k_{1}^{(n)}\!\ldots\,k_{n}^{(n)}\}=\{p_{1}\ldots p_{n}\}$, and at each branching in between, $\{k_{1}^{(t)}\!\ldots\,k_{t}^{(t)}\}\to\{k_{1}^{(t+1)}\!\ldots\,k_{t+1}^{(t+1)}\}$ involves a single $1\to 2$ momentum splitting. Junipr learns to compute the probability $P_{\mathcal{J}}$(jet) of the jet, meaning the probability that the corresponding set of final state momenta $\{p_{1}\ldots p_{n}\}$ would be found in the given sample. This probability can be factorized as a product over branching steps in the clustering tree:

[TABLE]

To learn these probability distributions, Junipr introduces a quantity $h^{(t)}$ as a representation of $\{k^{(t)}_{1}\!\ldots\,k^{(t)}_{t}\}$, i.e. the “state” of the jet at branching step $t$. Junipr learns to compute $h^{(t)}$ in training. In machine-learning language, $h^{(t)}$ is the autoregressive latent variable, which in our implementations is taken to be the latent state of a recurrent neural network. Then we can write, e.g.,

[TABLE]

where $P_{\text{end}}(\text{true}|h^{(n)})$ is the binary probability that the clustering tree ends at branching step $n$.

Junipr further factorizes the branching-probabilities of Eq. (1) into more intuitive probability distributions

[TABLE]

Here, $P_{\text{end}}(\text{false}|h^{(t)})$ is the binary probability that the clustering tree does not end at branching step $t$, $P_{\text{mother}}(m^{(t)}|h^{(t)})$ is the discrete probability that tree-momentum $k_{m}^{(t)}$ will participate in the $1\to 2$ branching at step $t$, and P_{\text{branch}}\big{(}k_{d_{1}}^{(t+1)}k_{d_{2}}^{(t+1)}\big{|}k_{m}^{(t)}\,h^{(t)}\big{)} is the distribution over daughters of the branching. Structuring the probabilistic model in terms of the product of these parts is essential to the interpretability of the model’s output, as each part has separate physical meaning.

The latent state $h^{(t)}$ has access to the global content of the jet at branching step $t$, i.e. to all the momenta $\{k_{1}^{(t)}\!\ldots\,k_{t}^{(t)}\}$. The factorization over branching steps is powerful, and useful to the extent that the $1\to 2$ branching dynamics encoded in P_{\text{branch}}\big{(}k_{d_{1}}^{(t+1)}k_{d_{2}}^{(t+1)}\big{|}k_{m}^{(t)}\,h^{(t)}\big{)} are local, depending only weakly on $h^{(t)}$. Even if there were no evidence for this factorization in the training data (as was explored with “printer jets” in JUNIPR ), Junipr would still learn the probability distributions, but physical interpretability would be lost.

In JUNIPR , Junipr was trained to model jet dynamics via unsupervised learning. In that approach, the probabilistic model is learned by maximising the log likelihood of $P_{\mathcal{J}}$ over the training data:

[TABLE]

where the sum is over jets $\{p_{1}\ldots p_{n}\}$ in the training set. We call this the “unary objective function”. Despite being unsupervised, this approach can be used to discriminate between two jet types, say $a$ and $b$. To accomplish this, one trains two separate Junipr models: $P_{\mathcal{J}}(\text{jet}|a)$ on a data set containing predominantly type-$a$ jets and $P_{\mathcal{J}}(\text{jet}|b)$ on predominantly type-$b$ jets. Discrimination between $a$ and $b$ is then achieved by thresholding the likelihood ratio $P_{\mathcal{J}}(\text{jet}|a)/P_{\mathcal{J}}(\text{jet}|b)$.

While discrimination by likelihood ratio is theoretically optimal in the perfect-model limit, it has been shown that deep neural networks classify out-of-distribution data poorly nguyen2015deep ; bendale2016towards . That is, e.g., the $P_{\mathcal{J}}(\text{jet}|a)$ model is not expected to behave well on type-$b$ jets. It is thus advantageous in practice to refine the training for discrimination. By training directly for discrimination, Junipr can also focus model capacity on learning the often-subtle differences between type-$a$ and type-$b$ jets. In fact, Junipr’s probabilistic nature makes supervised discrimination-learning very straightforward. Assuming a mixed sample of both jet types, the probability that a given jet drawn at random belongs to class $a$ is, through Bayes’ theorem, given by

[TABLE]

For binary discrimination, $P(a|\text{jet})+P(b|\text{jet})=1$, so

[TABLE]

Here $P(a)$ and $P(b)$ are simply the composition fractions $f_{a}$ and $f_{b}$ of the mixed sample, while $P(\text{jet}|a)$ and $P(\text{jet}|b)$ can be computed using two separate Junipr networks as laid out in the paragraphs above. This leads directly to the binary cross-entropy objective function one should use to train Junipr for discrimination:

[TABLE]

where the sums extend over type-$a$ and type-$b$ jets in the training data, respectively. We call training with this objective function “Binary Junipr”. Note that Binary Junipr still learns the probabilities for type-$a$ and type-$b$ jets and still trains the same neural-network functions; however, it uses a more effective objective function for discrimination applications. We also note that training can easily be generalized to multiclass classification.

As a test of the advantage that the binary objective function provides over its unary counterpart, we applied Binary Junipr to the discrimination of quark- and gluon-jets. We used a mixed sample of $10^{6}$ Pythia quark-jets and $10^{6}$ Pythia gluon-jets from the data set at energyflow.network Komiske:2018cqr ; komiske_patrick_2019_3164691 . We set aside $10^{5}$ jets of each type into a test set, $10^{5}$ for validation, and used the remaining 80% of the jets for training. For the Junipr models, $P_{\mathcal{J}}(\text{jet}|\text{quark})$ and $P_{\mathcal{J}}(\text{jet}|\text{gluon})$, we used an LSTM of dimension 30 to model $h^{(t)}$ and separate feed-forward networks, each with a single hidden layer of dimension 10, to model $P_{\text{end}}$, $P_{\text{mother}}$, and $P_{\text{branch}}$.111The Binary Junipr architecture is available at github.com/andersjohanandreassen/JUNIPR with example code.

We began by pre-training the two Junipr models using the original unary objective function of Eq. (4). We followed the same training schedule as in JUNIPR , but scaled down the number of epochs by a factor of 5 because this data set is larger than the one used there. Pre-training took about 5 hours on a 16-core CPU server for each model. After pre-training, we optimized the binary objective function of Eq. (7) using Adam with standard settings kingma2014adam and the following batch-size schedule:

[TABLE]

This segment of training took 12 hours on a 16-core CPU server. Binary Junipr parameters were decided upon by evaluating the AUC (area under the ROC curve) on the validation set 10 times per epoch and choosing the model that achieved the maximal AUC during the final 10 training epochs. Note that different hyperparameters might be appropriate for different applications.

In Fig. 1 we show the quark-versus-gluon Significance Improvement Curve Gallicchio:2010dq , (SIC), $\frac{\varepsilon_{S}}{\sqrt{\varepsilon_{B}}}$, achieved by Binary Junipr and compare it to recent results with previous state-of-the-art discriminants: a CNN approach based on jet images Komiske:2016rsd (with architecture from Komiske:2018cqr ) and Particle Flow Networks Komiske:2018cqr . One can see that Binary Junipr offers a small-but-significant advantage. Quantitatively, Binary Junipr achieves an AUC of $0.8986\pm 0.0004$, as compared to $0.8911\pm 0.0008$ for Particle Flow Networks, and $0.8799\pm 0.0008$ for the CNN. (Each reported number is the mean and semi-interquartile range over 10 trainings.) Unary Junipr, trained with Eq. (4), performs significantly worse than the other methods, achieving an AUC of $0.6968\pm 0.0008$. This demonstrates the importance of training Junipr with the binary objective function of Eq. (7) for classification.

As a second experiment, we trained and tested Binary Junipr for boosted top-jet identification. We used the same architecture and training schedule that were optimized for quark-versus-gluon discrimination. In doing so, we obtain a sense of the performance one might expect from Binary Junipr without specialized hyperparameter tuning. The training, validation, and test data for this experiment are taken from Butter:2017cot . We found that untuned Binary Junipr comes close to state-of-the-art top discrimination. Specifically, Junipr achieves an AUC of $0.9810\pm 0.0002$ as compared to $0.9819\pm 0.0001$ attained using Particle Flow Networks Komiske:2018cqr , and 0.9848 reported for ParticleNet Qu:2019gqs ; all significantly outperform traditional boosted top-tagging methods Kaplan:2008ie . For a recent overview of machine learning in top tagging, see Kasieczka:2019dbj .

Next we discuss the interpretability of Junipr models. As discussed below Eq. (3), each component of Junipr’s output has a well-defined physical meaning. Moreover, the output is structured along a physically-motivated binary tree, defined by clustering the momenta in a jet. One can thus decompose Junipr’s prediction, say $P_{\mathcal{J}}(\text{top}\,|\,\text{jet})$ as in Eq. (6), visually along the clustering tree. In Fig. 2, we show the clustering tree for an easily classifiable top jet drawn from the mixed top/QCD test set. In the figure, we label the $t$-th node with $P^{(t)}(\text{top}\,|\,\text{jet})$, i.e. the probability that the jet is top-type, given only the information present at branching step $t$; this is computed with Binary Junipr by substituting Eq. (3) into Eq. (6). One can see, for example, that the 3-prong structure characteristic of $t\to W^{+}b\to u\,\bar{d}\,b$ contributes to large $P_{\mathcal{J}}(\text{top}\,|\,\text{jet})$. Quantitatively, this results in the two hard branchings, with $P^{(t)}(\text{top}\,|\,\text{jet})=0.72$ and $0.71$, dominating the prediction. By analyzing such trees, one can develop intuition for which branchings are most decisive in classifying different types of jets.

To be concrete, let us return to the Binary Junipr model used to create Fig. 1, which learned to discriminate quark and gluon jets from Pythia. Much is already known about the difference between quark and gluon jets: gluon jets are known to be bigger, with larger multiplicity and larger shape parameters such as mass and width Gallicchio:2012ez ; Gallicchio:2011xq . Although many methods exist for quark/gluon discrimination, including other machine-learning approaches Komiske:2016rsd ; Fraser:2018ieu , it is not clear how well these methods will work on actual data. In particular, it is known that real gluon jets are more similar to real quark jets than Pythia leads us to believe Aad:2014gea . In particular, it is the modeling of gluon jets that seems most inaccurate. An alternative generator, Herwig, produces and gluon jets that are more similar to its quark jets Gras:2017jty . Thus, we also considered a secondary challenge: determine how Pythia and Herwig gluon jets differ. To explore their differences, we trained a second Binary Junipr model to discriminate Pythia 8.226 and Herwig 7.1.4 gluon jets using $10^{6}$ samples of each from komiske_patrick_2019_3164691 ; pathak_aditya_2019_3066475 .

Fig. 3 shows another visualization, complementary to Fig. 2, of exploring how Junipr discriminates. The top row of Fig. 3 shows how Junipr separates Pythia quarks from Pythia gluons, and the bottom row shows how Junipr separates Pythia gluons from Herwig gluons. In the middle column, the overall probability that Junipr uses for discrimination is decomposed into branching steps $t$, averaged over all jets of the given class. From this, we see that near $t=20$ – 50 there is roughly 3 times the quark/gluon discrimination power per branching step as for $t=1$ – 10. This echoes the well-known fact that multiplicity allows one to separate quark and gluon jets better than perturbatively calculable observables sensitive to only the first few splittings Gallicchio:2012ez . The lower-middle plot shows that differences in Pythia and Herwig gluon jets are more uniformly spread over branching steps.

Not only can Junipr break discrimination power down into branching steps; Junipr can further decompose classification probability into components at each branching. These components are displayed in the right column of Fig. 3; there are discrete components, such as whether branchings should end, as well as the energy $z$ and angles $\theta,\phi,\delta$ of the branching itself. While multiplicity ($P_{\text{end}}$) is the main driver of performance for quark/gluon discrimination, the angle $\theta$ also contributes significantly over a wide range of branchings, echoing the importance of jet width in this context. For the Pythia/Herwig task, both the angle $\theta$ and energy fraction $z$ play a significant role in discrimination on early branchings, and multiplicity becomes important on later branchings.

It is interesting that a significant fraction of the difference between Pythia and Herwig results from the way energy and angles are distributed early on in the clustering trees. Early branchings are controlled primarily by perturbative elements of the simulated parton showers. This suggests that a substantial portion of the difference between Pythia and Herwig gluon jets may be driven by the parton-shower implementations, rather than exclusively by the modelling of non-perturbative effects. To gain further insight into the importance of non-perturbative effects like hadronization in discrimination, Junipr could be upgraded to include quantum numbers of final state particles — a straightforward next step.

In JUNIPR , Junipr was introduced as a new framework for unsupervised machine learning in particle physics that prioritizes interpretability. Given a jet, i.e. a set of momenta, Junipr learns to compute the probability of that jet, i.e. how consistent the distribution of momenta is with the training data. In this paper, we used the same probabilistic framework as in JUNIPR , but we augmented the training to learn subtle differences between two samples, an enhancement we call Binary Junipr. We demonstrated both its effectiveness and interpretability, using quark/gluon jets, boosted top jets, and Monte-Carlo-generator-dependence as examples. It is satisfying that demanding interpretability does not lead to a loss in effectiveness: Binary Junipr discriminates at levels competitive with the best machine-learning methods available.

While these case studies were all simulation-based, there is a straightforward path to repeating these exercises on collider data. Although real data does not come with truth labels, there are established methods for working with mixed samples Komiske:2018oaa ; Metodiev:2017vrx which can be adapted to Junipr without much modification. Then one could use a data/simulation Binary Junipr model to understand deficiencies in simulations. One could also use insights derived from Binary Junipr trees to judge whether predictions should be trusted experimentally (was information below experimental resolution deemed important?) or to design new calculable observables (sensitive to previously-overlooked decisive branchings). Having interpretable methods opens the door to whole new approaches to understanding data from particle colliders.

We thank P. Komiske and A. Parthak for assistance with the samples. This work is supported in part by the U.S. Department of Energy under contract DE-SC0013607. Computations were performed on the CORI supercomputing resources at NERSC, Harvard’s Odyssey cluster, and the Faculty Platform for machine learning.