Interpretable Deep Learning for Two-Prong Jet Classification with Jet Spectra

TL;DR

This paper introduces an interpretable deep learning model for jet classification based on jet spectra, achieving comparable performance to CNNs while maintaining transparency through a truncated functional Taylor series.

Contribution

The authors develop a novel interpretable neural network using jet spectra and a truncated Taylor series, enhancing interpretability without sacrificing accuracy.

Findings

Performance comparable to CNNs on jet classification tasks

Uses jet spectrum $S_{2}(R)$ for interpretability

Simpler architecture than traditional CNN classifiers

Abstract

Classification of jets with deep learning has gained significant attention in recent times. However, the performance of deep neural networks is often achieved at the cost of interpretability. Here we propose an interpretable network trained on the jet spectrum $S_{2}(R)$ which is a two-point correlation function of the jet constituents. The spectrum can be derived from a functional Taylor series of an arbitrary jet classifier function of energy flows. An interpretable network can be obtained by truncating the series. The intermediate feature of the network is an infrared and collinear safe C-correlator which allows us to estimate the importance of a $S_{2}(R)$ deposit at an angular scale R in the classification. The performance of the architecture is comparable to that of a convolutional neural network (CNN) trained on jet images, although the number of inputs and complexity of architecture is significantly simpler than the CNN classifier. We consider two examples: one is the classification of two-prong jets which differ in color charge of the mother particle, and the other is a comparison between Pythia 8 and Herwig 7 generated jets.

Figures

Equations

$$S_2(R)$$

$$P_T(\vec{R})$$

$$S_2(R;\Delta R)$$

$$I_A(x)$$

$$\int d\vec{R}\delta P_T(\vec{R})w(\vec{R})$$

$$\int dR\delta S_2(R)w(R)=2\sum _{j\in \mathbf{J}}\left[\delta p_{T,i_0}+p_{T,i_1}(\delta {\vec{R}}_{i_1}\cdot {\nabla }_{\vec{R}})+p_{T,i_2}(\delta {\vec{R}}_{i_2}\cdot {\nabla }_{\vec{R}})+\cdots \right]p_{T,j}w(R_{i_{0}j}).$$

$$\int d{\vec{R}}_1\cdots d{\vec{R}}_n{P}_T({\vec{R}}_1)\cdots P_T({\vec{R}}_n)w({\vec{R}}_1,\cdots ,{\vec{R}}_n),$$

$$\int dR{S}_2(R)w(R)=\int d{\vec{R}}_{1}d{\vec{R}}_2{P}_T({\vec{R}}_1)P_T({\vec{R}}_2)w(R_{12}).$$

$${\mathbf{J}}_{\mathrm{trim}}=\bigcup _{\begin{array}{c}a\\ \frac{p_{T,{\mathbf{J}}_a}}{p_{T,\mathbf{J}}}\ge f_{\mathrm{trim}}\end{array}}{\mathbf{J}}_a.$$

$$S_{2,\mathrm{trim}}(R)$$

$$P_{T,\mathrm{trim}}(\vec{R})$$

$$S_{2,\mathrm{trim}}(R;\Delta R)$$

$$P_T(\vec{R})=\sum _a{P}_{T,a}(\vec{R}),P_{T,a}(\vec{R})=\sum _{i\in {\mathbf{J}}_a}{p}_{T,i}\delta (\vec{R}-{\vec{R}}_i).$$

$$S_2(R;\Delta R)$$

$$S_{2,\mathrm{soft}}(R;\Delta R)=S_2(R;\Delta R)-S_{2,\mathrm{trim}}(R;\Delta R).$$

$$S_{2,\mathrm{trim}}(R;\Delta R)$$

$$S_{2,\mathrm{soft}}(R;\Delta R)$$

$$h_i$$

$$\hat{y}$$

$${\phi }_{\mathrm{softmax},k}(\vec{z})=\frac{e^{z_k}}{{\sum }_{k=1}^K{e}^{z_k}}.$$

$$P_{T,a}(\vec{R})=\{\begin{array}{ll}{P}_{T,\mathrm{trim}}(\vec{R})& a=1,\\ P_T(\vec{R})-P_{T,\mathrm{trim}}(\vec{R})& a=2.\end{array}$$

$$h_i$$

$$\int d\vec{R}{P}_{T,1}(\vec{R})\simeq p_{T,\mathbf{J},\mathrm{trim}},\int d\vec{R}{P}_{T,2}(\vec{R})\simeq p_{T,\mathbf{J}}-p_{T,\mathbf{J},\mathrm{trim}}.$$

$$h_i$$

$$S_{2,ab}(R)$$

$$S_2(R)=\sum _{a,b}{S}_{2,ab}(R),S_{2,\mathrm{trim}}(R)=S_{2,11}(R).$$

$$h_i={\Psi }_i[S_{2,A};{\vec{x}}_{\mathrm{kin}}],$$

$$h_i$$

$$h_i=\frac{1}{2}\int dR{S}_{2,A}(R)w_{i,A}^{(2)}(R;{\vec{x}}_{\mathrm{kin}}).$$

$${\mathrm{EFP}}_{2,ab}^{(n)}$$

$${\mathrm{EFP}}_{2,ab}^{(n-1)}$$

$$D_2^{\beta }$$

$$e_2^{\beta }$$

$$e_3^{\beta }$$

Full text

††preprint: KEK-TH-2117††arxiv: 1904.02092

1 Introduction

Deep learning is gaining significant interest recently in the field of collider data analysis. One of the primary motivations is to extract the maximum information from the complex collision events. The deep learning in collider physics takes advantage of a large influx of data from experiments, more precise theoretical predictions, significant improvement in computing power, and ongoing progress in the field of machine learning itself. Such techniques offer advances in areas ranging from event selection to particle identification.

The large center-of-mass energy at the Large Hadron Collider (LHC) enables the production of boosted particles whose decay products are highly collimated. These collimated objects are reconstructed as a jet, and it is often misidentified as a QCD jet originated from light quarks or gluons. Many jet substructure techniques using the information of subjets Butterworth:2008iy ; Thaler:2008ju ; Kaplan:2008ie ; Plehn:2009rk ; Plehn:2010st ; Soper:2011cr ; Soper:2012pb ; Dasgupta:2013ihk ; Soper:2014rya ; Larkoski:2014wba and the distribution of jet constituents Gallicchio:2010sw ; Thaler:2010tr ; Gallicchio:2011xq ; Chien:2013kca ; Larkoski:2013eya ; Larkoski:2014gra ; Moult:2016cvt have been developed in order to improve the sensitivity of tagging and to classify these boosted particle jets. The deep learning methods Almeida:2015jua ; deOliveira:2015xxd ; Komiske:2016rsd ; Kasieczka:2017nvn ; Louppe:2017ipp ; Komiske:2017ubm ; Butter:2017cot ; Cheng:2017rdo ; Andreassen:2018apy ; Choi:2018dag ; Lim:2018toa ; Dreyer:2018nbf ; Lin:2018cin ; Komiske:2018cqr ; Martinez:2018fwc ; Kasieczka:2018lwf ; Qu:2019gqs have provided useful insight into the internal structure of the jets and, thereby, shown better performances than those jet substructure techniques.111For a review on the recent theoretical and machine learning developments in jet substructure techniques at the LHC, we refer Larkoski:2017jix ; Asquith:2018igt . The flexibility of deep learning also enables us to solve problems beyond supervised classifications, such as weakly supervised learning Dery:2017fap ; Cohen:2017exh ; Metodiev:2017vrx , adversarial learning to suppress learning from unwanted information Louppe:2016ylz ; Shimmin:2017mfk , and unsupervised learning for finding anomalous signatures Chakraborty:2017mbz ; Hajer:2018kqm ; Heimel:2018mkt ; Farina:2018fyg ; Cerri:2018anq ; Roy:2019jae . The neural network can also be useful to new physics searches with deep learning at the LHC Roxlo:2018adx ; Brehmer:2018kdj ; Brehmer:2018eca ; Guo:2018hbv ; Collins:2018epr ; DAgnolo:2018cun ; DeSimone:2018efk ; Englert:2019xhk ; Collins:2019jip .

The output of a neural network is, in general, a highly non-linear function of the inputs. A neural network classifier often acts like a “black box.” One may consider architectures with post-hoc interpretability DBLP:journals/corr/Lipton16a , which allows us to extract information other than its prediction from the learned model after training. A simple strategy is using a predefined functional form to restrict the representation power of the neural network Komiske:2018cqr ; Datta:2019ndh . Then the network is interpreted in terms of the functional form. The aim of this paper is also to construct an interpretable neural network architecture that allows us not only to interpret the predictions of the network but also to visualize it in terms of trained weights connected to physical variables.

In Lim:2018toa , a multilayer perceptron (MLP) trained on two-point correlation functions $S_{2}$ and $S_{2,\mathrm{trim}}$ of angular scale $R$ was introduced. The $S_{2}(R)$ and $S_{2,\mathrm{trim}}(R)$ spectra are constructed from the constituents of a jet before and after the trimming Krohn:2009th respectively. The angular scale $R$ is an important parameter for describing the kinematics of a decaying particle and parton shower (PS); hence, these spectra efficiently encode the radiation pattern inside a jet. The MLP trained on these inputs learns relevant features for the classification among the Higgs boson jet (Higgs jet) and QCD jet.

In this paper, we connect the spectra to energy flow functionals $P_{T}(\vec{R})$ Tkachov:1995kk , i.e., we consider transverse energy of a jet constituent as particle-specific information at $\vec{R}$ in the $\eta-\phi$ plane 7974879 . The spectra are basis vectors of infrared and collinear (IRC) safe variables called bilinear $C$-correlators Tkachov:1995kk whose angular weighting function depends only on the relative distance between two constituents. Those correlators naturally appear in the functional Taylor series of a classifier of $P_{T}(\vec{R})$, and the MLP can be considered as a subseries of the Taylor series. We show that the performance of the MLP and neural networks trained on jet images Cogan:2014oua ; Almeida:2015jua ; deOliveira:2015xxd ; Kasieczka:2017nvn are comparable. This strongly suggests that $S_{2}$ and $S_{2,\mathrm{trim}}$ contain sufficient information for jet classification. Encouraged by this feature, we construct an interpretable architecture by truncating the series. Namely, $\int dR\,S_{2}(R)w(R;\vec{x}_{\mathrm{kin}})$ can be implemented in a classifier after proper discretization in $R$, where $\vec{x}_{\mathrm{kin}}$ is a set of kinematic variables of the jet and $w$ is a smooth function. By reading the weights $w(R;\vec{x}_{\mathrm{kin}})$, we could quantify important features for the given classification problem.

Jet substructure studies often suffer from systematic uncertainties of soft activities. The soft radiations generated by a Monte Carlo program are strongly model dependent. While this mismodeling could be corrected by using real data, it is certainly useful to use input variables with less systematic uncertainties. When hard substructures are important for solving the problem, we may use jet grooming techniques Butterworth:2008iy ; Krohn:2009th ; Ellis:2009su ; Ellis:2009me ; Larkoski:2014wba to remove the soft activity. Instead of throwing this soft activity away, we encode it in $S_{2,\mathrm{soft}}(R)$, which is $S_{2}(R)-S_{2,\mathrm{trim}}(R)$. Then, the inputs $S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$ include hard and soft substructure information, respectively. The interpretable architecture trained on $S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$ is able to quantify these features. We study two classification problems: one is a classification of two-prong jets to understand their hard substructures and color coherence, and the other is a comparison of Pythia 8 Sjostrand:2014zea and Herwig 7 Bellm:2015jjp ; Bahr:2008pv events to quantify the differences.

The paper is organized as follows. In section 2, we review $S_{2}$ and $S_{2,\mathrm{trim}}$ and show its relation to energy flow and $C$-correlators. We also show $S_{2}$ and $S_{2,\mathrm{trim}}$ distributions of typical Higgs jet and QCD jet. A hypothetical color octet scalar particle, sgluon, decaying to $b\bar{b}$ is considered to study the color connection in two-prong jets. In section 3, we first discuss the capability of $S_{2}$ and $S_{2,\mathrm{trim}}$ for the classification of two-prong jets and show the result of an MLP trained on those inputs. The results are then compared with that of a CNN trained on jet images. In section 4, we introduce a two-level architecture consists of a softmax classifier and an MLP trained on $S_{2}$ and $S_{2,\mathrm{trim}}$. The intermediate feature of this architecture is the bilinear $C$-correlator whose basis vectors are $S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$, and the MLP generates its components. We visualize and interpret the weights of the given classification problem. Finally, the summary and outlook are given in section 5.

2 Two-Point Correlation Spectrum and Two-Prong Jets

2.1 Jet Spectra

In Lim:2018toa , we introduced a two-point correlation spectral function $S_{2}(R)$ which maps a jet to a function of angular scale $R$,

[TABLE]

where $\mathbf{J}$ is a set of jet constituents, $\vec{R}_{i}=(\eta_{i},\phi_{i})$ is the position of the $i$-th jet constituent in the pseudorapidity-azimuth plane, $R_{ij}=\sqrt{(\eta_{i}-\eta_{j})^{2}+(\phi_{i}-\phi_{j})^{2}}$ is the angular distance between the two jet constituents $i$ and $j$, and $P_{T}(\vec{R})$ is an energy flow functional Tkachov:1995kk of $\mathbf{J}$. For practical purpose, $S_{2}(R)$ is discretized as below,

[TABLE]

where $I_{A}(R_{ij})$ is an indicator function of the angular distance $R_{ij}$ of the domain $A$,

[TABLE]

The spectral function $S_{2}(R;\Delta R)$ is, therefore, the sum of the product of $p_{T}$’s of the two jet constituents with an angular distance $R_{ij}$ lying between $R$ and $R+\Delta R$.

We obtain IRC safe quantities by multiplying smooth functions222 Continuous functions are sufficient for the convergence and IRC safety Tkachov:1995kk , but we further restrict $w$’s to smooth functions for perturbative calculations.

$w(\vec{R})$ and $P_{T}(\vec{R})$ (or $S_{2}(R)$), and integrating over $\vec{R}$. To understand the IRC safety of $P_{T}(\vec{R})$, let us consider splitting of a given constituent $i_{0}$ in $\mathbf{J}$ into two constituents, $i_{0}\rightarrow i_{1}i_{2}$. The inner product of $w(\vec{R})$ and the difference of the energy flow before and after the splitting, $\delta P_{T}(\vec{R})$, is given as follows,

[TABLE]

where $\delta p_{T,i_{0}}=p_{T,i_{1}}+p_{T,i_{2}}-p_{T,i_{0}}$, and $\delta\vec{R}_{i_{1}(i_{2})}=\vec{R}_{i_{1}(i_{2})}-\vec{R}_{i_{0}}$. The soft limit, where $i_{2}$ carries a small momentum, corresponds to $\delta p_{T,i_{0}},\,\delta\vec{R}_{i_{1}},\;p_{T,i_{2}}\rightarrow 0$, while $\delta p_{T,i_{0}},\,\delta\vec{R}_{i_{1}},\;\delta\vec{R}_{i_{2}}\rightarrow 0$ in the collinear limit. The integral vanishes in these limits, namely the energy flow after parton splitting converges weakly Tkachov:1995kk to the one before splitting.

The spectrum $S_{2}(R)$ inherits the same property. The inner product of the smooth function $w(R)$ and the difference of the spectrum, $\delta S_{2}(R)$, before and after the splitting $i_{0}\rightarrow i_{1}i_{2}$ is given as follows,

[TABLE]

Again, this integral vanishes in the IRC limits. Note that the binned spectrum $S_{2}(R;\Delta R)$ is not completely IRC safe because of the discontinuity of the indicator function at the bin boundaries. Nevertheless, when the domain is discretized into small sections $[R_{i},R_{i}+\Delta R_{i})$, the IRC unsafe terms cancel in the sum, $\sum_{i}S_{2}(R_{i};\Delta R_{i})\,w(R_{i})$, and it is approximately IRC safe up to binning errors.

The resulting IRC safe observables belong to $C$-correlators Tkachov:1995kk , which are multilinear forms of the energy flow. An $n$-linear $C$-correlator is expressed as follows,

[TABLE]

where $w$ is a continuous function of $\vec{R}_{1},\cdots,\vec{R}_{n}$. For example, an inner product of $P_{T}(\vec{R})$ and $w(\vec{R})$ is a linear $C$-correlator, and an inner product of $S_{2}(R)$ and $w(R)$ is a bilinear $C$-correlator with $w$ depending only on the relative distance $R_{12}$,

[TABLE]

Many well-known jet observables belong to the $C$-correlator, for example, a jet transverse momentum $p_{T,\mathbf{J}}$ is a linear $C$-correlator with $w(\vec{R}_{1})\approx 1$, a jet mass $m_{\mathbf{J}}$ is a bilinear $C$-correlator with $w(\vec{R}_{1},\vec{R}_{2})\approx R_{12}^{2}/2$.

The $S_{2}(R)$ spectra use all the jet constituents, but it is useful to separate the correlations of constituents of the hard subjets; we consider jet trimming for this purpose. We recluster the constituents of a jet of a radius parameter $R_{\mathbf{J}}$ to subjets with a smaller radius parameter $R_{\mathrm{trim}}$. A subjet $\mathbf{J}_{a}$ is discarded if $p_{T,\mathbf{J}_{a}}<f_{\mathrm{trim}}\,p_{T,\mathbf{J}}$, where $p_{T,\mathbf{J}}$ and $p_{T,\mathbf{J}_{a}}$ are the transverse momenta of the jet and $a$-th subjet respectively. The trimmed jet $\mathbf{J}_{\mathrm{trim}}$ is defined as a union of the remaining subjets,

[TABLE]

The jet trimming is beneficial because it does not introduce additional angular scale parameters other than $R_{\mathrm{trim}}$. The trimmed spectrum is then calculated using the constituents of the trimmed jet. We denote it as $S_{2,\mathrm{trim}}(R)$ and its binned version $S_{2,\mathrm{trim}}(R;\Delta R)$, which are defined as follows:

[TABLE]

where $P_{T,\mathrm{trim}}(\vec{R})$ is the energy flow of $\mathbf{J}_{\mathrm{trim}}$.

In the limit of the constituents of each subjet $\mathbf{J}_{a}$ are localized, the energy flow and the jet spectrum can be approximated in terms of the subjet momenta. The energy flow of such a jet is decomposed into a sum of energy flows of all the subjets,

[TABLE]

The energy flow of each subjet converges weakly to $p_{T,\mathbf{J}_{a}}\delta(\vec{R}-\vec{R}_{\mathbf{J}_{a}})$. The $S_{2}(R)$ spectrum can be approximated by the momenta of the subjets, i.e.,

[TABLE]

The jet trimming also introduces a $p_{T}$ scale hierarchy among the subjets, and so their pairwise contributions to $S_{2}(R;\Delta R)$ can be classified by the scale. We define a quantity $S_{2,\mathrm{soft}}(R;\Delta R)$ where

[TABLE]

In the r.h.s. of the above equation, the correlations among the constituents of the hard subjets are canceled, and we have

[TABLE]

The dominant contributions to $S_{2,\mathrm{soft}}(R;\Delta R)$ (i.e., the $\mathcal{O}\left[f_{\mathrm{trim}}\right]$ terms) come from the correlations between a constituent in $\mathbf{J}_{\mathrm{trim}}$ and a constituent in $\mathbf{J}-\mathbf{J}_{\mathrm{trim}}$. The subleading $\mathcal{O}\left[f_{\mathrm{trim}}^{2}\right]$ terms denote the correlations among the constituents in $\mathbf{J}-\mathbf{J}_{\mathrm{trim}}$.

2.2 Derivation of Classifiers based on Energy Flows and Jet Spectra

We discuss the relation between $S_{2}(R)$ and neural network classifiers trained on the energy flow $P_{T}(\vec{R})$. A general softmax classifier that solves $K$-class jet classification problem can be expressed as a functional $\hat{\Psi}_{i}$ which maps the energy flow to real numbers $h_{i}$, i.e.,

[TABLE]

where $w^{\mathrm{(out)}}_{ki}$ and $b^{\mathrm{(out)}}_{k}$ are the weights and biases of the output layer, and $\hat{y}$ is the prediction of the classifier. Here the $\varphi_{\mathrm{softmax}}$ is the softmax function whose $k$-th component is expressed as follows,

[TABLE]

Many jet classifiers can be expressed in the form of eq. (17). For example, in the cut-based analysis, $\hat{\Psi}_{i}$ is a jet substructure variable, such as a ratio of $n$-subjettiness Thaler:2010tr , a ratio of energy correlation functions Larkoski:2014gra ; Moult:2016cvt , etc. The deep neural network classifiers, such as artificial neural network tagger Almeida:2015jua , convolutional neural network using pixelated jet images deOliveira:2015xxd , energy flow network Komiske:2018cqr , etc., are also described by eq. (17). The neural networks that are introduced in section 3 and section 4 also belong to this category.

The jet spectra $S_{2}$ and $S_{2,\mathrm{trim}}$ can be derived from eq. (17) using a functional Taylor expansion. The energy flow is decomposed by trimming as follows,

[TABLE]

One can express $\hat{\Psi}_{i}[P_{T,a}]$ as a functional series at a reference point $P_{T,a}(\vec{R})=0$,

[TABLE]

where $w^{(n)}_{i,a_{1}\cdots a_{n}}(\vec{R}_{1},\cdots,\vec{R}_{n})$ is the coefficient of $n$-th correlation function. The first order coefficient $w_{i,a}^{(1)}$ can be chosen as a constant if we are not interested in features depending on reference vectors, for example, jet axes, beam directions, etc. The linear term in $P_{T}(\vec{R})$ of eq. (21) is related to the jet momentum $p_{T,\mathbf{J}}$ and trimmed jet momentum $p_{T,\mathbf{J},\mathrm{trim}}$ as follows,

[TABLE]

The second order coefficient $w^{(2)}_{i,ab}$, the first non-trivial term of the series expansion, is a function of the relative distance of $\vec{R}_{1}$ and $\vec{R}_{2}$. The basis vectors of $w^{(2)}_{i,ab}$ are two-point correlation functions $S_{2,ab}(R)$,

[TABLE]

The spectra $S_{2}$ and $S_{2,\mathrm{trim}}$ are expressed in terms of $S_{2,ab}$ as follows,

[TABLE]

Instead of the energy flows, we consider a classifier of $S_{2,A}$ ($A$ = $\mathrm{trim}$, $\mathrm{soft}$),

[TABLE]

where $\vec{x}_{\mathrm{kin}}$ is a set of additional inputs to the classifier based on the kinematics of the jet. Similar to eq. (23), we expand eq. (26) around $S_{2,A}(R)=0$ as

[TABLE]

where $w_{i,A_{1}\cdots A_{\frac{n}{2}}}^{(n)}$ is the weight function corresponding to $w_{i}^{(n)}$ in eq. (26). One may further truncate the series to get a linear form,

[TABLE]

The above-mentioned linear setup has an advantage on the interpretability and visualization of the network predictions; we discuss more on this network in section 4.

2.3 Relation between Two-Point Correlation Spectra and Energy Flow Polynomials

Both the two-point correlation spectra and the energy flow polynomials Komiske:2017aww with two vertices span the set of bilinear $C$-correlators; therefore, there is a transformation rule between them. We first extend the definition of the energy flow polynomials to compare them to $S_{2,ab}$. Since $S_{2,ab}$ is a multivariate function of energy flows, we introduce a multivariate energy flow polynomial with two labeled vertices,

[TABLE]

This expression suggests that $R^{n}$ can be considered as an angular weighting function $w_{i,ab}^{(2)}(R)$ in eq. (23).

The resulting transformation from $S_{2,ab}(R)$ to ${\mathrm{EFP}}_{2,ab}^{(n)}$ is the Mellin transformation,

[TABLE]

The integral on the right-hand side is finite because $S_{2,ab}(R)$ vanishes on $R\gg 2R_{\mathbf{J}}$. The inverse transform is also well-defined if we allow the exponent $n$ in the angular weighting function of ${\mathrm{EFP}}_{2,ab}^{(n)}$ to be a complex number.

2.4 Spectra of Two-Prong Jets

The $S_{2}(R;\Delta R)$ spectrum is useful to identify the substructures of the jet and also to characterize the jet. Typically, $S_{2}(R;\Delta R)$ of QCD jet has a peak at $R=0$ with a long tail towards large $R$. The peak originates from the autocorrelation term $\sum_{i}p^{2}_{T,i}$ in eq. (3). On the other hand, if a jet originates from a Higgs boson decaying into $b\bar{b}$, the $b$-partons create two isolated cores inside the jet. The spectrum of the Higgs jet has a peak at the angular scale equal to the angle between the two clusters. In addition, $S_{2}(R;\Delta R)$ encodes the fragmentation pattern of $b$-partons.

At the LHC, boosted heavy objects such as top quark, gauge bosons and Higgs boson decaying into quarks can be studied by identifying jet substructures. Usually, these substructures are characterized by parameters such as $D_{2}$ defined as,

[TABLE]

where $\beta$ is the angular exponent. If a jet has a two-prong substructure, $D_{2}$ is much less than one. The jet spectrum $S_{2}(R;\Delta R)$ contains more information than $D_{2}$, and therefore, the analysis with $S_{2}(R;\Delta R)$ goes beyond the one using $D_{2}$. It was shown that a neural network trained on $S_{2}(R;\Delta R)$ distinguishes Higgs jet from QCD jet better than the one trained on $D_{2}$ Lim:2018toa .

To study the fragmentation pattern of the $b$-partons and their color connection to the mother particle, we introduce a color-octet scalar, sgluon ($\sigma$). We assume that the Higgs boson ($h$) and $\sigma$ decay into $b\bar{b}$ through the interaction,

[TABLE]

The Higgs boson is a color singlet particle, and the decay $h\rightarrow b\bar{b}$ is isolated in color flows. Therefore, $S_{2}(R;\Delta R)$ beyond the angle between the b-partons, i.e., $R_{b\bar{b}}$, is suppressed due to the color coherence. No such constraint on the angular scale exists for sgluon and QCD jets. Meanwhile, the Higgs jet and sgluon jet have the same two-prong substructure, unlike the QCD jet, as both are originating from a particle decaying into $b\bar{b}$ final states.

To study the spectra of two-prong jets, we simulate events as follows. We use Madgraph5 2.6.1 Alwall:2014hca to generate the events of $pp\rightarrow Zh$, $pp\rightarrow Z\sigma$, and $pp\rightarrow Zj$ processes with a collision energy of 13 TeV and the $Z$ boson decaying to a pair of neutrinos. These events are then passed to Pythia 8.226 Sjostrand:2014zea for the parton shower and hadronizations. To study the impact of the parton shower and hadronization schemes, we also pass those parton level events to Herwig 7.1.3 Bellm:2015jjp ; Bahr:2008pv . A color octet scalar UFO model Degrande:2014sta ; NLOModels generated by Feynrule 2.0 Alloul:2013bka ; Degrande:2011ua is used to simulate $pp\rightarrow Z\sigma$ process. The masses and widths of Higgs boson and sgluon are $m_{h}=m_{\sigma}=$ 125 GeV and $\Gamma_{\sigma}=\Gamma_{h}=6.4$ MeV. The detector response is simulated by Delphes 3.4.1 deFavereau:2013fsa with the default ATLAS detector configuration. We use FastJet 3.3.0 Cacciari:2011ma ; Cacciari:2005hq to reconstruct jets from the calorimeter towers using anti-$k_{T}$ algorithm Cacciari:2008gp with the radius parameter $R_{\mathbf{J}}=1.0$. For jet trimming, we use $R_{\mathrm{trim}}=0.2$ and $f_{\mathrm{trim}}=0.05$. We select the events with the leading jet transverse momentum $p_{T,\mathbf{J}}\in[300,400]$ GeV and its mass $m_{\mathbf{J}}\in[100,150]$ GeV. For Higgs jet and sgluon jet, we additionally require that the two $b$-partons originating from their decay are located within $R_{\mathbf{J}}$ from the jet axis. More details on our simulations are described in appendix A.

In figure 1, we show the pixelated jet image (left panel) and $S_{2}$ and $S_{2,\mathrm{trim}}$ spectra (right panel) of a QCD jet. There are high energy deposits in the jet image near the jet center along with a wide spray of soft activity. It also has a moderate amount of radiation at $(-0.4,0.0)$. As a result, $S_{2}(R;\Delta R)$ spectra has a long tail starting from $R=0.4$. The jet trimming eliminates a significant amount of soft particles and, therefore, the tail does not appear in $S_{2,\mathrm{trim}}(R;\Delta R)$. The remaining cross-correlations contributing to $S_{2,\mathrm{trim}}(R;\Delta R)$ are the ones between high and moderate energy deposits. Most of the energy deposits are concentrated at the center, and the peak intensity at $R=0.4$ is much lower than the intensity from autocorrelations at $R=0$.

In figure 2, we show $S_{2}(R;\Delta R)$ and $S_{2,\mathrm{trim}}(R;\Delta R)$ distributions of a Higgs jet. For this particular event, $S_{2}(R;\Delta R)$ distribution is similar to $S_{2,\mathrm{trim}}(R;\Delta R)$ distribution, and their difference $S_{2,\mathrm{soft}}(R;\Delta R)$ is hard to be seen. No significant activity has been observed beyond the peak at $R\sim 0.8$, mostly because the Higgs jet is very compact compared to the QCD jet. Correspondingly, there are two prominent subjets in the jet image, while most of the cells have no jets.

Finally, we show the $S_{2}(R;\Delta R)$ and $S_{2,\mathrm{trim}}(R;\Delta R)$ distributions of a sgluon jet in figure 3. The $S_{2}(R;\Delta R)$ distribution has a large peak at $R=0.6$ which is as significant as the one at $R=0$. This spectrum is qualitatively similar to the Higgs jet in figure 2. However, the $S_{2}(R;\Delta R)$ spectrum has a long tail beyond $R_{\mathbf{J}}$ as compared with that of a Higgs jet. The tail disappears after jet trimming, like the QCD jet in figure 1, that makes the $S_{2,\mathrm{trim}}(R;\Delta R)$ distribution more compact. From figure 1-3, we observe that $S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$ include useful complementary information.

In Lim:2018toa , it was shown that a neural network classifier trained on $S_{2}(R;\Delta R)$ and $S_{2,\mathrm{trim}}(R;\Delta R)$ spectra performs better than one without $S_{2,\mathrm{trim}}(R;\Delta R)$. The reason is that the hard and soft correlation terms in $S_{2}(R)$, i.e., $\mathcal{O}[1]$ terms in eq. (15) and $\mathcal{O}[f_{\mathrm{trim}}]+\mathcal{O}[f_{\mathrm{trim}}^{2}]$ in eq. (16) respectively, can be resolved by the jet trimming. Therefore, we use the orthogonal combinations, namely $S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$, throughout this paper.

The $S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$ spectra encode the important features of the parton shower and fragmentation, and, thus, may be regarded as a well-motivated prototype. The hard partons evolve by the parton splittings $i\rightarrow i_{1}i_{2}$, which are parameterized by the angle $R_{i_{1}i_{2}}$ and momentum fraction $z$ with $p_{T,i_{1}}=zp_{T,i}$ and $p_{T,i_{2}}=(1-z)p_{T,i}$. The splitting generates two-point correlation $z(1-z)p_{T,i}^{2}$ at $R_{i_{1}i_{2}}$. Therefore, $S_{2}$ spectra encode the parton splitting at any angular scale.

3 Classifying Higgs jet, sgluon jet, and QCD jet

In this section, we introduce a neural network trained on the jet spectra for classifying Higgs jet, sgluon jet, and QCD jet. We first discuss the basic kinematic features of these jets and then outline their dependence on the parton shower simulators. Afterward, we show the details of the neural network and then present our results in terms of the receiver operating characteristic (ROC) curves.

3.1 Basic Kinematics

In figure 4, we show $p_{T,\mathbf{J}}$ and $m_{\mathbf{J}}$ distributions for the Higgs boson, sgluon, and QCD jets. The solid and dashed lines correspond to Pythia 8 (PY8) and Herwig 7 (HW7) generated jets, respectively. The mild differences in the $p_{T}$ distribution are due to the difference in their matrix elements. The Higgs jet is produced via $s$-channel process, while the sgluon and QCD jet are produced via $t$-channel and $u$-channel processes; hence, $p_{T,\mathbf{J}}$ scalings are different. Not much difference is observed between the $p_{T,\mathbf{J}}$ distributions of PY8 and HW7 samples. This is because $p_{T,\mathbf{J}}$ is mostly determined by the matrix level $p_{T}$ of the leading parton and the jet algorithm with large radius parameter clusters most of the radiations from this parton into a single jet. However, the difference between $m_{\mathbf{J}}$ distributions is large. The peak at $m_{\sigma}$ of sgluon jet is significantly broader than that of Higgs jet because radiations of the b-partons from the Higgs boson decay are mostly confined due to the color coherence, but those of the sgluon are not. As a consequence, PY8 and HW7 generate different $m_{\mathbf{J}}$ distributions.

We assume that both Higgs boson and sgluon have narrow-widths although sgluon width can be large. An increase in the width will broaden $R_{b\bar{b}}$ distribution of $\sigma\rightarrow b\bar{b}$ that has a peak at the characteristic angular scale,

[TABLE]

For example, the variation of $R_{b\bar{b}}$ is only about 0.07 for $p_{T,\mathbf{J}}=300$ GeV, $\Gamma_{\sigma}=10$ GeV and 0.05 for $p_{T,\mathbf{J}}=400$ GeV, $\Gamma_{\sigma}=10$ GeV. Those variations are close to the calorimeter angular resolution $\sim$ 0.1 and do not affect the calorimeter level analysis.

We first make a quantitative estimate of the radiation pattern inside the jet. To do so, we define two quantities comparing $S_{2}(R)$ spectra around $\hat{R}_{b\bar{b}}$,

[TABLE]

with $a=0.75$, $a^{\prime}=1.25$ and $C=40$. The ratio $R_{\mathrm{sym}}$ compares energy deposits around $\hat{R}_{b\bar{b}}$ and in its surrounding angular scales Lim:2018toa . The ratio is sensitive to the correlation between the two hard substructures of the Higgs jet. On the other hand, The $R_{\mathrm{rad}}$ is sensitive to the color of mother particle as it compares energy deposits in the large angular scales.

We show the $R_{\mathrm{sym}}$ distributions in the left panel of figure 5. The distributions of the Higgs jet and sgluon jet are similar because both of the $S_{2}(R)$ spectra peak at ${R}_{b\bar{b}}$. Meanwhile, the two-point correlations for the QCD jet are not localized around the ${R}_{b\bar{b}}$ scale, so the $R_{\mathrm{sym}}$ is smaller than that of a Higgs jet and a sgluon jet. In the right panel of figure 5, we show the $R_{\mathrm{rad}}$ distributions. The $R_{\mathrm{rad}}$ of the sgluon jet and QCD jet are large on average, while $R_{\mathrm{rad}}$ is smaller for Higgs jet because large angle radiations are suppressed.

The difference in $R_{\mathrm{sym}}$ and $R_{\mathrm{rad}}$ distributions between PY8 and HW7 samples is small; however, there is an appreciable difference in the restricted phase space. In figure 6, we plot $R_{\mathrm{rad}}$ distributions after the selection, $R_{\mathrm{sym}}>0.85$, so that the jets always contain two hard subjets with similar transverse momenta. The PY8 (solid line) and HW7 (dashed line) samples have significantly different $R_{\mathrm{rad}}$ distributions for the Higgs jet. Such a difference is not observed for the QCD/sgluon jets. The observed deviation for the Higgs jets could be originating from the difference of the parton shower scheme. The angular-ordered shower is adopted in HW7. On the other hand, the $p_{T}$-ordered shower is the default shower algorithm for PY8 where angular ordering is enforced by hand. An artificial veto in $p_{T}$-ordered shower introduces the mismatch to the angular-ordered shower at double-leading log level Webber:2010vz ; Bhattacherjee:2015psa .

3.2 Multilayer Perceptron of Spectra

We introduce a neural network trained on the kinematic and spectrum ($S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$) variables to classify the jets. A schematic diagram of the architecture of the classifier is shown in figure 7. The following set of inputs is used,

[TABLE]

where $p_{T,\mathbf{J},\mathrm{trim}}$ and $m_{\mathbf{J},\mathrm{trim}}$ are the transverse momentum and mass of the trimmed jet, respectively. The discretized spectra $S_{2,\mathrm{trim}}^{k}$ and $S_{2,\mathrm{soft}}^{k}$ are used to analyze the radiation pattern of the jet,

[TABLE]

Here we take the bin width $\Delta R=0.1$, which is approximately the angular resolution of hadronic calorimeter of the ATLAS detector. Note that the maximum separation between any two constituents of the jet is $2R_{\mathbf{J}}$.

A multilayer perceptron (MLP) with $L$ layers is used to map the inputs to the class prediction. The following first-order recurrence relation between the layers describes an MLP,

[TABLE]

where $w^{(\ell)}_{ij}$ and $b^{(\ell)}_{i}$ are the weight and bias of the $\ell$-th layer. The activation function of the $\ell$-th layer, $\varphi^{(\ell)}:\mathbb{R}\rightarrow\mathbb{R}$, is a monotonic and nonlinear function. We use four hidden layers with 1000, 800, 400, and 200 nodes, respectively, with a rectified linear unit (ReLU), $\varphi_{\mathrm{ReLU}}(x)=\max(0,x)$, as the activation function. This MLP will identify important features of inputs for the classification after training. To make a class prediction, we provide the outputs of the MLP to a softmax classifier in eq. (18). The whole network architecture is illustrated in figure 8.

The MLP is trained by minimizing a loss function including categorical cross-entropy and $L_{2}$ weight regularization NIPS1991_563 ,

[TABLE]

where $N_{\mathrm{events}}$ is the total number of events in the training data set, $\lambda$ is a weight decay constant associated to the $L_{2}$ regularization. We choose $\lambda=0.01$. The $y_{i}$ ($\hat{y}_{i}$) denotes the components of the truth (predicted) label vector $\vec{y}$ ($\hat{y}$). The $L_{2}$ weight regularization reduces the over-fitting on the training data and also allows smooth extrapolation to the phase space that is not covered by the training sample. The minimization is done with ADAM optimizer DBLP:journals/corr/KingmaB14 . We stop training when the validation loss has stopped improving for 50 epochs. After the minimization of the loss function, the softmax layer provides scores of the classes of a given event. The truth label vectors are defined as follows,

[TABLE]

The unnecessary symmetries in the neural network are broken by using the Glorot uniform initialization method pmlr-v9-glorot10a . The weights in the hidden layers are initialized by assigning random numbers between $[-\sqrt{6/(N_{\mathrm{in}}+N_{\mathrm{out}})},\sqrt{6/(N_{\mathrm{in}}+N_{\mathrm{out}})}]$, where $N_{\mathrm{in}}$ and $N_{\mathrm{out}}$ are numbers of inputs and outputs of a layer, respectively. The biases are initialized to zero. All the inputs are standardized before training. The architecture is implemented in Keras chollet2015keras with backend TensorFlow tensorflow2015-whitepaper .

In figure 9, we show ternary plots of the predicted label vector $\hat{y}$. The three sides of the triangle (starting from the base of the triangle and then counterclockwise) are $\hat{y}_{1}$, $\hat{y}_{2}$, $\hat{y}_{3}$ axis; we denote them as $\hat{y}_{h}$, $\hat{y}_{\sigma}$ and $\hat{y}_{\rm QCD}$, respectively. The $\hat{y}$ distributions of the Higgs jet and QCD jet have high-density spots that do not overlap with each other. It means that the network has found the exclusive features of those two kinds of jets. The two-prong substructure of a Higgs jet and the one-prong structure of a QCD jet are the exclusive features. However, the two-prong substructure of a sgluon jet is more radiative and less exclusive, and therefore, there are no high-density spots in the $\hat{y}$ distribution of the sgluon jet.

Next, we show ROC curves of binary classifications in figure 10 with the red dotted lines. The following signal-background classifications are considered: Higgs-QCD, sgluon-QCD, and Higgs-sgluon. We assign the truth label vectors $\vec{y}=(1,0)$ for the signal and $\vec{y}=(0,1)$ for the background. The QCD jet mistag rates are comparable for both Higgs-QCD and sgluon-QCD classifications; however, the separation between the Higgs jet and sgluon jet is weaker.

We now compare these ROC curves with that of CNN trained on jet images.333The CNN setup is explained in detail in appendix C. The CNN classifier takes $20\times 20$ inputs of the jet images, while $2\times 20$ inputs of $S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$ spectra are used for the MLP. The solid blue lines in figure 10 denote the ROC curves of the CNN. Some improvement in the background mistag rates is observed compared with the MLP classifier. Quantitatively, it is only 0.2% ($=2.5\%-2.3\%$) at the signal acceptance of 20% for Higgs-QCD classification.

3.3 Event Generator Dependence

The classifier introduced in the previous subsection uses not only the information of hard subjets encoded in $S_{2,\mathrm{trim}}$ but also the soft activities captured in $S_{2,\mathrm{soft}}$ as well. This leads to concerns about the accuracy of the models of soft physics. Specifically, the performance of the classifier could be sensitive to the soft activities in the jet while the simulated soft activities may be significantly different from the truth.

In figure 11, we compare the ROC curves of the MLP trained with PY8 and HW7 samples. As these two event generators are based on different modeling of parton shower and hadronization scheme, the comparison would give us a reasonable estimate of the systematic uncertainty originating from the generator choice.

In the left panel of figure 11, we compare the ROC curves of the Higgs jet vs. QCD jet classification for different generator choices. By doing this exercise, we estimate a systematic uncertainty in the predictions of the classifier by comparing ROC(PY8, PY8) and ROC(HW7, HW7) curves, where the first and second entries in the parenthesis correspond to the generators used to simulate the training and test samples, respectively. On the other hand, ROC(HW7, PY8) and ROC(PY8, HW7) show the degradation of the performance of classifier trained on the “wrong sample” to analyze “real events.”

The performance of the classifier improves as we vary generator combinations in the following order: ROC(PY8, HW7), ROC(HW7, HW7), ROC(HW7, PY8), and ROC(PY8, PY8). We find that the classification performance is significantly better for PY8 test samples than that of HW7 samples. On the other hand, the classification performance for the same test samples hardly depends on the classifiers, namely ROC(PY8, HW7) $\sim$ ROC(HW7, HW7) and ROC(HW7, PY8) $\sim$ ROC(PY8, PY8). For the Higgs jet vs. QCD jet classification, the classifier mostly concentrates on the core substructures within the jet, and here both PY8 and HW7 provide similar kinematics and radiation spectra. Therefore, we do not observe any significant change in the ROC curves by varying training samples while keeping the test samples the same.

In the middle panel of figure 11, we compare the classifier performance for the sgluon jet vs. QCD jet classification. It improves in the following order: ROC(PY8, HW7), ROC(HW7, PY8), ROC(HW7, HW7), and ROC(PY8, PY8). The classifiers are indeed sensitive to the choice of generators. The network trained with PY8 (HW7) samples has failed to capture the features of HW7 (PY8) test samples. The networks have focused on different portions of the distribution of the fragmentation functions. In the right panel of figure 11, the ROC curves for the Higgs jet vs. sgluon jet classification show similar behavior.

4 Interpretable Two-level Architecture

A quantitative understanding of a neural network is not straightforward because the parameters and intermediate outputs of the neural network are less readable. In this section, we propose an architecture constructed from the truncated series in eq. (28) and try to explain quantitatively how this network classifies events. In the case of binary classifications, the discretized architecture is defined as follows,

[TABLE]

where $w^{k}_{A}$ is a trainable weight. We change the activation of the output layer to a sigmoid activation since the softmax function for binary classification is essentially a sigmoid with a scale factor on its argument. The loss function is the categorical cross-entropy as defined in eq. (43). This setup is effectively a logistic classifier on $S_{2,\mathrm{trim}}^{k}$ and $S_{2,\mathrm{soft}}^{k}$. After the training, the magnitude of $S_{2,\mathrm{trim}}^{k}w^{k}_{\mathrm{trim}}$ or $S_{2,\mathrm{soft}}^{k}w^{k}_{\mathrm{soft}}$ is high when the corresponding $S_{2,\mathrm{trim}}^{k}$ or $S_{2,\mathrm{soft}}^{k}$ is useful for the classification.

The logistic classifier does not take into account the $p_{T,\mathbf{J}}$ dependence of $\hat{R}_{b\bar{b}}$; therefore, we introduce a two-level architecture, which is a variant of the logistic classifier. The weights $w_{A}^{k}$ are calculated by a kinetic module $\Phi^{k}_{A}(\vec{x}_{\mathrm{kin}})$ of an MLP trained on $\vec{x}_{\mathrm{kin}}=(p_{T,\mathbf{J}},m_{\mathbf{J}})$,

[TABLE]

A schematic diagram of this setup is shown in figure 12. The inputs $\vec{x}_{\mathrm{kin}}$ are standardized before training, and $S_{2,\mathrm{trim}}^{k}$ and $S_{2,\mathrm{soft}}^{k}$ are divided by their maximum value of the training sample because standardizing the spectra reintroduce the zeroth order term of eq. (27). This architecture is similar to the self-explaining neural network NIPS2018_8003 . The $\Phi^{k}_{A}$ is modeled with an MLP of two hidden layers with exponential linear unit (ELU) activations DBLP:journals/corr/ClevertUH15 ,

[TABLE]

The nodes of the successive layers are configured as 400 $\mathrm{ELU}$, 200 $\mathrm{ELU}$, 2 $\times$ 20 linear respectively. We do not use $\mathrm{ReLU}$ in modeling $\Phi^{k}_{A}$ because dead $\mathrm{ReLU}$ nodes with a zero gradient kill $\vec{x}_{\mathrm{kin}}$ dependency of the weights. This is known as the dying $\mathrm{ReLU}$ problem. In this case, the architecture is reduced to the logistic classifier.

The vanishing gradient problem arises when the momentum range of the training sample is too small, which generates constant weights. The characteristic scale of $\hat{R}_{b\bar{b}}$ is $[0.625,0.833]$ for the $p_{T,\mathbf{J}}$ range [300, 400] GeV. The variation in $\hat{R}_{b\bar{b}}$ is about 0.2, which is not significantly large compared with the calorimeter resolution of 0.1. Therefore, we extend the $p_{T,\mathbf{J}}$ range of all the samples to $[300,\,600]$ GeV. In addition, we avoid vanishing gradient problem by using He uniform initializer He_2015_ICCV . The weights and biases are initialized by uniform random numbers in $[-\sqrt{6/N_{\mathrm{in}}},\sqrt{6/N_{\mathrm{in}}}]$ where $N_{\mathrm{in}}$ is the number of inputs to the layer. The advantage of using He initializer over Glorot initializer is that it generates random numbers in a wider range so that the neural network can start up from wider initial weights and gradient. The weight decay parameter $\lambda$ of the $L_{2}$ weight regularizer is set to 0.001 so that the weights do not vanish too early.

After the successful training, the performance of the two-level architecture is close to that of the MLP in section 3. The green dotted lines in figure 10 are the ROC curves of the classifier. The difference is smaller than the systematic uncertainty shown in figure 11. This makes a good reason to believe that the weights in eq. (47) capture the essential features of the MLP and CNN in section 3. The correlation between the output of the two-level architecture and the CNN model is shown in figure 13. We can see a positive correlation between them, but the correlation is slightly tilted towards the lower triangle (upper triangle) for small (large) $\hat{y}_{h}$ values because the CNN performs better than the two-level architecture.

In figure 14, we show the weight functions $w_{\mathrm{trim}}$ (left) and $w_{\mathrm{soft}}$ (right) of Higgs jet vs. QCD jet classifier trained with MG5+PY8+Delphes samples. Note that the weights are outputs of the neural network for the given $\vec{x}_{\mathrm{kin}}$ and not the output with the sample at the indicated $p_{T,\mathbf{J}}$ and $m_{\mathbf{J}}$. The weights in $R>R_{\mathbf{J}}$ and $w_{\mathrm{soft}}$ in $R<R_{\mathrm{trim}}$ are large compared to the weights in other angular scales. The $S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$ on these angular scales are typically smaller than that in other scales. Therefore, their weights become large to compensate for the energy difference when the corresponding value of the spectrum is useful for jet classification. The dotted lines in figure 14 denote the values of $\hat{R}_{b\bar{b}}$, defined in eq. (36), for different values of $p_{T,\mathbf{J}}$.

The $w_{\mathrm{trim}}$ around $\hat{R}_{b\bar{b}}$ is positive because the correlation at the scale is a characteristic feature of Higgs jet. If a Higgs boson is decaying to a pair of bottom quarks perpendicular to the boosted direction in its rest frame, the relative angular separation of the decay products is $\hat{R}_{b\bar{b}}$ in the lab frame. Due to the phase space of the decay, most of the events are distributed near $\theta=\pi/2$, where $\theta$ is Higgs decay angle relative to the boost direction. As $p_{T,\mathbf{J}}$ increases, $\hat{R}_{b\bar{b}}$ decreases, and the lower edge of $w_{\mathrm{trim}}>0$ moves toward smaller values. The region with $w_{\mathrm{trim}}>0$ also shifts towards smaller values of $R$. The weight $w_{\mathrm{trim}}$ on $R>\hat{R}_{b\bar{b}}$ is positive for capturing a Higgs jet whose $\theta$ is smaller than $\pi/2$. These events have $p_{T}$ asymmetric subjets, and the cross-correlation terms in $S_{2,\mathrm{trim}}(R)$ are smaller than that of $p_{T}$ symmetric case. These correlations are still useful for the classification because $S_{2}$ spectrum of QCD jet reduces much faster than that of Higgs jet. As a result, the $w_{\mathrm{trim}}$ is an increasing function in this region to compensate for the $S_{2,\mathrm{trim}}$ reduction.

For $R\gtrsim R_{\mathbf{J}}$, weight $w_{\mathrm{trim}}$ is negative. The score $\hat{y}_{h}$ decreases whenever there are any energy deposits at $R\gtrsim R_{\mathbf{J}}$. The crossover point from $w_{\mathrm{trim}}>0$ to $w_{\mathrm{trim}}<0$ shifts towards smaller values of $R$ with an increase of $p_{T,\mathbf{J}}$ because the Higgs decay products become more asymmetric with respect to the boost direction. In such a case, one of the subjets tends to be soft so that the two-point correlations are included in $S_{2,\mathrm{soft}}$, instead of $S_{2,\mathrm{trim}}$. These contributions to $S_{2,\mathrm{soft}}$ do not affect $w_{\mathrm{soft}}$, because $S_{2,\mathrm{soft}}$ spectra of QCD jets are overwhelming at large $R$.

The $S_{2,\mathrm{soft}}$ on $R>R_{\mathrm{trim}}$ always reduces $\hat{y}_{h}$, and there is no prominent structure around $\hat{R}_{b\bar{b}}$. Moreover, $|w_{\mathrm{soft}}|$ decreases as $p_{T,\mathbf{J}}$ increases. The reduction of $w_{\mathrm{soft}}$ compensates the increase of $S_{2,\mathrm{soft}}$, and the prediction is more or less $p_{T,\mathbf{J}}$ independent. On $R>R_{\mathbf{J}}$, $|w_{\mathrm{soft}}|$ increases with $R$ because activity in this region is a sign of QCD jet even though corresponding $S_{2,\mathrm{soft}}$ decreases due to suppressed large angle radiations.

The $w_{\mathrm{soft}}$ on $R\lesssim R_{\mathrm{trim}}$ is positive and $w_{\mathrm{soft}}$ has a break at $R\sim R_{\mathrm{trim}}$. Correlations between the constituents in a soft subjet contributes to the $S_{2,\mathrm{soft}}$ on $R<R_{\mathrm{trim}}$, i.e., $S_{2,\mathrm{soft}}(0;R_{\mathrm{trim}})\propto f^{2}_{\mathrm{trim}}$. Let us assume $\mathbf{J}_{a}$ is a single soft subjet, then $S_{2,\mathrm{soft}}(0;R_{\mathrm{trim}})\sim p_{T,\mathbf{J}_{a}}^{2}\sim(p_{T,\mathbf{J}}-p_{T,\mathbf{J},\mathrm{trim}})^{2}$. If there are multiple soft subjets, then $S_{2,\mathrm{soft}}(0;R_{\mathrm{trim}})\sim\sum_{a}p_{T,\mathbf{J}_{a}}^{2}<(\sum_{a}p_{T,\mathbf{J}_{a}})^{2}\sim(p_{T,\mathbf{J}}-p_{T,\mathbf{J},\mathrm{trim}})^{2}$. This triangular inequality suggests that the magnitude of $S_{2,\mathrm{soft}}(0;R_{\mathrm{trim}})$ is small for a jet with a given $p_{T,\mathbf{J}}-p_{T,\mathbf{J},\mathrm{trim}}$ when there are multiple soft jets. The positive $w_{2}$ on $R<R_{\mathrm{trim}}$ means that Higgs jet has less soft subjets than QCD jet. The $S_{2,\mathrm{soft}}^{0}$ consists of the autocorrelation of soft subjets, which has different energy scaling behavior compared with the other $S_{2,\mathrm{soft}}^{k}$. The $S_{2,\mathrm{soft}}$ on $R<R_{\mathrm{trim}}$ is $S_{2,22}\sim\mathcal{O}[f_{\mathrm{trim}}^{2}]$ in eq. (24). On the other hand, $S_{2,\mathrm{soft}}^{k}$ $(k\geq 1)$ is dominated by $S_{2,12}$. The $S_{2,12}$ on $R<R_{\mathrm{trim}}$ does not contribute to $S_{2,\mathrm{soft}}^{0}$ because it vanishes. Therefore, we may rewrite $h$ as follows,

[TABLE]

where $w_{\mathrm{soft}}^{\prime}$ and $w_{\mathrm{soft}}^{\prime\prime}$ are continuous functions with $w_{\mathrm{soft}}^{\prime}(R_{\mathrm{trim}})=0$ and $w_{\mathrm{soft}}(R)=w_{\mathrm{soft}}^{\prime}(R)+w_{\mathrm{soft}}^{\prime\prime}(R)$. The second term is essentially the same as $\int dR\,S_{2,22}(R)w_{2,22}^{(2)}(R)$ in eq. (23), and the last term is $\int dR\,S_{2,12}(R)w_{2,12}^{(2)}(R)+\int dR\,S_{2,21}(R)w_{2,21}^{(2)}(R)$.

The sudden changes of $w_{\mathrm{trim}}$ and $w_{\mathrm{soft}}$ at $R\simeq 1.8$ are due to the statistical fluctuations of the training sample. The $S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$ may have a non-zero value at large R if the jet has multiple large angle radiations with the opposite direction from the jet axis; however, the probability of such a radiation pattern is small. As a result, the number of events used for training the weights at large $R$ is not sufficient. The large weights also do not contribute much to the classification (see figure 16).

To estimate this statistical uncertainty, we use both training and test datasets. The merged dataset is divided into ten subsets and we train a network for each of them. This will decrease the number of events in each subset by a factor 5. We estimate the uncertainty of the fit by calculating the mean and variance of the $w_{\mathrm{trim}}^{k}$ and $w_{\mathrm{soft}}^{k}$ from the ten subsets. The green band in the bottom panels of figure 14 represents the estimated uncertainty. Note that $w_{\mathrm{trim}}^{k}$ and $w_{\mathrm{soft}}^{k}$ are not sensitive to the network initialization because the two-level network is effectively a logistic regression for a fixed $p_{T,\mathbf{J}}$ and $m_{\mathbf{J}}$ and the loss function is a convex function of them.

In the top panels of figure 15, we show the weights $w_{\mathrm{trim}}$ and $w_{\mathrm{soft}}$ for the sgluon jet vs. QCD jet classification with $p_{T,\mathbf{J}}=350$ GeV and $m_{\mathbf{J}}=115$ GeV. The $w_{\mathrm{trim}}$ distribution is similar to that of the Higgs jet vs. QCD jet classification. However, the $|w_{\mathrm{soft}}|$ is much smaller. This comes from the fact that $S_{2,\mathrm{soft}}$ of sgluon jet is similar to that of QCD jet and it is less important in the classification. No peak of $S_{2,\mathrm{soft}}$ around $R\lesssim R_{\mathrm{trim}}$ also indicates that the soft substructures of sgluon jet are as radiative as QCD jet. Additionally, there is no color coherence restriction of soft radiations for the sgluon jet. This leads to small $|w_{\mathrm{soft}}|$ for $R>R_{\mathbf{J}}$. In the bottom panels of figure 15, we show the weights for the Higgs jet vs. sgluon jet classification. The peak of $w_{\mathrm{trim}}$ around $R=\hat{R}_{b\bar{b}}$ is small as the hard substructures of Higgs jet and sgluon jet are (almost) the same. However, a sgluon is more radiative than a Higgs boson, and $w_{\mathrm{soft}}$ is negative in the entire region of $R<1.5$.

As described above, weights $w_{\mathrm{trim}}$ and $w_{\mathrm{soft}}$ may take large values, but it does not necessarily mean that the corresponding $S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$ contribute dominantly in the jet classification. The energy scaling factors on the $S_{2,\mathrm{trim}}$ ($S_{2,\mathrm{soft}}$) and its weight $w_{\mathrm{trim}}$ ($w_{\mathrm{soft}}$) cancel out in the quantity of our interest $h=\sum_{k}(S_{2,\mathrm{trim}}^{k}w_{\mathrm{trim}}^{k}+S_{2,\mathrm{soft}}^{k}w_{\mathrm{soft}}^{k})$. For example, $\mathcal{O}[1]$ terms in $S_{2,\mathrm{trim}}$ and $\mathcal{O}[f_{\mathrm{trim}}]$ terms in $S_{2,\mathrm{soft}}$ contribute equally to the classifier if $w_{\mathrm{soft}}$ is around $f_{\mathrm{trim}}\,w_{\mathrm{trim}}$. In the left panel of figure 16, we draw the mean values $\langle S_{2,\mathrm{trim}}^{k}w^{k}_{\mathrm{trim}}\rangle$ and $\langle S_{2,\mathrm{soft}}^{k}w^{k}_{\mathrm{soft}}\rangle$ of Higgs jet vs. QCD jet classification, which are more directly related to the jet classification. The solid and dashed red lines correspond to the distributions of the Higgs jet, while the solid and dashed blue lines are for the QCD jet. The regions where Higgs jet and QCD jet distributions differ significantly are important for the network predictions. We find $\langle S_{2,\mathrm{trim}}^{k}w^{k}_{\mathrm{trim}}\rangle$ around $R\sim\hat{R}_{b\bar{b}}$ and $\langle S_{2,\mathrm{soft}}^{k}w^{k}_{\mathrm{soft}}\rangle$ in the region $R<1.2$ mostly contribute to the jet classification.

The average distribution may not illustrate all the features of the classifier performance. The energy deposits in each bin fluctuate, and the bins with hits higher than the average value contribute more to the network decisions. For example, soft emissions outside the angle between the two hardest subjets are rare in the Higgs jet. Once there is large angle radiation outside the cone of hard subjets, the network is likely to identify the jet as a QCD jet. In the right panel of figure 16, we plot the $\langle S_{2,\mathrm{trim}}^{k}w_{\mathrm{trim}}^{k}\rangle$ and $\langle S_{2,\mathrm{soft}}^{k}w_{\mathrm{soft}}^{k}\rangle$ distributions of the Higgs jet (QCD jet) with $\hat{y}_{h}$ ($\hat{y}_{\mathrm{QCD}}$) higher than the 95th percentile. The distributions indicate that the selected Higgs jets are mostly classified because of the $S_{2,\mathrm{trim}}$ excess at $\hat{R}\sim R_{b\bar{b}}$, while the QCD jets are classified using $S_{2,\mathrm{soft}}$ excess above $R>0.2$.

We now use the two-level architecture to compare PY8 and HW7. As we have already shown in section 3, the performance of the classifier depends on the event generators significantly. We show the weights of the classifiers for $p_{T,\mathbf{J}}=350$ GeV and $m_{\mathbf{J}}=115$ GeV trained with Higgs jets in figure 17, sgluon jets in figure 18, and QCD jets in figure 19. For each plot, the signals are PY8 events, and the backgrounds are HW7 events. The $w_{\mathrm{trim}}$ in $R\lesssim R_{\mathbf{J}}$ is close to zero everywhere, representing $S_{2,\mathrm{trim}}$ spectra of PY8 and HW7 events are similar. It is not surprising because both PY8 and HW7 events from identical hard partons and these partons create the trimmed subjets inside the jet. On the other hand, the correlation involving constituents of the soft activities, $S_{2,\mathrm{soft}}$, is manifestly different and so $w_{\mathrm{soft}}$ is nonzero. For Higgs jet and sgluon jet, the $w_{\mathrm{soft}}$ distribution of PY8 events is significantly large (and positive) for $R\sim R_{\mathrm{trim}}$ and it decreases as $R$ increases. The weight $w_{\mathrm{soft}}$ is negative for $R>R_{\mathbf{J}}$, which means that HW7 events have more soft activity in the region $R\gg\hat{R}_{b\bar{b}}$.

For the case of QCD jet, the distribution of $w_{\mathrm{soft}}$ is always positive and flat for $R<R_{\mathbf{J}}$ and negative for $R>1.5$. It would be interesting to evaluate the weights for the classifiers trained with the experimental data and compare with the simulated results to tune the parameters of the event generators further.

5 Summary and Outlook

The classification of jets with deep learning has gained significant attention in recent times. Majority of these analyses take advantage of the significant development in computing power. These deep learning architectures utilize the complete event information in terms of low-level observables. These deep learning based strategies can be compared with the previous approaches for tagging jets, for example, mass drop tagger, $n$-subjettiness, energy correlation function where each of them has solid physics motivation.

In this paper, we introduce neural networks trained on “jet spectrum” $S_{2}(R)$, which are essentially two-point correlation functions of jet constituents. We also introduce $S_{2,\mathrm{trim}}(R)$, which is $S_{2}(R)$ calculated from the trimmed jet, to encode hard substructure of the jet. The difference $S_{2,\mathrm{soft}}(R)=S_{2}(R)-S_{2,\mathrm{trim}}(R)$ encodes the remaining correlations with soft radiations and is less affected by the correlations among the hard constituents. Our neural networks are trained on $S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$ integrated over certain bins. If the $S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$ spectra are multiplied by smooth functions and integrated over $R$, it forms an IRC safe C-correlator. This feature assures that the classifiers trained on binned $S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$ are approximately IRC safe.

The performance of MLP trained on $S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$ is compared to that of CNN trained on jet images. The CNN shows better performance than the MLP, but the difference is small. The key reason is efficient preprocessing of parton shower effects with less number of free parameters. Parton shower is the multiple splittings of the partons where each splitting is parametrized by the angular scale and momentum fraction of the partons. The binned $S_{2}(R)$ spectra collect the information of the parton splitting successfully. The spectra provide comparable jet classification performance with a fewer number of inputs. Furthermore, the MLP is computationally economical than the CNN because the MLP has smaller complexity than the CNN and takes only $\mathcal{O}(40)$ inputs.

The $S_{2,\mathrm{trim}}(R)$ and $S_{2,\mathrm{soft}}(R)$ spectra can be obtained from a functional Taylor series of an arbitrary classifier in energy flows. The spectra are basis vectors of the second order term in the expansion. In this context, the MLP trained on $S_{2,\mathrm{trim}}(R)$ and $S_{2,\mathrm{soft}}(R)$ can be considered as a sum of $2n$-linear $C$-correlators which can be reduced to products of the bilinear $C$-correlators in $S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$. The (mild) difference in the performance of the CNN and MLP comes from the remaining irreducible $n$-linear $C$-correlators.

The terms linear in $S_{2,\mathrm{trim}}(R)$ and $S_{2,\mathrm{soft}}(R)$ provide an opportunity to visualize and interpret the network predictions; therefore, we study a novel two-level architecture that involves an interpretable layer of a single node in the form of a $C$-correlator. The output is the sum of the product of the trained weights ($w_{\mathrm{trim}}$ and $w_{\mathrm{soft}}$) and jet spectra ($S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$). The absolute values of the weights signify the impact of the corresponding $S_{2,\mathrm{trim}}$ and $S_{2,\mathrm{soft}}$ bin values on the jet classification. In the context of classification between Higgs jet and QCD jet, the distribution of $w_{\mathrm{trim}}$ shows that $S_{2,\mathrm{trim}}$ spectrum around $\hat{R}_{b\bar{b}}=2m_{h}/p_{T,\mathbf{J}}$ increases the output, and the classifier regards the jet as a Higgs jet. We have also shown that the dependence of $w_{\mathrm{trim}}$ on jet $p_{T}$ can be qualitatively understood (at the parton level) from the decay of a boosted Higgs boson. In short, the network is using $S_{2,\mathrm{trim}}$ inputs to obtain the core substructure information inside the jet.

The soft activity is also useful for Higgs jet vs. QCD jet classification. The probability for assigning the given jet as a QCD jet increases with increase in $S_{2,\mathrm{soft}}$ on $R>R_{\mathrm{trim}}$. To study the impact of soft physics in jet classification, we also introduce sgluon, a hypothetical color octet scalar, and compare the classifier performance among Higgs jet, sgluon jet, and QCD jet. The network predictions for sgluon jet vs. QCD jet classification are primarily determined by the core substructure information as expected. However, the network uses the difference in the $S_{2,\mathrm{soft}}$ spectra arising from the different color structure of the decaying particle for Higgs jet vs. sgluon jet classification.

The non-trivial role of soft radiations in the predictions of the classifiers implies the results are highly sensitive to the choice of event generators. The weights associated with the $S_{2,\mathrm{trim}}$ are almost insensitive to the choice; however, the weights of the $S_{2,\mathrm{soft}}$ are strongly affected. This behavior is expected as modeling of soft physics is quite different in Pythia 8 and Herwig 7.

The two-point correlation spectra and the architectures introduced in this paper can be applied for solving other interesting problems, thanks to flexibility on designing neural network. For jets with more complex substructures, e.g., top jet, the higher order terms in the energy flow series expansion may be included. It would be worthwhile to study the classifier performances when the network is trained with the experimental data and compare with the predictions of event generators to tune their parameters to reduce the uncertainty in modeling the soft physics. It is also interesting to use this interpretable architecture as a model-agnostic interpreter for black box architectures 2016arXiv160605386T . We leave these possibilities for future works.

Acknowledgements

We thank the organizers of Beyond the BSM and Machine Learning for Jet Physics 2018 workshops. This work is supported by the Grant-in-Aid for Scientific Research on Scientific Research B (No.16H03991, 17H02878) and Innovative Areas (16H06492), and by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan.

Appendix A Event generation and reconstruction

The parton level event samples, namely $pp\rightarrow Zj$, $pp\rightarrow Zh$ and $pp\rightarrow Z\sigma$ events, are generated at the leading order in QCD using MadGraph5_aMC@NLO 2.6.1 Alwall:2014hca . We force the Higgs boson ($h$) and the sgluon ($\sigma$) to decay to a pair of bottom quarks, while $Z$ boson to decay invisibly. For sgluon, we use a UFO model in Degrande:2014sta ; NLOModels with the following interaction term for the decay,

[TABLE]

The parton distribution function (PDF) set NNPDF 2.3 LO at $\alpha_{S}(m_{Z})=0.130$ Ball:2012cx is used. To generate Higgs jets and Sgluon jets, we impose a parton level selection criterion on the $Z$ boson transverse momentum, $p_{T,Z}>250$ GeV. We simulate approximately 3 million events of $pp\rightarrow Zh$ and $pp\rightarrow Z\sigma$ processes and 18 million events of $pp\rightarrow Zj$.

We use two parton shower and hadronization simulators to compare the results. Namely, we use Pythia 8.226 Sjostrand:2014zea with Monash tune Skands:2014pea and Herwig 7.1.3 Bellm:2015jjp ; Bahr:2008pv with default tune herwigtune ; Gieseke:2012ft . The shower starting scale is $H_{T}/2$ for $pp\rightarrow Zh$ and $pp\rightarrow Z\sigma$ processes and $p_{T,j}$ for $pp\rightarrow Zj$ process, where $H_{T}$ is the transverse energy sum of the produced partons. The effects of underlying events and multi-parton interactions are taken into account, but we neglect the contaminations coming from the pile-ups. The PDF set for simulating all these effects are the same as that in the parton level simulation.

We use Delphes 3.4.1 deFavereau:2013fsa with its default ATLAS configuration for fast detector simulations. Jets are reconstructed from the calorimeter towers using FastJet 3.3.0 Cacciari:2011ma ; Cacciari:2005hq with anti-$k_{T}$ algorithm Cacciari:2008gp and jet radius parameter $R_{\mathbf{J}}=1$. The leading jet of each event with $p_{T,\mathbf{J}}\in[300,400]\,\mathrm{GeV}$ and $m_{\mathbf{J}}\in[100,150]\,\mathrm{GeV}$ is selected. Since the scale $H_{T}/2$ for $pp\rightarrow Zh$ and $pp\rightarrow Z\sigma$ is higher than $p_{T,h}$ and $p_{T,\sigma}$, respectively, there is a chance that the leading jet is from the initial state radiation rather than from the decay of Higgs boson or sgluon. To filter out such jets from the Higgs jet and sgluon jet samples, we require that $b$-partons produced from the decay are within $R_{\mathbf{J}}$ from the leading jet axis.

Appendix B Oversampling and $p_{T,\mathbf{J}}$-bias removal

The neural networks may learn the inherent $p_{T,\mathbf{J}}$ difference among Higgs jet, sgluon jet, and QCD jet in figure 4 to classify them instead of learning the difference in their substructures. To penalize the learning from the $p_{T,\mathbf{J}}$ distribution, we augment the training and validation samples by oversampling as follows,

The samples (of each class) are binned in $p_{T,\mathbf{J}}$ with bin-width 1 GeV with $b_{i}$ entries in the $i$-th bin and $b_{\max}=\max\{b_{i}\}$. 2. 2.

For each bin, oversample the events so that the number of the bin contents becomes a certain value $n_{\max}=c_{d}\,b_{\max}$. This oversampling is identical to repeating events in bin $i$ sequentially for

[TABLE]

times and stop the oversampling when the number of events in the bin reaches $n_{\max}$. It may introduce a small bias due to unequal oversampling among different bins. We choose $c_{d}=2$ and ignore the small residual bias.

In the upper (lower) panel of figure 20, we show the conditional probability density of the predicted class vector $\hat{y}_{h}$ ($\hat{y}_{\sigma}$) for a given $p_{T}$. The probability density has a mild dependence on $p_{T,\mathbf{J}}$ which is originating from the interplay of the phase-space selection and the jet radius parameter.

Appendix C Jet Image and Convolution Neural Network

We obtain and pre-process the jet image as follows,444This setup is similar to Kasieczka:2017nvn .

Recluster the jet constituents by $k_{T}$ algorithm with a jet radius parameter $R_{\mathbf{J}}=0.2$. 2. 2.

Set the center of $(\eta,\phi)$ coordinate to the leading (in $p_{T}$) subjet. 3. 3.

If a second leading subjet is found, rotate the jet constituents on $(\eta,\phi)$ plane about the jet center so that the sub-leading jet is on the positive $y$-axis. 4. 4.

If a third leading subjet is found, flip the image about $y$-axis when $x$ coordinate of the subjet is negative. 5. 5.

Select jet constituents within $[-1.5,1.5]\otimes[-1.5,1.5]$. 6. 6.

Finally, pixelate the jet constituents with pixel size $0.1\times 0.1$. The $(k,l)$-th pixel intensity $P_{T}^{k,l}$ is determined by the total transverse energy of the jet constituents present in a given pixel, i.e.,

[TABLE]

where $\mathrm{bin}_{k,l}$ is the region of $(k,l)$-th bin. 7. 7.

Standardize all the $P_{T}^{k,l}$.

This jet image is analyzed by a CNN which consists of two-dimensional convolutional layers (CONV) and max-pooling layers (figure 21). In particular, we use the following CNN setup,

• •

Layer 1: Convolutional layer with 64 filters with kernel size $3\times 3$, and $\mathrm{ReLU}$ activation,

• •

Layer 2: Max-pooling layer with pool size $2\times 2$,

• •

Layer 3: Convolutional layer with 32 filters with kernel size $4\times 4$, and $\mathrm{ReLU}$ activation,

• •

Layer 4: Max-pooling with pool size $2\times 2$.

The first convolutional layer deals with angular scale up to 0.3 to treat collinear radiations while the second convolutional layer operates up to 0.8. The outputs of Layer 4 are flattened into a one-dimensional array and concatenated with a set of kinematic inputs, $\{p_{T,\mathbf{J}},m_{\mathbf{J}},p_{T,\mathbf{J},\mathrm{trim}},m_{\mathbf{J},\mathrm{trim}}\}$. The flattened output array is fed into an MLP with two hidden layers with 300, 100 filters respectively, and $\mathrm{ReLU}$ activation function. The outputs of the MLP are fed into a softmax layer to make a prediction. The training setup is the same as in section 3.