aainstitutetext: High Energy Physics Division, Argonne National Laboratory,
Lemont, IL 60439bbinstitutetext: Department of Physics, University of Arizona,
Tucson, AZ 85721

# Explainable AI classification for parton density theory

Brandon Kriesten a,b , Jonathan Gomprecht a , T. J. Hobbs<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Quantitatively connecting properties of parton distribution functions (PDFs,
or parton densities) to the theoretical assumptions made within the QCD
analyses which produce them has been a longstanding problem in HEP
phenomenology. To confront this challenge, we introduce an ML-based
explainability framework, XAI4PDF, to classify PDFs by parton flavor or
underlying theoretical model using ResNet-like neural networks (NNs). By
leveraging the differentiable nature of ResNet models, this approach deploys
guided backpropagation to dissect relevant features of fitted PDFs,
identifying $x$-dependent signatures of PDFs important to the ML model
classifications. By applying our framework, we are able to sort PDFs according
to the analysis which produced them while constructing quantitative, human-
readable maps locating the $x$ regions most affected by the internal theory
assumptions going into each analysis. This technique expands the toolkit
available to PDF analysis and adjacent particle phenomenology while pointing
to promising generalizations.

###### Keywords:

PDFs; AI; Collider Phenomenology; Explainability; Machine Learning; QCD

††preprint: ANL-189807

##  1 Introduction

Parton distribution functions (PDFs) are an indispensable input informing
standard model (SM) theory predictions for high-energy experiments. In the
case of ongoing measurements at the Large Hadron Collider (LHC), its planned
high-luminosity upgrade (HL-LHC) Apollinari:2017lan ; Cepeda:2019klc , and the
upcoming Electron-Ion Collider (EIC) Accardi:2012qut , the collinear
unpolarized PDFs of the proton play a central role, given their importance in
determining standard model (SM) baselines for numerous processes while
carrying signatures of the still imperfectly understood dynamics of QCD. In
many instances, control over the PDFs limits measurements at colliders both in
their sensitivity to new physics beyond the standard model (BSM) as well as
the statistical rigor of associated precision tests of the SM ZEUS:2019cou ;
Carrazza:2019sec ; Greljo:2021kvv ; Madigan:2021uho ; Gao:2022srd ;
Iranipour:2022iak . This reality has motivated a worldwide effort
Alekhin:2017kpj ; Barry:2018ort ; Hou:2019efy ; Bailey:2020ooq ; NNPDF:2021njg
; PDF4LHCWorkingGroup:2022cjn ; Kotz:2023pbu  to model the PDFs in a
systematically-improving fashion while continuously updating extractions in
light of a steadily growing set of high-energy data recorded at the LHC and
other facilities.

By their comprehensive nature, modern PDF fits entail many internal pieces
Kovarik:2019xvh  — among them, aspects related to the deployment of
perturbative QCD and electroweak theory; a host of nonperturbative effects
essential to the modeling of highly sensitive data (e.g., the inclusion of
non-leading twist or nuclear corrections); the curation and treatment of
hadronic data and their accompanying (correlated) systematic and statistical
uncertainties; methodological choices related to analysis procedures and
workflows in practical PDF fits; and open issues in statistical theory which
become unavoidable as the parametric complexity of PDF analysis continues to
expand. Notably, collinear PDF analyses take place within a still larger
landscape of quantum correlation functions (QCFs) of higher dimensionality,
including transverse-momentum dependent (TMD) PDFs as well as generalized
parton distribution functions (GPDs) Ji:1996ek ; Radyushkin:1997ki ; these are
central to studies of the nucleon’s 3D structure, involving global data sets
including deeply virtual Compton scattering and related experiments planned at
Jefferson Lab, the EIC, and beyond Ji:1996nm ; Belitsky:2001ns ;
Kriesten:2019jep ; Kriesten:2020wcx ; Guo:2022cgq ; Accardi:2023chb ;
Qiu:2022pla ; Qiu:2022bpq . While there is some success in modeling these
functions Kriesten:2021sqc ; Bertone:2021wib ; Guo:2023ahv ; Moffat:2023svr ,
it remains an open question in how to perform a global analysis with full
consideration of all theoretical and experimental constraints.

QCD global analyses have progressed to a sufficient level that an era of PDF-
driven phenomenology with $\sim$1% precision is increasingly plausible
Amoroso:2022eow  for many high-energy observables of interest. At the same
time, however, increasing working accuracy in PDF fits and the accumulation of
high-energy data are such that variations in analysis assumptions can produce
fluctuations in likelihood functions and/or extracted PDFs which had
previously been minor relative to the comparatively larger uncertainties which
typified previous-generation analyses. For instance, improved control over
scale uncertainties achieved at NNLO relative to NLO is such that inclusion of
few-percent corrections in other sectors of the global analysis — for example,
related to light nuclear corrections relevant for key data sets
Accardi:2021ysh  — may induce $\chi^{2}$ variations of similar magnitude to
those obtained from shifts in the assumed NNLO formalism. As a consequence, it
is increasingly necessary to track a significant number of internally coupled
degrees-of-freedom in modern PDF fits, beyond selections of data sets or PDF
parametrization forms.

Extracting clear insights regarding how theory or analysis assumptions made in
a given fitting framework influence the behavior of the PDFs they ultimately
determine remains the subject of ongoing research. Historically, the most
straightforward approach to assessing the relationship between analysis
settings and fitted PDFs was the construction of arrays of fits in serial
fashion. While simple to execute, this method is both computationally costly
and carries the risk of obscuring subtleties owing to the relative smallness
of effects involved as well as the inherent difficulty of separating the
combined impacts of simultaneous variations in analysis settings; these latter
effects may interfere both constructively or destructively in surprising ways,
potentially conspiring to obscure the effect of varying one internal parameter
against another. In the meantime, the growing complexity of QCD analyses as
noted above greatly heightens the challenge of interpreting the origin of
specific features observed in the fitted PDFs they produce.

To ameliorate this problem, recent efforts have relied on detailed
benchmarking of PDF fits, with a strong focus on achieving a common basis on
which to compare global fits as in the PDF4LHC PDF4LHCWorkingGroup:2022cjn ;
Butterworth:2015oua ; Ball:2012wy ; Botje:2011sn  exercises. Similarly, a
number of methods based in classical statistical theory have been developed
based on metrics to quantify the statistical pulls on PDFs derived from
variations in the $\chi^{2}$ or benchmark comparisons to higher-level
calculated physics processes. For example, Ref. Hobbs:2019gob  proposed the
notion of $L_{2}$ sensitivity, discussed in slightly more detail in Sec. 5.2
and Ref. Jing:2023isu , which offers a quantitative metric to assess how
shifting PDFs within their $1\sigma$ uncertainties induces corresponding
shifts in the $\chi^{2}_{E}$ of key experiments constraining the global fit.
This method furnishes $x$-dependent information on the PDF-level pulls of the
most sensitive data sets, and builds upon the ($L_{1}$) PDF Sensitivity,
$|S_{f}|(x,Q)$, introduced in Ref. Wang:2018heo  to quantify and visualize the
projected PDF sensitivities of data sets within the context of a baseline
assumed fit. Both approaches are principally developed within the Hessian
global fit method and contributed to the analysis in the most recent PDF4LHC21
study PDF4LHCWorkingGroup:2022cjn .

Motivated in part by the growing size of global fits and use of neural
networks (NNs) in PDF parametrization and related analysis elements, a number
of recent studies have implemented ML methods, including generative models, to
understand properties of PDFs and their uncertainties; these have included
inverse mapper solutions with autoencoders 9534012 ; Kriesten:2023uoi ,
Gaussian mixture models Yan:2024yir , and GANs Carrazza:2021hny , which have
generally explored efficient representations or reconstructions of PDFs as
well as studies of their associated uncertainties. These approaches offer
parallel and complementary avenues to exploring the complicated correlations
between PDFs and their internal analysis elements. In Ref. Kriesten:2023uoi ,
for example, inverse problem solvers were used to reconstruct PDFs from their
Mellin-space behavior (i.e., their lattice-calculable moments); the resulting
interpretable framework could be used as a generative model to predict PDFs,
interpolate among encoded $x$-dependent shapes, or explore connections between
PDFs and higher-dimensional distributions like GPDs. In this study, the
initial input space of (variational) autoencoders was informed by $x$-gridded
arrays of PDFs obtained from toy-model parametrizations. A natural extension
of these approaches involves training feed-forward networks on
phenomenological PDFs themselves as in Ref. Liu:2022plj , an approach which
allows a variety of fast studies, including rapid exploration of combined PDF-
EFT theory spaces Gao:2022srd . In addition, some groups have implemented
machine learning (ML) methods Kumericki:2016ehc ; Cuic:2020iwt ;
Almaeen:2022imx ; Almaeen:2024guo  to understand how physics constraints such
as PDF limits can be incorporated into full fits of the GPDs.

Although these foregoing studies enlarge the toolkit of PDF analysis methods,
they have largely focused on the task of encapsulating the PDFs themselves in
NN parametrizations and/or training to a likelihood function, rather than
attempting to disentangle the role model assumptions play in fitted PDFs, let
alone doing so in a common setting which simultaneously admits collections of
global fits. Moreover, previous investigations were more immediately related
to tasks of PDF reconstruction or regression; meanwhile, an adjacent problem
is that of performing model discrimination among distinct families of parton
densities — a problem we recognize as inherently one of classification.
Motivated by this problem as well as the need to diagnose the PDF-level
impacts of assumed theory in an apple-to-apples framework, in the present work
we propose a new, generalizable scheme based on the ML modeling of fitted
phenomenological PDFs. Specifically, we concentrate on model discrimination
through classification111In a forthcoming study XAIreg , we revisit the
question of regression within an explainable AI framework. within a closed set
of PDF analyses, as this provides a means of cleanly distinguishing among
unique phenomenological fits (when statistically possible) and applying
explainability methods to trace the ability of ML models to identify
phenomenological fits based on unique features of the predicted PDFs. This
technique, which we call XAI4PDF, enables us in turn to dissect the dependence
of the fitted PDFs on elements of the assumed theory by leveraging recent
developments in explainable AI (XAI). In addition, by structuring our model
around the logic of QCD global fits and thereby integrating domain knowledge,
the resulting approach helps bridge the gap between ML-driven modeling and
traditional methods for QCD analysis.

Using a ResNet-like NN architecture optimized for classification, we identify
PDFs from the CTEQ-TEA (CT) family of global fits — specifically, the most
recent NNLO main release, CT18 Hou:2019efy . For our proof-of-principle
classification tasks, we concentrate on sorting parton densities according to
their associated parton flavor and/or the specific set of global analysis
assumptions which produced them. Saliency maps are constructed using gradients
calculated by a modified autodifferentiation algorithm called guided
backpropagation. These saliency maps trace the impact of theoretical
assumptions and the contribution of significant experiments to the $x$ and
flavor dependence of the PDFs. The procedure described here also introduces a
notion of analysis simultaneity — rather than serially comparing one pair of
analysis settings and corresponding fitted PDFs, it becomes possible to train
a space of PDF analyses of essentially arbitrary size while classifying the
results within a consistent (ML) model or framework.

Our paper is organized as follows: in Sec. 2, we discuss the phenomenological
and theoretical background that motivate this work, particularly related to
PDFs. We explain the guided backpropagation technique as related to
classification algorithms in Sec. 3. In Sec. 4 we describe the model
architecture and PDF data that we use to train. Sec. 5 is where we provide
details of our application of the XAI4PDF framework in two example cases for
PDFs: in Sec. 5.1 we use classifying CT18 NNLO PDF replicas by flavor, and in
Sec. 5.2 we classify PDF replicas by phenomenological fit with underlying
theory assumptions. Finally, we summarize conclusions and possible future
directions in Sec. 6.

##  2 Parton densities: background theory and phenomenology

PDFs are inherently probabilistic objects providing the likelihood of finding
a parton of flavor, $f$, carrying some fractional momentum of the proton, $x$,
at a given hard scale, $Q$. BSM searches now engage a variety of channels,
such as inclusive and coincident Higgs production; processes involving
$W^{\pm}$/$Z$ gauge bosons ATLAS:2016nqi ; and various top-production
Ablat:2023tiy  channels. In many cases, these experiments not only provide
fertile BSM search contexts, but also constraining data on hadronic structure
which may be actively included in PDF global analyses. Predictions for
collider processes are made on the basis of QCD factorization theorems
Collins:1985ue ; Collins:2011zzd ; CTEQ:1993hwr , which render the PDF
dependence of theoretical cross sections explicit. In a simple but ubiquitous
example, the single-inclusive production of an electroweak boson in $pp$
scattering at the LHC, $p_{1}+p_{2}\to W^{\pm}/Z+X$, is given schematically by
the convolution of a short-distance matrix element, $\hat{\sigma}_{ij}$, with
unpolarized proton PDFs as

${\sigma^{pp\to W/Z}\over dy_{W/Z}}\ =\ \sum_{i,j}\int
dx_{1}dx_{2}\,f_{i}(x_{1})f_{j}(x_{2})\,{d\hat{\sigma}^{W/Z}_{ij}\over
dy_{W/Z}}(x_{1},x_{2})\ .$ (1)

The parton-level arguments entering hard matrix elements and PDFs in the
expression above are typically evaluated via Born-level kinematical matchings,
such that, for the gauge-boson production scenario above, the relevant parton
fractions and hard scales are

$x_{1},x_{2}\sim{M_{W/Z}\over\sqrt{s}}\,\exp{\pm y_{W/Z}}\ ,\ \ \ \ \ Q\sim
M_{W/Z}\ ,$ (2)

respectively. We note that this matching is subject to higher-order
corrections motivating the inclusion of greater theoretical working accuracy
in both $\alpha_{s}$ and the electroweak interaction strength. In particular,
from Eq. (2) we conclude that rapidity-dependent cross sections,
$d\sigma/dy_{W/Z}$, involve an $x$ dependence spanning well beyond the
nominally expected low-$x$ region, with forward-boosted electroweak events
probing the high-$x$ behavior of the PDFs, in addition to possible suppressed
BSM contributions. This reality underscores the fundamental challenge of PDF
phenomenology for high-stakes HEP observables: precision control for
quantities like EW boson rapidity or invariant-mass distributions requires
detailed knowledge of the PDFs over a wide range of $x$ and a multiplicity of
parton-level channels, with combinations like
$\sim\\!\bar{u}(x_{1})u(x_{2}),\dots,\bar{s}(x_{1})s(x_{2})$ relevant for $Z$,
or $\sim\\!\bar{d}(x_{1})u(x_{2}),\dots,\bar{s}(x_{1})c(x_{2})$ implicated in
$W^{+}$ production. Moreover, these PDFs must be DGLAP-evolved Gribov:1972ri ;
Dokshitzer:1977sg ; Altarelli:1977zs  to the appropriate hard scale, $Q$, for
a given experiment, a process which itself introduces complicated
entanglements among the PDFs’ dependence on flavor, $x$, and $Q$. These
challenges are such that the only currently-available path to determine PDFs
for precision HEP phenomenology is performing comprehensive QCD global
analyses of the world’s hadronic data.

Phenomenological PDF analyses have an established history dating to the early
1990s, with a variety of methods blossoming in the ensuing decades. The logic
of QCD global analysis is commonly shared across these approaches as motivated
by the above discussion, with parameters fitted using data from a diverse set
of hadronic processes, including deep inelastic scattering (DIS), Drell-Yan
(dilepton) production, top-quark processes, and $W^{\pm}$/$Z$ gauge-boson
production as noted around Eq. (1). State-of-the-art implementations of
perturbative calculations are currently at the order of NNLO in the strong
coupling, $\alpha_{s}$, although there is ongoing work on incorporating
specific analysis elements at N3LO to produce approximate aN3LO fits
McGowan:2022nag ; Cooper-Sarkar:2024crx ; NNPDF:2024nan ; there is not yet
sufficient progress, however, to obtain a comprehensive fit at true N3LO.
Additionally, there are many nuanced theoretical considerations for a full
phenomenological fit which we will not address in this study.

Ultimately, despite the considerable variation in analysis settings implied
above by the wide range of data sets and theoretical choices available to QCD
fits, the great preponderance of global analyses falls into one of two main
categories: those involving explicit analytical (functional) forms to
parametrize the $x$ dependence of the PDFs, and methods involving (feed-
forward) neural networks (NNs) NNPDF:2021njg  as the basis for PDF
parametrization. Effectively, these parametrizations assume a common form
partly motivated by Regge theory, with flavor-specific PDF combinations given
by

$f_{i}(x,Q)=N_{i}x^{\alpha_{i}}\,(1-x)^{\beta_{i}}\ \mathcal{P}_{i}(x)\ ,$ (3)

and higher parametric flexibility attained by taking a polynomial of arbitrary
rank, $\mathcal{P}_{i}(x)\\!=\\!P^{j}_{i}(x)$, in the analytic case
Kotz:2023pbu , or a NN, $\mathcal{P}_{i}(x)\\!=\\!\mathrm{NN}_{i}(x)$
NNPDF:2021njg , in which the PDF shapes are encoded by the activation
functions and weights of the network. Alongside a chosen parametrization
basis, global analyses must also adopt a prescription for the comprehensive
determination of the associated PDF uncertainties. In this respect, the
commonly used paradigms are the Hessian formalism Hou:2019efy ; Bailey:2020ooq
; Alekhin:2017kpj  typically employed alongside analytical PDF
parametrizations, and the method of Monte Carlo (MC) resampling
Costantini:2024wby . In the former approach, a Hessian matrix is diagonalized
at the global minimum of the analysis likelihood function, $\chi^{2}$,
allowing a basis of eigendirections to be defined; with these, it is then
possible to saturate a tolerance criterion through finite displacements along
each eigendirection, such that the parametric dependence of the likelihood
function can be enfolded into the uncertainty via a finite collection of PDF
eigensets. The approach of MC resampling is commonly used in both NN
parametrizations of the PDFs NNPDF:2021njg  as well as fits based on
analytical forms Hunt-Smith:2022ugn ; in this scenario, statistical ensembles
of PDF error replicas are generated by optimizing model parameters to a set of
pseudodata sampled from Gaussian distributions about the fitted data in the
analysis. The resulting MC replicas then collectively encapsulate the PDF
uncertainty. In the present study, given the need to generate sufficiently
plentiful instances of PDF behavior to train the ResNet-like model introduced
in Sec. 1, we assume the validity of MC resampling for the purpose of
quantifying the true uncertainty, following the examples of Refs.
Kriesten:2023uoi ; Liu:2022plj , and generate sizable sets of replicas as
shown in Fig. 1.

This is possible due to the fact that systematic conversions between the
Hessian and MC approaches have been developed Hou:2016sho . For this work, PDF
MC replicas, labeled by index $(k)$ in the following equations, are created
from the Hessian error sets using the mcgen code Hou:2016sho  wherein the
central fitted PDF, $X_{0}$, is randomly displaced by a replica-specific
value, $d^{(k)}$, and global shift, $\Delta$:

$\displaystyle X^{(k)}=X_{0}+d^{(k)}-\Delta\ .$ (4)

The random displacements are given by positive and negative variations along
successive Hessian eigendirections, $X_{+i}$ and $X_{-i}$, about the central-
fitted PDF, $X_{0}$, according to a random Gaussian sampling,
$R_{i}^{(k)}\\!\sim\\!\mathcal{N}(0,1)_{i}^{(k)}$; these are finally summed
over each eigendirection, $i$, among the full set of $D$ eigendirections in
the phenomenological fit:

$\displaystyle d^{(k)}=\sum_{i=1}^{D}$ $\displaystyle\Bigg{(}{1\over
2}[X_{+i}-X_{-i}]\,R_{i}^{(k)}+{1\over
2}[X_{+i}+X_{-i}-2X_{0}]\,\left(R_{i}^{(k)}\right)^{2}\Bigg{)}\ .$ (5)

The global shift is the expectation value of the random displacements,

$\displaystyle\Delta$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{k=1}^{N}d^{(k)}.$ (6)

Figure 1: Monte Carlo replicas generated for the CT18 NNLO fit, plotted as
ratios to the central-fitted PDF; we take a perturbative scale $Q^{2}=100$
GeV2 and display the range $x=[10^{-4}-0.4]$ over which the ResNet model
inputs (training data) are sampled.

In the case of the MC replica method, the results of a PDF analysis are
encapsulated by the resulting statistical ensemble, with the mean over the
full set corresponding to the optimal fit222In practice, small deviations can
exist between the set of PDFs corresponding to the global minimum of
$\chi^{2}$ and the statistical mean over the MC ensemble; in the present work,
we explicitly ensure that any such deviations are zero.; the PDF uncertainties
can then be calculated as statistical deviations from this mean. In
particular, the average (central) prediction in the MC method for a PDF-
dependent quantity, $X$, is denoted as $\langle X\rangle$ and given by the
discrete expectation over replicas; namely,

$\displaystyle\langle X\rangle\approx X_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{k=1}^{N}X^{(k)}\ .$ (7)

While the symmetric PDF uncertainty for $X$ can then be obtained as the
standard deviation of replica-by-replica differences from this mean, the more
general case of the asymmetric uncertainty may be determined from the ensemble
of MC replicas at the $68\%$ confidence level (c.l.) as

$\displaystyle\delta^{MC>}X$ $\displaystyle=$
$\displaystyle+\sqrt{\langle(X-\langle X\rangle)^{2}\rangle_{X>\langle
X\rangle}}$ $\displaystyle\delta^{MC<}X$ $\displaystyle=$
$\displaystyle-\sqrt{\langle(X-\langle X\rangle)^{2}\rangle_{X<\langle
X\rangle}}\ ,$ (8)

where $>$ denotes the $68\%$ c.l. upward shift from the mean value and $<$
denotes the corresponding shift below. We note that Eq. (8) is applied by only
admitting those replicas predicting values of $X$ greater than $\langle
X\rangle$ in the calculation of the upward variation, and similarly in the
calculation of the downward asymmetric uncertainty. Taking $X=f_{i}(x,Q^{2})$,
this procedure may be applied on a point-by-point basis for a given value of
$x$ for each PDF flavor to obtain asymmetric PDF uncertainty bands about the
mean (central) PDF. We apply this procedure to generate $10^{4}$-member MC
replica sets for each of the CT18 fits explored in this work; we show an
example of central-set-normalized PDF ratio replicas for CT18 NNLO in Fig. 1.

PDF | Factorization scale | ATLAS 7 TeV $W/Z$ | CDHSW $F_{2}^{p,d}$ | Pole charm
---|---|---|---|---
fits | in DIS | data included? | data included? | mass, GeV
CT18 NNLO | $\mu_{F,DIS}^{2}=Q^{2}$ | No | Yes | 1.3
CT18A NNLO | $\mu_{F,DIS}^{2}=Q^{2}$ | Yes | Yes | 1.3
CT18X NNLO | $\mu_{F,DIS}^{2}=0.8^{2}\left(Q^{2}+\frac{0.3\mbox{ GeV}^{2}}{x_{B}^{0.3}}\right)$ | No | Yes | 1.3
CT18Z NNLO | $\mu_{F,DIS}^{2}=0.8^{2}\left(Q^{2}+\frac{0.3\mbox{ GeV}^{2}}{x_{B}^{0.3}}\right)$ | Yes | No | 1.4

Table 1:  The four PDF fits analyzed in this work and the differences in their
theory assumptions. Table obtained from Hou:2019efy .

To train our classifier network on a theoretically meaningful and self-
contained discrimination task, we take the set of global fits accompanying the
CT18 NNLO main release Hou:2019efy . This collection consisted of four QCD
analyses spanning several variations in the chosen theory and/or global data
set: CT18 as well as the three CT18A/X/Z alternative fits, where the different
analysis assumptions informing this collection are summarized in Table 1.
Taken together, the four NNLO fits of CT18 represent the most recent global
fit after the CT14 Dulat:2015mca  and CT14HERA2 Hou:2016nqm  releases,
following the inclusion of a variety of data from Run-1 of the LHC, such as
inclusive jet production, $W^{\pm}$/$Z$ boson production, and $t\bar{t}$ pair
production from the ATLAS, CMS, and LHCb. Among these, the CT18 NNLO main fit
constitutes a baseline fitting scenario to which the CT18A/X/Z alternative
fits were comparative variations allowing exploration of the relevant physics.
CT18A included an additional data set on 7 TeV $W^{\pm}$/$Z$ production from
ATLAS, discussed in detail in App. A of Ref. Hou:2019efy . Meanwhile, CT18X
implemented a modified scaling prescription for low-$x$ DIS, mimicking the
effects of low-$x$ resummation. Lastly, CT18Z combined the assumptions of the
CT18A and CT18X NNLO fits while also removing the CDHSW $F_{2}^{p,d}$ data and
shifting the charm mass to $m_{c}\\!=\\!1.4$ GeV to agree with HERA inclusive
DIS and LHC vector boson measurements Hou:2019efy . The CT18 and CT18Z fits
are the most distinct from each other in theory space due to their differing
underlying theory assumptions. We stress that the 4-member collection of fit
variations provided by CT18 provides a ideal testbed to demonstrate the XAI
methods developed in Sec. 5; this owes to its self-contained nature of this
simple ensemble based on a discrete set of analysis variations; these features
collectively make for a training task well-disposed to straightforward ML
labeling and classification.

Our data set for training is composed of MC replicas in the 6-fold flavor
basis of the CT18 fits, spanned by the $g$-, $s$-, $u$-, $d$-, $\bar{u}$-, and
$\bar{d}$-PDFs; given its proximity to the electroweak and Higgs scales, we
take each of these at a common physical scale of $Q^{2}\\!=\\!100$ GeV2, where
$Q=\mu_{F}=\mu_{R}$. In practice, we generate MC replicas from the available
LHAPDF grids Buckley_2015  associated with the original Hessian error sets
using default settings in the mcgen package Hou:2016sho . We use a split
$x$-scale that is logarithmic in the low-$x$ region, $(10^{-4}\\!-\\!10^{-1})$
and linear from $0.1\\!-\\!0.4$, with approximately $100$ points in each range
for a total of $200$ combined points covering the sampled range in $x$. We
note that the large- and small-$x$ PDFs are comparatively difficult to
constrain experimentally, leading to a sparsity of sensitive data relevant to
those regions. Consequently, the PDF shapes at high and low $x$ are
significantly determined by the assumed parametrization, with corresponding
increases in the PDF uncertainty. To mitigate the impact of the
parametrization on the classifications made by the NNs, we elect in this study
to train on self-normalized ratios, dividing replicas by their corresponding
central PDFs ($f$/$f_{0}$) to produce the information in Fig. 1.

As argued in Sec. 1, it would be beneficial to comprehensively connect PDF
analysis theory settings to the $x$-dependence of the associated PDFs in a
common and generalizable framework. We perform this task by classifying PDFs
across the CT18 series of phenomenological fits using the XAI approach
discussed in Sec. 3 below.

##  3 Explainable AI for PDF classification

Explainable AI (XAI) refers to a collection of techniques used to interpret a
ML model’s predictions in a human-readable format (see Refs.
gohel2021explainable ; xaioverview21  and references therein). Typically,
these techniques come in the form of saliency maps which highlight relevant
areas of the input feature space that were crucial during the prediction
process. The attribution scores calculated with XAI methods come in many forms
such as gradients, classification scores, or other numerical ratings. Given
our desire to trace the effects of theory assumptions on the $x$-dependent
behavior of individual instances of the fitted PDFs, in this study we
principally explore local methods; these explain the decisions of a NN with
respect to an individual sample of data. This is contrasted to global
attribution methods which explain decisions of a NN using features which were
significant across all inputs. Given the size and complexity of most NN
models, global methods are often not feasible.

Saliency maps are often used in computer vision to define regions where pixels
are substantially different from their surroundings. In classification tasks,
they highlight the regions of an image which most significantly impact the
classification score. We can translate this concept from computer vision to
QCD analysis to disentangle important physics that is learned by the NN during
classification. There are many different types of saliency methods that use
gradients simonyan2014deep ; sundararajan2017axiomatic ; shrikumar2017just ;
springenberg2015striving ; simonyan2014deep ; smilkov2017smoothgrad , class-
activated maps Selvaraju_2019 , occlusion zeiler2013visualizing , and even
edge detectors/image segmentation minaee2020image ; 4767851 . In this work, we
focus on a particular gradient-based saliency map called guided
backpropagation simonyan2014deep .

We start with a training data set denoted as
$\mathcal{D}=\\{\mathbf{x}^{(i)},y^{(i)}\\}$ which contains a set of input
features, $\mathbf{x}^{(i)}$ (here, $x$-dependent values of the PDFs), and
associated categorical labels, $y^{(i)}$ (corresponding parton flavors or PDF
models). Samples from the data set are drawn from the joint probability
distribution $p(\mathbf{x},y)$ where $\mathbf{x}$ and $y$ are random variables
of the input and label spaces, respectively. For a classification problem with
$C$ total classes, the classifier model can be functionally described by the
mapping
$f(\mathbf{x},\theta):\mathbb{R}^{\mathcal{D}}\rightarrow\mathbb{R}^{C}$
through model parameters, $\theta$. The classification model outputs
parameters of a categorical distribution represented as a vector of class
probabilities, $\hat{\pi}=\sigma(f(\mathbf{x},\theta))$ where the ML model
learns to approximate a categorical multinomial distribution:

$\displaystyle\text{Cat}(y|\mathbf{\hat{\pi}})$ $\displaystyle=$
$\displaystyle p(y=c|\mathbf{x}^{(i)},\theta)$ (9) $\displaystyle=$
$\displaystyle\prod_{c=1}^{C}\pi_{c}^{\mathbb{1}(y^{(i)}=c)}\ .$

Here, $\pi_{c}$ represent the model predictions for the probability that the
input belongs to the particular class $c$, and the function
$\mathbb{1}(\cdot)$ is $1$ when the condition is satisfied and $0$ otherwise.
The index $i$ runs through the full collection of training examples, of number
$N$, of feature/label pairs in the data set $\mathcal{D}$. The set of
predicted output scores is numerically unbounded,
$\\{f_{c}(\mathbf{x}^{(i)},\theta)=S_{c}^{(i)}\\}_{c=1}^{C}=S^{(i)}$;
therefore, in order to predict probabilities, the outputs of this function are
passed through a SoftMax activation, $\sigma$, defined as:

$\displaystyle\sigma\big{[}S^{(i)}\big{]}=\left[\frac{e^{S_{c=1}^{(i)}}}{\sum_{c=1}^{C}e^{S_{c}^{(i)}}},\dots,\frac{e^{S_{c=C}^{(i)}}}{\sum_{c=1}^{C}e^{S_{c}^{(i)}}}\right]\
.$ (10)

We train the model through maximum likelihood estimation which can be written
in terms of predicted labels,
$\displaystyle\hat{y}=\operatorname*{arg\,max}_{c}(\hat{\pi})$, as:

$\displaystyle\theta^{*}=\operatorname*{arg\,max}_{\theta}\left(\prod_{i=1}^{N}\prod_{c=1}^{C}p(\hat{y}=c|\mathbf{x}^{(i)},\theta)^{\mathbb{1}(y^{(i)}=c)}\right)\
.$ (11)

It is equivalent, and computationally more tractable, to minimize the negative
log likelihood, as doing so stabilizes the training regimen:

$\displaystyle\theta^{*}=\operatorname*{arg\,min}_{\theta}\Big{[}$
$\displaystyle\mathcal{L}(\theta)\Big{]}\ $
$\displaystyle\mathcal{L}(\theta)=-\frac{1}{N}\sum_{i=1}^{N}\sum_{c=1}^{C}{\mathbb{1}(y^{(i)}=c)}\,\ln{\left(p(\hat{y}=c|\mathbf{x}^{(i)},\theta)\right)}$
(12)

The latter expression above specifies the categorical cross-entropy loss, and
is most often used while training multi-class classification problems.

To further illustrate the mathematical structure of classifier models and
coincident use of backpropagation, we consider the functional form of a
shallow NN:

$\displaystyle
f^{\text{out}}=\sigma^{\ell}\Big{[}\omega^{\ell}\Big{(}\sigma^{\ell-1}\big{(}\omega^{\ell-1}\left(\dots\sigma^{0}\left(\omega^{0}x+b^{0}\right)+\dots\right)+b^{\ell-1}\big{)}+b^{\ell}\Big{)}\Big{]}\
,$ (13)

where $b^{\ell}$ and $\omega^{\ell}$ are the respective biases and weights of
layer $\ell$; $\sigma^{0}$, $\dots$, $\sigma^{\ell-1}$ represent ReLU
activation functions, $\text{ReLU}(x)=\text{max}(0,x)$, and $\sigma^{\ell}$ is
the SoftMax activation function defined in Eq. (10) above. During training, we
use autodifferentiation to efficiently compute gradients of the (categorical
cross-entropy) loss function, $\mathcal{L}$ of Eq. (12), with respect to the
NN’s parameters (i.e., its weights, $\omega$, and biases, $b$)
nielsen2015neural . Schematically, we can write these gradients in terms of a
specific layer $n$, where $n$ exists within a larger set of layers,
$\\{\ell,\ell-1,\dots,n+1,n,n-1,\dots,1,0\\}$:

$\displaystyle\frac{\partial\mathcal{L}}{\partial b^{n}}$ $\displaystyle=$
$\displaystyle[(\omega^{n+1})^{T}\delta^{n+1}]\odot\frac{\partial\sigma}{\partial
z}(z^{n})$ (14) $\displaystyle\frac{\partial\mathcal{L}}{\partial\omega^{n}}$
$\displaystyle=$ $\displaystyle
f^{n-1}[(\omega^{n+1})^{T}\delta^{n+1}]\odot\frac{\partial\sigma}{\partial
z}(z^{n})\ ,$ (15)

where $f^{n}=\sigma^{n}(z^{n})$ are the outputs of the $n$th layer,
$z^{n}=\omega^{n}f^{n-1}+b^{n}$, and the errors, $\delta$, are backpropagated
through the full NN from the final layer, $\ell$, to the initial layer:

$\displaystyle\delta^{\ell}$ $\displaystyle=$
$\displaystyle\nabla_{f}\mathcal{L}\odot\frac{\partial\sigma}{\partial
z}(z^{\ell})$ (16) $\displaystyle\delta^{n}$ $\displaystyle=$
$\displaystyle[(\omega^{n+1})^{T}\delta^{n+1}]\odot\frac{\partial\sigma}{\partial
z}(z^{n})\ .$ (17)

In the above notation, the symbol $\odot$ is the Hadamard product which is the
element-wise multiplication of two matrices. The backpropagation algorithm
begins with a forward pass from input feature space to the calculation of
$f^{out}$. At the final layer, $\ell$, the error vector $\delta^{\ell}$ of Eq.
(16) is computed and propagated back through to the previous layers using Eq.
(17). The gradients can then be subsequently calculated at each layer using
Eqs. (14) and (15). The optimizer uses the gradients to update the network
weights and biases. Autodifferentiation automates the application of the chain
rule through the functional mapping, tracking the mathematical operations and
their derivatives during the forward pass at each layer to facilitate
efficient gradient computation during the backward pass.

Once the model is trained, we can use autodifferentiation to compute the
gradients of the trained model’s pre-SoftMax classification score,
$S_{c}^{(i)}$, with respect to the input feature space in a single backward
pass Lerma_2023 . The guided backpropagation algorithm modifies how this
differentiation is performed over the rectified linear unit (ReLU) activations
by applying a double filter such that only positive gradients in the current
layer and positive inputs from the previous layer are propagated during a
single backward pass. This double masking is accomplished by replacing the
ReLU activation layers with custom layers which satisfy the following
condition:

$\frac{\partial f^{\text{out}}}{\partial
f_{i}^{n}}=\left(f_{i}^{n}>0\right)\cdot\left(\frac{\partial
f^{\text{out}}}{\partial f_{i}^{n+1}}>0\right)\cdot\frac{\partial
f^{\text{out}}}{\partial f_{i}^{n+1}}\ .$ (18)

Here, $f_{i}^{n}$ refer to the outputs of the $n^{\mathit{th}}$ layer and
$f^{\text{out}}$ is the output of the final pre-SoftMax layer. By double
masking, we constrain the autodifferentiation procedure to propagate only
those gradients that positively influence the classification score. This
approach ensures that we capture only the essential information necessary for
the model’s prediction, resulting in precision maps capable of fine detail.

##  4 ML model architecture and optimization

Figure 2: The architecture of the ResNet-like NN used to classify PDFs. The
model contains an input layer, four dense ResModules, followed by a final
fully-connected (FC) dense layer to predict classification scores,
$S_{c}^{(i)}$, and a SoftMax layer which renders the model outputs in
probabilistic terms.

We construct a ResNet-like he2015deep  NN, depicted in Fig. 2 to perform the
parton-level classification task(s) of our study. The model is constructed
from dense ResModules which contain a fully-connected (FC) dense layer,
followed by an activation layer and a dropout layer with skip connections that
alternate layers. We implement four such ResModules with a final output dense
layer, for a total of 9 dense layers. The output FC dense layer outputs the
classification scores, $S^{(i)}_{c}$, for each example and each class. The
final layer is a SoftMax layer which converts the output scores into a vector
of probabilities. We employ the ReLU activation function in the ResModules
which prevents negative inputs from propagating through the NN activations,
while dropout layers are used for regularization to ensure that the ML model
does not overfit.

As we are primarily interested in the regions significant to the
classification of all PDFs of a given class (per flavor and per
phenomenological fit), we must use a scoring method less specific than the
gradients corresponding to a single PDF MC replica. However, a global
explainability technique would be too broad, as it would not discriminate by
member of the error set. For each ratio, we therefore compute standard
deviations333Alternative prescriptions, for instance, based on computing
simple averages of the backpropagated classification-score gradients over all
members of the test set, produce very similar results. of the gradients of all
PDFs belonging to a specific class at each $x$ point. This allows us to
address ML-specific issues that arise from saliency-based interpretations of
ML predictions such as those discussed in Refs. kindermans2017unreliability ;
ghorbani2018interpretation ; adebayo2020sanity . The standard deviations tend
to be largest in those regions of the PDFs’ $x$ dependence corresponding to
the most significant shape changes, implying that these regions provide the
most sensitive discrimination basis for the PDF classifier. Additionally, we
smooth the gradients in $x$-space by evaluating a moving average, further
increasing the readability of the saliency maps without compromising their
accuracy. Finally, we obtain human-readable saliency maps by plotting the
values of our smoothed and normalized standard deviations as colors between
the uncertainty bounds of the ensemble of PDF MC replicas.

Along with the saliency maps, we report the confusion matrices for the ML
model’s predictions of the various PDF classifications. The entries of the
matrix are reported as percentages of correct predictions and are normalized
across the true labels, of which there are 1500 for each class.

##  5 XAI4PDF classification results

In this section, we summarize our main findings regarding the explainable
classification of PDFs. First, we report the saliency maps and confusion
matrices for the parton-flavor classification task within the CT18 NNLO
baseline fit. This first study allows us to identify which regions in
$x$-space the NN identifies as most important for accurately classifying the
parton flavor of the associated PDF. We then perform a similar analysis but
with the classification task of discriminating among the CT18(A/X/Z)
phenomenological fits. The XAI4PDF method identifies regions of importance
that correlate to specific theory choices, highlighting how the NN identifies
the manner in which specific choices in the underlying theory change the shape
of the downstream fitted PDFs.

Figure 3: Confusion matrix of the output classification scores for a PDF
flavor classifier tested on the CT18 NNLO PDF fit. Depicted are the true PDF
ratio labels versus predicted PDF ratio labels, normalized to the true labels,
and scores are given in percentages. The color bar indicates scores as well in
percentages for presentation clarity.

###  5.1 Parton flavor identification

We first demonstrate XAI techniques within a single phenomenological analysis,
CT18 NNLO, by identifying prominent $x$-dependent features informing the
accurate classification of the parton flavor associated with each PDF. In
particular, the XAI model learns which region(s) in $x$ makes a PDF of
specific flavor unique amongst the others. We train the ResNet-like
architecture described in Sec. 4 to distinguish among parton flavors by
correctly classifying them, specifically looking at $u(x)$, $d(x)$, $s(x)$,
$g(x)$, $\bar{u}(x)$, and $\bar{d}(x)$. The model was trained on the ratios of
PDF MC error sets to the central fit, with the aim of eliminating the
dependence on artifacts of the PDFs’ absolute normalizations and instead
emphasizing the uncertainty bands. We note that the XAI4PDF approach is of
sufficient generality that alternative training and classification strategies
based on some admixture of PDFs’ uncertainty bands and absolute [i.e.,
$\sim\\!\\!f_{i}(x,Q)$] behavior may also be well-motivated. We defer such
detailed studies to future work.

In Fig. 3, we show the results of our trained classification model on the PDF
ratio data for the CT18 NNLO fit, reported through a confusion matrix. The
model predictions are evaluated on a holdout test set which remains unseen
during the training process. The confusion matrix is normalized across rows
such that each row sums to $100\%$, with the columns depicting the percent
that the predicted label matched the true label. The confusion matrix shows
that the ratios for the gluon ($g$/$g_{0}$), strange ($s$/$s_{0}$), and up-
quark ($u$/$u_{0}$) PDFs are most accurately classified, with positive
identification rates of $99.5\%$, $95.8\%$, and $96\%$, respectively.
Conversely, the $d$-quark ratio is confused for that of the $u$-PDF by over
$\sim\\!10\%$, while the $\bar{d}$ ratio remains unseparated from the
$\bar{u}$ ratio by more than $\sim\\!17\%$. This confusion between the
$\bar{u}$ and $\bar{d}$ ratios is related to statistical indistinguishability
arising from flavor-separation challenges in phenomenological fits of the sea-
quark densities; this highlights the importance of further experimental
measurements of the $\bar{d}$/$\bar{u}$ asymmetry in experiments such as
SeaQuest Hou:2022sdf ; Guzzi:2021fre ; Hou:2022ajg .

Figure 4: Results of the guided backpropagation algorithm on the trained
ResNet-like architecture for PDF flavor classification. The gradients are
ensembled and we report a normalized variance. The gradients are superimposed
over the uncertainty bands of each PDF at $Q^{2}=100$ GeV2 over which the
model was tested. Plotted colors represent the normalized magnitudes of the
backpropagated gradient standard deviations as given by the scale at right.

The results of our guided backpropagation method on the trained ResNet-like
architecture for flavor classification are shown in Fig. 4. We report the
normalized variances for each ensemble of gradients per parton flavor as a
color superimposed over the PDF uncertainty band on which the model was
tested. We note that each flavor-dependent PDF band has a unique color
striping depicting the regions in $x$-space distinguishing that PDF as unique
relative to the others. The large variances in the gradients (depicted by red
bands) occur in regions where there are more significant shape changes
amalgamated over the set of MC replicas as compared to the other PDF flavors.
For example, at the $Q\\!=\\!100$ GeV scale of Fig. 4, the moderate-to-
high-$x$ region of $x\\!\in\\![10^{-2}\\!-\\!10^{-1}]$ is influential for
separating the $u$-, $d$-, $\bar{u}$-, and $\bar{d}$-PDFs as indicated by
gradient highlighting over that $x$-range. The characteristic region used by
the ML model to distinguish the gluon PDF is located at larger
$x\\!\sim\\!0.3$ (of interest in connection with jet and dijet production
Hou:2019efy ), while the strange PDF ratios are most discriminated by their
$x$ dependence near $x\\!=\\!10^{-1}$. By visual inspection of the uncertainty
band, the regions where the gradients vary the most are regions where there
are shape changes that are unique to that particular PDF. We deduce this to be
the likely reason the NN model occasionally confuses $\bar{d}$-PDF replicas
with those of the $\bar{u}$-PDF, but not reciprocally. The $\bar{d}$-PDF
contains replicas which encapsulate more of the features of the $\bar{u}$-PDF
replicas than the other way around. This pattern suggests that the
$\bar{d}$-PDF may be slightly less constrained than the corresponding
$\bar{u}$-PDF for the purposes of the PDF classifier.

###  5.2 Disentangling analysis elements among phenomenological fits

Many theory settings which enter phenomenological QCD fits can be actively
adjusted prior to optimizing the shape parameters for the PDFs. As argued
above, we aim to identify kinematical regions of the PDFs (in particular in
$x$-space) where these theory assumptions affect the fitted shape(s) of the
parton densities. We base our demonstration on the CT18(A/X/Z) analysis series
discussed in Sec. 2, which entailed a closed group of systematic variations in
the PDF analysis settings suitable for classifier labeling as summarized in
Table 1. Much as discussed for the parton-flavor classification task of Sec.
5.1, we now train our ResNet-like NN to identify PDF combinations according to
which of the four CT18 fits produced them.

Figure 5: The confusion matrix for the ResNet-like model trained to classify
members of the 4-fold collection of CT18(A/X/Z) fits. The scores are given as
percentages normalized across the rows according to the true label. The color
gradient indicates the score as well for clarity of presentation.

In Fig. 5, we again report the results of our trained ResNet-like model in the
form of a confusion matrix based on a holdout test set. The model generally
performs quite well, accurately classifying PDFs belonging to each of the
phenomenological fits, with the lowest performance at $83.1\%$. This most
common confusion occurs between the CT18 and CT18A replicas; as noted in Table
1, the latter of these differs from the CT18 baseline only in the inclusion of
the ATLAS 7 TeV data set. The observed behavior of the confusion matrix is
thus consistent with expectation, as the modified DIS factorization scale used
in both CT18X and CT18Z has a more apparent impact on the PDF uncertainties
over a range of flavors than the inclusion of the ATLAS 7 TeV data in the $x$
regions analyzed Hou:2019efy , especially given that the latter variation
primarily affected the strange PDF. Additionally, the CT18Z replicas are
confused with CT18X at a rate of $5.5\%$. This relatively mild rate of
confusion likely reflects the fact that, although CT18Z differs from CT18X in
every analysis setting toggled in Table 1 except the factorization scale, this
latter variation of the DIS scale is a dominant effect for the classifier and
somewhat overwhelms the variations related to the ATLAS and CDHSW sets and
$m_{c}$ shift. It is also interesting to confirm that the two analyses which
are “furthest” from each other based on the choices of Table 1, CT18 vs.
CT18Z, are also the least confused, confirming that the cumulative shift in
analysis assumptions drives the fits’ statistical distinguishability as
inferred by the XAI calculation.

Figure 6: Results of the XAI4PDF guided backpropagation calculation on the
trained ResNet-like architecture for classifying across phenomenological fits.
The gradients are ensembled and we report a normalized variance. As before,
the gradients overlay the uncertainty bands of each PDF at $Q^{2}=100$ GeV2 on
which the model was tested.

We show in Fig. 6 the results of the guided backpropagation saliency map
applied to the ResNet-like model trained on the PDF replicas from the CT18
series of fits. As in Fig. 4, colors indicate variances in the normalized
gradients, highlighting regions with prominent feature changes connected to
the ML model’s fit classifications. In general, the model identified the
gluon- and strange-PDF ratios as being most significant in distinguishing
among the phenomenological fits, as evidenced by the gradient patterns. In
contrast, the other PDF ratios which we have not explicitly plotted (i.e., the
scored ratios for $u$, $d$, $\bar{u}$, $\bar{d}$) showed very little gradient
structure, indicating that the model did not identify significant feature
changes among the various fits for these other PDF replicas. Consequently, we
report only the gluon and strange saliency maps as corresponding to the PDF
replicas most important in the guided backpropagation model classifications.
We find this observation to align closely with the findings in Ref.
Hou:2019efy , which explored the significant PDF impacts of the CT18 fit
variations with special attention on the gluon and strange parton densities.
The ability of the model to pick out this behavior and emphasize these
variations suggests that it is effectively learning the flavor-dependent
effects of the varying theory assumptions in the fits.

The results for the strange- and gluon-PDF ratios shown in Fig. 6 demonstrate
additional interesting behavior. For example, it is noteworthy that the gluon
replicas play a dominant role in the classifications among the CT18 series, as
reflected by the large peaked values of the gradient deviations in the
neighborhood of $x\\!\sim\\!10^{-3}$. The strange-PDF replicas have smoother,
less peaked gradients, indicating that their relevance in the classification
is important but perhaps not as strong for the ML model. The gradients tend to
peak in some cases at $x\\!\sim\\!10^{-2}$ but also at more moderate-$x$
values such as $x\\!\sim\\!0.2$. The $x$-regions which are highlighted by the
gradients are those where the uncertainty changes most between fits. These
results are supported by Ref. Hou:2019efy , where the $s$ and $g$ PDF
uncertainties are found to be affected more dramatically in those regions,
particularly at low-$x$.

Finally, it is natural to ask whether meaningful correlations exist which
might further reinforce the interpretability of the classification model,
e.g., between the NN classification and initial PDF; such quantities could
yield additional insights into the meaning of the saliency maps, while also
possessing calculable analogues in the input feature spaces on which the
classification is trained. Borrowing an analogy to Hessian analyses, we note
that the $L_{2}$ sensitivity method Hobbs:2019gob ; Jing:2023isu  can be
defined in the parametric vicinity of the global $\chi^{2}$ minimum. This
leads to the definition,

$S_{f}=\vec{\nabla}\chi^{2}_{E}\cdot{\vec{\nabla}f\over|\vec{\nabla}f|}=\Delta[\chi^{2}]\
\mathrm{Corr}_{\mathrm{H}}(\chi^{2},f)\ .$ (19)

On the logic that classification scores can be computed on the basis of MC
replicas generated from an original Hessian analysis and informed by notions
of statistical separability in that progenitor analysis, an approximate
analogue expression for the pre-SoftMax classification score, $S_{c}$, may
assume the form

$\Sigma_{c}\equiv\vec{\nabla}S_{c}\cdot{\vec{\nabla}f\over|\vec{\nabla}f|}\sim\Delta[S_{c}]\
\mathrm{Corr}_{\mathrm{MC}}(S_{c},f)\ ;$ (20)

this would be reasonable for scenarios in which classifications can be
directly projected onto the Hessian fits by which they were produced.
Ultimately, developing a more systematic matching between ML-based
classification methods like those embodied by the XAI4PDF framework could thus
proceed on the basis of such statistically-matched calculations, an
undertaking we expect to pursue in future work.

##  6 Conclusion

In this study, we have explored for the first time NN explainability
techniques for PDF phenomenology through a novel framework, XAI4PDF, applied
to the task of PDF classification and model discrimination. Using this
approach, we classified PDF replicas within a particular phenomenological fit
by flavor and identified key regions in the $x$-space that distinguished them.
Additionally, we classified PDF replicas across phenomenological fits with
incrementally shifted analysis internals, uncovering specific $x$-space
regions corresponding to these changes. We calculated confusion matrices for
each fit to ensure that our models achieved peak accuracy, and to understand
the specific errors made during classification and how those errors related to
aspects of the phenomenological fit.

We conclude with thoughts on the generalizability of the calculations explored
above as they relate to full-scale fitting frameworks. With sufficient
sampling, it may be possible to discriminate the downstream realization of
very subtle theoretical corrections or effects within a full phenomenological
fit, including those related to possible BSM scenarios. The use of the guided
backpropagation technique with its double masking procedure during
autodifferentiation allows for highly expressive and fine-grained gradient
structures, as one would expect from potential BSM contaminations.

In this work, we utilized a single explainability technique in the form of
guided backpropagation that fit the parameters of our specific problem;
however, many other explainability techniques exist that can be employed for
different but complementary problems in PDF theory. A further opportunity for
improvement pertains to developing additional quantitative checks of the
method 2023NatSR..1316887S  against which to correlate the saliency maps, for
instance, along the lines of the discussion at the end of Sec. 5.2. Extending
XAI4PDF with quantifiable measures of uncertainty would allow for a deeper and
more comprehensive analysis of the impact of variations in the PDF fitting
procedure.

This approach could also aid in the construction of PDF parametrizations by
visualizing the impact of different theoretical choices implemented during
fitting. Explainability techniques, as demonstrated in this paper, could serve
as valuable tools to analyze how and where varying the PDF parameters affects
the downstream uncertainty. A goal that is particularly attractive for the
advancement of high-precision SM tests and new physics discovery at future
machines where uncertainty quantification is crucial. By leveraging novel ML
explainability techniques, this work suggests options for more precise studies
of parton densities. Moreover, refinements of the approach used in this study
would be further applicable within analyses of the type leading to the
combined fits of the PDF4LHC variety, which are based on constructing MC
samples from a series of phenomenological analyses before statistically
enfolding these into a single amalgamated set. The XAI4PDF method would permit
an investigation of the statistical compatibility of input sets in the apples-
to-apples MC basis.

Code availability. We release a working version of XAI4PDF for example cases
of both PDF flavor and phenomenological fit classification at a dedicated
GitHub repository, https://github.com/ML4HEP-Theory/XAI4PDF.

##  7 Acknowledgments

The work of BK and TJH at Argonne National Laboratory was supported by the
U.S. Department of Energy under contract DE-AC02-06CH11357. The work of JG was
supported by the U.S. Department of Energy, Office of Science, Office of
Workforce Development for Teachers and Scientists (WDTS) under the Science
Undergraduate Laboratory Internship (SULI) program.

## References

  * (1) “High-Luminosity Large Hadron Collider (HL-LHC): Technical Design Report V. 0.1,” 2017.
  * (2) M. Cepeda et al., “Report from Working Group 2: Higgs Physics at the HL-LHC and HE-LHC,” CERN Yellow Rep. Monogr. 7 (2019) 221–584, arXiv:1902.00134 [hep-ph].
  * (3) A. Accardi et al., “Electron Ion Collider: The Next QCD Frontier: Understanding the glue that binds us all,” Eur. Phys. J. A 52 no. 9, (2016) 268, arXiv:1212.1701 [nucl-ex].
  * (4) ZEUS Collaboration, H. Abramowicz et al., “Limits on contact interactions and leptoquarks at HERA,” Phys. Rev. D 99 no. 9, (2019) 092006, arXiv:1902.03048 [hep-ex].
  * (5) S. Carrazza, C. Degrande, S. Iranipour, J. Rojo, and M. Ubiali, “Can New Physics hide inside the proton?,” Phys. Rev. Lett. 123 no. 13, (2019) 132001, arXiv:1905.05215 [hep-ph].
  * (6) A. Greljo, S. Iranipour, Z. Kassabov, M. Madigan, J. Moore, J. Rojo, M. Ubiali, and C. Voisey, “Parton distributions in the SMEFT from high-energy Drell-Yan tails,” JHEP 07 (2021) 122, arXiv:2104.02723 [hep-ph].
  * (7) M. Madigan and J. Moore, “Parton Distributions in the SMEFT from high-energy Drell-Yan tails,” PoS EPS-HEP2021 (2022) 424, arXiv:2110.13204 [hep-ph].
  * (8) J. Gao, M. Gao, T. J. Hobbs, D. Liu, and X. Shen, “Simultaneous CTEQ-TEA extraction of PDFs and SMEFT parameters from jet and $t\overline{t}$ data,” JHEP 05 (2023) 003, arXiv:2211.01094 [hep-ph].
  * (9) S. Iranipour and M. Ubiali, “A new generation of simultaneous fits to LHC data using deep learning,” JHEP 05 (2022) 032, arXiv:2201.07240 [hep-ph].
  * (10) S. Alekhin, J. Blümlein, S. Moch, and R. Placakyte, “Parton distribution functions, $\alpha_{s}$, and heavy-quark masses for LHC Run II,” Phys. Rev. D 96 no. 1, (2017) 014011, arXiv:1701.05838 [hep-ph].
  * (11) P. C. Barry, N. Sato, W. Melnitchouk, and C.-R. Ji, “First Monte Carlo Global QCD Analysis of Pion Parton Distributions,” Phys. Rev. Lett. 121 no. 15, (2018) 152001, arXiv:1804.01965 [hep-ph].
  * (12) T.-J. Hou et al., “New CTEQ global analysis of quantum chromodynamics with high-precision data from the LHC,” Phys. Rev. D 103 no. 1, (2021) 014013, arXiv:1912.10053 [hep-ph].
  * (13) S. Bailey, T. Cridge, L. A. Harland-Lang, A. D. Martin, and R. S. Thorne, “Parton distributions from LHC, HERA, Tevatron and fixed target data: MSHT20 PDFs,” Eur. Phys. J. C 81 no. 4, (2021) 341, arXiv:2012.04684 [hep-ph].
  * (14) NNPDF Collaboration, R. D. Ball et al., “The path to proton structure at 1% accuracy,” Eur. Phys. J. C 82 no. 5, (2022) 428, arXiv:2109.02653 [hep-ph].
  * (15) PDF4LHC Working Group Collaboration, R. D. Ball et al., “The PDF4LHC21 combination of global PDF fits for the LHC Run III,” J. Phys. G 49 no. 8, (2022) 080501, arXiv:2203.05506 [hep-ph].
  * (16) L. Kotz, A. Courtoy, P. Nadolsky, F. Olness, and M. Ponce-Chavez, “Analysis of parton distributions in a pion with Bézier parametrizations,” Phys. Rev. D 109 no. 7, (2024) 074027, arXiv:2311.08447 [hep-ph].
  * (17) K. Kovařík, P. M. Nadolsky, and D. E. Soper, “Hadronic structure in high-energy collisions,” Rev. Mod. Phys. 92 no. 4, (2020) 045003, arXiv:1905.06957 [hep-ph].
  * (18) X.-D. Ji, “Gauge-Invariant Decomposition of Nucleon Spin,” Phys. Rev. Lett. 78 (1997) 610–613, arXiv:hep-ph/9603249.
  * (19) A. V. Radyushkin, “Nonforward parton distributions,” Phys. Rev. D 56 (1997) 5524–5557, arXiv:hep-ph/9704207.
  * (20) X.-D. Ji, “Deeply virtual Compton scattering,” Phys. Rev. D 55 (1997) 7114–7125, arXiv:hep-ph/9609381.
  * (21) A. V. Belitsky, D. Mueller, and A. Kirchner, “Theory of deeply virtual Compton scattering on the nucleon,” Nucl. Phys. B 629 (2002) 323–392, arXiv:hep-ph/0112108.
  * (22) B. Kriesten, S. Liuti, L. Calero-Diaz, D. Keller, A. Meyer, G. R. Goldstein, and J. Osvaldo Gonzalez-Hernandez, “Extraction of generalized parton distribution observables from deeply virtual electron proton scattering experiments,” Phys. Rev. D 101 no. 5, (2020) 054021, arXiv:1903.05742 [hep-ph].
  * (23) B. Kriesten and S. Liuti, “Theory of deeply virtual Compton scattering off the unpolarized proton,” Phys. Rev. D 105 no. 1, (2022) 016015, arXiv:2004.08890 [hep-ph].
  * (24) Y. Guo, X. Ji, B. Kriesten, and K. Shiells, “Twist-three cross-sections in deeply virtual Compton scattering,” JHEP 06 (2022) 096, arXiv:2202.11114 [hep-ph].
  * (25) A. Accardi et al., “Strong Interaction Physics at the Luminosity Frontier with 22 GeV Electrons at Jefferson Lab,” arXiv:2306.09360 [nucl-ex].
  * (26) J.-W. Qiu and Z. Yu, “Single diffractive hard exclusive processes for the study of generalized parton distributions,” Phys. Rev. D 107 no. 1, (2023) 014007, arXiv:2210.07995 [hep-ph].
  * (27) J.-W. Qiu and Z. Yu, “Exclusive production of a pair of high transverse momentum photons in pion-nucleon collisions for extracting generalized parton distributions,” JHEP 08 (2022) 103, arXiv:2205.07846 [hep-ph].
  * (28) B. Kriesten, P. Velie, E. Yeats, F. Y. Lopez, and S. Liuti, “Parametrization of quark and gluon generalized parton distributions in a dynamical framework,” Phys. Rev. D 105 no. 5, (2022) 056022, arXiv:2101.01826 [hep-ph].
  * (29) V. Bertone, H. Dutrieux, C. Mezrag, H. Moutarde, and P. Sznajder, “Shadow generalized parton distributions: a practical approach to the deconvolution problem of DVCS,” SciPost Phys. Proc. 8 (2022) 107, arXiv:2107.11312 [hep-ph].
  * (30) Y. Guo, X. Ji, M. G. Santiago, K. Shiells, and J. Yang, “Generalized parton distributions through universal moment parameterization: non-zero skewness case,” JHEP 05 (2023) 150, arXiv:2302.07279 [hep-ph].
  * (31) E. Moffat, A. Freese, I. Cloët, T. Donohoe, L. Gamberg, W. Melnitchouk, A. Metz, A. Prokudin, and N. Sato, “Shedding light on shadow generalized parton distributions,” Phys. Rev. D 108 no. 3, (2023) 036027, arXiv:2303.12006 [hep-ph].
  * (32) S. Amoroso et al., “Snowmass 2021 Whitepaper: Proton Structure at the Precision Frontier,” Acta Phys. Polon. B 53 no. 12, (2022) 12–A1, arXiv:2203.13923 [hep-ph].
  * (33) A. Accardi, T. J. Hobbs, X. Jing, and P. M. Nadolsky, “Deuterium scattering experiments in CTEQ global QCD analyses: a comparative investigation,” Eur. Phys. J. C 81 no. 7, (2021) 603, arXiv:2102.01107 [hep-ph].
  * (34) J. Butterworth et al., “PDF4LHC recommendations for LHC Run II,” J. Phys. G 43 (2016) 023001, arXiv:1510.03865 [hep-ph].
  * (35) R. D. Ball et al., “Parton Distribution Benchmarking with LHC Data,” JHEP 04 (2013) 125, arXiv:1211.5142 [hep-ph].
  * (36) M. Botje et al., “The PDF4LHC Working Group Interim Recommendations,” arXiv:1101.0538 [hep-ph].
  * (37) T. J. Hobbs, B.-T. Wang, P. M. Nadolsky, and F. I. Olness, “Charting the coming synergy between lattice QCD and high-energy phenomenology,” Phys. Rev. D 100 no. 9, (2019) 094040, arXiv:1904.00022 [hep-ph].
  * (38) X. Jing et al., “Quantifying the interplay of experimental constraints in analyses of parton distributions,” Phys. Rev. D 108 no. 3, (2023) 034029, arXiv:2306.03918 [hep-ph].
  * (39) B.-T. Wang, T. J. Hobbs, S. Doyle, J. Gao, T.-J. Hou, P. M. Nadolsky, and F. I. Olness, “Mapping the sensitivity of hadronic experiments to nucleon structure,” Phys. Rev. D 98 no. 9, (2018) 094030, arXiv:1803.02777 [hep-ph].
  * (40) M. Almaeen, Y. Alanazi, N. Sato, W. Melnitchouk, M. P. Kuchera, and Y. Li, “Variational autoencoder inverse mapper: An end-to-end deep learning framework for inverse problems,” in 2021 International Joint Conference on Neural Networks (IJCNN), pp. 1–8.  2021\.
  * (41) B. Kriesten and T. J. Hobbs, “Learning PDFs through Interpretable Latent Representations in Mellin Space,” arXiv:2312.02278 [hep-ph].
  * (42) M. Yan, T.-J. Hou, Z. Li, K. Mohan, and C. P. Yuan, “A generalized statistical model for fits to parton distributions,” arXiv:2406.01664 [hep-ph].
  * (43) S. Carrazza, J. M. Cruz-Martinez, and T. R. Rabemananjara, “Compressing PDF sets using generative adversarial networks,” Eur. Phys. J. C 81 no. 6, (2021) 530, arXiv:2104.04535 [hep-ph].
  * (44) D. Liu, C. Sun, and J. Gao, “Machine learning of log-likelihood functions in global analysis of parton distributions,” JHEP 08 (2022) 088, arXiv:2201.06586 [hep-ph].
  * (45) K. Kumericki, S. Liuti, and H. Moutarde, “GPD phenomenology and DVCS fitting: Entering the high-precision era,” Eur. Phys. J. A 52 no. 6, (2016) 157, arXiv:1602.02763 [hep-ph].
  * (46) M. Čuić, K. Kumerički, and A. Schäfer, “Separation of Quark Flavors Using Deeply Virtual Compton Scattering Data,” Phys. Rev. Lett. 125 no. 23, (2020) 232005, arXiv:2007.00029 [hep-ph].
  * (47) M. Almaeen, J. Grigsby, J. Hoskins, B. Kriesten, Y. Li, H.-W. Lin, and S. Liuti, “Benchmarks for a Global Extraction of Information from Deeply Virtual Exclusive Scattering,” arXiv:2207.10766 [hep-ph].
  * (48) M. Almaeen, T. Alghamdi, B. Kriesten, D. Adams, Y. Li, and H.-W. L. a. S. Liuti, “VAIM-CFF: A variational autoencoder inverse mapper solution to Compton form factor extraction from deeply virtual exclusive reactions,” arXiv:2405.05826 [hep-ph].
  * (49) B. Kriesten and T. J. Hobbs, 2024.  In preparation.
  * (50) ATLAS Collaboration, M. Aaboud et al., “Precision measurement and interpretation of inclusive $W^{+}$ , $W^{-}$ and $Z/\gamma^{*}$ production cross sections with the ATLAS detector,” Eur. Phys. J. C 77 no. 6, (2017) 367, arXiv:1612.03016 [hep-ex].
  * (51) A. Ablat, M. Guzzi, K. Xie, S. Dulat, T.-J. Hou, I. Sitiwaldi, and C. P. Yuan, “Exploring the impact of high-precision top-quark pair production data on the structure of the proton at the LHC,” Phys. Rev. D 109 no. 5, (2024) 054027, arXiv:2307.11153 [hep-ph].
  * (52) J. C. Collins, D. E. Soper, and G. F. Sterman, “Factorization for Short Distance Hadron - Hadron Scattering,” Nucl. Phys. B 261 (1985) 104–142.
  * (53) J. Collins, Foundations of perturbative QCD, vol. 32.  Cambridge University Press, 11, 2013.
  * (54) CTEQ Collaboration, R. Brock et al., “Handbook of perturbative QCD: Version 1.0” Rev. Mod. Phys. 67 (1995) 157–248.
  * (55) V. N. Gribov and L. N. Lipatov, “Deep inelastic e p scattering in perturbation theory,” Sov. J. Nucl. Phys. 15 (1972) 438–450.
  * (56) Y. L. Dokshitzer, “Calculation of the Structure Functions for Deep Inelastic Scattering and e+ e- Annihilation by Perturbation Theory in Quantum Chromodynamics.,” Sov. Phys. JETP 46 (1977) 641–653.
  * (57) G. Altarelli and G. Parisi, “Asymptotic Freedom in Parton Language,” Nucl. Phys. B 126 (1977) 298–318.
  * (58) J. McGowan, T. Cridge, L. A. Harland-Lang, and R. S. Thorne, “Approximate N3LO parton distribution functions with theoretical uncertainties: MSHT20aN3LO PDFs,” Eur. Phys. J. C 83 no. 3, (2023) 185, arXiv:2207.04739 [hep-ph]. [Erratum: Eur.Phys.J.C 83, 302 (2023)].
  * (59) A. Cooper-Sarkar, T. Cridge, F. Giuli, L. A. Harland-Lang, F. Hekhorn, J. Huston, G. Magni, S. Moch, and R. S. Thorne, “A Benchmarking of QCD Evolution at Approximate $N^{3}LO$,” arXiv:2406.16188 [hep-ph].
  * (60) NNPDF Collaboration, R. D. Ball et al., “The Path to N3LO Parton Distributions,” arXiv:2402.18635 [hep-ph].
  * (61) M. N. Costantini, M. Madigan, L. Mantani, and J. M. Moore, “A critical study of the Monte Carlo replica method,” arXiv:2404.10056 [hep-ph].
  * (62) N. T. Hunt-Smith, A. Accardi, W. Melnitchouk, N. Sato, A. W. Thomas, and M. J. White, “Determination of uncertainties in parton densities,” Phys. Rev. D 106 no. 3, (2022) 036003, arXiv:2206.10782 [hep-ph].
  * (63) T.-J. Hou et al., “Reconstruction of Monte Carlo replicas from Hessian parton distributions,” JHEP 03 (2017) 099, arXiv:1607.06066 [hep-ph].
  * (64) S. Dulat, T.-J. Hou, J. Gao, M. Guzzi, J. Huston, P. Nadolsky, J. Pumplin, C. Schmidt, D. Stump, and C. P. Yuan, “New parton distribution functions from a global analysis of quantum chromodynamics,” Phys. Rev. D 93 no. 3, (2016) 033006, arXiv:1506.07443 [hep-ph].
  * (65) T.-J. Hou, S. Dulat, J. Gao, M. Guzzi, J. Huston, P. Nadolsky, J. Pumplin, C. Schmidt, D. Stump, and C. P. Yuan, “CTEQ-TEA parton distribution functions and HERA Run I and II combined data,” Phys. Rev. D 95 no. 3, (2017) 034003, arXiv:1609.07968 [hep-ph].
  * (66) A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page, M. Rüfenacht, M. Schönherr, and G. Watt, “Lhapdf6: parton density access in the lhc precision era,” The European Physical Journal C 75 no. 3, (Mar., 2015) . http://dx.doi.org/10.1140/epjc/s10052-015-3318-8.
  * (67) P. Gohel, P. Singh, and M. Mohanty, “Explainable ai: current status and future directions,” 2021.
  * (68) P. P. Angelov, E. A. Soares, R. Jiang, N. I. Arnold, and P. M. Atkinson, “Explainable artificial intelligence: an analytical review,” WIREs Data Mining and Knowledge Discovery 11 no. 5, (2021) e1424, https://wires.onlinelibrary.wiley.com/doi/pdf/10.1002/widm.1424. https://wires.onlinelibrary.wiley.com/doi/abs/10.1002/widm.1424.
  * (69) K. Simonyan, A. Vedaldi, and A. Zisserman, “Deep inside convolutional networks: Visualising image classification models and saliency maps,” 2014.
  * (70) M. Sundararajan, A. Taly, and Q. Yan, “Axiomatic attribution for deep networks,” 2017.
  * (71) A. Shrikumar, P. Greenside, A. Shcherbina, and A. Kundaje, “Not just a black box: Learning important features through propagating activation differences,” 2017.
  * (72) J. T. Springenberg, A. Dosovitskiy, T. Brox, and M. Riedmiller, “Striving for simplicity: The all convolutional net,” 2015.
  * (73) D. Smilkov, N. Thorat, B. Kim, F. Viégas, and M. Wattenberg, “Smoothgrad: removing noise by adding noise,” 2017.
  * (74) R. R. Selvaraju, M. Cogswell, A. Das, R. Vedantam, D. Parikh, and D. Batra, “Grad-cam: Visual explanations from deep networks via gradient-based localization,” International Journal of Computer Vision 128 no. 2, (Oct., 2019) 336–359. http://dx.doi.org/10.1007/s11263-019-01228-7.
  * (75) M. D. Zeiler and R. Fergus, “Visualizing and understanding convolutional networks,” 2013.
  * (76) S. Minaee, Y. Boykov, F. Porikli, A. Plaza, N. Kehtarnavaz, and D. Terzopoulos, “Image segmentation using deep learning: A survey,” 2020.
  * (77) J. Canny, “A computational approach to edge detection,” IEEE Transactions on Pattern Analysis and Machine Intelligence PAMI-8 no. 6, (1986) 679–698.
  * (78) M. A. Nielsen, Neural networks and deep learning, vol. 25.  Determination press San Francisco, CA, USA, 2015.
  * (79) M. Lerma and M. Lucas, “Pre or post-softmax scores in gradient-based attribution methods, what is best?” in 2023 IEEE 13th International Conference on Pattern Recognition Systems (ICPRS).  IEEE, July, 2023.  http://dx.doi.org/10.1109/ICPRS58416.2023.10179032.
  * (80) K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” 2015.
  * (81) P.-J. Kindermans, S. Hooker, J. Adebayo, M. Alber, K. T. Schütt, S. Dähne, D. Erhan, and B. Kim, “The (un)reliability of saliency methods,” 2017.
  * (82) A. Ghorbani, A. Abid, and J. Zou, “Interpretation of neural networks is fragile,” 2018.
  * (83) J. Adebayo, J. Gilmer, M. Muelly, I. Goodfellow, M. Hardt, and B. Kim, “Sanity checks for saliency maps,” 2020.
  * (84) T.-J. Hou, H.-W. Lin, M. Yan, and C. P. Yuan, “Impact of lattice $s(x)-\bar{s}(x)$ data in the CTEQ-TEA global analysis,” in Snowmass 2021.  4, 2022.  arXiv:2204.07944 [hep-ph].
  * (85) M. Guzzi et al., “NNLO constraints on proton PDFs from the SeaQuest and STAR experiments and other developments in the CTEQ-TEA global analysis,” SciPost Phys. Proc. 8 (2022) 005, arXiv:2108.06596 [hep-ph].
  * (86) T.-J. Hou, M. Yan, J. Liang, K.-F. Liu, and C. P. Yuan, “Connected and disconnected sea partons from the CT18 parametrization of PDFs,” Phys. Rev. D 106 no. 9, (2022) 096008, arXiv:2206.02431 [hep-ph].
  * (87) K. Szczepankiewicz, A. Popowicz, K. Charkiewicz, K. Nałecz-Charkiewicz, M. Szczepankiewicz, S. Lasota, P. Zawistowski, and K. Radlak, “Ground truth based comparison of saliency maps algorithms,” Scientific Reports 13 (Oct., 2023) 16887.