# Overview of intermittency analysis in heavy-ion collisions

Zhiming Li$\star$ Key Laboratory of Quark and Lepton Physics (MOE) and
Institute of Particle Physics,
Central China Normal University, Wuhan 430079, China
*<EMAIL_ADDRESS>
###### Abstract

In this paper, a search for power-law fluctuations with fractality and
intermittency analysis to explore the QCD phase diagram and the critical point
is summarized. Experimental data on self-similar correlations and fluctuations
with respect to the size of phase space volume in various high energy heavy-
ion collisions are presented, with special emphasis on background subtraction
and efficiency correction of the measurement. Phenomenological modelling and
theoretical work on the subject are discussed. Finally, we highlight possible
directions for future research.

###### keywords:

Power-law fluctuations; intermittency; QCD phase diagram; critical point;
heavy-ion collisions.

Received (Day Month Year)Revised (Day Month Year)

##  1 Introduction

One of the major goals in current heavy-ion collisions is to locate the
critical point (CP) in the phase diagram of strongly interacting matter
predicted by the Quantum Chromodynamics (QCD) [1, 2, 3, 4]. Lattice QCD has
shown that there occurs a smooth crossover from hadronic phase to the Quark
Gluon Plasma (QGP) phase at $\mu_{B}=0$ [5, 6, 7]. QCD based model
calculations indicate that the transition could be a first-order at large
$\mu_{B}$ [8, 9]. The endpoint of the first-order phase transition to the
crossover is referred to as the critical point [10, 11, 12, 13]. However,
Lattice QCD calculations at finite $\mu_{B}$ face numerical challenges in
computing. The location of the CP is highly theoretically uncertain.
Investigations from heavy-ion experiments, such as the Beam Energy Scan (BES)
II program at the Relativistic Heavy Ion Collider (RHIC) [14, 15], and
phenomenological modeling are required.

Fluctuations of conserved quantities, which behave differently between the
hadronic and the QGP phases, are normally considered to be promising
signatures for the QCD phase transition [16, 17, 18, 19, 20, 21, 22, 23, 24,
25, 26, 3, 4, 27]. The singularity at the CP, at which the transition is
believed to be the second-order, may cause enhancement of fluctuations if
fireballs created in the collision pass the CP region during the time
evolution [28, 29, 30, 31]. In heavy-ion collisions, in analogy to the
critical opalescence phenomenon observed in Quantum Electrodynamics [32],
large density fluctuation of late-stage baryon numbers is expected to develop
for the system near the CP due to a rapid increase of correlation length [1,
33, 34]. Strong baryon density fluctuations should be observed in experiments
at kinetic freeze-out if such fluctuations can survive during final-state
interactions and hadronic evolution of the system.

It is suggested that the expected large baryon density fluctuations near the
CP can be expressed in terms of density-density correlator of baryon numbers,
which exhibiting a fractal and scale invariant behavior [34, 35, 36]. These
fluctuations are supposed to form a unique pattern of self-similar or
intermittency behavior of finite state particles in high energy collisions
[34, 37, 38, 39]. Intermittency refers to the property of scale invariance,
fractality and a stochastic nature of the underlying scaling law [36]. It can
be measured by calculation of a scaled factorial moment (SFM) at various
system scales in transverse momentum space [34, 40]. The advantage of using
SFM is that it can characterize the non-statistical fluctuations in momentum
spectra, which is supposed to be connected with the dynamics of particle
production.

The last decade has witnessed remarkably exciting experimental explorations of
intermittency and fractality in high energy collisions. A systematic search
for QCD critical fluctuations has been performed by NA49 and NA61/SHINE
collaborations at CERN SPS with measurements of intermittency in A+A
collisions [41, 40, 42]. At RHIC energies, the studies on both charged
particles [43] and strangeness [44] in Au+Au collisions at BES-I energies are
reported. With the LHC at CERN ushering in a new era of TeV-scale high energy
physics, the CMS collaboration has measured SFM and intermittency in $pp$
collisions up to 8 TeV [45]. It is expected that fractality or self-
similarity, analogous to that encountered in complex non-linear systems, might
open a new way leading towards deeper insight into the CP and phase transition
of QCD.

At the same time, different phenomenological models have been used to study
the unique behavior of SFMs and intermittency under various underlying
mechanisms [46, 47, 48, 49, 50, 51]. With a Critical Monte-Carlo (CMC) model
belonging to the 3D Ising universality class, it is found that the self-
similar or intermittency nature of particle correlations is closely related to
the large baryon density fluctuation associated with the critical point [46].
By including hadronic mean-field potentials [47] or hydrodynamical
descriptions [48] into the data sample of a transport UrQMD model, the result
exhibits a clear self-similar behavior. It infers that intermittency is
associated with the mechanism of evolution of the medium produced in heavy-ion
collisions.

Intermittency is measured by calculations of SFM in momentum space at various
scales. Whereas, the single-particle multiplicity spectra of late-stage
particles in heavy-ion collisions are largely influenced by background effects
[52, 53, 54, 55]. In experiments, insufficient detector efficiency [56, 57]
will also modify the value of SFM. To understand the underlying physics of
this measurement, a careful study on the non-critical effects to get a clean
signal is needed. After that, we can compare results among different
experiments and with those from theoretical or phenomenological predictions.

This paper contains a review of the present study on intermittency in heavy-
ion collisions over the last years. In Section 2, we introduce the method of
intermittency analysis and the necessary formalism. Section 3 discusses the
estimation and subtraction of background contributions in the calculation of
SFMs. In section 4, the efficiency correction formula on SFM is derived.
Section 5 gives an overview of experimental measurements on intermittency in
current heavy-ion collisions. Section 6 is devoted to phenomenological model
results on the subject and the search for power-laws. Conclusions and outlooks
are summarized in Section 7.

##  2 Method of Analysis

In high-energy collisions, a power-law density fluctuation is proposed to be
detectable in momentum space through intermittency analysis. The observables
of our interests in this analysis are chosen to be sensitive to the power-law
singularity of the density-density correlation function. The SFM of
multiplicity distributions of final state particles is suggested as one of the
most suitable quantities [34, 40]. For this purpose, a selected D-dimensional
momentum space is partitioned into equal-sized cells. The $q$th-order SFM is
defined as:

$F_{q}(M)=\frac{\langle\frac{1}{M^{D}}\sum_{i=1}^{M^{D}}n_{i}(n_{i}-1)\cdots(n_{i}-q+1)\rangle}{\langle\frac{1}{M^{D}}\sum_{i=1}^{M^{D}}n_{i}\rangle^{q}},$
(1)

where $M$ is the number of cells in one dimension, $n_{i}$ is the measured
number of particles in the $i$th cell, and the angular bracket denotes an
average over the event sample.

For the system near the critical point, it is expected to observe a scaling or
power-law behavior of the SFM on the partitioned number of cells $M^{D}$,
given $M$ is large enough:

$F_{q}(M)\sim(M^{D})^{\phi_{q}},M\rightarrow\infty.$ (2)

Here $\phi_{q}$ is called intermittency index which specifies the strength of
the self-similar property. The method of using SFM to search for the critical
fluctuation was firstly proposed in Refs [52, 58] several years ago. By using
the effective action of a 3D Ising system belonging to the same universality
class of the QCD critical point [59, 60, 61], the critical $\phi_{2}$ is
predicted to be $\frac{5}{6}$ for baryon density [34] and $\frac{2}{3}$ for
pion density [62].

Besides the power-law behavior of SFMs on the partitioned numbers as described
above, another promising relation between SFMs of different orders is
suggested [63, 64, 65]:

$F_{q}(M)\propto F_{2}(M)^{\beta q},$ (3)

where $\beta_{q}=\phi_{q}/\phi_{2}$. According to Ginzburg-Landau (GL) theory,
the power-law behavior of $F_{q}(M)\sim M^{D}$ might not be observed near the
QCD critical point since $\phi_{q}$ depends on particular critical parameters
which would vary with temperatures of the system and are unknown in nuclear
collisions [63, 64]. However, $\beta_{q}$ is independent of these parameters
and thus the power-law behavior of $F_{q}(M)$ on $F_{2}(M)$ is feasible to
measure in experiments of nuclear-nuclear collisions.

The scaling exponent $\nu$ specifies the power-law behavior of $F_{q}(M)$ on
$F_{2}(M)$ and quantitatively describes the values of $\beta_{q}$ [63, 64,
65]:

$\beta_{q}\propto(q-1)^{\nu}.$ (4)

This critical exponent $\nu$, which is supposed to be independent on the
chosen values of specific critical parameters, can be used to investigate the
presence of intermittency in the QCD phase transition [63, 64, 65, 66].
Theoretical prediction for the critical value of $\nu$ is equal to 1.304 in
the full phase space based on Ginzburg-Landau (GL) theory [63], and 1.0 from
calculations of a two-dimensional Ising model [64].

##  3 Background Subtraction

In the calculations of SFMs, it is important to estimate and subtract trivial
contributions from background. It has been found that the values of the
measured SFMs can be significantly modified by adding uncorrelated particles
as a background to the event samples containing self-similar signals. [54].
These background effects should be eliminated in the intermittency analysis.

###  3.1 Mixed Event Method

By analyzing the NA49 and NA61 results, it is shown that the experimentally
measured scaling behavior of SFMs and the intermittency index can be
reproduced by adding more than $90\%$ uncorrelated random tracks into the
event sample generated by the CMC model [40, 67]. In this purpose, they
propose to use the mixed event method to estimate background contributions
[40, 68]. Mixed events are constructed by randomly selecting particles from
original events while reproducing the same multiplicity distributions. The
correlations between pairs of particles which exist in the original event, are
eliminated in the mixed event samples since each particle now is chosen from
different events. By assuming that multiplicity distributions of measured
particles in each partitioned cell can be simply divided into background and
critical contributions, the correlator $\Delta F_{q}(M)$, supposed to contain
only critical contributions, is defined as:

$\Delta F_{q}(M)=F_{q}^{data}(M)-F_{q}^{mix}(M).$ (5)

After subtracting background contributions, the intermittency index $\phi_{q}$
can be obtained from $\Delta F_{q}(M)$ instead of $F_{q}(M)$ by using Eq. (2).

###  3.2 Cumulative Variable Method

The cumulative variable method was proposed for intermittency analysis in Refs
[69, 70] several years ago. It has been proved to be able to effectively
reduce the distortions of a non-uniform single-particle spectrum [69, 71].

Assuming the variable $x$, e.g., momentum $p_{x}$ or $p_{y}$, is measured by a
density distribution function $\rho(x)$. The cumulative variable $X(x)$ is
defined as the transformation of [69, 70]

$\large
X(x)=\frac{\int_{x_{min}}^{x}\rho(x)dx}{\int_{x_{min}}^{x_{max}}\rho(x)dx},$
(6)

where $x_{min}$ and $x_{max}$ are minimum and maximum values of the selected
variable $x$.

The advantage of the cumulative variable $X(x)$ is that the value of the new
variable does not depend on particular variable $x$ but on its density
distribution function $\rho(x)$, which provides a way to compare results among
different experiments. The other benefit is that the probability distribution
of $X(x)$ is uniform which can remove the dependence of intermittency on
density distribution of $\rho(x)$ [69].

In the calculation of $F_{q}(M)$ of cumulative variable in a 2D momentum
region, the transverse momentum space of $p_{x}p_{y}$ is transfer to the one
of $p_{X}p_{Y}$. To be clear, we denote the $F_{q}(M)$ after cumulative
transformation as $CF_{q}(M)$. The process of fitting intermittency index
$\phi_{q}^{c}$ from $CF_{q}(M)$ is the same as $\phi_{2}$ from $F_{2}(M)$ by
Eq. (2).

Figure 1: (a) The second-order $F_{2}(M)$ (black stars) and $CF_{2}(M)$ (blue
circles) as a function of number of partitioned bins in the CMC model. (b)
$F_{2}(M)$ and $CF_{2}(M)$ for the same CMC event samples but being
contaminated by uncorrelated particles with a statistical Gaussian
distribution.

To verify the validity of the this method, the CMC model [34, 46] belonging to
the 3D Ising universality is applied to generate event samples incorporating
critically correlated particles in momentum space. The CMC model involves
large density fluctuations risen from the self-similar correlations among
particles and leads to the intermittency index of $\phi_{2}=\frac{5}{6}$ [34,
46]. Figure 1 (a) shows the $F_{2}(M)$ and $CF_{2}(M)$ as a function of
$M^{2}$ from the CMC model in a double-logarithmic scale. It is clearly
observed that both $F_{2}(M)$ (black stars) and $CF_{2}(M)$ (blue circles)
obey a good power-law dependence on the number of division bins. The fitted
$\phi_{2}$ is equal to $\phi_{2}^{c}$ within statistical uncertainty. It
implies that the intermittency behavior keeps unchanged by applying the
cumulative variable method to a critical event sample with pure self-similar
fluctuations.

We further test the method by adding background fluctuations into the signal.
The CMC samples are contaminated by particles with a Gaussian distribution as
background by a mixed ratio of $\lambda=95\%$. In Fig. 1 (b), one observes
that the directly calculated $F_{2}(M)$ deviates substantially from the power-
law dependence on $M^{2}$ and can not be fitted by the relation of Eq. (2).
However, $CF_{2}(M)$ still follows the same scaling behavior with increasing
$M^{2}$ as that in Fig. 1 (a). The fitted $\phi_{2}^{c}$ is found to be equal
to the one in the pure CMC event sample. It confirms that the cumulative
variable method is able to effectively erase background effects in the
intermittency analysis.

##  4 Efficiency Correction

Besides background subtraction, efficiency correction for SFMs is also an
important aspect in the calculations of intermittency. In experimental
measurement, some particles are missing due to a limited capacity of the
detector, which leads to the measured multiplicity distribution to be
different from the originally produced one in the momentum space. Therefore,
the calculated SFM will be changed accordingly. To recover the true SFM from
the experimentally measured one, a cell-by-cell efficiency correction method
is suggested in Ref. [71].

The probability function of $p(n)$ is related to the one of $p(N)$ with [72,
73]:

$p(n)=\sum_{N}W(n|N)p(N).$ (7)

Here $W(n|N)$ is the probability of experimentally measured $n$ particles,
given $N$ originally produced particles in the event. In general, $W(n|N)$ can
be approximated by a binomial distribution as [72, 73, 74]:

$W(n|N)=B(n,N;\epsilon)=\frac{N!}{n!(N-n)!}\epsilon^{n}(1-\epsilon)^{N-n},$
(8)

where $\epsilon$ is the particle detection efficiency.

Based on Eq. (7) and (8), the true factorial moment $f_{q}^{true}=\langle
N(N-1)...(N-q+1)\rangle$ will be restored by dividing the measured moment
$f_{q}^{measured}=\langle n(n-1)...(n-q+1)\rangle$, with certain power of the
experimental detector efficiency $\epsilon$ [74, 72, 75]:

$f_{q}^{corrected}=\frac{f_{q}^{measured}}{\epsilon^{q}}=\frac{\langle
n(n-1)...(n-q+1)\rangle}{\epsilon^{q}}.$ (9)

This method has been extensively used for efficiency corrections on the high-
moment studies in experiments [74, 26, 75, 25].

We use the same strategy in the intermittency analysis. The $q$-th moment
$\langle n_{i}(n_{i}-1)...(n_{i}-q+1)\rangle$ of SFM in each bin is corrected
according to Eq. (9). Then, the SFM defined in Eq. (1) can be corrected as:

$\displaystyle
F_{q}^{corrected}(M)=\frac{\langle\frac{1}{M^{2}}\sum_{i=1}^{M^{2}}\frac{n_{i}(n_{i}-1)\cdots(n_{i}-q+1)}{\bar{\epsilon_{i}}^{q}}\rangle}{\langle\frac{1}{M^{2}}\sum_{i=1}^{M^{2}}\frac{n_{i}}{\bar{\epsilon_{i}}}\rangle^{q}}.$
(10)

Where $n_{i}$ is the number of particles located in the $i$-th cell.
$\bar{\epsilon_{i}}$ represents the event average of the mean efficiency of
all the particles in the $i$-th cell, which is calculated by
$\langle\frac{\sum_{j=1}^{n_{i}}\epsilon_{i}^{j}}{n_{i}}\rangle$. The
efficiency correction of Eq. (10) is known as the cell-by-cell method [71].

Figure 2: (a) The STAR experimental detection efficiency of protons as a
function of $p_{T}$ in the TPC detector in Au+Au collisions at
$\sqrt{s_{\mathrm{NN}}}$ = 19.6 GeV. (b) The original true $F_{2}(M)$ (open
squares) as a function of number of bins in the transverse momentum region,
together with the measured $F_{2}(M)$ (solid stars) after discarding particles
and efficiency corrected SFM (solid crosses) by using the proposed cell-by-
cell method. Figure 3: (a) The STAR experimental tracking efficiency as a
function of $p_{T}$ in the TPC+TOF detectors in Au+Au collisions at
$\sqrt{s_{\mathrm{NN}}}$ = 19.6 GeV. (b) The original true $F_{2}(M)$ (open
squares) as a function of $M^{2}$, together with the measured $F_{2}(M)$
(solid stars) after discarding particles and the efficiency corrected
$F_{2}(M)$ (solid crosses).

In order to test the accuracy of the cell-by-cell method, we apply the
tracking efficiencies in the real experiment to the event samples generated by
the transport UrQMD model [71]. The experimental detection efficiency of
tracks can be gained from Monte-Carlo embedding technique which considering
effects including detector efficiency, acceptance, decays and interaction
losses in the collisions [56, 76]. In the STAR experiment at RHIC, the
efficiency of the Time Projection Chamber (TPC) or the Time-of-Flight (TOF)
detector is a function of transverse momentum [56, 76]. Figure 2 (a) and 3 (a)
show the detection efficiencies of protons as a function of $p_{T}$ in the TPC
detector and in the TPC+TOF detectors, respectively.

The tracking efficiency is then employed into the UrQMD event samples by
keeping or discarding a particle based on the probability obtained from the
$p_{T}$-dependent efficiencies. Figure 2 (b) shows the original true
$F_{2}(M)$ (open squares), measured $F_{2}(M)$ (solid stars) after discarding
particles and efficiency corrected SFM (solid cross) calculated by the cell-
by-cell method, respectively. We observe that the values of the measured
$F_{2}(M)$ are smaller than those of the original true ones, especially in the
region of large $M^{2}$. However, the efficiency corrected $F_{2}(M)$ agree
well with the original true ones. Fig. 3 (b) shows the calculations for the
case of TPC+TOF detection efficiency. The corrected SFMs are confirmed to be
consistent with the original true ones.

Therefore, the proposed cell-by-cell method provides an accurate and effective
way of efficiency corrections in the measurement of SFMs in heavy-ion
collisions. This method has been adopted recently in the intermittency
analysis in the STAR experiment [43].

##  5 Search for Intermittency in Experiments

The exploration of experimental evidence of creation of the QGP and location
of the CP is one of the major objectives in current heavy-ion collisions. The
goal is pursued by measurements of correlations and fluctuations including
intermittency, which are supposed to be promising signatures of the QCD phase
transition.

###  5.1 Results from NA49 and NA61 Collaborations at SPS

In the experiments at SPS, the search for the QCD critical point is attempted
by NA49 experiment through changing system sizes of colliding nuclei (p+p,
C+C, Si+Si, Pb+Pb) at 158$A$ GeV/$c$ and NA61/SHINE experiment by varying
energies in p+p,p+Pb, Be+Be, Ar+Sc and Xe+La collisions.

Figure 4: The measured $F_{2}(M)$ as a function of $M^{2}$ of proton density
in the most central collisions at $\sqrt{s_{\mathrm{NN}}}$ = 17.3 GeV for (a)
C+C, (b) Si+Si, and (c) Pb+Pb collisions. Figure taken from Ref [40].

Figure 4 presents $F_{2}(M)$ as a function of $M^{2}$ for data and mixed
events at $\sqrt{s_{\mathrm{NN}}}$ = 17.3 GeV in C+C, Si+Si and Pb+Pb
collisions measured in NA49 experiment. It is observed that $F_{2}(M)$ in the
Si+Si collisions are larger than those calculated from mixed events when
$M^{2}$ is large. The correlator $\Delta F_{2}(M)$ of Si+Si collision
calculated by Eq. 5 is found to obey a good power-law dependence on $M^{2}$
with $\phi_{2}=0.96\pm 0.16$ which approaches theoretical prediction [34],
indicating a typical property of intermittency in this collision. However, the
$F_{2}(M)$ are almost overlap with those of mixed events in the C+C and Pb+Pb
systems, suggesting that intermittency is not visible in these collisions. The
reason could be that the critical self-similar fluctuations can not develop in
the small size of C+C system. As for the Pb+Pb collision, the signal might be
diluted during the longer evolution of the hadronic phase [40].

Figure 5: The measured $F_{2}(M)$ as a function of $M^{2}$ of proton density
in 0-20% central Ar+Sc collisions at 150$A$ GeV/$c$ (left pad) and in 0-10%
central Pb+Pb collisions at 30$A$ GeV/$c$ (right pad). Figure taken from Ref
[77].

Figure 5 illustrates the preliminary results from NA61 on the second-order
SFMs of proton numbers in 0-20% central Ar+Sc collisions at 150$A$ GeV/$c$ and
in 0-10% central Pb+Pb collisions at 30$A$ GeV/$c$ [77]. The cumulative
variable method is used here to subtract background effects of the
measurement. It is observed that $F_{2}(M)$ is nearly flat with increasing
$M^{2}$, indicating the absence of an intermittency behavior in Ar+Sc or in
Pb+Pb collisions.

###  5.2 Results from STAR Experiment at RHIC

To explore the QCD phase diagram, the STAR Collaboration has measured a few
experimental observables which may signature the QCD phase transition and CP
at RHIC BES-I energies. In particularly, some appealing non-monotonic
behaviors of the measurements have been observed in Au+Au collisions at
energies around $20<\sqrt{s_{NN}}<30$ GeV [25, 80, 81, 79, 78].

Figure 6: $F_{q}(M)$ (q=2-6) as function of $M^{2}$ for charged particles in
the 0-5% most central Au+Au collisions at $\sqrt{s_{\mathrm{NN}}}$ = 7.7-200
GeV in double-logarithmic scale. Red (black) marks represent $F_{q}(M)$ of
data (mixed event) as a function of $M^{2}$. Figure taken from Ref [43].

Recently, the STAR Collaboration reports the preliminary results on the first
measurements on intermittency of charged particles in Au+Au collisions at
$\sqrt{s_{\mathrm{NN}}}$= 7.7-200 GeV [43]. The SFMs of $p$, $\bar{p}$,
$K^{\pm}$ and $\pi^{\pm}$ have been measured within $|\eta|<0.5$ in a 2D
transverse momentum space. The SFMs are corrected for the finite tracking
reconstruction efficiencies by the cell-by-cell method. Figure 6 shows the
$F_{q}(M)$ of both the data (solid symbols) and the mixed events (open
symbols) for identified charged particles in the $0-5\%$ most central Au+Au
collision at BES-I energies. The $F_{q}(M)$ can be calculated up to the sixth-
order based on the event statistics. It is observed that $F_{q}(M)^{data}$ is
larger than $F_{q}(M)^{mix}$ at the region of large $M^{2}$ values at various
energies and thus a deviation of $\Delta F_{q}(M)$ from zero is found in these
collisions.

Figure 7: Energy dependence of the scaling exponent $\nu$ for charged
particles in Au+Au collisions at $\sqrt{s_{\mathrm{NN}}}$ = 7.7-200 GeV.
Figure taken from Ref [43].

After subtracting the background contributions, the critical related exponent
$\nu$, which extracted from the $\Delta F_{q}(M)/\Delta F_{2}(M)$ scaling
according to Eq. (3) and (4), exhibits a non-monotonic behavior on energy with
a dip located between 20 GeV and 30 GeV in the most central Au+Au collisions,
as shown in Fig. 7. The non-monotonic energy dependence of $\nu$ needs to be
understood with more theoretical inputs, especially a reliable theoretical
calculation in a 2D transverse momentum space with the same acceptance cuts
used in the experimental analysis.

##  6 Intermittency Analysis in Models

###  6.1 Scaled Factorial Moments of Protons in the Transport UrQMD Model

In this section, we illustrate the SFMs of proton density calculated in
different energies and centralities by using the UrQMD model. The UrQMD (Ultra
relativistic Quantum Molecular Dynamics) is a transport model [82, 83] which
is widely used to simulate relativistic heavy-ion collisions with a wide
energy coverage ranging from SIS energies to top of the RHIC energy. Since the
model does not include any correlations or fluctuations related to the QCD
phase transition or CP, it is feasible to be used to study background effects
or other trivial fluctuations in the search of the QCD CP in heavy-ion
collisions.

The cascade mode of the UrQMD model with the version 2.3 is applied to produce
events in Au+Au collision at $\sqrt{s_{\mathrm{NN}}}$ = 7.7, 11.5, 19.6, 27,
39, 62.4 and 200 GeV. In the model calculation, we apply the same kinematic
cuts and analysis methods as those used in the STAR experiment [76]. Protons
are selected at mid-rapidity($|y|<0.5$) in transverse momentum space
$0.4<p_{T}<2.0$. Centrality is defined by charged pions and kaons in pseudo-
rapidity region $|\eta|<1.0$.

Figure 8: The second-order SFM (black squares) as a function of $M^{2}$ in the
0-5% most central Au+Au collisions at $\sqrt{s_{\mathrm{NN}}}$ = 7.7-200 GeV.
The black solid lines are the power-law fitting according to Eq. (2). The
corresponding red ones are those calculated by the cumulative variable method.
Figure 9: The second-order SFM (black squares) as a function of $M^{2}$
measured in various collision centralities in Au+Au collisions at
$\sqrt{s_{\mathrm{NN}}}$ = 19.6 GeV. The corresponding red ones are those
calculated by the cumulative variable method.

Figure 8 shows the directly measured second-order SFMs (black squares) and the
ones with background subtraction by the cumulative variable method (red
stars), as a function of the number of partitioned cells for protons in the
0-5% central Au+Au collisions at $\sqrt{s_{\mathrm{NN}}}$ = 7.7-200 GeV. It is
found that directly measured $F_{2}(M)$ increases slowly with increasing
$M^{2}$. By using the cumulative method to subtract the background
contribution, the $CF_{2}(M)$ shows a nearly flat trend with $M^{2}$. The
value of $\phi_{2}^{c}$, which is fitted from SFMs by the cumulative variable
method, is found to be near zero at the measured energies.

In Fig. 9, we show the $F_{2}(M)$ and $CF_{2}(M)$ at $\sqrt{s_{\mathrm{NN}}}$
= 19.6 GeV in various collision centralities. The values of $\phi_{2}$
calculated from directly measured SFMs slightly rise from the most central to
the most peripheral collisions. However, $\phi_{2}^{c}$ calculated from
$CF_{2}(M)$ show a flat centrality dependence with all the values near zero.
The results verify that the background of the non-critical contributions can
be efficiently wiped off by using the proposed cumulative variable method in
the UrQMD model.

###  6.2 Intermittency in the CMC Model belonging to the 3D Ising Universality
Class

Based on an effective action of the three-dimensional Ising universality
class, there exists a power-law singularity of density-density correlator for
small momentum transfer $\vec{k}$ [34]:

$\lim\limits_{|\vec{k}|\rightarrow
0}\langle\rho_{\vec{k}}\rho_{\vec{k}}^{*}\rangle\sim|\vec{k}|^{-d_{F}}.$ (11)

Here, $\rho_{\vec{k}}$ represents the baryon-number density in momentum space
by Fourier transforms from the one in coordinate space.
$\langle\rho_{\vec{k}}\rho_{\vec{k}}^{*}\rangle$ is the density-density
correlator of particle pairs in momentum space. Eq. (11) shows a fractal
geometry with a power-law pattern in momentum space for critical systems
belonging to the 3D Ising universality class.

Figure 10: $F_{2}(M)$ as a function of $M^{2}$ from the CMC events (black
circles) and from the UrQMD model (red triangles) in Au+Au collisions at
$\sqrt{s_{\mathrm{NN}}}$ = 19.6 GeV in a double-logarithmic scale. Figure
taken from Ref [46].

In the study on intermittency in heavy-ion collisions, the CMC model [34] is
widely used to simulate critical event samples incorporating self-similar
fluctuations. The Monte-Carlo simulations involving critical density
fluctuations related to Eq. (11) can be simulated by the algorithm of Levy
random walk, which requiring the probability density $\rho(p)$ between two
adjacent walks follows:

$\rho(p)=\frac{\nu p_{\rm min}^{\nu}}{1-(p_{\rm min}/p_{\rm
max})^{\nu}}p^{-1-\nu}.$ (12)

Here $\nu$ is the Levy exponent directly connected to the intermittency index,
$p$ represents the momentum distance of two adjacent particles, $p_{\rm min}$
and $p_{\rm max}$ are the minimum and maximum values of $p$. The parameters of
Levy function are set to be $\nu=1/6$ and $p_{\rm min}/p_{\rm max}=10^{-7}$
for critical events leading to the second-order intermittency
$\phi_{2}=\frac{5}{6}$. The precise rules of the Levy algorithm and more
details about the CMC model could be found in Refs [84, 34].

In Fig. 10, the black circles show $F_{2}(M)$ as a function of $M^{2}$ and the
solid black line is the fitting according to Eq. (2). The results are
calculated from a generated sample of 600 critical events. The Levy algorithm
is parameterized to produce $p_{x}$ and $p_{y}$ distributions of particles in
transverse momentum space and the multiplicity is set to obey a Poisson
distribution with the mean value $\langle N\rangle=20$. It is found that
$F_{2}(M)$ from CMC events does obey a good power-law dependence on $M^{2}$
with the fitting slope $\phi_{2}=0.834\pm 0.001$. It infers that the CMC model
could excellently reproduce the intermittency behavior related to the self-
similar correlations of Eq. (11). The open red triangles are the results
calculated from the UrQMD event sample by the same mean multiplicity as used
in the CMC model. It gives a flat trend of the dependence on $M^{2}$. The
reason is that this model does not incorporate any self-similar fluctuations
in the mechanism of particle production process.

The definition of the relative density fluctuation of baryons, $\Delta n$, is
[30, 31]

$\Delta n=\frac{\langle(\delta n)^{2}\rangle}{\langle
n\rangle^{2}}=\frac{\langle n^{2}\rangle-\langle n\rangle^{2}}{\langle
n\rangle^{2}}.$ (13)

The angle bracket represents the average over the whole event sample.

Figure 11: The second-order intermittency index as a function of relative
density fluctuation in the CMC model. The dash lines represent the
experimental measured $\Delta n$ at $\sqrt{s_{\mathrm{NN}}}$ = 7.7-200 GeV
from the STAR experiment [79]. Figure taken from Ref [46].

Measurements of SFMs to get intermittency index and calculations of $\Delta n$
can be performed simultaneously in the same CMC event sample [46]. The solid
black line in Fig. 11 displays the second-order intermittency index as a
function of the relative density fluctuation. We observe a monotonically
increase of $\phi_{2}$ on $\Delta n$. It reveals a fact that great
intermittency can be achieved if large baryon density fluctuations are
developed for the system near the QCD critical point. It supplies another
experimentally measurable variable to obtain density fluctuations besides
measurement of light nuclei productions in the coalescence model [30, 31]. The
colored dash lines in the same figure show the STAR measured $\Delta n$ in the
0-10% central Au+Au collisions at $\sqrt{s_{\mathrm{NN}}}$ = 7.7, 11.5, 14.5,
19.6, 27, 39, 62.4 and 200 GeV, respectively [79].

Figure 12: The collision energy dependence of the intermittency index
$\phi_{2}$ in various collision systems. The solid blue and green symbols
represent the values from the NA49 [41, 40] and NA61 experiments [67, 68]. The
red stars are the $\phi_{2}$ indirectly extracted from the measured light
nuclei productions in the STAR experiment [79]. The black circles show the
intermittency index of protons from UrQMD model. The black arrow is the
theoretic expectation for a critical model [34].

Once the connection between $\phi_{2}$ and $\Delta n$ is gained, the
intermittency index can be obtained by measuring the density fluctuations
within the same event sample, or vice versa. In this way, we can calculate
$\phi_{2}$ indirectly by mapping the experimental measured $\Delta n$ into the
relation. In Fig. 12, the red stars show the energy dependence of indirectly
calculated $\phi_{2}$ in the 0-10% central Au+Au collisions. We observe a non-
monotonic dependence of $\phi_{2}$ on collision energy with a peak structure
at around 20-30 GeV. It indicates that the strength of intermittency is
strongest in this energy region. For comparison, data from the NA49 [41, 40]
(blue symbols) and the NA61/SHINE [67, 68] (green crosses) experiments are
also plotted in the same figure. In Si+Si collisions, the $\phi_{2}$ is found
to be close to the value of theoretical expectation which is illustrated as
the black arrow. The black circles represent $\phi_{2}$ from the UrQMD
calculations, which give a flat energy dependence with all the values around
zeros. It is because that no critical self-similar mechanism is implemented in
this model.

###  6.3 Intermittency Analysis in the UrQMD Model with Hadronic Potentials

As has been shown in Section 6.1, there is no intermittency behavior in the
original UrQMD model after subtracting background contributions. Recently,
attempts have been made to give obvious correlations of proton pairs and
intermittency by introducing hadronic potentials into the UrQMD model with the
mean-field mode [47].

By treating potentials of both the formed and performed hadrons from string
fragmentation by the same means, the density-dependent potentials can be
written as [47, 85]:

$U=\alpha(\frac{\rho_{h}}{\rho_{0}})+\beta(\frac{\rho_{h}}{\rho_{0}})^{\gamma},$
(14)

where $\rho_{h}$ represents the hadronic density, and $\rho_{0}$ is the
nuclear matter saturation density.

The momentum-dependent term of hadronic potential is defined as:

$U_{md}=\sum_{k=1,2}\frac{t^{k}_{md}}{\rho_{0}}\int
d\bm{p^{\prime}}\frac{f({\bm{r},\bm{p^{\prime}}})}{1+[(\bm{p}-\bm{p^{\prime}})/a_{md}^{k}]^{2}}.$
(15)

Here $t_{md}$ and $a_{md}$ are parameters. More details about the
implementation of hadronic potentials into the model could be found in Refs
[47, 85].

Figure 13: (a) The second-order SFM as a function of $M^{2}$ of proton density
in the 5-10% central Ar+Sc collisions at 40$A$ GeV/$c$. The solid squares
(circles) are $F_{2}(M)$ of real events (mixed events) from the UrQMD model
with hadronic potentials (UrQMD/M). The corresponding open squares and circles
represent the measurements without hadronic potentials (UrQMD/C). (b) $\Delta
F_{2}(M)$ as a function of $M^{2}$ of protons. Figure taken from Ref [47].

Figure 13 illustrates the second-order $F_{2}(M)$ and the correlator $\Delta
F_{2}(M)$ calculated in the UrQMD data samples both with and without hadronic
potentials in Ar+Sc collisions at 40A GeV/$c$. As shown in the left figure of
the UrQMD model with hadronic potential (UrQMD/M), $F_{2}(M)$ of original
events (solid blue squares) are larger than those of mixed events (solid blue
circles) at the large $M^{2}$ regions. In Fig. 13 (b), the correlator $\Delta
F_{2}(M)$ of hadronic potentials shows a power-law behavior with $M^{2}$. The
extracted $\phi_{2}=0.32\pm 0.03$ is similar to the one measured in the NA61
experiment [86, 68]. As for the UrQMD cascade model (UrQMD/C), the values of
$F_{2}(M)$ of real events are almost overlapped with those of mixed events and
thus $\Delta F_{2}(M)$ is around 0 in this case. It could infer that the
power-law character of $\Delta F_{2}(M)$ in the UrQMD/M model is introduced by
the hadronic interaction, particular nuclear potentials which cause
enhancements of proton pairs with small relative momenta [47].

Figure 14: The extracted intermittency index $\phi_{2}$ at different
centralities for Ar+Sc collisions in the UrQMD/M and UrQMD/C models. Figure
taken from Ref [47].

Figure 14 demonstrates centrality dependence of $\phi_{2}$ in 0-5%, 5-10%, and
10-15% central Ar+Sc collisions, respectively. The calculations from the UrQMD
model with (solid blue squares) and without (open pink squares) hadronic
potential are plotted along with the NA61 preliminary experimental results
(solid red stars). The second-order intermittency index $\phi_{2}$ calculated
in the UrQMD/C model is found to be nearly zero. For the UrQMD/M model, the
$\phi_{2}$ slightly increases from the most central to the mid-central
collisions, in coincide with the experimental result for the cases of purity
$>85\%$ [86, 68]. This centrality dependence could be explained by a shortened
hadronic freeze-out phase with decreasing size of the system [47].

###  6.4 Observation of Intermittency in the Hybrid UrQMD-hydro Model

Apart from the introduce of hadronic potentials into the transport UrQMD model
to obtain intermittency, another attempt is the hybrid UrQMD-hydro model which
incorporating both transport and hydrodynamical descriptions of heavy-ion
collisions [48, 87]. Hydrodynamic description has been suggested as an
effective tool to express the hot and dense stage of collision reactions and
model phase transitions [87]. In this study, the intermediate hydrodynamic
calculation is applied to the microscopic transport computation in initial
condition and freeze-out process [87, 83].

Figure 15: $\ln F_{q}(M)$ vs $\ln M^{2}$ of charged particles in Au+Au
collisions at 10A GeV/$c$ from the UrQMD-hydro model with hadronic EoS (red
circles) and chiral one (blue squares). Figure taken from Ref [48].

In Fig. 15, the values of $\ln F_{q}(M)$ are shown as a function of $\ln
M^{2}$ of charged particles in Au+Au collisions at 10A GeV/$c$ from the UrQMD-
hydro model with two different equations of state (EoS). The $F_{q}(M)$ of
various orders are carried out in $\eta-\phi$ space and the background has
been subtracted by the cumulative variable method. The value of $\ln F_{q}(M)$
is found to linearly develop with increasing $\ln M^{2}$ for both hadronic EoS
(red circles) and chiral one (blue squares) in the (0–5%) most central
collisions. It implies that the observed scaling behavior or intermittency is
closely related to the evolution of the medium produced in the collision. In
addition, the intermittency indices of chiral EoS are found to be larger than
those of hadronic one. It could be account for the production of cascading
particles in partonic level and the hydrodynamic evolution in the chiral EoS
case [48].

##  7 Summary and Outlook

In this review, we have summarized current status of the exploration of the
QCD critical point via intermittency analysis in heavy-ion collisions. The
main results from both experimental measurements and phenomenological model
calculations are presented. We highlighted two important issues, i.e.
background subtraction and efficiency correction, in this analysis. It
provides a way to obtain a clean signal in the calculations of intermittency.

Although the main contributions from background have been eliminated by using
the proposed cumulative variable method, the intermittency index calculated in
the transport UrQMD model which does not incorporate any critical self-similar
fluctuations is not exact zero, as shown in Fig. 8 and 9. It might attribute
to some effects which still remained in the results, such as correlation of
protons caused by Coulomb repulsion and Fermi–Dirac statistics [40], or width
of the experimental momentum resolution [54]. More investigations on these
effects should help to obtain a clear and convincing result.

On the search of self-similar and intermittency behavior in current heavy-ion
experiments, the NA49 Collaboration has found a clear signature in Si+Si
collisions at $\sqrt{s_{\mathrm{NN}}}$ = 17.3 GeV [40], with the measured
second-order intermittency index touching the theoretic expectation calculated
in a critical model [34]. Preliminary results from NA61/SHINE [77, 88] of
proton numbers in Ar+Sc collisions at 150A GeV/$c$ and in the most central
Pb+Pb collisions at 30A GeV/$c$ exhibit no evidence of power-law increase with
increasing number of division cells. However, the RHIC/STAR measurement [43]
at the near energy region illustrates a non-monotonic dependence of a scaling
exponent on collision energy in the most central Au+Au collisions with a dip
at around 20-30 GeV. Further investigations should be done to compare results
coming from both the STAR and the NA61 heavy-ion experimental data despite the
differences between two detectors.

The RHIC experiment has updated a few detectors and finished taking the second
phase of Beam Energy Scan program during 2018-2021 [14, 15]. It could be
exciting and anticipated if intermittency would be measured by the STAR
Collaboration on the new data with high statistics and improved particle
identification to explore the CP in the QCD phase diagram.

In addition to the benefits from the upgraded facilities discussed above, the
large data samples that collected in the high luminosity phase at BES-II may
bring computing difficulties since the calculation of SFMs at very small
scales requires a huge computational effort which may prevent its
implementation in experimental analysis. A new computational technique,
compared to conventional methods, is suggested to be much more efficient and
prominent in the calculations of the second-order SFM with increase of the
dimension of space [89]. The technique is also supposed to be a useful tool
for recognizing weak signals which may be hidden in the background with strong
noises.

With the development of modern computer hardware and artificial intelligence,
machine learning (ML) [90, 91], as a data-driven method, may be another
intriguing direction for future investigations. Based on a recent study [92],
a dynamical edge convolution plus point cloud neural network shows a strong
pattern recognition ability in identifying events with self-similar
fluctuations. It can figure out the most majority of signal particles which
could be used for decision-making in Monte-Carlo events. It is interesting to
investigate whether this cutting-edge ML method would be also helpful to pick
out possibly existing weak intermittency signals associated with critical
phenomena, which encountered in the NA61 experiment.

## Acknowledgments

I am grateful to many colleagues for conversations that have shaped my
understanding of self-similar property and intermittency over the years. In
particularly, I would like to express my sincere appreciations to Prof.
Lianshou Liu, Yuanfang Wu, Feng liu, Nu Xu, W. Kittel, A. Bialas, W. Metzger,
Jinghua Fu, Fuming Liu, Xiaofeng Luo and Lizhu Chen for collaborations and
instructive suggestions. I further thank the STAR and the NA61 Collaborations
for valuable discussions and comments. This work is supported by the Ministry
of Science and Technology (MoST) under grant No. 2016YFE0104800.

## References

  * [1] M. A. Stephanov, K. Rajagopal and E. V. Shuryak, Phys. Rev. Lett. 81, 4816 (1998), arXiv:hep-ph/9806219.
  * [2] STAR Collaboration, J. Adams et al., Nucl. Phys. A 757, 102 (2005), arXiv:nucl-ex/0501009.
  * [3] M. Asakawa, U. W. Heinz and B. Muller, Phys. Rev. Lett. 85, 2072 (2000), arXiv:hep-ph/0003169.
  * [4] V. Koch, A. Majumder and J. Randrup, Phys. Rev. Lett. 95, 182301 (2005), arXiv:nucl-th/0505052.
  * [5] Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz and K. K. Szabo, Nature 443, 675 (2006), arXiv:hep-lat/0611014.
  * [6] Y. Aoki, Z. Fodor, S. D. Katz and K. K. Szabo, Phys. Lett. B 643, 46 (2006), arXiv:hep-lat/0609068.
  * [7] M. Cheng et al., Phys. Rev. D 74, 054507 (2006), arXiv:hep-lat/0608013.
  * [8] E. S. Bowman and J. I. Kapusta, Phys. Rev. C 79, 015202 (2009), arXiv:0810.0042 [nucl-th].
  * [9] S. Ejiri, Phys. Rev. D 78, 074507 (2008), arXiv:0804.3227 [hep-lat].
  * [10] M. A. Stephanov, Prog. Theor. Phys. Suppl. 153, 139 (2004), arXiv:hep-ph/0402115.
  * [11] Z. Fodor and S. D. Katz, JHEP 04, 050 (2004), arXiv:hep-lat/0402006.
  * [12] R. V. Gavai and S. Gupta, Phys. Rev. D 78, 114503 (2008), arXiv:0806.2233 [hep-lat].
  * [13] R. V. Gavai and S. Gupta, Phys. Rev. D 71, 114014 (2005), arXiv:hep-lat/0412035.
  * [14] STAR Note 0598: BES-II whitepaper. https://drupal.star.bnl.gov/STAR/starnotes/public/sn0598.
  * [15] STAR Collaboration, C. Yang, Nucl. Phys. A 967, 800 (2017).
  * [16] M. A. Stephanov, Phys. Rev. Lett. 102, 032301 (2009), arXiv:0809.3450 [hep-ph].
  * [17] M. A. Stephanov, Phys. Rev. Lett. 107, 052301 (2011), arXiv:1104.1627 [hep-ph].
  * [18] B. J. Schaefer and M. Wagner, Phys. Rev. D 85, 034027 (2012), arXiv:1111.6871 [hep-ph].
  * [19] S. Gupta, X. Luo, B. Mohanty, H. G. Ritter and N. Xu, Science 332, 1525 (2011), arXiv:1105.3934 [hep-ph].
  * [20] STAR Collaboration, M. M. Aggarwal et al., Phys. Rev. Lett. 105, 022302 (2010), arXiv:1004.4959 [nucl-ex].
  * [21] STAR Collaboration, L. Adamczyk et al., Phys. Rev. Lett. 112, 032302 (2014), arXiv:1309.5681 [nucl-ex].
  * [22] STAR Collaboration, L. Adamczyk et al., Phys. Rev. Lett. 113, 092301 (2014), arXiv:1402.1558 [nucl-ex].
  * [23] STAR Collaboration, L. Adamczyk et al., Phys. Lett. B 785, 551 (2018), arXiv:1709.00773 [nucl-ex].
  * [24] PHENIX Collaboration, A. Adare et al., Phys. Rev. C 93, 011901 (2016), arXiv:1506.07834 [nucl-ex].
  * [25] STAR Collaboration, J. Adam et al., Phys. Rev. Lett. 126 (2021), arXiv:2001.02852 [nucl-ex].
  * [26] X. Luo and N. Xu, Nucl. Sci. Tech. 28, 112 (2017), arXiv:1701.02105 [nucl-ex].
  * [27] S. Jeon and V. Koch, Phys. Rev. Lett. 85, 2076 (2000), arXiv:hep-ph/0003168.
  * [28] Y. Hatta and T. Ikeda, Phys. Rev. D 67, 014028 (2003), arXiv:hep-ph/0210284.
  * [29] F. Li and C. M. Ko, Phys. Rev. C 93, 035205 (2016), arXiv:1601.00026 [nucl-th].
  * [30] K.-J. Sun, L.-W. Chen, C. M. Ko and Z. Xu, Phys. Lett. B 774, 103 (2017), arXiv:1702.07620 [nucl-th].
  * [31] K.-J. Sun, L.-W. Chen, C. M. Ko, J. Pu and Z. Xu, Phys. Lett. B 781, 499 (2018), arXiv:1801.09382 [nucl-th].
  * [32] H. Brumberger, N. G. Alexandropoulos and W. Claffey, Phys. Rev. Lett. 19, 555 (Sep 1967).
  * [33] B. Berdnikov and K. Rajagopal, Phys. Rev. D 61, 105017 (2000), arXiv:hep-ph/9912274.
  * [34] N. G. Antoniou et al., Phys. Rev. Lett. 97, 032002 (2006).
  * [35] H. Satz, Nucl. Phys. B 326, 613 (1989).
  * [36] E. A. De Wolf et al., Phys. Rept. 270, 1 (1996).
  * [37] N. G. Antoniou, F. K. Diakonos and A. S. Kapoyannis, Phys. Rev. C 81, 011901 (2010), arXiv:0809.0685 [hep-ph].
  * [38] N. G. Antoniou, N. Davis and F. K. Diakonos, Phys. Rev. C 93, 014908 (2016), arXiv:1510.03120 [hep-ph].
  * [39] N. G. Antoniou, F. K. Diakonos, X. N. Maintas and C. E. Tsagkarakis, Phys. Rev. D 97, 034015 (2018), arXiv:1705.09124 [hep-ph].
  * [40] NA49 Collaboration, Anticic et al., Eur. Phys. J. C 75, 587 (2015), arXiv:1208.5292 [nucl-ex].
  * [41] NA49 Collaboration, Anticic et al., Phys. Rev. C 81, 064907 (2010), arXiv:0912.4198 [nucl-ex].
  * [42] NA61/SHINE Collaboration, M. Maćkowiak-Pawłowska, Nucl. Phys. A 1005, 121753 (2021), arXiv:2002.04847 [nucl-ex].
  * [43] STAR Collaboration, J. Wu, Measurement of Intermittency for Charged Particles in Au + Au Collisions at $\sqrt{s_{\mathrm{NN}}}$ = 7.7-200 GeV from STAR, in 50th International Symposium on Multiparticle Dynamics, (10 2021).  arXiv:2110.09794 [nucl-ex].
  * [44] M. V. Tokarev and I. Zborovsky (9 2021), arXiv:2109.01524 [hep-ph].
  * [45] Z. Ong, P. Agarwal, H. W. Ang, A. H. Chan and C. H. Oh, Intermittency in $pp$ collisions at $\sqrt{s}=$ 0.9, 7 and 8 TeV from the CMS collaboration, in 50th International Symposium on Multiparticle Dynamics, (10 2021).  arXiv:2110.13087 [hep-ph].
  * [46] J. Wu, Y. Lin and Y. W. Li, Zhiming, Phys. Lett. B 801, 135186 (2020), arXiv:1901.11193 [nucl-th].
  * [47] P. Li et al., Phys. Lett. B 818, 136393 (2021), arXiv:2104.08483 [nucl-th].
  * [48] S. Gope and B. Bhattacharjee, Eur. Phys. J. A 57, 44 (2021), arXiv:2008.01943 [nucl-th].
  * [49] S. Ahmad, M. M. Khan, S. Khan, A. Khatun and M. Irfan, Int. J. Mod. Phys. E 23, 1450065 (2014), arXiv:1510.03176 [nucl-ex].
  * [50] Y.-L. Xie, G. Chen, J.-L. Wang, Z.-H. Liu and M.-J. Wang, Nucl. Phys. A 920, 33 (2013), arXiv:1303.2183 [nucl-th].
  * [51] R. Gupta and S. K. Malik, Adv. High Energy Phys. 2020, 5073042 (2020), arXiv:1911.13111 [hep-ex], [Erratum: Adv.High Energy Phys. 2020, 7319894 (2020)].
  * [52] A. Bialas and R. B. Peschanski, Nucl. Phys. B 273, 703 (1986).
  * [53] K. Fialkowski, B. Wosiek and J. Wosiek, Acta Phys. Polon. B 20, 639 (1989).
  * [54] S. Samanta, T. Czopowicz and M. Gazdzicki, J. Phys. G 48, 105106 (2021), arXiv:2105.01763 [nucl-th].
  * [55] S. Samanta, T. Czopowicz and M. Gazdzicki, Nucl. Phys. A 1015, 122299 (2021), arXiv:2105.00344 [nucl-th].
  * [56] STAR Collaboration, L. Adamczyk et al., Phys. Rev. C 96, 044904 (2017), arXiv:1701.07065 [nucl-ex].
  * [57] X. Luo, Phys. Rev. C 91, 034907 (2015), arXiv:1410.3914 [physics.data-an], [Erratum: Phys.Rev.C 94, 059901 (2016)].
  * [58] A. Bialas and R. B. Peschanski, Nucl. Phys. B 308, 857 (1988).
  * [59] S. Gavin, A. Gocksch and R. D. Pisarski, Phys. Rev. D 49, R3079 (1994), arXiv:hep-ph/9311350.
  * [60] A. M. Halasz, A. D. Jackson, R. E. Shrock, M. A. Stephanov and J. J. M. Verbaarschot, Phys. Rev. D 58, 096007 (1998), arXiv:hep-ph/9804290.
  * [61] F. Karsch, E. Laermann and C. Schmidt, Phys. Lett. B 520, 41 (2001), arXiv:hep-lat/0107020.
  * [62] N. G. Antoniou, Y. F. Contoyiannis, F. K. Diakonos and G. Mavromanolakis, Nucl. Phys. A 761, 149 (2005), arXiv:hep-ph/0505185.
  * [63] R. C. Hwa and M. T. Nazirov, Phys. Rev. Lett. 69, 741 (1992).
  * [64] R. C. Hwa, Phys. Rev. D 47, 2773 (1993).
  * [65] R. C. Hwa and C. B. Yang, Phys. Rev. C 85, 044914 (2012), arXiv:1111.6651 [nucl-th].
  * [66] A. Kamal, N. Ahmad and M. M. Khan, Acta Phys. Polon. B 46, 1549 (2015).
  * [67] N. G. Antoniou et al., Nucl. Phys. A 1003, 122018 (2020), arXiv:2004.12133 [nucl-th].
  * [68] NA61/SHINE Collaboration, D. Prokhorova and N. Davis, Universe 5, 103 (2019).
  * [69] A. Bialas and M. Gazdzicki, Phys. Lett. B 252, 483 (1990).
  * [70] W. Ochs, Z. Phys. C 50, 339 (1991).
  * [71] J. Wu, Y. Lin, Z. Li, X. Luo and Y. Wu, Phys. Rev. C 104, 034902 (2021), arXiv:2104.11524 [nucl-th].
  * [72] A. Bzdak and V. Koch, Phys. Rev. C 91, 027901 (2015), arXiv:1312.4574 [nucl-th].
  * [73] V. Vovchenko and V. Koch, Nucl. Phys. A 1010, 122179 (2021), arXiv:2101.02182 [nucl-th].
  * [74] X. Luo, Phys. Rev. C 91, 034907 (2015), arXiv:1410.3914 [physics.data-an].
  * [75] HADES Collaboration, J. Adamczewski-Musch et al., Phys. Rev. C 102, 024914 (2020), arXiv:2002.08701 [nucl-ex].
  * [76] STAR Collaboration, M. Abdallah et al., Phys. Rev. C 104 (2021), arXiv:2101.12413 [nucl-ex].
  * [77] NA61/SHINE Collaboration, M. Kuich, EPJ Web Conf. 259, 01001 (2022), arXiv:2112.07987 [nucl-ex].
  * [78] STAR Collaboration, L. Adamczyk et al., Phys. Rev. Lett. 112, 162301 (2014), arXiv:1401.3043 [nucl-ex].
  * [79] STAR Collaboration, D. Zhang, Nucl. Phys. A 1005, 121825 (2021), arXiv:2002.10677 [nucl-ex].
  * [80] STAR Collaboration, L. Adamczyk et al., Phys. Rev. C 92, 014904 (2015), arXiv:1403.4972 [nucl-ex].
  * [81] R. A. Lacey, Phys. Rev. Lett. 114, 142301 (2015), arXiv:1411.7931 [nucl-ex].
  * [82] S. A. Bass et al., Prog. Part. Nucl. Phys. 41, 255 (1998), arXiv:nucl-th/9803035.
  * [83] M. Bleicher et al., J. Phys. G 25, 1859 (1999), arXiv:hep-ph/9909407.
  * [84] P. A. Alemany and D. H. Zanette, Phys. Rev. E 49, R956 (Feb 1994).
  * [85] Q. Li, M. Bleicher and H. Stocker, Phys. Lett. B 659, 525 (2008), arXiv:0709.1409 [nucl-th].
  * [86] Na61/Shine Collaboration, N. Davis, N. Antoniou and F. K. Diakonos, PoS CORFU2018, 154 (2019).
  * [87] H. Petersen, J. Steinheimer, G. Burau, M. Bleicher and H. Stöcker, Phys. Rev. C 78, 044901 (2008), arXiv:0806.1695 [nucl-th].
  * [88] NA61/SHINE Collaboration, M. Mackowiak-Pawlowska, Study of the phase diagram of strongly interacting matter in the NA61/SHINE experiment, in Particles and Nuclei International Conference , (12 2021).  arXiv:2112.01877 [nucl-ex].
  * [89] F. K. Diakonos and A. S. Kapoyannis (9 2021), arXiv:2109.12571 [hep-ph].
  * [90] G. E. HINTON and R. R. SALAKHUTDINOV, Science 313, 504 (Jul 2006).
  * [91] Y. B. Yann LeCun and G. Hinton, Nature 521, 436 (May 2015).
  * [92] Y.-G. Huang, L.-G. Pang, X. Luo and X.-N. Wang (Jul 2021), arXiv:2107.11828 [nucl-th].


# Extracting Explanations, Justification, and Uncertainty
from Black-Box Deep Neural Networks

Paul Ardis  GE Aerospace Research, 1 Research Circle, Niskayuna, NY 12309, USA
Arjuna Flenner  GE Aerospace, 3290 Patterson Avenue SE, Grand Rapids, MI
49512, USA

###### Abstract

Deep Neural Networks (DNNs) do not inherently compute or exhibit empirically-
justified task confidence. In mission critical applications, it is important
to both understand associated DNN reasoning and its supporting evidence. In
this paper, we propose a novel Bayesian approach to extract explanations,
justifications, and uncertainty estimates from DNNs. Our approach is efficient
both in terms of memory and computation, and can be applied to any black box
DNN without any retraining, including applications to anomaly detection and
out-of-distribution detection tasks. We validate our approach on the CIFAR-10
dataset, and show that it can significantly improve the interpretability and
reliability of DNNs.

## Introduction

Deep neural networks (DNNs) are powerful machine learning models that can
learn complex patterns from data. However, their complexity can make it
difficult to understand why a DNN makes a particular prediction purely from
their mathematical construction. We build upon the prior Explainable AI (XAI)
work from Virani et al. [1] (see Figure 1) which effectively explains if there
is proper training data to support the decisions. However, their memory and
computational footprint prevents the approach from being effective for very
large data sets and on edge devices. We exploit the concept of Sparse Gaussian
processes [2, 3, 4, 5] to overcome these two computational challenges while
maintaining the their method’s accuracy and explainability. Our approach is a
computationally efficient XAI method that extracts example-based
justifications and uncertainty estimates that can be applied to any pre-
trained DNN.

Explainable artificial intelligence addresses the black-box problem in many
ways: it decomposes the model construction for more intuitive human
understanding; it supplements models with tools that explain their decision-
making; or a structured combination of the previous two. By providing insights
into how AI systems make decisions, XAI can help humans to identify and
correct potential errors before they cause harm. As such, XAI is a critical
tool for ensuring the safety and human trust in mission-critical tasks such as
jet engine maintenance and airport safety. In the context of jet engine
maintenance, XAI can be used to provide detailed reasoning about the evidence
suggesting potential problems before they can cause failures, recommend
corrective actions to prevent potential issues, and monitor the performance of
engines to ensure that they continue to operate safely. In the context of
airport safety, XAI can be used to detect and explain potential hazards on the
airport grounds, expedite human understanding of the evidence behind security
threats, and monitor the movement of aircraft to ensure that they follow safe
takeoff and landing directions.

Figure 1: An illustration of the example-based XAI approach of Virani et al.
[1]. Their approach uses the transformation spaces from intermittent layers of
a pre-trained DNN to evaluate the model’s justification and builds a support
neighborhood around each test sample’s transformed point.

While there are a number of competing standards and meta-studies of XAI [6],
primary seminal papers explore XAI using one of three main approaches. The
first approach is to use gradient-based methods to identify the input features
that are most important for a particular prediction [7]. This can be done by
calculating the gradient of the loss function with respect to the input
features - the features with the largest gradients are the ones that have the
most influence on the prediction. The second approach to explainability is to
train the DNN to output text that explains its reasoning. This can be done by
using a technique called attention, which allows the DNN to focus on specific
input features when making a prediction [8]. The DNN can then be trained to
output text that describes the features that it focused on. The third approach
to explainability is to train the DNN using metric learning [9]. This involves
training the DNN to learn a distance metric that can be used to measure the
similarity between input features. The metric can then be used to explicitly
identify training set samples that are similar to any given test sample. The
DNN can also then be used to make predictions for the testing sample based
solely on the predictions that it made for the similar training samples.

Each of these approaches has its own advantages and disadvantages. Gradient-
based methods are relatively simple to implement, but they can be difficult to
interpret in their raw form. Text-based explanations can be more informative,
but they can be more difficult to generate without adding a complex model of
relevant human language. Metric learning can be used to generate both simple
and informative explanations, but it can be more computationally expensive as
it requires retraining the neural network.

We address the limitations of existing methods by presenting an approach to
extract example-based justifications and uncertainty estimates from pre-
trained DNNs. Our example-based justifications can be used to determine the
training set samples that were the most relevant to any given test sample. In
scenarios where there are an insufficient number of relevant training set
samples, we can conclude that the neural network is extrapolating and has a
high potential of misclassification. In addition, our uncertainty estimates
provide explicit information about the model’s prediction confidence based
upon local exemplar density and coherence, and this uncertainty can be
thresholded to limit operation to avoid potential misclassifications. Our
contributions are as follows:

  * •

We present an approach that uses sparse Gaussian processes to take the latent
embeddings from a pre-trained DNN and predict the model output. Such an
approach will allow us to estimate the prediction uncertainty while also
leveraging the DNN’s representation power.

  * •

We obtain example-based justifications for our model predictions by selecting
the training set samples that are highly correlated with the test samples. We
use the SGP’s kernel functions to determine how the samples are correlated.

### Advantages of our approach

  * •

Bayesian uncertainty quantification using sparse Gaussian process (SGPs [10,
3])

    * –

SGPs have $\mathcal{O}(nm^{2}+m^{3})$ computation cost [4], where n is the
number of training samples and m is the number of SGP inducing points

    * –

SGPs can be trained on large datasets using stochastic gradient descent [11]

    * –

SGPs support non-Gaussian likelihoods (e.g., softmax likelihood)

    * –

SGPs can be defined on Riemannian manifolds and leverage geodesics to
determine correlations [12]

  * •

Prediction justification from example-based explanations

    * –

Reduce neighborhood operators’ computation cost by leveraging SGP inducing
points

## Background: Gaussian Processes

Gaussian processes (GP), [13] can be used as priors in a non-parametric
Bayesian approach for regression and classification problems. Let
$\mathbf{x}\in\mathbb{R}^{d}$, then a Gaussian process is uniquely defined
using a mean function $m:\mathbb{R}^{d}\rightarrow\mathbb{R}$ and covariance
function $k:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$. The
covariance function $k$ must be a positive-definite kernel function which may
depend on a set of parameters $\Theta$. Given these two functions, a
stochastic process $f\in\mathbb{R}$ indexed by $\mathbf{x}\in\mathbb{R}^{d}$
such that

$\displaystyle m(\mathbf{x})$ $\displaystyle=$
$\displaystyle\mathbb{E}[f(\mathbf{x})],$ $\displaystyle
k(\mathbf{x}_{1},\mathbf{x}_{2})$ $\displaystyle=$
$\displaystyle\mathbb{E}[(f(\mathbf{x}_{1})-m(\mathbf{x}_{1}))(f(\mathbf{x}_{2})-m(\mathbf{x}_{2}))]$

is a GP. If $f$ is a Gaussian process with mean $m(\cdot)$ and covariance
kernel $k(\cdot,\cdot)$ write

$\displaystyle f(\mathbf{x})$ $\displaystyle\sim$
$\displaystyle\mathcal{GP}(f|m(\mathbf{x}),k(\mathbf{x},\mathbf{x})).$

For regression problems, the GP is often used to define a generative model of
the form:

$\displaystyle y_{i}$ $\displaystyle=$ $\displaystyle
f(\mathbf{x}_{i})+\epsilon_{i},$ $\displaystyle\epsilon_{i}$
$\displaystyle\sim$
$\displaystyle\mathcal{N}(\epsilon|0,\sigma_{\text{noise}}^{2}),$
$\displaystyle f(\mathbf{x}_{i})$ $\displaystyle\sim$
$\displaystyle\mathcal{GP}(f|m(\mathbf{x},k(\mathbf{x},\mathbf{x})).$

An essential concept in the following discussion is the number of variables
$\mathbf{x}$ used in the calculations. The following notation is used to
highlight the set sizes with the goal of illuminating the computational
complexity. Let
$\mathbf{X}_{n}=[\mathbf{x}_{1},\mathbf{x}_{2},...,\mathbf{x}_{n}]^{\top}$ be
a set of locations. Define $1\leq i,j\leq n$ and let $\mathbf{K}_{nn}$ be a
matrix with elements $k(\mathbf{x}_{i},\mathbf{x}_{j})$. Note that the
subscript $nn$ refers to the number of rows and number of columns and not a
row or column index. Define the random vector
$\mathbf{f}_{n}=[f(\mathbf{x}_{1}),f(\mathbf{x}_{2}),...,f(\mathbf{x}_{n})]^{\top}$.
Given the location $\mathbf{x}$, define
$\mathbf{k}_{\mathbf{x}n}=[k(\mathbf{x},\mathbf{x}_{1}),\dots,k(\mathbf{x},\mathbf{x}_{n})]^{\top}$.

Bayesian methods assume a prior data model and calculate a posterior
distribution once data is observed. Consider an observed data set
$\mathcal{D}_{n}=\\{(\mathbf{x}_{i},y_{i})\;|\;1\leq i\leq
n,\;\mathbf{x}_{i}\in\mathbb{R}^{n},\;y_{i}\in\mathbb{R}\\}$. Define the
vector $\mathbf{y}=[y_{1},\dots,y_{n}]^{\top}$ and assume that the mean
function $m(\cdot)=0$. With this mean function, recall that the definition of
a GP designates
$p(\mathbf{f}|\mathbf{X}_{n})=\mathcal{N}(\mathbf{f}|\mathbf{0},\mathbf{K}_{nn})$.
Since all the random variables in the GP regression model are Gaussian
distributed, the posterior distribution is also a GP with mean function
$m_{\mathbf{y}}(\cdot)$ and covariance kernel $k_{\mathbf{y}}(\cdot,\cdot)$
given by

$\displaystyle m_{\mathbf{y}}(\mathbf{x})$
$\displaystyle=\mathbf{k}_{\mathbf{x}n}^{\top}(\mathbf{K}_{nn}+\sigma_{\text{noise}}^{2}I)^{-1}\mathbf{y}\,,$
(1) $\displaystyle k_{\mathbf{y}}(\mathbf{x},\mathbf{x}^{\prime})$
$\displaystyle=k(\mathbf{x},\mathbf{x}^{\prime})-\mathbf{k}_{\mathbf{x}n}^{\top}(\mathbf{K}_{nn}+\sigma_{\text{noise}}^{2}I)^{-1}\mathbf{k}_{\mathbf{x}^{\prime}n}.$

The parameters $\\{\Theta,\sigma_{\text{noise}}\\}$ can be estimated by
maximizing the log marginal likelihood

$\displaystyle\log
p(\mathbf{y})=\log[\mathcal{N}(\mathbf{y}|\mathbf{0},\mathbf{K}_{nn}+\sigma^{2}_{\text{noise}}I)].$

Inspection of the posterior parameters illustrates a computational bottleneck.
Specifically, computation is via an inversion of a matrix of size $n\times n$,
where $n$ is the cardinality of the data set, which is a $\mathcal{O}(n^{3})$
operation. For large data sets, this is computationally infeasible. This
observation has motivated the development of Sparse Gaussian processes.

Sparse Gaussian processes (SGPs) [2, 3, 4, 5] are designed to overcome the
cubic nature of calculating GP posteriors by approximating the GP using
another Gaussian process supported with $m\ll n$ inducing points. This reduces
the posterior calculation to a matrix inversion using a $m\times m$ matrix;
SGPs reduce the computations to $\mathcal{O}(m^{3})$.

This paper uses the sparse variational GP (SVGP) [4, 14, 15]. This approach
implements a variational Bayesian method that approximates the true posterior
distribution with a variational distribution $q$ defined using inducing points
$\mathbf{X}_{m}$, mean parameter $\boldsymbol{\mu}$, and covariance parameter
$\mathbf{A}$ as variational parameters.

Once these parameters are calculated, they can be used to define an
approximate GP with mean parameter $m_{\mathbf{y}}^{q}(\cdot)$ and covaraince
function $k_{\mathbf{y}}^{q}(\cdot,\cdot)$ given by

$\displaystyle m_{\mathbf{y}}^{q}(\mathbf{x})$
$\displaystyle=\mathbf{k}_{\mathbf{x}m}^{\top}\mathbf{K}_{mm}^{-1}\boldsymbol{\mu}\,,$
(2) $\displaystyle k_{\mathbf{y}}^{q}(\mathbf{x},\mathbf{x}^{\prime})$
$\displaystyle=k(\mathbf{x},\mathbf{x}^{\prime})-\mathbf{k}_{\mathbf{x}m}^{\top}\mathbf{K}_{mm}^{-1}\mathbf{k}_{m\mathbf{x}^{\prime}}$
$\displaystyle\quad+\mathbf{k}_{\mathbf{x}m}^{\top}\mathbf{K}_{mm}^{-1}\mathbf{A}\mathbf{K}_{mm}^{-1}\mathbf{k}_{m\mathbf{x}^{\prime}}\,.$

Given the inducing points $\mathbf{X}_{m}$, the parameters $\boldsymbol{\mu}$
and $\mathbf{A}$ can be calcualted as

$\displaystyle\boldsymbol{\mu}$ $\displaystyle=$
$\displaystyle\sigma_{\text{noise}}^{-2}\mathbf{K}_{mm}\boldsymbol{\Sigma}\mathbf{K}_{mn}\mathbf{y},$
$\displaystyle\mathbf{A}$ $\displaystyle=$
$\displaystyle\mathbf{K}_{mm}\boldsymbol{\Sigma}\mathbf{K}_{mm},$
$\displaystyle\boldsymbol{\Sigma}$ $\displaystyle=$
$\displaystyle(\mathbf{K}_{mm}+\sigma_{\text{noise}}^{-2}\mathbf{K}_{mn}\mathbf{K}_{nm})^{-1}.$

There are many techniques to determine the inducing points $\mathbf{X}_{m}$;
See Titsias et al.[4] for more details.

We point out that the subscripts of $\mathbf{K}_{mn}$ highlights if the
inducing points or the data correspond to the first and second component of
the covariance function. For example, $\mathbf{K}_{mn}$ varies the inducing
point across the rows and varies the data points across the columns.

Bauer et al. [15] provides an in-depth analysis of the SVGP’s lower bound
which illuminates role of the inducing points. The variational approximation
maximizes the evidence lower bound (ELBO) $\mathcal{F}$ with regard to the
variational distribution:

$\displaystyle\mathcal{F}$
$\displaystyle=\underbrace{\frac{n}{2}\log(2\pi)}_{\text{constant}}+\underbrace{\frac{1}{2}\mathbf{y}^{\top}(\mathbf{Q}_{nn}+\sigma_{\text{noise}}^{2}I)^{-1}\mathbf{y}}_{\text{data
fit}}$ (3)
$\displaystyle+\underbrace{\frac{1}{2}\log|\mathbf{Q}_{nn}+\sigma_{\text{noise}}^{2}I|}_{\text{complexity}}-\underbrace{\frac{1}{2\sigma_{\text{noise}}^{2}}Tr(\mathbf{K}_{nn}-\mathbf{Q}_{nn})}_{\text{trace}}\,,$

where $\mathbf{Q}_{nn}=\mathbf{K}_{nm}\mathbf{K}_{mm}^{-1}\mathbf{K}_{mn}$ is
the Nyström approximation to $\mathbf{K}_{nn}$. The separate terms have the
following interpretations. The data fit term penalizes inaccuracies in
outputs. The complexity term penalizes local density in inducing point
selection. When the trace term approaches zero, the $m$ inducing points become
a sufficient statistic for the $n$ training samples, meaning that SGP with
inducing points $\mathbf{X}_{m}$ has the same distribution on the
$\mathbf{X}_{n}$ samples points.

## Method

Our approach replaces the support neighborhoods of Virani et al. [1] with a
SGP; thus, reducing the computational cost and memory footprint from $n$ data
points to $m$ inducing points. This greatly reduces the required computational
resources when $m\ll n$, with $m$ tunable based upon validation performance.

A review of the concepts in Virani et al. [1] motivates our SGP approach. They
first calculate the activations of each of the layers for the training data
points $\mathbf{X}_{n}$. Given a new data point, if there are enough training
point activations within an $\epsilon$-neighborhood around the new data point,
then there is local evidence to support the final classification decision. If
the number of points within the ball is defined as a decision supporting set,
then there is support-set-sufficiency to trust the output. At wrost, though,
this requires a comparison of all new operating points to all training points.

Instead of using $\epsilon$-neighborhoods to calculate sufficient evidence,
our approach includes a probability based measure for sufficient evidence.
Specifically, if the training set activation density is estimated using a
Gaussian process, then the Gaussian process’s training index set and
covariance provides an estimate of the local supporting evidence. However, the
computations required for a Gaussian process requires a large memory and
computational load. We overcome this computational burden using SGPs.

The SGP’s inducing points $\mathbf{X}_{m}$ are a set of locations that can be
used to effectively approximates the Gaussian distribution of the training
data set $\mathcal{D}$. The SGP calculates the covariance using only these
inducing points and a test sample exploits this computational efficient
covariance to predict the support-set-sufficiency.

The covariances determines the limits of the support neighborhood; it gives
example-based justification for the model’s predictions using the inducing
points. Compared with GPs, SGPs reduce the computation cost from
$\mathcal{O}(n^{3})$ to $\mathcal{O}(nm^{2}+m^{3})$; Thus, our support
neighborhood determination is significantly faster than the prior approach for
practical use. Figure 2 illustrates our approach.

Figure 2: An illustration of the proposed example-based XAI approach.

To define the specific calculations, consider $r$ observations
$\mathbf{X}_{r}$, define $a=\min(\mathbf{K}_{mm})$, $b=\max(\mathbf{K}_{mm})$,
and $\mathbf{X}=\text{concat}[\mathbf{X}_{r},\mathbf{X}_{m}]$. We estimate the
distance to each inducing point, the covariance and posterior likelihood of
generation corresponding to the embedding around that inducing point, and a
covariance-adjusted distance:

$\displaystyle\mathbf{D}_{rm}$ $\displaystyle=\text{pairwise-
distance}(\mathbf{X}_{r},\mathbf{X}_{m}),$ (4) $\displaystyle\mathbf{K}$
$\displaystyle=k_{\mathbf{y}}^{q}(\mathbf{X},\mathbf{X}),$
$\displaystyle\mathbf{K}_{rm}$
$\displaystyle=\frac{\mathbf{K}^{-1}[:r,r:]-a}{b-a},$
$\displaystyle\mathbf{P}_{rm}$
$\displaystyle=\text{clip}(\mathbf{K}_{rm},0,1)$
$\displaystyle\mathbf{D}_{rm}^{\text{cov}}$
$\displaystyle=\mathbf{D}_{rm}+\lambda\mathbf{P}_{rm}$

With $\epsilon$ as a threshold parameter, for each new sample in
$\mathbf{X}_{r}$ counting the number of inducing points such that
$\mathbf{D}_{rm}<\epsilon$ or $\mathbf{D}_{rm}^{\text{cov}}<\epsilon$ provide
a means to define $\epsilon$-neighborhoods and support sets. These
neighborhoods give local evidence to support the final classification
decision. As in Virani et al., this procedure defines epistemic operating
points based label-coherent support from the estimated local evidence. At a
threshold level, epistemic uncertainty becomes apparent in the presence of
label-conflicting inclusions. Moreover, assuming a sufficient number of
samples, for inducing points with conflicting labels the function
$\mathbf{P}_{rm}$ directly estimate label prediction uncertainty. We showcase
these uses below.

## Experiments

We compare the efficacy of our approach to Virani et al.’s performance. Using
the CIFAR-10 dataset [16], we applied Virani et al.’s method that uses all $n$
data points against our method that uses SGPs and $m<<n$ inducing points using
three metrics: final labeling accuracy; new label confidence; and computation
time. The new label confidence requires explanation.

Figure 3: Label accuracy for CIFAR-10 Figure 4: Epistemic operation for
CIFAR-10 Figure 5: Inference time for CIFAR-10

New label confidence measures the supporting evidence of a new data point’s
label and is measured using tunable
$\epsilon$-neighborhoods. For Virani et al., the $\epsilon$-neighborhood was
an $\epsilon$-ball around the new points. For our case, the
$\epsilon$-neighborhood is defined by counting the number of inducing points
with $D^{\text{cov}}_{rm}$ or $D_{rm}$ in Equation 4. For both methods, if
there are significant (thresholded) number of label-coherent points within the
respectively defined $\epsilon$-neighborhoods, then there is sufficient
training samples to support the decision and the sample point passes the ”I
Know” (IK) requirement.

The metrics are shown in three figures. The accuracy of sample points that
pass the IK requirement is in Figure 3, the fraction of points that pass the
IK requirement is in figure 4, and execution time is shown in Figure 5. In
each case, performance compared with the Virani et al. ”Baseline” is shown via
a black dashed line. The ”SGP” method uses $\mathbf{D}_{rm}$ and ”Cov-SGP”
uses $\mathbf{D}_{rm}^{\text{cov}}$, and, instead of using variational methods
to find the inducing points, ”Random Subset” uses a randomly selected set of
inducing points. The addition of a random subset highlights the utility of
including the inducing points as variational parameters in the Sparse Gaussian
Process versus less nuanced embeddings. Performance variance over 10 random
restarts is shown as a shaded region around each operating curve to help
clarify confidence in the application of the approach even with varied
instantiations. We observed that there was no significant measurable benefit
if the number of inducing points was incrased past 2500. For ease of visual
comparison, this plateau effect of epistemic inclusion thresholding is
indicated with a dashed red line.

As shown, our SGP or Cov-SGP solution provides assured label performance
slightly improving upon the baseline regardless of the number of inducing
points selected. Using the IK requirements as an epistemic selection criteria
restricts the CIFAR-10 dataset to label 79% of all validation points vs. 82%
for Virani et al. Of most importance, however, is the reduction of time
required for inference by two orders of magnitude; this goes hand-in-hand with
the reduction of data points which need to be retained from the original
50,000 by at least an order of magnitude.

## Conclusions

Based upon our findings, we believe that there is strong potential in
deploying epistemic uncertainty-aware systems for safe practice even in
extremely confined size, weight, and power environments. The efficacy of SGPs
for use in embedding training experience suggests an initial path to high-
confidence use, while we intend to further explore alternate embeddings and
other improvements in future work.

## References

  * [1] Virani, N., Iyer, N., and Yang, Z., “Justification-based reliability in machine learning,” in [Proceedings of the AAAI Conference on Artificial Intelligence ], 34(04), 6078–6085 (2020).
  * [2] Snelson, E. and Ghahramani, Z., “Sparse Gaussian Processes using Pseudo-inputs,” in [Advances in Neural Information Processing Systems ], Weiss, Y., Schölkopf, B., and Platt, J., eds., 18, MIT Press (2006).
  * [3] Bui, T. D., Yan, J., and Turner, R. E., “A Unifying Framework for Gaussian Process Pseudo-Point Approximations Using Power Expectation Propagation,” Journal of Machine Learning Research 18(104), 1–72 (2017).
  * [4] Titsias, M., “Variational Learning of Inducing Variables in Sparse Gaussian Processes,” in [Proceedings of the Twelth International Conference on Artificial Intelligence and Statistics ], van Dyk, D. and Welling, M., eds., 567–574, PMLR, Florida, USA (2009).
  * [5] Hoang, T. N., Hoang, Q. M., and Low, B. K. H., “A Unifying Framework of Anytime Sparse Gaussian Process Regression Models with Stochastic Variational Inference for Big Data,” in [Proceedings of the 32nd International Conference on Machine Learning ], Bach, F. and Blei, D., eds., 37, 569–578, PMLR, Lille, France (2015).
  * [6] Schwalbe, G. and Finzel, B., “A comprehensive taxonomy for explainable artificial intelligence: a systematic survey of surveys on methods and concepts,” Data Mining and Knowledge Discovery (2023).
  * [7] Ancona, M., Ceolini, E., Öztireli, C., and Gross, M., [Gradient-Based Attribution Methods ], 169–191, Springer International Publishing, Cham (2019).
  * [8] Neely, M., Schouten, S. F., Bleeker, M. J. R., and Lucic, A., “Order in the court: Explainable AI methods prone to disagreement,” CoRR abs/2105.03287 (2021).
  * [9] Suarez, J. L., Garcia, S., and Herrera, F., “Ordinal regression with explainable distance metric,” Machine Learning 110, 2729–2762 (2021).
  * [10] Quinonero-Candela, J., Rasmussen, C. E., and Williams, C. K. I., “Approximation Methods for Gaussian Process Regression,” in [Large-Scale Kernel Machines ], 203–223, MIT Press (2007).
  * [11] van der Wilk, M., Dutordoir, V., John, S., Artemev, A., Adam, V., and Hensman, J., “A Framework for Interdomain and Multioutput Gaussian Processes,” ArXiv (2020).
  * [12] Borovitskiy, V., Terenin, A., Mostowsky, P., and Deisenroth (he/him), M., “Matérn Gaussian processes on riemannian manifolds,” in [Advances in Neural Information Processing Systems ], Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H., eds., 33, 12426–12437, Curran Associates, Inc. (2020).
  * [13] Rasmussen, C. E. and Williams, C. K. I., [Gaussian Processes for Machine Learning ], MIT Press, Cambridge, USA (2005).
  * [14] Burt, D., Rasmussen, C. E., and Van Der Wilk, M., “Rates of Convergence for Sparse Variational Gaussian Process Regression,” in [Proceedings of the 36th International Conference on Machine Learning ], Chaudhuri, K. and Salakhutdinov, R., eds., 97, 862–871, PMLR (Jun 2019).
  * [15] Bauer, M., van der Wilk, M., and Rasmussen, C. E., “Understanding Probabilistic Sparse Gaussian Process Approximations,” in [Advances in Neural Information Processing Systems ], 1533–1541 (2016).
  * [16] Krizhevsky, A., “Learning multiple layers of features from tiny images,” (2009).