# A fast map-making preconditioner for regular scanning patterns

Sigurd K. Næss<EMAIL_ADDRESS>Thibaut Louis
<EMAIL_ADDRESS>Sub-department of Astrophysics, University of
Oxford, Keble Road, Oxford, OX1 3RH, UK

###### Abstract

High-resolution Maximum Likelihood map-making of the Cosmic Microwave
Background is usually performed using Conjugate Gradients with a
preconditioner that ignores noise correlations. We here present a new
preconditioner that approximates the map noise covariance as circulant, and
show that this results in a speedup of up to 400% for a realistic scanning
pattern from the Atacama Cosmology Telescope. The improvement is especially
large for polarized maps.

##  I Introduction

As the resolution and sensitivity of Cosmic Microwave Background (CMB)
experiments increase, so do the computational resources needed to analyze
their data. Because modern detectors are background-limited, the only way to
significantly increase sensitivity is to increase the number of detectors. The
last decades have seen an increase from tens of detectors to thousands of
detectors in experiments like ACT Dünner et al. (2013); Niemack et al. (2010),
SPT Carlstrom et al. (2011); Austermann et al. (2012), POLARBEAR Kermish et
al. (2012) and Keck Kernasovskiy et al. (2012), and plans already exist for
experiments with $10^{5}-10^{6}$ detectors Abazajian et al. (2013). Reducing
the data from all these detectors into a coherent map of the sky presents a
significant computational challenge, and already with 1000-detector-class
experiments this step is the most important bottleneck of the data analysis
pipeline Dünner et al. (2013). It is therefore important to investigate ways
to speed up this process.

Three main classes of map-makers are in popular use: _Maximum likelihood map-
makers_ Tegmark (1997); Cantalupo et al. (2010); de Gasperis et al. (2005);
Traficante et al. (2011); Dünner et al. (2013); QUIET Collaboration et al.
(2012), which are slow, but produce unbiased, optimally noise-weighted maps;
faster but slightly less accurate _destripers_ Keihänen et al. (2005); Sutton
et al. (2010); and biased and sub-optimal but very fast _naive map-makers_
QUIET Collaboration et al. (2012); Schaffer et al. (2011). The topic of this
paper is a method for significantly speeding up maximum likelihood map-makers.

Assuming a linear detector response, we can model the time-ordered data $d$
via the the linear system

$\displaystyle d=Pm+n,$ (1)

where $m$ is the pixelized map of the sky, the pointing matrix $P$ is a
sparse111The pointing matrix will be sparse if we solve for a beam-convolved
map. For variable beams or asymmetric beams, one may want to reconvolve to a
standard beam as part of map-making. This can be done using $P$, at the cost
of some of its sparsity. We do not consider this case here. matrix mapping
from pixels to samples, and $n$ is the time-domain noise, which we assume to
be gaussian with covariance matrix $N$. The maximum likelihood solution for
$m$ is given by the map-making equation Tegmark (1997),

$\displaystyle(P^{T}N^{-1}P)m$ $\displaystyle=P^{T}N^{-1}d.$ (2)

This is of the form $Ax=b$, and while the matrices involved are usually too
large to solve by direct inversion, the system is amenable to solution by
Preconditioned Conjugate Gradients (PCG) Press et al. (2007) provided a good
preconditioner can be found222Without a preconditioner, the number of
iterations needed for conjugate gradients is proportional to the condition
number of the matrix $A$. By applying a preconditioner $M$, one is effectively
solving the system $MAx=Mb$. The goal is then to choose $M$ such that $MA$ is
as well-conditioned as possible. This can be acchevied if $M\approx A^{-1}$..
The most commonly used map-making preconditioners are the binned Natoli et al.
(2001); Ashdown et al. (2007); Cantalupo et al. (2010) and Jacobi Doré et al.
(2001); de Gasperis et al. (2005); Cantalupo et al. (2010); Dünner et al.
(2013) preconditioners. The binned preconditioner $M_{\textrm{B}}$
approximates the time covariance matrix $N$ as diagonal (i.e. it ignores
correlations), which for pointing matrices where only one pixel is hit per
sample results in a diagonal pixel-space covariance matrix.

$\displaystyle N^{-1}_{ij}$ $\displaystyle\approx
N^{-1}_{ii}\delta_{ij}\textrm{ and }P_{ti}P_{tj}\propto\delta_{ij}\Rightarrow$
(3) $\displaystyle A_{ij}$ $\displaystyle\equiv
P_{ti}N^{-1}_{tt^{\prime}}P_{t^{\prime}j}\approx
P_{ti}^{2}N^{-1}_{tt}\delta_{ij}\equiv{M_{\textrm{B}}}^{-1}_{ij}$ (4)

The Jacobi preconditioner $M_{\textrm{J}}$ simplifies one step further, and
assumes that every detector has the same variance, $N=aI$, resulting in

$\displaystyle A_{ij}$ $\displaystyle\approx
aP_{ti}^{2}\delta_{ij}\equiv{M_{\textrm{J}}}^{-1}_{ij}.$ (5)

The proportionality factor $a$ is usually set to 1, as PCG is insensitive to
an overall scaling of the preconditioner.

$-10^{-4}$ $10^{-4}$

Figure 1: A $4.2^{\circ}$ by $3.2^{\circ}$ subset of a row from the pixel-
space inverse correlation matrix from a patch from ACT (each row in the matrix
corresponds to a two-dimensional map). In order to highlight the correlation
structure, the color scale is capped at $\pm 10^{-4}$. The significantly
correlated area has a quite complicated shape, which is driven by the scanning
pattern and focal plane layout.

$-10^{-3}$ $10^{-3}$

Figure 2: A mosaic of the local correlation structure for a $12^{\circ}$ by
$9^{\circ}$ subset of an ACT patch. A set of evenly spaced pixels were chosen
as reference points, and for each the correlation in a neighborhood around it
was calculated. Each pixel in the map is colored according to its correlation
relative to the nearest reference point. The correlation structure is very
uniform in the whole central region of the map.

The assumption of independent noise is quite inaccurate. All realistic
experiments have at least some time-correlation in the noise, and usually also
correlations between different detectors. Additionally, filters will also
generally introduce correlations.

While computing the full, exact $A$ is often too expensive, single rows of it
can be computed at the same expense as one CG iteration:

$\displaystyle A_{ij}$ $\displaystyle=A_{ik}\delta_{kj}=(A\vec{e}_{i})_{j}$
(6)

Here $(\vec{e}_{i})_{j}=\delta_{ij}$ is the pixel-space basis vector
corresponding to pixel $i$. An example of what such a row looks like for the
Atacama Cosmology Telescope (ACT) Dünner et al. (2013) can be seen in figure
1.

It is clear that the approximation of no correlation is quite inaccurate.
However, since the correlation structure is driven by the scanning pattern and
relative position of the detectors in the focalplane, the correlation
structure should be the same for all pixels which are scanned the same way.
ACT, which used long-duration, small-amplitude drift scans, pixels at the same
declination but different right ascension will be hit by the same phase of the
same scanning motion, and should therefore have the same correlation
structure.333 Circulant correlation is not a good approximation for every
experiment. It is suitable for constant elevation dift scans, as employed by
ACT, SPT and POLARBear, but we expect it to work poorly for full-sky scanning
patterns.

Indeed, that is what measurements show (see figure 2). It may therefore be a
good approximation to assume that every point in the map has the same relative
correlation structure, i.e. that the correlation between two points on the sky
only depends on their relative position. If this is the case, then it is
possible to choose a pixelization where the correlation only depends on the
difference between pixel numbers, and hence that the pixel correlation matrix
is _circulant_ 444 For example, if
$\textrm{corr}(\vec{x}_{1},\vec{x}_{2})=f(\vec{x}_{1}-\vec{x}_{2})$, where
$\vec{x}$ are coordinates, then a pixelization scheme $\vec{x}=G(p)$, where
$p$ is a pixel index and G is a linear function will fulfill
$\textrm{corr}(p_{1},p_{2})=f(G(p_{1}-p_{2}))$, resulting in a circulant
correlation matrix..

A circulant matrix has the nice property of being diagonal in the frequency
domain, which means that it can be computed, stored and applied cheaply, at
the cost of a few FFTs. Hence, the constant correlation approximation is
promising as a preconditioner for solving the map-making equation.

##  II Implementation

The inverse pixel covariance matrix $A$ can be decomposed into variance and
correlation such that $A=\Sigma^{T}U\Sigma$. Here $\Sigma$ is diagonal (block-
diagonal in the case of polarization) in pixel space, and corresponds to a map
of the inverse standard deviation per pixel. As per the binned preconditioner,
this can be approximated as

$\displaystyle(\Sigma^{T}\Sigma)_{ij}$
$\displaystyle=P_{ti}^{2}N^{-1}_{tt}\delta_{ij}.$ (7)

The correlation matrix $U$ is in general a dense matrix, but as noted above,
it can often be approximated as circulant555The matrix will be circulant
provided that the correlation structure is position-independent, that a
constant offset in each coordinate corresponds to a constant pixel offset, and
provided that indices wrap around at the edges.. Therefore, the constant
correlation preconditioner replaces $U$ with a circulant matrix $Q$, such that

$\displaystyle A$
$\displaystyle\approx\Sigma^{T}Q\Sigma\equiv{M_{\textrm{C}}}^{-1}$ (8)

This relation can be inverted to give us an expression for $Q$ in terms of
$A$,

$\displaystyle Q=$ $\displaystyle{\Sigma^{T}}^{-1}A\Sigma^{-1}\Rightarrow$ (9)
$\displaystyle Q_{ij}=$
$\displaystyle{\Sigma^{T}}^{-1}_{ii}(A\vec{e}_{i})_{j}\Sigma^{-1}_{jj}.$ (10)

Since $Q$ is circulant, i.e. $Q_{ij}=Q_{0,j-i}=q_{j-i}$, we have

$\displaystyle(FQF^{-1})_{ff^{\prime}}$
$\displaystyle=(Fq)_{-f}\delta_{ff^{\prime}}=(Fq)^{*}_{f}\delta_{ff^{\prime}}$
(11)

using forward and backwards Fourier transforms
$F_{fj}\equiv\mathrm{e}^{-\frac{2\pi ijf}{N}}$ and
$F^{-1}_{jf}\equiv\frac{1}{N}\mathrm{e}^{\frac{2\pi ijf}{N}}$, where $N$ is
the number of rows in the matrix.

With this in hand, the preconditioner can be applied as

$\displaystyle M_{\textrm{C}}b$
$\displaystyle=\Sigma^{-1}Q^{-1}{\Sigma^{-1}}^{T}b$
$\displaystyle=\Sigma^{-1}F^{-1}{(Fq)^{*}}^{-1}F{\Sigma^{-1}}^{T}b.$ (12)

$\Sigma^{-1}$ and ${(Fq)^{*}}^{-1}$ can be precomputed, so the cost of
applying the preconditioner is simply that of two FFTs and three diagonal
matrix multiplications.

The choice of the reference pixel at which the correlation is measured is
somewhat arbitrary. We used the pixel nearest the center of the map, but other
locations not too close to the edge of the map should also work.

Figure 3: Example TOD noise power spectrum from the simulation, based on noise
behavior from ACT. Low frequencies are dominated by atmospheric noise, while
the increase at high frequencies is caused by a Butterworth filter.

##  III Regularization

While constant correlation is a good approximation for relatively short-scale
correlations, it works less well for long-distance correlations, and regions
near the edges666Near the edges the telescope must decelerate in order to
reverse the scanning direction, which makes the correlation structure
different there than in the center. Applying the preconditioner as described
above to realistic cases results in the appearance of large scale modes which
change extremely slowly during the subsequent CG iteration.

A way around this is to artificially limit the range of the correlations that
are modeled, by multiplying $q$ by a Gaussian. For ACT, a standard deviation
of 20 pixels was found to be effective, but this will depend on the scanning
pattern, and some experimentation may be needed to find the optimal number.

Figure 4: The residual A-norm Strakos & Tichy (2008) as a function of
iteration number for a simple binned preconditioner and the constant
correlation preconditioner. The latter converges roughly twice as fast as the
former according to this criterion. The horizontal blue line indicates the
level at which the maps have mostly stopped changing visually.

##  IV Polarization

The previous discussion assumed that each pixel only had a single degree of
freedom, e.g. temperature-only maps of the sky. In the case of polarization,
each pixel has several correlated components, typically the Stokes parameters
T, Q and U Zaldarriaga & Seljak (1997), but this can instead be expressed as a
larger number of block-correlated single-component pixels. This results in
$\Sigma$ being block-diagonal with e.g. one (T,Q,U)-block per physical pixel,
while $Fq$ becomes a vector of similar blocks. And instead of a single row of
$A$ needing to be measured, all the rows corresponding to a given physical
pixel now need to be computed (i.e. $\vec{e}_{i\alpha}$ for all components,
where Greek indices indicate polarization components). Hence, eq. (10) becomes

$\displaystyle Q_{i\alpha j\beta}=$ $\displaystyle{\Sigma^{T}}^{-1}_{i\alpha
i\gamma}(A\vec{e}_{i\gamma})_{j\delta}\Sigma^{-1}_{j\delta j\beta}$ (13)

Aside from that, everything works the same.

##  V Test setup

| Step 5 | Step 15 | Step 45 | Step 115 | Step 340
---|---|---|---|---|---
Binned |  |  |  |  |
Const. corr. |  |  |  |  |
Figure 5: Example temperature maps from the CG solution process for the binned (top) and constant correlation (bottom) preconditioners. The rows correspond to steps 5, 15, 45, 115 and 340 from left to right. The steps are chosen such that the binned map in column $n$ is as similar as possible as the constant correlation map in column $n-1$. We see that the constant correlation preconditioner visually converges about 3 times faster than the binned one. The maps have been cropped for compactness of presentation. Step 15 | Step 70 | Step 250
---|---|---
|  |
Figure 6: Comparison of the input (blue), binned (red) and constant
correlation (green) power spectra at CG steps 15 (left), 70 (middle) and 250
(right). Unlike the map-space plots in figure 5, where only the large-scale
convergence is visible, we can here clearly see the convergence at all scales.
Small-scale convergence is much slower for EE and BB than for TT, and form the
bottleneck for the CG solver if we ignore the $\ell<500$ modes. The constant
correlation preconditioner is 3-5 times faster than the binned one here. This
is quantified more precisely in figure 9.

We tested the preconditioners on a simulated time-ordered data (TOD) based on
the scanning pattern for 64 detectors from each of 417 15-minute scans of a
subset of ACT’s southern patch centered at
$\alpha=56^{\circ},\delta=-53^{\circ}$. Each scan was a constant elevation
drift scan with amplitude of $3.5^{\circ}$ in azimuth, and the scan centers
were spread over 10 steps in elevation, covering a patch of about $11^{\circ}$
by $8^{\circ}$ degrees. Odd steps in elevation scanned while rising and the
even ones when setting. This resulted in most pixels being hit from two
directions, and hence the x-shaped correlation pattern seen in figures 1-2. An
example of a noise power spectrum used in the simulation can be seen in figure
3.

The simulated detectors were polarization-sensitive, with each detector
measuring a linear combination $T+\cos(2\psi)Q+\sin(2\psi)U$ of the local
radiation field, with each detector having a different, randomly chosen
detector angle $\psi$. While an ACT-like noise model, including the effects of
atmosphere and inter-detector noise correlations was assumed in the map-making
step, no noise was added to the simulated TOD in order to allow the
convergence to be studied all the way to the highest multipoles777 This is
valid since the convergence rate of PCG is mostly independent of the noise
level of the right-hand side after the first few iterationsNatoli et al.
(2001). However, with higher noise, higher CG errors also become acceptable,
so the number of iterations needed for CG errors to be subdominant will be
smaller for realistic noise levels.. For the same reason, the simulated input
CMB did not include a beam, and was pixelated at the same resolution as the
output map, in order to avoid subpixel noise.

We then solved the map-making equation for this data set using PCG, first
using the binned preconditioner described in equation (4), and then the
constant correlation approximation described in this paper. Each was run for
600 CG iterations, with intermediate maps being output for every 5 steps.

##  VI Results

The constant correlation preconditioner visually converges roughly 3 times
faster than the binned one, as shown for the temperature map in figure 5.
Likewise, the residual A-norm Strakos & Tichy (2008) from the CG solver (shown
in figure 4) also shows a significant improvement in convergence: roughly a
factor 2 according to this metric888The A-norm $||x||_{A}$ of a vector $x$ is
defined as $\sqrt{x^{T}Ax}$, where $A$ is $P^{T}N^{-1}P$ in our case. The
error A-norm after $i$ CG steps is $||x_{i}-x||_{A}$, where $x_{i}$ is our
estimate after $i$ steps, and $x$ is the true map. When the true $x$ is
unknown, $||x_{i}-x||_{A}$ can be estimated as
$||x_{i}-x||_{A}\approx\sum_{j=i}^{i+d-1}\gamma_{j}||r_{j}||^{2}$, where
$r_{j}$ and $\gamma_{j}$ are two internal variables in the CG algorithm at
step $j$, and $d$ is is an integer that controls the accuracy of the estimate
(4 in our case).

Figure 7: An example of the convergence of a single bin in the power spectrum,
in this case that centered on $\ell=3200$. This shows the absolute difference
between the recovered spectrum $C_{\ell}$ and the spectrum $\hat{C}_{\ell}$ of
the input map in units of cosmic variance $\sigma_{\ell}^{\textrm{cv}}$, as a
function of the CG step in cosmic variance units. Shown are curves for TT
(red), EE (green) and BB (blue) for the binned (solid) and constant
correlation (dashed) preconditioners. The trend lines are based on binned
averages of the errors for many individual CG steps. This binning reduces the
jitter, and makes the trend lines suitable for measuring the time needed to
converge to a given level.

However, neither of these tests take into account the fact that not all scales
in the map are equally interesting. To remedy this, figure 6 compares the
binned and constant correlation power spectra with that of the simulated input
map. These spectra were computed using the method described in Louis et al.
(2013). On large scales ($\ell<500$), this tells the same story as the maps
did: The larger the scale, the more slowly it converges, with the constant
correlation preconditioner being about 3 times faster than the binned one.
Somewhat surprisingly, a similar phenomenon occurs at the small scales. For
$\ell>2000$, higher $\ell$ results in slower convergence, and this is
especially prominent for the EE and BB power spectra. On all scales, however,
the constant correlation preconditioner appears to converge several times
faster than the binned one.

In order to quantify the convergence more precisely, we consider the time at
which the absolute error in a given multipole-bin reaches 0.1 times cosmic
variance in that bin. While somewhat arbitrary, this choice ensures that CG
errors are guaranteed to be sub-dominant in the power spectrum, regardless of
the noise properties of the actual experiment. Figure 7 shows the convergence
of TT, EE and BB for both preconditioners for a typical multipole-bin. The
overall trend for each component is an initial rapid fall followed by a slower
decay, with both being significantly faster for the new preconditioner,
particularly for the polarization spectra.

Figure 8: The number of CG steps needed for each multipole-bin to converge to
0.1 times cosmic variance, for each of TT (red), EE (green) and BB (blue) for
the binned (solid) and constant correlation (dashed) preconditioners.

We found that binning the errors in bins of $\sim 50$ conjugate gradients
steps and using linear interpolation between these bins resulted in a robust
estimate of when each spectrum bin reaches the convergence criterion. The
resulting convergence times can be seen in figure 8, and confirm our earlier
finding that the largest and smallest scales converge more slowly. The figure
also highlights how much trouble the binned preconditioner has with the EE and
especially BB spectra, where it performs much more poorly relative the
constant correlation preconditioner than we see for the TT spectrum. We
speculate that this is due to the X-shaped correlation structure introduced by
our scanning pattern. In binned maps, which ignore the correlations, this
introduces spurious X-shaped patterns in both Q and U, corresponding to
spurious signal in both E and B. With the constant correlation preconditioner,
these are partially corrected because some of the correlation structure is
taken into account.

Figure 9: The ratio of the convergence times to 0.1 times cosmic variance for
the binned and constant correlation preconditioners, per multipole bin for
each of TT, EE and BB. For TT, the speedup is typically between 20% and 200%,
while for BB the speedup ranges from 200% at large scales to about 400% at
small scales, with EE being intermediate. Large points here indicate cases
where the binned preconditioner did not converge to 0.1 times cosmic variance
in time. For these points, the comparison was performed at the lowest
threshold where both converged.

We summarize the performance characteristics of the new preconditioner in
figure 9, which shows its relative speed gain vs. the baseline binned
preconditioner. Depending on $\ell$, we have a speedup ranging from 20% to
200% for TT, from 150% to 300% for EE, and from 200% to 400% for BB, with the
greatest relative improvement happening at the scales that converge most
slowly.

##  VII Summary

The structure of a CMB map’s pixel covariance matrix is determined by the
noise properties of the time-ordered data and the scanning pattern of the
telescope. For constant elevation drift scans like those employed by ACT, SPT
and POLARBEAR, this results in an approximately circulant covariance. We have
developed a new preconditioner for conjugate gradient solutions of the map-
making equation which exploit this property by deconvolving the correlations
in harmonic space, an operation which is very cheap due to the Fourier
representation of a circulant matrices being diagonal.

For a realistic scanning pattern and noise model the preconditioner results a
speedup of 20% to 200% for temperature and 150% to 400% for polarization
compared to a binned preconditioner.

Convergence speed might potentially be further improved by allowing the
correlation pattern to change slowly across the map, for example by
decomposing the map into overlapping tiles, and applying the constant
correlation preconditioner separately to each tile, followed by a merging
operation. Our preliminary attempts at such a tiled preconditioner have
however not been able to beat the simple constant correlation approximation
presented here.

###### Acknowledgements.

The authors would like to thank Jo Dunkley and Johannes Noller for useful
discussion and suggestions, and Jon Sievers for testing the preconditioner in
another map-maker. We also thank the ACT collaboration for access to internal
ACT data used in the simulations. Computations were performed on the gpc
supercomputer at the SciNet HPC Consortium. SciNet is funded by: the Canada
Foundation for Innovation under the auspices of Compute Canada; the Government
of Ontario; Ontario Research Fund - Research Excellence; and the University of
Toronto. SN and TL are supported by ERC grant 259505.

## References

  * Abazajian et al. (2013) Abazajian, K. N., Arnold, K., Austermann, J., et al. 2013, ArXiv e-prints, arXiv:1309.5383
  * Ashdown et al. (2007) Ashdown, M. A. J., Baccigalupi, C., Balbi, A., et al. 2007, A&A, 467, 761
  * Austermann et al. (2012) Austermann, J. E., Aird, K. A., Beall, J. A., et al. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8452, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series
  * Cantalupo et al. (2010) Cantalupo, C. M., Borrill, J. D., Jaffe, A. H., Kisner, T. S., & Stompor, R. 2010, ApJS, 187, 212
  * Carlstrom et al. (2011) Carlstrom, J. E., Ade, P. A. R., Aird, K. A., et al. 2011, PASP, 123, 568
  * de Gasperis et al. (2005) de Gasperis, G., Balbi, A., Cabella, P., Natoli, P., & Vittorio, N. 2005, A&A, 436, 1159
  * Doré et al. (2001) Doré, O., Teyssier, R., Bouchet, F. R., Vibert, D., & Prunet, S. 2001, A&A, 374, 358
  * Dünner et al. (2013) Dünner, R., Hasselfield, M., Marriage, T. A., et al. 2013, ApJ, 762, 10
  * Keihänen et al. (2005) Keihänen, E., Kurki-Suonio, H., & Poutanen, T. 2005, MNRAS, 360, 390
  * Kermish et al. (2012) Kermish, Z. D., Ade, P., Anthony, A., et al. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8452, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series
  * Kernasovskiy et al. (2012) Kernasovskiy, S., Ade, P. A. R., Aikin, R. W., et al. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8452, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series
  * Louis et al. (2013) Louis, T., Næss, S., Das, S., Dunkley, J., & Sherwin, B. 2013, ArXiv e-prints, arXiv:1306.6692
  * Natoli et al. (2001) Natoli, P., de Gasperis, G., Gheller, C., & Vittorio, N. 2001, A&A, 372, 346
  * Niemack et al. (2010) Niemack, M. D., Ade, P. A. R., Aguirre, J., et al. 2010, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7741, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series
  * Press et al. (2007) Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2007, Numerical Recipes 3rd Edition: The Art of Scientific Computing, 3rd edn. (Cambridge University Press)
  * QUIET Collaboration et al. (2012) QUIET Collaboration, Araujo, D., Bischoff, C., et al. 2012, ApJ, 760, 145
  * Schaffer et al. (2011) Schaffer, K. K., Crawford, T. M., Aird, K. A., et al. 2011, Astrophys. J. , 743, 90
  * Strakos & Tichy (2008) Strakos, Z., & Tichy, P. 2008, ETNA, 13, 56
  * Sutton et al. (2010) Sutton, D., Zuntz, J. A., Ferreira, P. G., et al. 2010, MNRAS, 407, 1387
  * Tegmark (1997) Tegmark, M. 1997, ApJL, 480, L87
  * Traficante et al. (2011) Traficante, A., Calzoletti, L., Veneziani, M., et al. 2011, MNRAS, 416, 2932
  * Zaldarriaga & Seljak (1997) Zaldarriaga, M., & Seljak, U. 1997, Phys. Rev. D, 55, 1830


∎

11institutetext: Emanuel Pfeffer1 (corresponding author) 22institutetext:
22email<EMAIL_ADDRESS>33institutetext: Michael Waßmer1
44institutetext: 44email<EMAIL_ADDRESS>55institutetext: Yee-Ying
Cung 66institutetext: 66email<EMAIL_ADDRESS>77institutetext: Roger
Wolf1 88institutetext: 88email<EMAIL_ADDRESS>99institutetext: Ulrich
Husemann1 1010institutetext: 1010email<EMAIL_ADDRESS>1111institutetext: 1 Karlsruhe Institute of Technology, Institute of
Experimental Particle Physics, Karlsruhe, Germany

# A case study of sending graph neural networks back to the test bench for
applications in high-energy particle physics

Emanuel Pfeffer Michael Waßmer Yee-Ying Cung  Roger Wolf Ulrich Husemann

###### Abstract

In high-energy particle collisions, the primary collision products usually
decay further resulting in tree-like, hierarchical structures with a priori
unknown multiplicity. At the stable-particle level all decay products of a
collision form permutation invariant sets of final state objects. The analogy
to mathematical graphs gives rise to the idea that graph neural networks
(GNNs), which naturally resemble these properties, should be best-suited to
address many tasks related to high-energy particle physics. In this paper we
describe a benchmark test of a typical GNN against neural networks of the
well-established deep fully-connected feed-forward architecture. We aim at
performing this comparison maximally unbiased in terms of nodes, hidden
layers, or trainable parameters of the neural networks under study. As physics
case we use the classification of the final state $\mathrm{X}$ produced in
association with top quark-antiquark pairs in proton-proton collisions at the
Large Hadron Collider at CERN, where $\mathrm{X}$ stands for a bottom quark-
antiquark pair produced either non-resonantly or through the decay of an
intermediately produced $\mathrm{Z}$ or Higgs boson.

###### Keywords:

Graph Neural Networks Deep Neural Networks High-Energy Particle Physics LHC

##  1 Introduction

The continuous rise and flourish of deep learning has significantly impacted
also the community of high-energy particle physics, where modern algorithms of
deep learning —mostly in the form of various neural network (NN)
architectures— find applications as automation tools, for (multiclass)
classification, parameter regression, or universal function approximation. The
Large Hadron Collider (LHC) at CERN offers a unique test environment for such
algorithms providing a large amount of independent identically distributed
(i.i.d.) data from proton-proton ($\mathrm{pp}$) collisions under well
controlled laboratory conditions. These data feature a rich hierarchical
structure, optimally suited for the application of all kinds of general
methods of statistical data analysis. Moreover, the underlying physics laws
and statistical models, which have emerged over many decades of research, are
scrutinized to a level that allows the reliable estimation of particle
properties with a relative accuracy ranging far below the per-mille level, in
rare cases even below $10^{-10}$ Muong-2:2023cdq . This circumstance offers a
toolbox for generating a large amount of perfectly known, complex, synthetic
data, with a high relation to experimental observations, through the
application of Monte Carlo (MC) methods MC-method-1 ; MC-method-0 . These data
are usually obtained as samples from an intractable though well-known
likelihood function $\mathcal{L}$. This setup provides a unique opportunity to
thoroughly benchmark any kinds of machine learning (ML) algorithms under
complex, real-life laboratory conditions.

At the LHC, data analysts strive for the application of more and more
sophisticated ML-models with more and more not further processed —and in this
sense “raw”— input data. This strategy is fed by the belief that automated
algorithms might find ways of extracting information of interest to the
analyst, which are superior to selection strategies that are vulnerable to the
bias of human prejudice. On the other hand, ML-algorithms should not be forced
to learn already known and well-established physics principles, like
symmetries inherent to the presented task. While such information can only be
insufficiently passed through the necessarily finite training samples, it can
be intrinsically incorporated either into the loss functions used for
training, or in the NN architectures.

At the large multi-purpose LHC experiments, ATLAS Aad:2008zzm  and CMS
Chatrchyan:2008aa , $\mathrm{pp}$ collisions at a center-of-mass energy of,
e.g., $13\,\mathrm{TeV}$ result in the creation of thousands of collision
products to be recorded by the experiments. Primary collision products might
decay further resulting in tree-like, hierarchical structures with a priori
unknown multiplicity. The collision process can be described in a factorized
approach:

During the hard scattering process, the fundamental constituents of the
protons, i.e., the quarks and gluons which are also collectively referred to
as partons, interact via the fundamental interactions under investigation. We
refer to the result of these interactions as the partonic final state. It
cannot be observed directly in an experiment. Rather, each parton undergoes a
series of theoretically well-known processes, setting in at lower energy
scales, resulting in stable particles. The inverse problem usually subject of
high-energy particle physics is to infer the presence and properties of the
stable particles and eventually the partonic final state from their observable
energy deposits in the detectors.

At the stable particle level all decay products of a collision form
permutation invariant sets of final state objects, which may emerge from the
collision in the form of collimated particle jets Salam:2009jx , forming well-
suited proxies for strongly-interacting final state partons, or individual,
spatially isolated particles, like leptons. From the preparation of the
collision’s initial state and energy and momentum conservation, physicists may
infer the presence of non- or weakly-interacting particles, like neutrinos, in
the collision’s final state, through the principle of missing transverse
energy (MET) CMS:2019ctu . A natural representation of this richly structured
data is in the form of mathematical graphs $\mathcal{G}$, which are indeed
also the basis of theoretical amplitude calculations of the quantum-mechanical
wave function in the form of Feynman graphs PhysRev.76.749 .

Within the high-energy particle physics community, this observation has lead
to an increased interest in NNs based on mathematical graphs (GNNs)
GNN_original_1 ; GNN_original_2 ; GNN_original_3 ; GNN_original_4 ;
GNN_original_5 , where nodes are usually identified by particles and edges
potentially by relations across particles. A comprehensive review can be found
in Ref. Shlomi:2020gdn . However, on closer inspection, several features of
GNNs that count for seemingly obvious advantages reveal subtle challenges:

  * •

The mathematical model to match physics entities like stable particles to the
graph $\mathcal{G}$ of the GNN, which is subject to mathematical operations,
bears ambiguities. The representation of particles by the nodes of
$\mathcal{G}$ appears as an obvious choice. This choice raises the question of
the (potentially physical) meaning to the connecting edges, which could
represent mother-daughter relations, or proximity in an arbitrarily defined
space. This issue is emphasized once mathematical weights or even trainable
parameters (TPs) used during node aggregation are assigned to the edges.

  * •

Training and application imply the potentially complex and computationally
time-consuming task of building $\mathcal{G}$, based on the physics inputs.

  * •

In general, the more complex structure of the GNN compared to other NN
architectures complicates the users’ comprehension of how the GNN arrives at
its prediction.

In this paper we describe a comparison of typical GNN architectures with NN
models based on the deep fully-connected feed-forward architecture (DNN),
which has been studied intensely in the past. We aim at comparisons that are
to best effort unbiased in terms of expressiveness and information provided to
the NN models to solve a given task. In Section 2 we give an introduction to
the task that serves as benchmark for this comparison. In Section 3 the
architectures and training setups of the NNs under study are described. In
Section 4 we present the results of the comparison. We conclude the paper with
a summary in Section 5.

##  2 Neural network task

###  2.1 Physics processes

As benchmark for the comparison we use the classification of the final state
$\mathrm{X}$ produced in association with top ($\mathrm{t}$) quark-antiquark
pairs ($\mathrm{t\bar{t}}$) in $\mathrm{pp}$ collisions at the LHC, where
$\mathrm{X}$ stands for a bottom ($\mathrm{b}$) quark-antiquark pair
($\mathrm{b\bar{b}}$) produced either through non-resonant gluon exchange or
through the decay of a massive $\mathrm{Z}$ or Higgs ($\mathrm{H}$) boson as
intermediate particle, as discussed, e.g., in Ref. CMS:2018hnq . Under
realistic conditions, the collision of interest might be overlaid by several
tens of additional collisions, referred to as pileup. The complete detectable
final state of a collision of interest, including pileup, is referred to as an
event, whose feature vector $\mathbf{x}$ would be presented to the NN. An
arbitrary number of such events may be generated synthetically by evaluating
$\mathcal{L}$ of the full process via the MC method. In this study, we focus
on the classification of the underlying hard process neglecting the effects of
pileup.

The interest in the chosen classification task arose from studies of
$\mathrm{H}$ production in association with $\mathrm{t\bar{t}}$ in the
subsequent $\mathrm{H}\to\mathrm{b}\bar{\mathrm{b}}$ decay
($\mathrm{ttH(bb)}$), for which $\mathrm{t\bar{t}}$ associated $\mathrm{Z}$
boson production ($\mathrm{ttZ}$) in the
$\mathrm{Z}\to\mathrm{b}\bar{\mathrm{b}}$ decay channel ($\mathrm{ttZ(bb)}$)
and non-resonant $\mathrm{b\bar{b}}$ production in association with
$\mathrm{t\bar{t}}$ ($\mathrm{ttbb}$) are important background processes.
Exemplary Feynman diagrams of these processes in leading-order (LO) of
perturbation theory are shown in Figure 1.

Figure 1: Exemplary Feynman diagrams for the processes of interest to this
study: (left) $\mathrm{ttbb}$, (middle) $\mathrm{ttH(bb)}$, and (right)
$\mathrm{ttZ(bb)}$.

The decay products, resembled by the outgoing lines in the diagrams, shown in
Figure 1, represent the partonic final state of interest to this study, which
is the same for all processes. Therefore, the processes can only be
distinguished by the kinematic properties of the particles, in particular
$\mathrm{b}$ quarks. This situation is complicated by the fact that the
$\mathrm{t}$ quark also decays into $\mathrm{b}$ quarks radiating a quasi-real
$\mathrm{W}$ boson with a branching fraction of nearly 100% Workman:2022ynf .
The $\mathrm{W}$ boson subsequently decays either into quarks, which further
on form jets in the detector, or leptons. For the presented study the semi-
leptonic $\mathrm{t\bar{t}}$ final state has been chosen, in which the
$\mathrm{W}$ boson of one $\mathrm{t}$ decays into an electron or muon,
further on referred to as $\mathrm{\ell}$, and a corresponding neutrino
$\mathrm{\nu_{\ell}}$. The other $\mathrm{W}$ boson decays into quarks. Due to
the radiation of additional gluons and the splitting of gluons into quark-
antiquark pairs additional colored particles and consequently jets might
emerge from the process. This constellation implies a richly structured final
state of an event with at least four $\mathrm{b}$ quark- and two predominantly
light-quark-induced jets; an $\mathrm{\ell}$, which is spatially isolated from
any other activity originating from the hard scattering process in the
detector; and MET, due to the emitted $\mathrm{\nu_{\ell}}$. The $\mathrm{b}$
quark-induced jets, referred to as $\mathrm{b}$ jets in the following, may be
identified experimentally with a finite purity and efficiency, as, e.g.,
described in Ref. CMS:2017wtu ; the methods of how to achieve this are not
subject of this paper.

###  2.2 Sample preparation

Samples for all processes in question have been generated synthetically from a
corresponding likelihood $\mathcal{L}$ using the MC technique. The tools used
for event generation are the matrix-element generator MadGraph5_aMC@NLO
Alwall:2014hca ; Alwall:2011uj  in version 2.9.9 interfaced with the Pythia
event generator Sjostrand:2014zea  in version 8.306 to map the partonic final
state to the stable-particle level. All processes in consideration have been
generated at LO in perturbative quantum chromodynamics (QCD), in the four-
flavor scheme. The same setup has been used for the generation of all samples
to avoid spurious differences due to the use of different generation tools.

All generated events have been passed through a simplified simulation of the
CMS detector as configured during the LHC Run-2 data-taking period in the
years 2016–2018, using the DELPHES simulation package deFavereau:2013fsa . For
this study only reconstructed leptons, jets, and MET are considered. All
detected and reconstructed final-state objects have been selected to fulfill a
set of selection criteria typically used for the analysis of data collected by
the CMS experiment, as summarized in Table 1. These selection criteria
comprise the following observables:

  * •

The transverse momentum $p_{\mathrm{T}}$ and pseudorapidity $\eta$ of the
reconstructed $\mathrm{\ell}$ and jets.

  * •

A variable $I_{\mathrm{rel}}^{\Delta R}$ that corresponds to the scalar sum of
energy deposits $E_{i}$ detected within the radius $\Delta
R=\sqrt{\Delta\eta^{2}+\Delta\phi^{2}}$ around $\mathrm{\ell}$, divided by the
$p_{\mathrm{T}}$ of $\mathrm{\ell}$, where $\Delta\eta$ refers to the
difference between $E_{i}$ and $\mathrm{\ell}$ in $\eta$ and $\Delta\phi$ to
the corresponding difference in azimuthal angle $\phi$, based on the
coordinate system deployed by CMS Chatrchyan:2008aa . For $\mathrm{\ell}$
originating from $\mathrm{W}\to\mathrm{\ell}\nu_{\mathrm{\ell}}$ low values of
$I_{\mathrm{rel}}^{\Delta R}$ are expected.

  * •

The output of a specific $\mathrm{b}$ jet identification algorithm represented
by the discrete observable $\beta\in\\{\mathrm{0,1,2,3}\\}$ indicating whether
an object has been identified as a $\mathrm{b}$ jet under a specific working
point $\alpha$ ($\beta\geq\alpha$). The value of $\alpha$ represents
($\alpha=1$) loose, ($\alpha=2$) medium, and ($\alpha=3$) tight selection
criteria, corresponding to a rate of non-$\mathrm{b}$ jets wrongly identified
as a $\mathrm{b}$ jet (false-positive rate), of approximately 10%, 1%, and
0.1%, respectively.

Table 1: Selection requirements on the reconstructed final-state objects. The quantity $I_{\mathrm{rel}}^{\Delta R}$ corresponds to the scalar sum of energy deposits detected within the radius $\Delta R$ around $\mathrm{\ell}$ in $\eta$-$\phi$ divided by the $p_{\mathrm{T}}$ of $\mathrm{\ell}$, as defined in the text. A lower value of $I_{\mathrm{rel}}^{\Delta R}$ implies less activity in the spacial vicinity of $\mathrm{\ell}$, indicating that $\mathrm{\ell}$ originates from $\mathrm{W}\to\mathrm{\ell}\nu_{\mathrm{\ell}}$. The variable $\beta$ refers to the working points of a specific $\mathrm{b}$ jet identification algorithm, as described in the text. Object | $p_{\mathrm{T}}\,(\,\mathrm{GeV})$ | $\left|\eta\right|$ | $I_{\mathrm{rel}}^{\Delta R}$ | $\beta$
---|---|---|---|---
Electron | $\geq 25$ | $<2.5$ |  $<0.12$ ($\Delta R=0.3$) | $-$
Muon | $\geq 25$ | $<2.4$ |  $<0.25$ ($\Delta R=0.4$) | $-$
Jet (anti-$k_{\text{T}}$, $R_{0}$=0.4 Cacciari:2008gp ) | $\geq 20$ | $<2.4$ | $-$ | $-$
$\mathrm{b}$ jet | $\geq 20$ | $<2.4$ | $-$ | $\beta\geq 2$

All events have been selected to exhibit at least six and not more than eight
jets, at least four of which are assumed to be correctly identified as
$\mathrm{b}$ jets according to $\beta\geq 2$, and exactly one $\mathrm{\ell}$,
matching all selection criteria. The selection of one $\mathrm{\ell}$ and six
jets, of which at least four are identified as $\mathrm{b}$ jets, is motivated
by the partonic final state under study, as depicted in Figure 1. The
selection of up to two additional jets, increases the chance that the complete
partonic final state can be matched to the selected jets.

For each reconstructed jet the attempt is made to assign the initiating
particle of the partonic final state, based on the distance $\Delta R$ between
the jet and the corresponding parton. Only partonic final state objects with
transverse momentum of $p_{\mathrm{T}}>20\,\mathrm{GeV}$ and
$\left|\eta\right|<2.4$ are considered for this assignment. From the
assignment five mutually exclusive jet classes are build, referring to the
(ADDB) $\mathrm{b}$ quarks not originating from any $\mathrm{t}$ decay; the
$\mathrm{b}$ quark originating from the (HTB) hadronic and (LTB) leptonic
$\mathrm{t}$ decay, (HTQ) quarks originating from the hadronic $\mathrm{W}$
decay; and (NA) jets not assigned to the partonic final state.

The assignment of the partonic final state to the reconstructed jets may be
incomplete, for a given event. Events for which no or only one jet is assigned
to the ADDB class are discarded from the training and test samples. In all
other cases, if the assignment by $\Delta R$ did not result in one jet of
class LTB, one jet of class HTB, and two jets of class HTQ, the remaining not
associated jets of class NA are ordered by decreasing $p_{\mathrm{T}}$ and the
leading jets in $p_{\mathrm{T}}$ are re-assigned to these classes in the order
of LTB, HTB, HTQ. A summary of all jet classes is given in Table 2

Table 2: Requirements for the jet-class definition, according to the matching to the partonic final state, where $\mathrm{b_{add}}$ corresponds to a $\mathrm{b}$ quark not originating from a $\mathrm{t}$ decay, $\mathrm{b_{t_{had}}}$ ($\mathrm{b_{t_{lep}}}$) to a $\mathrm{b}$ quark originating from $\mathrm{t}\to\mathrm{b}\mathrm{W}(\mathrm{q}\mathrm{q}^{\prime})$ ($\mathrm{t}\to\mathrm{b}\mathrm{W}(\ell\nu_{\ell})$), and $\mathrm{q_{W_{had}}}$ to a quark originating from $\mathrm{W}\to\mathrm{q}\mathrm{q}^{\prime}$. We note that the classes ADDB and HTQ comprise two jets. Jets not assigned to any other class are assigned to the NA class and sorted by decreasing $p_{\mathrm{T}}$. If no jet is found to fulfill the corresponding $\Delta R$ criterion, the leading jet from the NA class is re-assigned to the classes HTB, HTL, HTQ, in that order. Class label | Assignment | Description
---|---|---
ADDB | $\Delta R(\mathrm{jet},\mathrm{b_{add}})<0.4$ |  $\mathrm{b}$ jets not from $\mathrm{t}$ decays
HTB | $\Delta R(\mathrm{jet},\mathrm{b_{t_{had}}})<0.4^{\dagger}$ |  $\mathrm{b}$ jet from $\mathrm{t_{had}}$
LTB | $\Delta R(\mathrm{jet},\mathrm{b_{t_{lep}}})<0.4^{\dagger}$ |  $\mathrm{b}$ jet from $\mathrm{t_{lep}}$
HTQ | $\Delta R(\mathrm{jet},\mathrm{q_{W_{had}}})<0.4^{\dagger}$ |  $\mathrm{q}$ jet from $\mathrm{W}$
NA | No match | Additional jet
† If no assignment by $\Delta R$ is found, the leading jet from NA is
assigned.

###  2.3 Task definition

The NN models are supposed to perform a classification task, in which a given
event should be assigned either to $\mathrm{ttbb}$, $\mathrm{ttH(bb)}$, or
$\mathrm{ttZ(bb)}$. The physics process from which the event was generated
constitutes the ground truth information for this task.

To simplify the task, the reconstructed jets are assumed to perfectly match
the initiating partons, whenever possible, through the parton association
algorithm, as described in Section 2.2. The ADDB class contains jets which
stem directly from the intermediate particle to distinguish the processes
under study. The features of these jets are assumed to provide the most
decisive contribution to the classification.

###  2.4 Training setup

All NN models are subject to a supervised training. The generated samples used
for this are split into a training, validation, and test sample, containing
60, 20, and 20% of the generated events, respectively. Event numbers, split by
process and sample, are given in Table 3.

Table 3: Numbers of events for each process, in the training, validation, and test samples. Process | Training | Validation | Test | Sum
---|---|---|---|---
$\mathrm{ttbb}$ | 41650 | 13883 | 13803 | $\hphantom{0}69336$
$\mathrm{ttH(bb)}$ | 99695 | 33329 | 33229 | 166253
$\mathrm{ttZ(bb)}$ | 88107 | 29272 | 29453 | 146832

Each training is performed for an ensemble of ten statistically independent
repetitions to obtain a rough estimate of the statistical spread of the
trained models, due to random choices in the training setup. Performance
measures of each model are reported as sample means $\mu$, of which the
corresponding uncertainty $\Delta\mu$ is estimated from the square root of the
sample variance. All repetitions are based on the same training, validation,
and test samples. Due to the large number of events in these samples,
randomization through data shuffling is assumed not to change the conclusions
of the studies significantly.

Each training is performed on CPUs through a distributed computing
infrastructure, where each training is performed on a dedicated CPU. We ensure
that each step affected by random choices is based on a different random seed.
The library used to build the GNNs is PyTorch Geometric v2.0.3
cite:PyTorch_Geometric  based on the PyTorch library cite:PyTorch . The same
library is used, to construct the corresponding DNNs under study. Further
parameter choices of the training setups are given in the upper part of Table
4. They are the same for both NN architectures.

Table 4: Common parameter choices for the NN architectures under study and their training setup. Parameter | Setting(s)
---|---
Loss function | Binary cross-entropy
Optimizer | Adam Kingma:2014vow  ($\gamma=0.01$)
Mini-batch size | 200
Maximum number of epochs | 200
Early-stopping |  $\Delta\text{epochs}=15$, $\Delta\text{loss}=0.001$
Use of weights and biases | Yes
Number of outputs | 1 (binary)
Activation function (for hidden layers) | ReLu
Activation function (for output layer) | Sigmoid

##  3 Neural networks under study

###  3.1 Architectures

For the GNN models, a graph representation of the final state is obtained from
the reconstructed jets, $\mathrm{\ell}$, and MET in an event. Each of these
objects is represented by a node $i$ in $\mathcal{G}$. This set of nodes may
be complemented by nodes for up to two more jets in the selection.
Accordingly, $\mathcal{G}$ has 8–10 nodes. For each $i$, a vector
$\mathbf{a}_{i}$ of attributes with $\dim(\mathbf{a}_{i})=n,\,\forall i$ is
defined, forming the embedding space. At initialization time, the
$\mathbf{a}_{i}$ are initialized by the feature vectors $\mathbf{x}_{i}$ with
$\dim(\mathbf{x}_{i})=n_{\mathrm{feat}},\,\forall i$. The following studies
also comprise configurations with $n>n_{\mathrm{feat}}$. In these cases,
attributes in the $\mathbf{a}_{i}$ without correspondence in $\mathbf{x}_{i}$
are set to zero (zero padding), at initialization time.

All nodes are connected with edges, resulting in a fully-connected non-
directed graph without self-loops. Relational information between nodes $i$
and $j$ may be assigned, in form of edge weights $\omega_{ij}$. For this
purpose three physics motivated choices are made: (i) the invariant mass $m$;
(ii) the distance $\Delta R$; and (iii) the reciprocal distance $\Delta
R^{-1}$ of the connected final-state objects. In addition, the cases of a
(one) constant, (rnd) random, and (zero) no edge connection at all are
studied, resulting in a total of six variants of edge connections. An
illustration of a resulting graph for an event with eight nodes and no
additional jets is given in Figure 2.

Figure 2: Translation of an (upper part) selected event (without jets in the
NA jet-class in this case) into a (middle part) graph $\mathcal{G}$ and
finally into the (lower part) GNN model. For the indication of the partonic
final state we do not distinguish particles from anti-particles. The
individual object-classes are indicated by different colors. The nodes of
$\mathcal{G}$ are labelled by $i=1\ldots 8$ and colored the same way as the
object-classes. The boxes next to the nodes indicate the embedding space of
the GNN model. The GNN output is indicated by $\hat{y}$.

The GNN algorithm to process the graph data is based on the layered GraphConv
operation, as introduced and described in Ref. GraphConv , using the sum over
all $i$ as aggregation function. After initialization, $k$ GraphConv
operations are applied, after which the resulting $\mathbf{a}_{i}$ are
transformed into a single vector of length $n$, averaging over all $i$. A
linear combination of the components of this vector, which is scaled to values
between 1 (indicating $\mathrm{ttbb}$ as signal) and 0, for binary
classification, eventually forms the output $\hat{y}$ of the GNN. A graphical
illustration of this model is given in the lower part of Figure 2.

In a first study, the GNN are compared with corresponding DNN models with $k$
hidden layers, containing $n$ hidden nodes, each. The values of $k$, $n$, and
the choice of weights, steering the exploitation of relational information by
the GNN models are varied, resulting in 36 variants of parameter choices, as
summarized in Table 5. We note that here, as in the following, $n$ represents
a tuple of length $k$. For a GNN this tuple indicates the dimension of the
embedding space per GraphConv operation; for a DNN it represents the number of
nodes per hidden layer.

Table 5: Parameters varied for the comparison of GNN with corresponding DNN models, where $n$ corresponds to the dimension of the embedding space during a GraphConv operation (number of nodes in a hidden layer) and $k$ to the number of GraphConv operations (hidden layers) in the GNN (DNN) case. The choices of $n$ are motivated by the size of the input vector $n_{\mathrm{feat}}$ to the GNN, as described in Section 3.2. The choices of one, rnd, and zero for the use of edge information in GNN models are compared to DNN models without relational information between individual objects. These DNN models are indicated by the label none, in corresponding figures. Parameter | Setting(s)
---|---
$k$ | 1, 2
$n$ | 13, 26, 39
Edge weights ($\omega_{ij}$) |  $m$, $\Delta R$, $\Delta R^{-1}$, one, rnd, zero

Other parameter choices related to the NN architectures under study or the
setup of the NN training are made common and summarized in Table 4. Special
care is taken to compare the GNN with the corresponding DNN models on the same
footing, especially in terms of information about the feature space presented
to them, as discussed in the following section.

###  3.2 Presentation of the feature space

Primary features passed to the NNs are the invariant mass ($m$), energy ($E$),
$\eta$, and $\phi$ of each reconstructed final state object. The reconstructed
final state objects comprise $\mathrm{\ell}$, MET, at least four $\mathrm{b}$
jets, two additional jets, all of which are associated with the partonic final
state, and potentially two more non-associated jets from the event selection.
For MET the features $E$ and $m$ are set to zero. The primary features are
complemented by $\beta$. The selection requires values of $\beta\geq 2$ for
identified $\mathrm{b}$ jets. For $\mathrm{\ell}$ and MET $\beta$ is set to
zero.

We note that both NN architectures may profit from additional information,
which is not passed explicitly through $\mathbf{x}_{i}$, but implicitly
through the way the features are presented to the NNs. An obvious difference
between the NN architectures arises from the fact that the GNN naturally
supports processing of events with arbitrary object multiplicities. The
$\mathbf{x}_{i}$ are transferred to the GNN through the $\mathbf{a}_{i}$,
during initialization. During the GNN processing the information per $i$ is
aggregated over all nodes. In the DNN case such an aggregation step is absent.
Instead, the $\mathbf{x}_{i}$ are concatenated into an enlarged feature vector
$\mathbf{x}^{\mathrm{DNN}}$ of length $10\times n_{\mathrm{feat}}$, comprising
the $\mathbf{x}_{i}$ of the eight reconstructed objects, of which all jets
have been associated to the partonic final state, plus potentially two
additional selected jets. The order in which these objects are concatenated
has been chosen to follow the association to the partonic final state. For
events that contain fewer than two jets in addition to those that have been
matched to the partonic final state, the corresponding entries in
$\mathbf{x}^{\mathrm{DNN}}$ are filled with zeros. We point out two
subtleties, which are related to these choices:

One subtlety, in favor of the GNN, arises from incorporating relational
information through $\omega_{ij}$. This advantage is compensated for by
appending equivalent information to $\mathbf{x}^{\mathrm{DNN}}$. For up to ten
selected objects in an event this results in up to 45 additional features. The
choices of one, rnd, and zero for the use of edge information in GNN models
are compared to DNN models without relational information between individual
objects. These DNN models are indicated by the label none, in corresponding
figures.

Another subtlety, related to the same fundamental difference, but this time in
favor of the DNN, arises from the fact that through the concatenation of the
$\mathbf{x}_{i}$ into $\mathbf{x}^{\mathrm{DNN}}$, according to the
association to the partonic final state, the DNN receives extra information
through the positions of the $\mathbf{x}_{i}$ in $\mathbf{x}^{\mathrm{DNN}}$
which is not accessible to the GNN. This advantage is compensated for by
adding information about the association of the $i$ to the partonic final
state via one-hot encoding. For the five jet-classes defined in Table 2 plus
one label (LEP) for $\mathrm{\ell}$ and one label (MET) for MET this extension
increases $n_{\mathrm{feat}}$ by seven, leading to the dimension of
$\mathbf{x}_{i},\,\forall i$ of $n_{\mathrm{feat}}=13$.

For $\mathbf{x}^{\mathrm{DNN}}$ the $\mathbf{x}_{i}$ are concatenated for all
$i$ assuming two more jets in the NA class. For events with fewer than two
jets in the NA class the foreseen features are initialized with zero. Together
with the relational information between all potential objects in an event,
this results in a dimension of $\mathbf{x}^{\mathrm{DNN}}$ of 175 features, of
which up to 47 features might potentially be filled with zeros.

In this configuration the information about the association to the partonic
final state is presented in the form of one-hot encoding to both NN
architectures. To confirm to what extent the DNN may infer this information
already from the position of the $\mathbf{x}_{i}$ within
$\mathbf{x}^{\mathrm{DNN}}$ we also investigate configurations of the DNNs
without this information in form of on-hot encoding, resulting in a reduced
input vector $\mathbf{x}^{\mathrm{DNN}}_{\mathrm{red}}$ with dimension 105.
All input features in use are listed in Table 6.

Table 6: Input features used for the studies described in the text. Columns 2–4 indicate whether a given feature is continuous, discrete, or presented via one-hot encoding. The given features form the feature vectors $\mathbf{x}_{i}$ per object $i$. For the DNNs the $\mathbf{x}_{i}$ for potentially ten selected objects are concatenated into an extended feature vector $\mathbf{x}^{\mathrm{DNN}}$, according to their association to the partonic final state; for events with fewer than ten objects the $\mathbf{x}_{i}$ of the missing objects are filled with zeros. In addition, relational information between all potential objects is added to $\mathbf{x}^{\mathrm{DNN}}$. In a reduced configuration, the one-hot encoded information about the association of the objects to the partonic final state is omitted to form a reduced input vector $\mathbf{x}^{\mathrm{DNN}}_{\mathrm{red}}$ to the DNNs. Input feature | Continuous | Discrete | One-hot | Comment
---|---|---|---|---
$m$ | $\checkmark$ | $-$ | $-$ |  }6*[Primary features]
$E$ | $\checkmark$ | $-$ | $-$ |
$p_{\mathrm{T}}$ | $\checkmark$ | $-$ | $-$ |
$\phi$ | $\checkmark$ | $-$ | $-$ |
$\eta$ | $\checkmark$ | $-$ | $-$ |
$\beta$ | $-$ | $\checkmark$ | $-$ |
LEP | $-$ | $-$ | $\checkmark$ | Assoc. to $\mathrm{\ell}$
MET | $-$ | $-$ | $\checkmark$ | Assoc. to MET
ADDB | $-$ | $-$ | $\checkmark$ | Assoc. to additional $\mathrm{b}$ quarks
HTB | $-$ | $-$ | $\checkmark$ | Assoc. to $\mathrm{t}_{\mathrm{b}_{\mathrm{had}}}$
LTB | $-$ | $-$ | $\checkmark$ | Assoc. to $\mathrm{t}_{\mathrm{b}_{\mathrm{lep}}}$
HTQ | $-$ | $-$ | $\checkmark$ | Assoc. to $\mathrm{q}_{\mathrm{W}_{\mathrm{had}}}$
NA | $-$ | $-$ | $\checkmark$ | Not associated

All non-integer features are standardized to a mean of zero and a standard
deviation of one.

##  4 Results

To be able to draw fair conclusions from a comparison of different NN
architectures (of potentially different complexity) special care has to be
taken for this comparison to be based on the same ground. For this study we
have focused on a common choice of non-tunable (hyper-)parameters, which are
not subject to the NN training, as well as on an equal level of information
primarily passed to the NNs, through training conditions and input features.

An inevitable difference remains in the organization and layering of hidden
nodes, which when kept similar, may well lead to a different number of TPs and
therefore a priori different expressiveness of the NN models. Vice versa,
keeping the number of TPs similar, implies differences in the layering and
organization of hidden nodes. Since differences of one or the other kind may
not be overcome, both configurations, (i) similar layering of hidden nodes;
and (ii) similar number of TPs, are studied. In any case, the enlarged size of
$\mathbf{x}^{\mathrm{DNN}}$ ($\mathbf{x}^{\mathrm{DNN}}_{\mathrm{red}}$) with
respect to $\mathbf{x}$ will give higher emphasis to the first DNN layer
compared to the corresponding GNN architecture. In addition, the DNN
architecture features the less complex pre-processing of the inputs, since it
does not imply the creation of graphs. On the other hand, the potentially more
constrained GNN may have advantages over the DNN architecture in terms of
convergence properties of the training.

###  4.1 Comparison of neural networks with similar layering of hidden nodes

####  4.1.1 Neural networks with one layer of hidden nodes

A first comparison of GNN with corresponding DNN models, based on a similar
layering of hidden nodes, is shown in Figure 3.

Figure 3: Mean ROC-AUC $\mu_{\mathrm{AUC}}$ as obtained for 18 different
configurations of GNN and 24 corresponding configurations of DNN models with
$k=1$. The labels in brackets on the vertical axis indicate the use of
relational information, as discussed in Section 3.2, the numbers in
parentheses correspond to the choices of $n$. The circles refer to GNN and the
upward (downward) pointing triangles to DNN models with a default (reduced)
set of input features $\mathbf{x}^{\mathrm{DNN}}$
($\mathbf{x}^{\mathrm{DNN}}_{\mathrm{red}}$), as discussed in Section 3.2. For
better readability, markers of the same configuration are shifted vertically
along the y-axis. The bars are obtained from the sample variance of an
ensemble, as described in Section 2.4. Those NN architectures which belong to
the same choices of varying parameters are spatially grouped and shown in the
same color. Open markers indicate that significant outliers of the
corresponding distribution of ROC-AUC values have been removed from the
calculation of $\mu_{\mathrm{AUC}}$ and its variance, as described in the
text.

The metric by which to judge the success of an NN to fulfill the task is
chosen to be the mean of the ROC-AUC $\mu_{\mathrm{AUC}}$ based on the
training setup, as described in Section 2.4. The results are presented for the
variation of parameters as summarized in Table 5. For the presentation in
Figure 3 a simple architecture with one hidden layer for the DNN models and
one GraphConv operation of the GNN models ($k=1$) has been chosen. For the GNN
architecture this implies that there is only one exchange of information
across adjacent nodes of $\mathcal{G}$, i.e., each node receives information
only of its nearest and not the next to nearest neighbor in $\mathcal{G}$.
Since $\mathcal{G}$ has been chosen to be fully connected, there is no strict
suppression of information that way, in the sense that each node $i$ receives
information from any other node $j\neq i$ in $\mathcal{G}$. The input features
presented to each corresponding NN architecture are chosen, as described in
Section 3.2 and summarized in Table 6.

The results for the GNN models are represented by circles, the results for the
DNN models with $\mathbf{x}^{\mathrm{DNN}}$
($\mathbf{x}^{\mathrm{DNN}}_{\mathrm{red}}$) as input vector by upward
(downward) pointing triangles. The bars associated with the points indicate
the uncertainty $\Delta\mu_{\mathrm{AUC}}$ in $\mu_{\mathrm{AUC}}$ due to
random choices in the training, as described in Section 2.4. Open markers
indicate that at least one training repetition in an ensemble has been removed
as an outlier from the calculations of $\mu_{\mathrm{AUC}}$ and
$\Delta\mu_{\mathrm{AUC}}$. Candidates for outliers have been identified by
their $\mu_{\mathrm{AUC}}$ values exceeding $1.5\,\Delta\mu_{\mathrm{AUC}}$,
as obtained from the full ensemble. An outlier candidate has been definitely
removed and the ensemble size reduced accordingly, if doing so changed
$\Delta\mu_{\mathrm{AUC}}$ by a value of at least 0.0025. Following this
procedure, 54 outliers have been removed from a total of 1840, which
corresponds to a rate of 2.9%. Split by NN architectures, it corresponds to 17
(37) removed outliers for GNN (DNN) models from a total of 800 (1040). In no
case more than two outliers have been removed from the original ten training
repetitions of any individually model.

On the $x$-axis of the figure the corresponding values of
$\mu_{\mathrm{AUC}}$, ranging from 0.67 to 0.83, are shown. On the $y$-axis
the individual NN models are labelled, such that brackets indicate, what
relational information between final state objects has been used, and the
values in parentheses indicate the choices of $n$.

We conclude that all training setups have succeeded in the sense that all NN
models result in values of $\mu_{\mathrm{AUC}}>0.5$. The worst separation of
signal from background we observe for the GNN models for which no relational
information is exploited, indicated by three groups of NN architectures shown
in pink, purple and red colors in the lower part of the figure. It is noted
that the feature set zero refers to the case where the node convolution in the
GraphConv operation is forcefully suppressed, and deliberately no information
across nodes is exchanged at all. We anticipate that this approach counteracts
the whole GNN idea. We still keep this configuration as part of the study, to
gauge the effect and importance of the GraphConv operation itself. Compared to
the feature set zero, the feature sets one and rnd single out cases in which
information exchange across nodes takes place, but no real relational
information is associated with it. Instead, the embedding spaces of the
individual nodes are just mixed without particular prevalence. We note that
even when allowing node convolution the GNN architecture falls significantly
behind the comparable DNN architecture, even with a reduced input vector
$\mathbf{x}^{\mathrm{DNN}}_{\mathrm{red}}$, as long as no mindful relational
information across adjacent nodes $i$ and $j$ is provided according to
$\omega_{ij}$. This is true for assigning (one) the same or (rnd) random
weights to each edge, irrespective of the expressiveness of the NNs, indicated
by $n$. The superiority of the DNN architecture in this case cannot be
attributed to the additional information about the jet parton association,
since this information is available to all NN architectures under study. In
particular it is passed on to the GNNs, in the explicit form of one-hot
encoding, which is not even the case for the DNNs with
$\mathbf{x}^{\mathrm{DNN}}_{\mathrm{red}}$ as input vector. At this occasion,
we note that the additional (and in fact in this case redundant) information
of the parton association in form of one-hot encoding to the DNN does not lead
to a significant increase of $\mu_{\mathrm{AUC}}$, compared to the implicit
knowledge already provided by the positional information in
$\mathbf{x}^{\mathrm{DNN}}_{\mathrm{red}}$. Hence, if the way this information
is presented to the DNN were of influence, this influence is not significant
in the scope of our study. We also observe that $\Delta\mu_{\mathrm{AUC}}$ is
considerably larger for the DNNs.

In conclusion, if the advantage of the GNN over the DNN architecture were that
potentially excessive degrees of freedom in the DNN architecture are replaced
by built-in constraints, the GNN architecture appears too confined, until
these constraints are introduced mindfully. In turn the additional degrees of
freedom of the DNN architecture result in a larger spread
$\Delta\mu_{\mathrm{AUC}}$ due to random choices in the training setup.

The upper part of Figure 3 reveals that, as soon as a domain-knowledge
motivated ranking of information exchange across neighboring nodes is
introduced, the GNN architecture significantly gains in separation power. Also
here, this gain comes with an increase in $\Delta\mu_{\mathrm{AUC}}$. The
choices of (green) $\Delta R^{-1}$, (orange) $\Delta R$, and (blue) $m$ as
weights leads to an increase in separation power in the given order, where for
each choice the values of $\mu_{\mathrm{AUC}}$ can be grouped with a
corresponding internal spread. On the other hand, a significant gain in
$\mu_{\mathrm{AUC}}$, when increasing the expressiveness of an NN within a
certain configuration group in terms of $n$, for the utilized NN models, is
not observed.

We note that, as in the case of the positional encoding of parton association
in $\mathbf{x}^{\mathrm{DNN}}_{\mathrm{red}}$, all choices of relational
information are intrinsic to the training sample and implicitly accessible to
all NNs through their feature vectors. For $\Delta R$ this is, e.g., the case
through $\phi$ and $\eta$ in the primary features of each reconstructed
object. However, the information of $\Delta\phi$ and $\Delta\eta$ between
pairs of reconstructed objects appears too subtle in the high-dimensional
feature space, so that none of the chosen architectures could grasp it without
the assistance of an accordingly conditioned representation of $\mathbf{x}$
and $\mathbf{x}^{\mathrm{DNN}}$, even from a training sample with more than
200,000 events.

We further note that when turning the edge-weights of the GNN structure into
TPs we did not obtain a separation of signal from background better than the
domain-knowledge supported use of $m$. At the same time we observed a
significantly increased spread in the achieved separation power based on
random choices of the training setup.

Our physics prior assigns more physical meaning to the choice of $\Delta R$
over $\Delta R^{-1}$, since due to causality we expect a closer relation
between objects with smaller than larger spatial distance in $\Delta R$. The
observation that both choices of relational information lead to nearly similar
results in $\mu_{\mathrm{AUC}}$ we explain by a special characteristic of NN-
based classification tasks in the given setup. For the NN decision,
downgrading information from further-away objects is equivalent to upgrading
close-by objects. The fact that corresponding DNN architectures follow the
trends of the GNN architectures, as long as equipped with the same
information, supports the assumption that it is this additional relational
information between objects, rather than the GNN-specific operation of mixing
features across nodes that leads to the increase of $\mu_{\mathrm{AUC}}$.

We note that, consistently for all architectures, the highest values of
$\mu_{\mathrm{AUC}}$ are obtained with an energy-weighted distance measure
like $m$, which again follows our prior physics intuition. It is noteworthy
that in this case the DNN architecture with
$\mathbf{x}^{\mathrm{DNN}}_{\mathrm{red}}$ seems to significantly lose in
separation power, compared to the other architectures. In fact
$\mu_{\mathrm{AUC}}$ even decreases for increasing values of $n$. Also
$\Delta\mu_{\mathrm{AUC}}$ appears consistently higher for all DNN compared to
the GNN architectures. These observations may be interpreted as indications of
the advantage of careful guidance of the NN training- and model-setup over
just confronting a highly expressive NN architecture, represented by a large
number of TPs, with an even excessively large training sample. This guidance
may be provided through the choice and representation of input features, as
well as through the choice of a more constrained NN architecture. We note that
the highest value of $\mu_{\mathrm{AUC}}$ with the smallest spread
$\Delta\mu_{\mathrm{AUC}}$ is indeed obtained from the GNN architecture with
highest expressiveness, given for $n=39$.

####  4.1.2 Neural networks with two layers of hidden nodes

Moving on to an NN architecture with $k=2$ introduces the ambiguity of how to
choose $n$ for each individual layer. To prevent any kinds of potential
selection biases, all ways of allocating the tested values of $n=13,\,26,\,39$
to the individual layers/embedding spaces are shown in Figure 4.

Figure 4: Mean ROC-AUC $\mu_{\mathrm{AUC}}$ as obtained for 54 different GNN
and 72 corresponding DNN models with $k=2$. The labels in brackets on the
vertical axis indicate the use of relational information, as discussed in
Section 3.2, the numbers in parentheses correspond to $n$. The circles refer
to GNN and the upward (downward) pointing triangles to DNN models with a
default (reduced) set of input features $\mathbf{x}^{\mathrm{DNN}}$
($\mathbf{x}^{\mathrm{DNN}}_{\mathrm{red}}$), as discussed in Section 3.2. For
better readability, markers of the same configuration are shifted vertically
along the y-axis. The bars are obtained from the sample variance of an
ensemble, as described in Section 2.4. NN architectures which belong to the
same choices of varying parameters are spatially grouped and shown in the same
color. Open markers indicate that significant outliers of the corresponding
distribution of ROC-AUC values have been removed from the calculation of
$\mu_{\mathrm{AUC}}$ and its variance, as described in the text.

For the case of completely suppressed relational information for the GNN
architectures (zero) $\mu_{\mathrm{AUC}}$ remains lowest and unaffected by the
choice of $k$, as expected for a setup in which any information transfer
through a GaphConv operation is deliberately suppressed. At the same time we
generally observe that no superior choice of allocating $n$ across layers can
be pointed to, throughout all tested architecture configurations. Especially
the gain of using a configuration with $n=(39,\,39)$ for both hidden
layers/GraphConv operations over a configuration with $n=(13,\,13)$ within a
given architecture appears marginal.

However, the increase in $k$ consistently mitigates the previously observed,
clearly inferior separation of signal from background, of the GNN compared to
the corresponding DNN architectures, for randomly (rnd) and unweighted (one)
relational information. The $\mu_{\mathrm{AUC}}$ values of these groups of GNN
models start to clearly supersede the $\mu_{\mathrm{AUC}}$ values of the one-
layered DNN models with unweighted relational information, labeled by none in
Figure 3, even slightly taking the lead over the two-layered DNNs of the same
kind, in terms of $\mu_{\mathrm{AUC}}$.

We note that all two-layered NN architectures without use of relational
information still result in lower values of $\mu_{\mathrm{AUC}}$ than all
tested one-layered NN models that profit from the use of relational
information, as presented in Figure 3. At the same time, they are subject to
an increased spread $\Delta\mu_{\mathrm{AUC}}$ compared to their one-layered
counterparts, in most cases. In this sense, the wise choice of relational
information outweighs the presumable advantage in expressiveness provided by
$k=2$, irrespective of the allocated values of $n$.

For the NN architectures including relational information, we observe no
further, dramatic gains, with respect to their one-layered counterparts in
$\mu_{\mathrm{AUC}}$, apart from a slight advantage of the GNN over the
corresponding DNN architectures that seems to become more manifest. While this
advantage is below the 1%-level it is still significant compared to
$\Delta\mu_{\mathrm{AUC}}$.

Figure 5:  Summary of the achieved values of $\mu_{\mathrm{AUC}}$ for the GNN
models with (upper half) one and (lower half) two GraphConv operations, with
different use of relational information. For this summary, the associations of
$n$ with the highest values of $\mu_{\mathrm{AUC}}$ in each group of GNN
models have been used. The value of $\mu_{\mathrm{AUC}}$ is displayed on the
$x$-axis. Improvements relative to the least separating GNN with no relational
information at all (zero) is given in numbers to the right of the bars. The
use of relational information, as defined in Table 5, is indicated in
brackets, on the $y$-axis.  Figure 6:  Summary of the achieved values of
$\mu_{\mathrm{AUC}}$ for the DNN models with (upper half) one and (lower half)
two hidden layers, with different use of relational information. For this
summary, the associations of $n$ with the highest values of
$\mu_{\mathrm{AUC}}$ in each group of DNN models have been used. The DNN
configurations with $\mathbf{x}^{\mathrm{DNN}}$ and
$\mathbf{x}^{\mathrm{DNN}}_{\mathrm{red}}$ are shown separately. The values of
$\mu_{\mathrm{AUC}}$ are displayed on the $x$-axis. Improvements relative to
the least separating DNN with no relational information (none) is given in
numbers to the right of the bars. The use of relational information, as
defined in Table 5, is indicated in brackets, on the $y$-axis.

A summary of the achieved values of $\mu_{\mathrm{AUC}}$ for the one- and two-
layered GNN models is shown in Figure 5. An equivalent summary for the
corresponding DNN models is shown in Figure 6.

From the study we conclude that the external information of the $\omega_{ij}$
seems to give slight advantages to the GNNs with $k=2$. We note that two
subsequent GraphConv operations indeed convey more information than a DNN
model with two hidden layers. Viewing $\Delta R^{-1}$, $\Delta R$, and $m$ as
distance measures, the first GraphConv operation conveys information about the
nearest neighborhood of each $i$. The second GraphConv operation conveys
information about the nearest neighborhood of the nearest neighbors, which is
not the same as in the case of the first operation. This information is indeed
not primarily accessible to the DNN architectures, but it emerges from the
definition of the GraphConv operation. Along this line, once again, we
conclude that not the mixing of features across neighboring nodes $i$ during
the GraphConv operation, but the additional (implicit) information accessible
to the GNN through this operation is the source for the slight gain in
$\mu_{\mathrm{AUC}}$. This interpretation is supported by the observation that
an increase in $n$ does not lead to any significant improvements in
$\mu_{\mathrm{AUC}}$ despite the increase in expressiveness of the models.
This also explains why the GNN models with constant (one) or randomly
associated (rnd) weights suffer in their performance: Through both ways of
assigning weights across the nodes $i$ the information about any kind of
distance measure across nodes, in what ever space, is omitted.

A summary of the GNN and DNN configurations with the highest values of
$\mu_{\mathrm{AUC}}$ is given in Table 7.

Table 7:  Summary of $n$, number of TPs ($N_{\mathrm{TP}}$), and $\mu_{\mathrm{AUC}}$ of the (upper part) GNN and (lower part) DNN models with the highest results in $\mu_{\mathrm{AUC}}$. For the DNN models a number of effective TPs $N^{\mathrm{eff}}_{\mathrm{TP}}$, as defined in the text, is also given in parentheses. Corresponding summaries are given for the cases of $k=1$ and 2. In all cases $m$ has been used as relational information between nodes/physics objects. Also shown are the configurations of three additional NN models discussed in Section 4.2: the $\mathrm{GNN}^{\uparrow}$ model with $N_{\mathrm{TP}}$ comparable within 1% with $N_{\mathrm{TP}}(\mathrm{DNN}_{k=2})$ and the $\mathrm{DNN}^{\downarrow}$ ($\mathrm{DNN}^{\downarrow}_{\rm{eff}}$) model with $N_{\mathrm{TP}}$ ($N^{\mathrm{eff}}_{\mathrm{TP}}$) comparable within 1% with $N_{\mathrm{TP}}(\mathrm{GNN}_{k=2})$.  NN arch. | $n$ |  $N_{\mathrm{TP}}$ ($N^{\mathrm{eff}}_{\mathrm{TP}}$) | $\mu_{\mathrm{AUC}}$ | Label
---|---|---|---|---
{3*[GNN ] | (39) | 1093 |  $0.838\,701\,956\,542\,409\,4$$\pm$ | $0.000\,598\,750\,425\,9$ | $\mathrm{GNN}_{k=1}$
| (26, 39) | 2809 |  $0.848\,356\,936\,914\,419\,4$$\pm$ | $0.000\,772\,640\,704\,1$ | $\mathrm{GNN}_{k=2}$
| (29, 28, 29, 29) | 5829 |  $0.854\,602\,560\,655\,694$$\pm$ | $0.000\,577\,756\,711\,8$ | $\mathrm{GNN}^{\uparrow}$
{4*[DNN ] | (26) | 5773 (4753) |  $0.830\,037\,944\,444\,558\,9$$\pm$ | $0.001\,116\,829\,909$ | $\mathrm{DNN}_{k=1}$
| (26, 39) | 6839 (5819) |  $0.838\,765\,664\,034\,957\,2$$\pm$ | $0.000\,728\,252\,091\,4$ | $\mathrm{DNN}_{k=2}$
| (13, 14, 14) | 3294 (2784) |  $0.840\,039\,549\,541\,018\,7$$\pm$ | $0.000\,580\,934\,236\,8$ | $\mathrm{DNN}^{\downarrow}_{\rm{eff}}$
| (11, 10, 11, 11) | 2816 (2385) |  $0.838\,641\,371\,508\,374\,3$$\pm$ | $0.000\,741\,198\,589\,1$ | $\mathrm{DNN}^{\downarrow}$

###  4.2 Neural networks with comparable numbers of trainable parameters

As stated before, two principally different NN architectures lack
comparability in the sense that it might be more natural to pick up certain
information from the training sample through one or the other architecture. As
long as it addresses an intrinsic property of one or the other architecture,
such a difference is part of the benchmark comparison. If, on the other hand,
such an inequality results from withholding primary information from one or
the other architecture, or potentially inappropriate advantages in terms of
expressiveness, an effort should be made to study and estimate the effect of
it.

We have identified and noted such an inequality, in the beginning of Section
4, in terms of the number of TPs ($N_{\mathrm{TP}}$), which turns out to be
naturally higher for the DNN compared to the GNN architecture, due to the
usually much larger input layer. We therefore complete our study by three
additional setups, for which we drop the restrictions on the layered structure
of the trained NN models in favor of comparable numbers of TPs.

As shown in Table 7, the GNN model with the highest value of
$\mu_{\mathrm{AUC}}$ (labelled as $\mathrm{GNN}_{k=2}$) is based on $k=2$ and
$n=(26,\,39)$, with $N_{\mathrm{TP}}(\mathrm{GNN}_{k=2})=2809$ and a value of
$\mu_{\mathrm{AUC}}=0.8484\pm 0.0008$. The DNN with the highest result of
$\mu_{\mathrm{AUC}}$ (labelled as $\mathrm{DNN}_{k=2}$) is based on the same
configuration in terms of $k$ and $n$, with
$N_{\mathrm{TP}}(\mathrm{DNN}_{k=2})=6839$ and a value of
$\mu_{\mathrm{AUC}}=0.8388\pm 0.0007$. One may argue that for the DNN model,
not the full set of TPs is really actively contributing to the solution of the
task, since part of the input space is regularly filled with zeros, e.g., if
fewer than eight jets are selected in an event. Therefore, we estimate, in
addition to $N_{\mathrm{TP}}$, a number of effective TPs
($N^{\mathrm{eff}}_{\mathrm{TP}}$) from the product of $N_{\mathrm{TP}}$ with
the average number of nonzero input nodes in $\mathbf{x}^{\mathrm{DNN}}$,
evaluated on the training dataset. This results in a value of
$N^{\mathrm{eff}}_{\mathrm{TP}}(\mathrm{DNN}_{k=2})=5819$.

In a first approach we survey varying DNN structures with $N_{\mathrm{TP}}$
($N^{\mathrm{eff}}_{\mathrm{TP}}$) comparable to
$N_{\mathrm{TP}}(\mathrm{GNN}_{k=2})$. We do this based on the following
algorithm: We allow $k\leq 4$ and any number of nodes per hidden layer ($n$).
Of all DNN configurations for which $N_{\mathrm{TP}}(\mathrm{GNN}_{k=2})$ is
matched by $N_{\mathrm{TP}}$ ($N^{\mathrm{eff}}_{\mathrm{TP}}$) within a
margin of 1%, the model with the smallest spread of $n$ across layers is
selected. If no DNN configuration with $N_{\mathrm{TP}}$
($N^{\mathrm{eff}}_{\mathrm{TP}}$) within a 1% margin of the target value can
be found the closest possible model is chosen. This situation occurs only in
models with one hidden layer.

This procedure ensures a homogeneous structure of hidden layers. As a result,
e.g., a DNN configuration with $k=4$ and $n=(11,\,10,\,11,\,11)$, with
$N_{\mathrm{TP}}=2816$, further on referred to as $\mathrm{DNN}^{\downarrow}$,
is preferred over a model with $k=4$ and $n=(10,\,8,\,14,\,24)$, even though
the latter results in an exact match with
$N_{\mathrm{TP}}(\mathrm{GNN}_{k=2})$. The $\mathrm{DNN}^{\downarrow}$ model
achieves a value of $\mu_{\mathrm{AUC}}=0.8386\pm 0.0007$. A second DNN
configuration with $k=3$ and $n=(13,\,14,\,14)$, with
$N^{\mathrm{eff}}_{\mathrm{TP}}=2784$ within the 1% margin of
$N_{\mathrm{TP}}(\mathrm{GNN}_{k=2})$ is also considered and further on
referred to as $\mathrm{DNN}^{\downarrow}_{\rm{eff}}$. This model achieves a
value of $\mu_{\mathrm{AUC}}=0.8400\pm 0.0006$. We observe that, although the
$\mathrm{DNN}^{\downarrow}$ ($\mathrm{DNN}^{\downarrow}_{\rm{eff}}$) model
uses only 41% (48%) of $N_{\mathrm{TP}}$ ($N^{\mathrm{eff}}_{\mathrm{TP}}$) of
$\mathrm{DNN}_{k=2}$, this does not result in any significant loss in
separation power, after training.

In a second approach we survey varying GNN structures with $N_{\mathrm{TP}}$
comparable to $N^{\mathrm{eff}}_{\mathrm{TP}}(\mathrm{DNN}_{k=2})$. For this
purpose we exploit the same algorithm as described above, resulting in a GNN
with $k=4$ and $n=(29,\,28,\,29,\,29)$ with $N_{\mathrm{TP}}=5829$, further on
referred to as $\mathrm{GNN}^{\uparrow}$. This model achieves a value of
$\mu_{\mathrm{AUC}}=0.8546\pm 0.0006$. It reveals the highest value of
$\mu_{\mathrm{AUC}}$ across all tested models. The difference in
$\mu_{\mathrm{AUC}}$ with respect to other NN configurations is only at the
1%-level, but it is still significant in terms of $\Delta\mu_{\mathrm{AUC}}$.

From this finding we conclude that the GNN model with the same expressiveness
as a maximally comparable DNN must have intrinsic advantages over the DNN
model in extracting additional information from the given, large training
sample. For the benchmark setup in use, this advantage is small but
significant in the scope of the study. It emerges after external augmentation
with an energy-weighted distance measure like $m$ between the input
objects/nodes, and more clearly manifests itself in the study for $k>1$. We
anticipate that this gain originates from the hierarchically structured
information about nearest neighbors and the nearest neighborhood of nearest
neighbors of node $i$, when viewing the edge weights $\omega_{ij}$ as a
distance measure. This information is an intrinsic property of the GNN model
and not easily accessible through the more simplistic DNN structure. An
increase of TPs of the simpler DNN structure does not compensate for this
informational advantage. In this interpretation the gain of
$\mathrm{GNN}^{\uparrow}$ over all other configurations should mostly be
attributed to the increase in $k$ over the association of $n$ per GraphConv
operation. The parameter choices of the $\mathrm{DNN}^{\downarrow}$,
$\mathrm{DNN}^{\downarrow}_{\rm{eff}}$, and $\mathrm{GNN}^{\uparrow}$ models
and correspondingly achieved values of $\mu_{\mathrm{AUC}}$ are also given in
Table 7.

###  4.3 Convergence behavior

In this study we have investigated the capacities of GNN and DNN models to
fulfill their primary target, i.e., to provide the best possible solutions to
the classification task defined in Section 2.3. We anticipate that, in
particular in practical life, the properties of an NN architecture may be
evaluated in other terms, viz. the mean of the training speed
$\mu_{\mathrm{Train}}$, which we evaluate as the inverse of the epoch with the
highest value of the ROC-AUC on the validation dataset and the mean of the
empirical risk obtained from the test dataset $\mu_{R}$, which we take as a
measure of the generalization property of the given NN model under study. To
conclude our studies we provide a visualization of these properties and all
other properties of the NN models that have been discussed throughout the
paper so far, in Figure 7.

Figure 7:  Visualization of GNN and DNN models with the highest values in
$\mu_{\mathrm{AUC}}$ and $k=1$ ($\mathrm{GNN}_{k=1}$, $\mathrm{DNN}_{k=1}$)
and $k=2$ ($\mathrm{GNN}_{k=2}$, $\mathrm{DNN}_{k=2}$). Also shown are the DNN
models $\mathrm{DNN}^{\downarrow}$ and $\mathrm{DNN}^{\downarrow}_{\rm{eff}}$
with a comparable number of (effective) TPs as the $\mathrm{GNN}_{k=2}$ model,
and the GNN model $\mathrm{GNN}^{\uparrow}$ with a comparable number of
effective TPs as the $\mathrm{DNN}_{k=2}$ model. Five different metrics to
evaluate the properties of an NN model are shown: $N_{\mathrm{TP}}$
($N^{\mathrm{eff}}_{\mathrm{TP}}$), the mean convergence rate
$\mu_{\mathrm{Train}}$, the mean empirical risk evaluated on the test sample
$\mu_{R}$, and $\mu_{\mathrm{AUC}}$. The spanned area in the figures indicates
the capability of an NN model to fulfill the task.

In Figure 7 (left) $\mu_{\mathrm{Train}}$, $\mu_{R}$, $\mu_{\mathrm{AUC}}$,
$N_{\mathrm{TP}}$, and $N^{\mathrm{eff}}_{\mathrm{TP}}$ for the
$\mathrm{DNN}_{k=1}$, $\mathrm{DNN}_{k=2}$, $\mathrm{DNN}^{\downarrow}$, and
$\mathrm{DNN}^{\downarrow}_{\rm{eff}}$ models are shown, on five independent
axes. The axes are defined such that values closer to the common origin of the
figure are disfavored. This is in particular true for $N_{\mathrm{TP}}$,
$N^{\mathrm{eff}}_{\mathrm{TP}}$, and $\mu_{R}$, where the values are given in
descending order when moving away from the origin. In this sense a larger size
of the correspondingly colored area indicates the larger capability of a given
NN model to adequately solve the presented task. In Figure 7 (right) the same
quantities are shown for the $\mathrm{GNN}_{k=1}$, $\mathrm{GNN}_{k=2}$, and
$\mathrm{GNN}^{\uparrow}$ models. The axes ranges are kept the same to ease
comparison between both architectures.

In terms of $\mu_{\mathrm{Train}}$ the DNN models usually fall behind their
corresponding GNN counterparts. This finding, as well as the observation that
the GNN models usually achieve a comparable or slightly larger value of
$\mu_{\mathrm{AUC}}$, after training, indicates the assumed effect of guiding
the convergence by constraints, which are built-in to the GNN architecture.
This property allows the GNN architecture not only to converge to a solution
that is equally good or slightly better than the solution found by
corresponding DNN architectures, but also to converge faster and typically
with a smaller number of TPs. We note though that a clear correlation between
$\mu_{\mathrm{Train}}$ and $N_{\mathrm{TP}}$
($N^{\mathrm{eff}}_{\mathrm{TP}}$) cannot be deduced from our study. We
understand this situation such that a less expressive NN model, with fewer
degrees of freedom for the minimization process, may well lead to a more
pronounced landscape of the expected risk and thus reduced
$\mu_{\mathrm{Train}}$. In addition we note that also $\mu_{R}$ for the DNN
models falls behind, compared to their GNN counterparts. Here we observe the
benefit of a regularizing effect of the built-in constraints, which correlates
with reduced numbers of (effective) TPs. It is obvious that an NN model with
more TPs reveals a higher vulnerability to specific properties of the training
sample. A summary of all NN configurations that have been discussed in this
section is shown in Figure 8.

Figure 8:  Overview of the GNN and DNN models with $k=1$ and 2, and the
highest achieved values of $\mu_{\mathrm{AUC}}$, as well as a GNN model with
an increased number of TPs ($\mathrm{GNN}^{\uparrow}$) and a DNN model
$\mathrm{DNN}^{\downarrow}$ ($\mathrm{DNN}^{\downarrow}_{\rm{eff}}$) with a
restricted (effective) number of TPs. The GNN models are indicated by circles
and the DNN models by triangles. The models are shown in a three-dimensional
space built from $\mu_{\mathrm{AUC}}$, $N_{\mathrm{TP}}$ and
$\mu_{\mathrm{Train}}$. The bars in $\mu_{\mathrm{AUC}}$ are obtained from the
sample variance of a training ensemble, as described in Section 2.4. The
quantity $\mu_{\mathrm{Train}}$, indicated by the color code of the points
corresponds to the inverse of the epoch with the highest value of the ROC-AUC
on the validation dataset.

##  5 Summary

With this paper we have made an effort to put the comparison of graph neural
network (GNN) and equivalent fully-connected feed-forward neural network (DNN)
architectures on a maximally fair ground. Under the laboratory conditions of a
high-energy physics process of interest, at the CERN LHC, we have controlled
the definition of the task, choice of non-tunable (hyper-)parameters of the
models, which are not subject to the training, and amount of information
primarily presented to the neural network models, through their input feature
vectors.

Within the scope of the study we have demonstrated clear evidence that the
presumable advantage of the more complex GNN over an equivalent DNN structure
does not originate from an uncontrolled mixture of features in the embedding
space of the graph nodes, but from the extra relational information between
nodes, that we have added based on domain-knowledge. Without this extra
knowledge the GNN models fall behind equivalent DNN models in terms of their
capability to separate signal from background. Neither are the GNN models
superior in terms of their separation power to equivalent DNN models, as long
as these are equipped with the same information in the input space.

Both, the built-in permutation invariance and the circumstance that the GNN
architecture a priori is not bound to a fixed number of nodes might be viewed
as advantages. They might also have positive influence on the convergence
behavior of the training. On the other hand it cannot be deduced that either
of these properties significantly contributes to an increase, e.g., in the
power to separate a given signal from background. Any advantage of the GNN
over the DNN architecture that we observed in our studies could be traced back
to the access to more information, which, when given to the other architecture
lead to the same performance also for the other architecture.

The real advantage of the GNN over the DNN structure emerges as soon as more
than one GraphConv operations are applied, during which the GNN structure
naturally accesses more relational information between nodes than a DNN has
access to. In the configurations investigated in our study this gain is tied
to the use of relational information that can be interpreted as a distance
measure to define proximity between two nodes. Apart from that we observe
significant advantages of the GNN over the DNN architecture in terms of
convergence and generalizability that we attribute to a level of built-in
implicit constraints to the GNN model resulting in a better ratio of
accessible information over trainable parameters of the model. We anticipate
that these advantages might be more pronounced the more hierarchical the
training data are. We assume that this property is the basis for the success
of GNN structures when applied to particle physics jets, which are highly
hierarchical objects. In conclusion we expect the highest gain of a GNN over a
DNN structure for tasks based on hierarchically structured data, ideally based
a known distance measure.

## ACKNOWLEDGEMENTS

This research was supported by the German Federal Ministry of Education and
Research (BMBF) under grant 05H21VKCCB.

## References

  * [1] D. Aguillard, T. Albahri, et al.  Measurement of the Positive Muon Anomalous Magnetic Moment to 0.20 ppm.  Phys. Rev. Lett., 131(16):161802, 2023.
  * [2] N. Metropolis, A. Rosenbluth, et al.  Equation of State Calculations by Fast Computing Machines.  The Journal of Chemical Physics, 21(6):1087, 2004.
  * [3] M. Rosenbluth.  Genesis of the Monte Carlo Algorithm for Statistical Mechanics.  AIP Conference Proceedings, 690(1):22, 2003.
  * [4] G. Aad, E. Abat, et al.  The ATLAS Experiment at the CERN Large Hadron Collider.  JINST, 3:S08003, 2008.
  * [5] S. Chatrchyan, G. Hmayakyan, et al.  The CMS Experiment at the CERN LHC.  JINST, 3:S08004, 2008.
  * [6] G. Salam.  Towards Jetography.  Eur.Phys.J., C67:637–686, 2010.
  * [7] A. Sirunyan, A. Tumasyan, et al.  Performance of missing transverse momentum reconstruction in proton-proton collisions at $\sqrt{s}=$ 13 TeV using the CMS detector.  JINST, 14(07):P07004, 2019.
  * [8] R. Feynman.  The theory of positrons.  Phys. Rev., 76:749–759, Sep 1949.
  * [9] F. Scarselli, M. Gori, et al.  The graph neural network model.  IEEE Transactions on Neural Networks, 20(1):61–80, 2009.
  * [10] A. Micheli.  Neural network for graphs: A contextual constructive approach.  IEEE Transactions on Neural Networks, 20(3):498–511, 2009.
  * [11] T. Kipf and M. Welling.  Semi-Supervised Classification with Graph Convolutional Networks.  arXiv e-prints, page arXiv:1609.02907, September 2016.
  * [12] W. Hamilton, R. Ying, et al.  Inductive representation learning on large graphs.  In Proceedings of the 31st International Conference on Neural Information Processing Systems, NIPS’17, page 1025–1035, Red Hook, NY, USA, 2017. Curran Associates Inc.
  * [13] P. Veličković, G. Cucurull, et al.  Graph Attention Networks.  International Conference on Learning Representations, 2018.
  * [14] J. Shlomi, P. Battaglia, et al.  Graph Neural Networks in Particle Physics.  7 2020.
  * [15] A. Sirunyan, A. Tumasyan, et al.  Search for $\mathrm{t}\overline{\mathrm{t}}\mathrm{H}$ production in the $\mathrm{H}\to\mathrm{b}\overline{\mathrm{b}}$ decay channel with leptonic $\mathrm{t}\overline{\mathrm{t}}$ decays in proton-proton collisions at $\sqrt{s}=13$ TeV.  JHEP, 03:026, 2019.
  * [16] R. Workman, V. Burkert, et al.  Review of Particle Physics.  PTEP, 2022:083C01, 2022.
  * [17] A. Sirunyan, A. Tumasyan, et al.  Identification of heavy-flavour jets with the CMS detector in pp collisions at 13 TeV.  JINST, 13(05):P05011, 2018.
  * [18] J. Alwall, R. Frederix, et al.  The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations.  JHEP, 07:079, 2014.
  * [19] J. Alwall, M. Herquet, et al.  MadGraph 5 : Going Beyond.  JHEP, 06:128, 2011.
  * [20] T. Sjöstrand, S. Ask, et al.  An introduction to PYTHIA 8.2.  Comput. Phys. Commun., 191:159–177, 2015.
  * [21] J. de Favereau, C. Delaere, et al.  DELPHES 3, A modular framework for fast simulation of a generic collider experiment.  JHEP, 02:057, 2014.
  * [22] M. Cacciari, G. Salam, et al.  The anti-$k_{t}$ jet clustering algorithm.  JHEP, 04:063, 2008.
  * [23] M. Fey and J. Lenssen.  Fast graph representation learning with pytorch geometric, 2019.
  * [24] A. Paszke, S. Gross, et al.  Pytorch: An imperative style, high-performance deep learning library, 2019\.
  * [25] D. Kingma and J. Ba.  Adam: A Method for Stochastic Optimization.  12 2014.
  * [26] C. Morris, M. Ritzert, et al.  Weisfeiler and leman go neural: Higher-order graph neural networks.  Proceedings of the AAAI Conference on Artificial Intelligence, 33(01):4602–4609, Jul. 2019.