# An algebraic approach to the Kuramoto model

Lyle Muller  Ján Mináč  Tung T. Nguyen  Department of Mathematics, Western
University
London, Ontario, Canada N6A 3K7

###### Abstract

We study the Kuramoto model with attractive sine coupling. We introduce a
complex-valued matrix formulation whose argument coincides with the original
Kuramoto dynamics. We derive an exact solution for the complex-valued model,
which permits analytical insight into individual realizations of the Kuramoto
model. The existence of a complex-valued form of the Kuramoto model provides a
key demonstration that, in some cases, re-formulations of nonlinear dynamics
in higher-order number fields may provide tractable analytical approaches.

###### pacs:

05.45.Xt,02.10.Ox

The dynamics of networks with many nodes and connections poses difficulties
for mathematical treatment. While linear dynamics on a network of interacting
units is straightforward to handle with matrix algebra, once nonlinearity is
introduced at individual units, the dynamics of the system becomes difficult
to study analytically. Here, we consider the Kuramoto model (KM), a
paradigmatic example of nonlinear network dynamics describing the
synchronization of coupled oscillators Acebron05 . This model provides a
canonical description of synchronization in nature, where populations of
interacting units (from neurons to Josephson junctions and fireflies
Strogatz03 ) coordinate the timing of their behavior in the absence of a
central coordinator Strogatz01 . Because of its wide applicability throughout
many systems, the KM is a central consideration in the study of nonlinear
dynamics.

Analytical approaches to the KM started with the original work of Kuramoto,
who introduced the model in Kuramoto75 . By changing to a rotating coordinate
frame and passing to a continuum description, Kuramoto provided a description
of the population dynamics in the infinite limit and the transition to
synchrony at a critical coupling strength. Since this pivotal work, much
research has gone on to analyze the dynamics of KM. In recent years, several
studies have focused on introducing complex versions of the KM Roberts08 ;
vanMieghem09 , with the goal of providing some analytical insight; however,
finding a complex KM that permits an exact solution for the evolution of this
system has proven evasive.

In this Letter, we provide a complex-valued matrix formulation of the KM whose
argument corresponds to the original KM. We derive an explicit analytical
solution for the dynamics using the eigenspectrum of the adjacency matrix. We
then use this solution to study synchronization in individual realizations of
the model. The existence of a complex-valued version of the KM permitting an
analytical solution provides a key demonstration for potential analytical
approaches to complex nonlinear dynamics at finite scales. The complex-valued
KM and its exact solution provide an example that, in some cases, re-
formulations of nonlinear dynamics in different number fields may provide
opportunities for analytical descriptions of nonlinear systems.

We start with the standard Kuramoto model (KM) on a general network of $N$
nodes:

$\dot{\theta}_{i}=\omega_{i}+\kappa\sum_{j=1}^{N}a_{ij}\sin(\theta_{j}-\theta_{i})\,,$
(1)

where $\theta_{i}(t)\in[-\pi,\pi]$ is the state of oscillator $i\in[1,N]$ at
time $t$, $\omega_{i}$ is the intrinsic angular frequency, $\kappa$ scales the
coupling strength, and element $a_{ij}\in\\{0,1\\}$ represents the connection
between oscillators $i$ and $j$. The standard sine coupling in the interaction
term causes phases of two connected oscillators $i$ and $j$ to attract, to an
extent depending on the homogeneous coupling strength $\kappa$. Here, we first
consider the KM defined on an undirected ring graph $\mathfrak{G}_{RG}$, where
nodes are arranged on a one-dimensional ring with periodic boundary conditions
and connected to their $k$ neighbors in each direction (later considering
Erdős-Rényi and Watts-Strogatz graphs). We also consider first the case of
homogeneous intrinsic frequencies $\omega_{i}=\omega\,\forall\,i\in[1,N]$
(with inhomogeneous case in Supplement Supplement ). For $k=\lfloor
N/2\rfloor$, $\mathfrak{G}_{RG}$ corresponds to the fully connected graph on
$N$ nodes, $K_{N}$. It is important to note that the approach described here
is general to the specific graph model we consider.

Figure 1:  Synchronization dynamics in a three-node Kuramoto model (KM). In
this example, $N=3$, $\kappa=1$, $\omega/2\pi=10\,\text{Hz}$, and the network
is fully connected (top left), so the adjacency matrix $\bm{A}$ is the
complete graph $K_{3}$ (top middle). Compared are the numerical result (dotted
gray line) and the analytical expression (solid black line) for the state
variables $\theta_{1}$ (top), $\theta_{2}$ (middle), and $\theta_{3}$
(bottom). Plotted is one second of simulation ($dt=0.001\,\text{s}$). The
difference between the numerical result and analytical prediction was verified
to remain less than $\pi/16$ in simulations as long as 10 seconds.

Starting with Equation (1), we can change to a rotating coordinate frame
Strogatz88  and set $\omega=0$ without loss of generality. We then subtract an
additional imaginary component in the interaction term:

$\dot{\theta}_{i}=\gamma\sum_{j=1}^{N}a_{ij}\big{[}\sin(\theta_{j}-\theta_{i})-\mathrm{i}\cos(\theta_{j}-\theta_{i})\big{]}\,.$
(2)

We note this expression now implies $\theta_{i}\in\mathbb{C}$ and requires a
scaling of the coupling strength ($\gamma=2\kappa/\pi$); as we will show, the
resulting complex-valued system admits a solution whose argument agrees with
the original KM.

Multiplying both sides of Equation (2) by $\mathrm{i}$, we have

$\mathrm{i}\dot{\theta}_{i}=\gamma\sum_{j=1}^{N}a_{ij}\big{[}\mathrm{i}\sin(\theta_{j}-\theta_{i})+\cos(\theta_{j}-\theta_{i})\big{]}.$
(3)

Applying Euler’s formula, we can rewrite the system as:

$\mathrm{i}\dot{\theta}_{i}=\gamma\sum_{j=1}^{N}a_{ij}e^{\mathrm{i}(\theta_{j}-\theta_{i})}=\gamma
e^{-\mathrm{i}\theta_{i}}\sum_{j=1}^{N}a_{ij}e^{\mathrm{i}\theta_{j}}.$ (4)

This equation now results in the following matrix form:

$\bm{\dot{\theta}}=\frac{\gamma}{\mathrm{i}}\,\text{diag}[e^{-\mathrm{i}\bm{\theta}}]\,\bm{A}e^{\mathrm{i}\bm{\theta}}\,,$
(5)

where $\bm{\theta}\in\mathbb{R}^{N}$ is the state vector across nodes, and
$\bm{A}\in\\{0,1\\}^{N\times N}$ is the adjacency matrix encoding the
connections between nodes in the network. We can then utilize the fact that
the inverse of the matrix $\text{diag}[e^{-\mathrm{i}\bm{\theta}}]$ is
$\text{diag}[e^{\mathrm{i}\bm{\theta}}]$ to arrive at

$\text{diag}[e^{\mathrm{i}\bm{\theta}}]\bm{\dot{\theta}}=\frac{\gamma}{\mathrm{i}}\bm{A}e^{\mathrm{i}\bm{\theta}}\,.$
(6)

Now, using the fact that

$\text{diag}[e^{\mathrm{i}\bm{\theta}}]\bm{\dot{\theta}}=\frac{1}{\mathrm{i}}\frac{d}{dt}e^{\mathrm{i}\bm{\theta}}\,,$
(7)

mentioned in vanMieghem09 , we have that

$\frac{d}{dt}e^{\mathrm{i}\bm{\theta}}=\gamma\bm{A}e^{\mathrm{i}\bm{\theta}}$
(8)

and letting $\bm{x}=e^{\mathrm{i}\bm{\theta}}$, we have

$\bm{\dot{x}}=\gamma\bm{A}\bm{x}$ (9)

whose general solution is

$\bm{x}(t)=e^{\gamma t\bm{A}}\bm{x}(0)\,.$ (10)

We can now utilize a diagonalization of the adjacency matrix,
$\bm{A}=VDV^{-1}$, to rewrite the solution as

$\bm{x}(t)=\bm{V}e^{\gamma t\bm{D}}\bm{V^{-1}}\bm{x}(0)\,,$ (11)

and the problem of computing the matrix exponential is now reduced to
computing the standard exponential function on the diagonal entries of
$\bm{D}$. We note that the adjacency matrix $\bm{A}$ is always diagonalizable
for the undirected graphs considered here.

Finally, as noted above, Equation (10) can be directly related to the original
KM. Let
$\bm{\theta}=\bm{\theta}_{{\text{re}}}+\mathrm{i}\bm{\theta}_{{\text{im}}}$ be
the decomposition of $\bm{\theta}$ into the real and imaginary parts. Then we
have

$\bm{x}=e^{\mathrm{i}\bm{\theta}_{{\text{re}}}-\bm{\theta}_{{\text{im}}}}=e^{-\bm{\theta_{{\text{im}}}}}e^{\mathrm{i}\bm{\theta}_{{\text{re}}}}$
(12)

$\bm{\theta_{{\text{re}}}}$ is thus the argument of the analytical solution
$\bm{x}$. In particular, we can take $\bm{\theta_{{\text{re}}}}\in[-\pi,\pi]$.
We will use this solution to compare with the numerical integration of
Equation (1).

Figure 2: Numerical and analytical realizations of the KM. Plotted are the
state variables ($\theta_{i}\,\forall\,i\in[1,N]$) produced by the model
($N=200$, $k=100$, resulting in the fully connected complete graph $K_{N}$),
with different nodes on the horizontal axis and time on the vertical. Dark
colors indicate values close to $-\pi$, and light colors indicate values close
to $\pi$.

What remains in deriving an expression for the resulting temporal dynamics is
to study the eigenspectrum of $\bm{A}$. To do this, we note the adjacency
matrix of the ring graph $\mathfrak{G}_{RG}$ is by definition a circulant
matrix Davis79 , and we can use the circulant diagonalization theorem (CDT) to
obtain its eigenspectrum analytically. The CDT states that all circulants
$c_{ij}=\text{circ}(c_{j})$, where $\text{circ}(c_{j})$ is the circulant
matrix constructed from the generating vector $c_{j}$, are diagonalized by the
same unitary matrix $\bm{U}$ with components

$u_{rs}=\frac{1}{\sqrt{N}}\exp\left[-\frac{2\pi\mathrm{i}}{N}(r-1)(s-1)\right],$
(13)

$r,s\in[1,N]$ and that the $N$ eigenvalues are given by

$E_{r}(\mathbf{C})=\sum\limits_{j=1}^{N}c_{j}\exp\left[-\frac{2\pi\mathrm{i}}{N}(r-1)(j-1)\right]\,.$
(14)

Using these expressions, we can then evaluate Equation (10) in terms of the
eigenspectrum of $\bm{A}$ to obtain a fully analytical evaluation of
$\bm{x}(t)$, which we can then compare to numerical integration of the
original KM.

We can now compare the argument of the analytical expression Equation (10)
with the result obtained by numerical integration of Equation (1). As a first
example, we considered a small network of $N=3$ nodes and $k=1$, resulting in
the complete graph $K_{3}$. Figure 1 (top left) shows the considered network
along with the state variables $\theta_{i}$ at each node. Initial conditions
$\theta_{i}(0)$ were selected randomly with uniform distribution
$\mathcal{U}(-\pi,\pi)$. Starting from these random initial conditions and a
homogeneous coupling strength $\kappa=1$, this small network synchronizes over
the course of the one second simulation (cf. dotted blue vertical lines across
bottom three panels, Fig. 1). Equations were integrated using a forward Euler
method with fine temporal precision ($dt=0.001\,\text{s}$) and were compared
with results from a Runge-Kutta method Dormand80 . Additional simulations were
run at very high temporal precision (up to $dt=10^{-6}\,\text{s}$) to ensure
the numerical validity of these results. Simulation code to reproduce all
figures in this work will be made availabe at our GitHub site
(http://mullerlab.github.io).

We next compare the analytical and numerical results for larger networks of
$N=200$ nodes on a ring graph with $k=100$. To create a synchronization
transition during a one second simulation interval, we scaled the coupling
strengths accordingly ($\kappa=6/N$). Figure 2 shows the results of this
comparison between the numerical (top) and analytical (bottom) evaluations.
The numerical and analytical evaluations exhibit similar dynamics, from the
macroscopic synchronization to the specific trajectories of individual
oscillators. Interestingly, a small fraction of the nodes in the numerical
simulation remain counterphase to the rest of the population, while this
behavior does not occur in the analytical version.

Figure 3: Synchronization in numerical and analytical realizations of the KM
($N=200$, $k=100$). The modulus of the synchronziation order parameter $r$ is
plotted as a function of the coupling strength $\kappa$ for the numerical
(black line) and analytical (red line) evaluations of the model. Solid lines
represent the average $|r|$ over 10 1-second simulations (with the same random
initial conditions for the numerical and analytical evaluation in each
realization) at a value of $\kappa$ on a logarithmic space of 1,000 points
between $10^{-3}$ and $10^{1}$. Shaded areas represent standard deviation over
realizations.

To understand this point further, we systematically studied synchronization in
the KM. To quantify the extent of synchronization, we use a standard approach
to study the sum over oscillators:

$r(t)=\frac{1}{N}\sum_{j=1}^{N}e^{\mathrm{i}\theta_{j}(t)}$ (15)

where $r(t)\in\mathbb{C}$, $\lvert r(t)\rvert\in[0,1]$ is the order parameter
at time $t$, and $\theta_{j}(t)\in[1,N]$ is the phase of oscillator $j$ at
time $t$ (n.b. $\theta_{j,{\text{re}}}(t)$ in the case of the analytical
expression Equation 10). Figure 3 shows the time-average order parameter over
1-second simulations. As $\kappa$ increases in the KM, the order parameter
begins at a low value (desynchronized state) and increases until approaching
unity (synchronized state) at values of $\kappa$ ranging from 1 to 10. As
observed in Figure 2, the numerical and analytical versions of the KM exhibit
very similar macroscopic synchronization dynamics.

Lastly, we considered the numerical and analytical solutions of the KM on
undirected Erdős-Rényi and Watts-Strogatz random graphs. For each realization
of a random graph, we obtain a numerical estimate of the eigenspectrum of its
adjacency matrix $\bm{A}$ to use in the analytical expression Equation (10).
We first considered the KM on an Erdős-Rényi random graph
($\mathfrak{G}_{ER}$), which displays synchronization dynamics similar to
those previously observed (Fig. 4, top). We then considered the KM on a Watts-
Strogatz network ($\mathfrak{G}_{WS}$), which is defined as a ring graph
$\mathfrak{G}_{RG}$ where each node is first connected to its $k$ neighbors in
each direction and each edge is rewired to another node with uniform
probability $q$ WattsStrogatz98 . The KM on $\mathfrak{G}_{WS}$ displays non-
trivial spatiotemporal dynamics before converging to the synchronized state
(Fig. 4, bottom middle). Importantly, these spatiotemporal dynamics were also
well described by the analytical expression introduced here (Fig. 4, bottom
right).

In this work, we have introduced a complex-valued formulation of the KM whose
argument corresponds to the original Kuramoto dynamics. This formulation
permits an exact analytical solution for individual realizations of the KM.
Here, we have first considered the case of homogeneous intrinsic frequency;
however, this approach generalizes to the inhomogeneous case (Supplement
Supplement ). We then compared the analytical version to numerical integration
of the KM equations and found the analytical version displays similar
dynamics, from the macroscopic synchronization behavior to the specific
trajectories of individual oscillators.

Previous research, including an inspiring technical report written by van
Mieghem vanMieghem09 , has studied complex-valued formulations of the KM
Roberts08 ; Cadieu10 ; however, no exact analytical expression was previously
obtained. In particular, the expression derived in vanMieghem09  was noted to
hold only for the repulsive cosine variant of the KM, and the linear
reformulation in Roberts08  requires tuning a parameter to create the
correspondence to the original KM. Thus, the results reported in this Letter,
while representing only an initial study utilizing the analytical expression
in Equation (10), to the best of our knowledge represent the first analytical
version of the Kuramoto model.

Figure 4: Numerical and analytical realizations on the Erdős-Rényi and Watts-
Strogatz random graphs. (top left) The Erdős-Rényi random graph
$\mathfrak{G}_{ER}$ is plotted with nodes (black dots) and edges (blue lines)
($N=200,p=0.2$). (top right) The numerical integration of the KM (Eq. 1) and
the analytical evaluation of Eq. (10) are plotted ($\kappa=50/N$). As
previously, dark colors indicate values close to $-\pi$, and light colors
indicate values close to $\pi$. (bottom left) The Watts-Strogatz (WS) graph
$\mathfrak{G}_{WS}$ ($N=200$, $k=10$, and rewiring probability $q=0.1$).
(bottom right) Numerical and analytical evaluations of the KM are plotted as
above.

Importantly, we emphasize that the analytical version introduced here is valid
at finite scales and for individual realizations of the KM. This analytical
version allows future mathematical study of synchronization dynamics in
networks with many nodes and connections, potentially using new tools from
spectral graph theory, in addition to allowing one to obtain the future state
of the system at an arbitrary future moment without numerical integration of
the differential equations defined on the network. Because the KM has been
extensively studied both as a model for neural dynamics Strogatz88 ; Mirollo90
; Abbott93  and as a fundamental model for neural computation Kleinfeld01 ,
these results open up several possibilities for studying the connections
between network structure, nonlinear dynamics, and computation. Recurrent
connections have previously been shown to produce powerful computations
through nonlinear interactions Sinha98 . The approach introduced here may have
applications in understanding such recurrent interactions, which have been
increasingly acknowledged to play an important and unexplained role in visual
processing in the brain Kietzmann19 . Understanding more clearly the
connection between networks and computation thus may have implications for
fields such as neuroscience and beyond.

In this work, we have specifically studied the KM defined on a ring graph
$\mathfrak{G}_{RG}$, whose highly regular structure permits analytical study
of the eigenspectrum of its adjacency matrix. The regular structure of this
graph means that its adjacency matrix belongs to the class of circulant
matrices, whose eigenspectrum can be calculated analytically. Further, when
$k$ is maximal, such that the number of neighbors to which a node is connected
equals the rest of the graph, $\mathfrak{G}_{RG}$ corresponds to $K_{N}$, the
complete graph on $N$ nodes. The KM defined on $K_{N}$, in turn, corresponds
to the case of all-to-all connectivity first considered by Kuramoto Kuramoto75
. For these reasons, we chose to focus on the KM defined on
$\mathfrak{G}_{RG}$ in this work. In previous work RudolphMuller14 , however,
we have introduced an operator-based approach to the structure of random
graphs. In future work, we aim to extend the present results to understand the
connection between nonlinear dynamics and random graphs with various
structural features.

The existence of a complex-valued formulation of the Kuramoto model that
permits an analytical solution raises an important example in nonlinear
dynamics. While Equation (1), which includes the sine coupling interaction
term and the network adjacency matrix, appears analytically intractable, in
this case a re-formulation in the complex domain provides an algebraic
approach to the Kuramoto dynamics which serve as a canonical description for
synchronization phenomena in nature. While this re-formulation may not
generalize beyond the Kuramoto model, the ability of this approach to provide
insight in this case suggests that representations in number fields beyond
$\mathbb{R}$ may represent opportunities for future insight.

###### Acknowledgements.

The authors wish to thank M Ly, E Bienenstock, and TP Coleman for insightful
discussions. L.M. was supported by BrainsCAN at Western University through the
Canada First Research Excellence Fund (CFREF) and by the NSF through a
NeuroNex award (#2015276). J.M. was supported by the Natural Sciences and
Engineering Research Council of Canada (NSERC) grant R0370A01.

## References

  * (1) J.A. Acebrón, L.L. Bonilla, C.J. Vicente, F. Ritort, R. Spigler, Rev Mod Phys 77, 2005.
  * (2) S.H. Strogatz, Sync: The Emerging Science of Spontaneous Order, Hyperion, 2003.
  * (3) S.H. Strogatz, Nature 410, 2001.
  * (4) Y. Kuramoto, Int Symp Math Theor Phys 30, 420, 1975.
  * (5) D.C. Roberts, Phys Rev E 77, 2008.
  * (6) P. van Mieghem, Delft U Technical Report, 2009.
  * (7) See Supplemental Material at [URL will be inserted by publisher] for additional calculations on the case of inhomogeneous intrinsic frequencies.
  * (8) S.H. Strogatz and R.E. Mirollo, J Phys A 21, L699-L705, 1988.
  * (9) P.J. Davis, Circulant Matrices, Wiley, 1979.
  * (10) J.R. Dormand and P.J. Prince, J. Comp. Appl Math, 6, 1980.
  * (11) D.J. Watts, S.H. Strogatz, Nature 393, 440 (1998).
  * (12) C.F. Cadieu and K. Koepsell, Neural Comp 22, 3107-3126, 2010.
  * (13) R.E. Mirollo, S.H. Strogatz, SIAM J Appl Math, 50, 1990.
  * (14) L.F. Abbott, C. van Vreeswijk, Phys Rev E, 48, 1993.
  * (15) G.B. Ermentrout, D. Kleinfeld, Neuron, 29, 2001.
  * (16) S. Sinha, W.L. Ditto, Phys Rev Lett, 81, 1998.
  * (17) T.C. Kietzmann, et al., PNAS, 116, 2009.
  * (18) M. Rudolph-Lilith and L. Muller, Phys Rev E 89, 012812, 2014.


# QUCE: The Minimisation and Quantification of Path-Based Uncertainty for
Generative Counterfactual Explanations

Jamie Duell  Hsuan Fu  Monika Seisenberger  Xiuyi Fan Corresponding Author.
Email<EMAIL_ADDRESS>School of Mathematics and Computer Science, Swansea
University  Department of Finance, Insurance and Real Estate, Université Laval
Lee Kong Chian School of Medicine, Nanyang Technological University

###### Abstract

Deep Neural Networks (DNNs) stand out as one of the most prominent approaches
within the Machine Learning (ML) domain. The efficacy of DNNs has surged
alongside recent increases in computational capacity, allowing these
approaches to scale to significant complexities for addressing predictive
challenges in big data. However, as the complexity of DNN models rises,
interpretability diminishes. In response to this challenge, explainable models
such as Adversarial Gradient Integration (AGI) leverage path-based gradients
provided by DNNs to elucidate their decisions. Yet the performance of path-
based explainers can be compromised when gradients exhibit irregularities
during out-of-distribution path traversal. In this context, we introduce
Quantified Uncertainty Counterfactual Explanations (QUCE), a method designed
to mitigate out-of-distribution traversal by minimizing path uncertainty. QUCE
not only quantifies uncertainty when presenting explanations but also
generates more certain counterfactual examples. We showcase the performance of
the QUCE method by comparing it with competing methods for both path-based
explanations and generative counterfactual examples. The code repository for
the QUCE method is available at: https://github.com/jamie-duell/QUCE.

123

##  1 Introduction

Given the prevelance of big data and increased computability, the application
of Deep Neural Network (DNN) methods are a commonality. However, the
intricacies and depth of DNN architectures lead to results that lack inherent
interpretability. In pivotal domains such as healthcare and finance,
interpretability is crucial and thus the application of eXplainable Artificial
Intelligence (XAI) to extract valuable insights from the DNN models is
widespread [6, 8].

The Path-Integrated Gradients (Path-IG) [30] formulation presents axiomatic
properties that are upheld solely by path-based explanation methods. The Out-
of-Distribution (OoD) problem is prevalent in the application of path-based
explanation methods [11]; here the intuition is that traveling along a
straight line path can incur irregular gradients and thus provide noisy
attribution values [17]. Another known limitation of many Integrated Gradient
(IG) [30] based approaches is the selection of a baseline reference; thus the
Adversarial Gradient Integration (AGI) [24] method relaxes this constraint by
generating baselines in adversarial classes. We note that AGI utilizes the
path-based approach for generating counterfactual examples, and for this
reason will be a primary baseline for our proposed method throughout this
paper.

Counterfactual explanations [13] are often presented in the form of
counterfactual examples [32, 23]; here the goal is to provide a counterfactual
example belonging to an alternative class with respect to a reference example.
Counterfactual approaches aim to answer the question:

> “Given an instance, what changes can be made to change the outcome for that
> instance?"

Naturally, this allows for empirical observation as to which changes could
provide an alternative outcome. The argument for using counterfactual methods
is often developed from a causal lens [15, 25]. It follows that to better
evaluate this causal relationship, a promising avenue is to unify feature
attribution with counterfactual examples, as demonstrated by the Diverse
Counterfactual Explanations (DiCE) [23] method. Naturally, given quantitative
approaches to feature attribution calculation such as these, ideally feature
attribution methods should adhere to desirable axioms across XAI literature
[30, 2]. Thus, we aim to utilize state-of-the-art feature attribution
assignment as to satisfy key axioms in our model development. Another concern
with production of counterfactual examples is the production of realistic
paths to successfully create a counterfactual example; therefore we shall be
exploring uncertainty.

Uncertainty quantification is not often considered when producing
explanations, although some approaches have explored this. Examples include
[29] where post-hoc model-agnostic approaches such as Local Interpretable
Model-Agnostic Explanations (LIME) [26] and kernel SHapley Additive
exPlanations (SHAP) [20] are adapted into a Bayesian framework to model the
uncertainty of the explanations produced. Since path-based methods are
implementation invariant with respect to the model, the explanations will be
consistent and thus there will be no variance in the explanations produced. In
this way uncertainty quantification in the form of repeated runs of the XAI
algorithm as elucidated in [22], while applicable to post-hoc approximation
XAI methods, will not suffice for implementation-invariant models.
Autoencoder-based frameworks have also been used to measure uncertainty for
both machine learning predictions and explanations [4], with integration of
uncertainty seen in the production of counterfactual examples [21, 14].
Inherently, autoencoder approaches provide a more suitable basis for
attribution settings by evaluating uncertainty in explanations with respect to
the uncertainty inherent in the data. The standard autoencoder approach
evaluates the reconstruction error, which is often utilized in work
surrounding anomaly detection [31, 3]; instead, we explore the use of a
variational autoencoder (VAE) for variational inference, and thus investigate
counterfactuals generated with respect to our approximation of the true data
distribution.

To address the above constraints, we propose the Quantified Uncertainty (Path-
Based) Counterfactual Explanations (QUCE) method. The focus of the proposed
method is three-fold. We aim to

  * •

minimize uncertainty and thus maximize the extent to which the generated paths
and counterfactual examples are within distribution;

  * •

relax the straight-line path constraints of Integrated Gradients;

  * •

provide uncertainty quantification for counterfactual paths and counterfactual
feature attribution.

In this work, we focus on the minimization of uncertain paths for
counterfactual generation with quantifiable uncertainty measures on the
generated counterfactual. QUCE’s learning process relaxes IG’s straight-line
path restrictions as part of the generative process.

Intuitively, it is unclear in many scenarios if one single best path toward an
alternative outcome exists; for example a patient’s treatment path may be
unclear [5], or there may be many viable paths to achieve the same outcome
[7]. Therefore, QUCE utilises both a single and multiple-paths approaches, so
we can observe a generalized explanation over all paths for an instance in
obtaining a desired class and likewise inspect many example paths. From the
multiple-paths approach, we are able to gauge a general approximation over
many piecewise linear paths of the most important features in obtaining a
desired counterfactual outcome, each path independently aiming to minimize the
uncertainty in its generative process and thus provide greater in-distribution
interpolation. Similarly, we present the optimisation over the key metrics –
_proximity_ [13], _validity_ [13] and _uncertainty_ [27]. We provide an
overview of some generative counterfactual methods and evaluate which account
for the aforementioned metrics in table 2.

Table 1: An overview of generative counterfactual methods and their consideration of key metrics. Here we observe of the three metrics QUCE is the method that accounts for all three. Properties | Proximity | Validity | Uncertainty
---|---|---|---
QUCE | ✓ | ✓ | ✓
DiCE | ✓ | ✓ | ✗
AGI | ✓ | ✓ | ✗

##  2 Explainable AI and Counterfactuals

Counterfactual explanations can be presented both in the form of
counterfactual examples [23, 32] and also via counterfactual feature
attribution [23].

######  Definition 1 (Counterfactual Example).

Given a probabilistic classifier $f:\mathbbm{R}^{J}\rightarrow\\{0,1\\}$,
differentiable probabilistic function $F:\mathbbm{R}^{J}\rightarrow[0,1]$, and
class $\tau=0\oplus 1$111The $\oplus$ symbol is an operator following the
standard definition of the logical XOR. an instance $\mathbf{x}=\langle
x^{1},\ldots,x^{J}\rangle\in X$ where $X\in\mathbbm{R}^{N\times J}$ and a
classification threshold $\vartheta\in[0,1]$ such that

$\displaystyle f(\mathbf{x})=\begin{cases}\neg\tau,\mbox{if
}F(\mathbf{x})\leq\vartheta\\\ \tau,\mbox{otherwise.}\end{cases}$ (1)

Then a counterfactual example of $\mathbf{x}$ is some $\mathbf{x}^{c}$ where
$f(\mathbf{x}^{c})\neq f(\mathbf{x})$.

Counterfactual examples are often produced through the use of a learned
generative function, examples of this can be seen in the work of [21].
Extending this, we define a _counterfactual generator_.

######  Definition 2 (Counterfactual Generator).

Given an instance $\mathbf{x}\in X$ and a classifier $f$, a counterfactual
generator is a function
$\mathcal{G}:\mathbbm{R}^{J}\rightarrow\mathbbm{R}^{J}$ that takes an instance
$\mathbf{x}$ and returns a counterfactual example $\mathbf{x}^{c}\notin X$.

Feature attribution is a common form of XAI, and the feature attribution
approach is seen in many methods [20, 18, 19, 9]. We continue by defining
_feature attribution_ :

######  Definition 3 (Feature Attribution).

Given an instance $\mathbf{x}$, a feature attribution method is a function
$\Phi:\mathbbm{R}^{J}\rightarrow\mathbbm{R}^{J}$, where $\Phi$ takes an
instance $\mathbf{x}$ and returns a vector of feature attribution values
$\langle\phi^{1},\ldots,\phi^{J}\rangle$.

Concatenating prior feature attribution methods and counterfactual examples
leads to _counterfactual feature attribution_ (CFA) which we define as
follows:

######  Definition 4 (Counterfactual Feature Attribution).

A counterfactual feature attribution method is a function
$\Phi_{CF}:\mathbbm{R}^{J}\times\mathbbm{R}^{J}\rightarrow\mathbbm{R}^{J}$
that takes an instance $\mathbf{x}$ and a counterfactual example
$\mathbf{x}^{c}$ and returns a vector of feature attribution values
$\langle\phi^{1},\ldots,\phi^{J}\rangle$.

##  3 Axioms for Path-Based Explainers

Here we informally introduce axioms used in XAI literature. In the seminal
work of [30], the authors introduce a set of axioms which play a foundational
role for the development of path-based feature attribution methods.
Informally, these encompass:

  * •

Completeness: The difference in prediction between the baseline and input
should be equal to the sum of feature attribution values.

  * •

Sensitivity(a): For every input and baseline that differ in one feature and
for which the subsequent prediction is different, feature attribution should
only be given to that one feature.

  * •

Sensitivity(b): If the neural network is not mathematically dependent on one
feature, the feature attribution assigned to that feature should be 0.

  * •

Implementation Invariance: Two functionally equivalent neural networks should
produce the same feature attribution as an explanation.

The above axiomatic guarantees are for path-based explainers and for this
reason we utilise the path-based formulation in our work.

##  4 Proposed Model: QUCE

###  4.1 Generating Counterfactuals with QUCE

To generate counterfactuals we propose a three-part objective function with a
composite weighting vector
$\Lambda=\langle\lambda_{1},\lambda_{2},\lambda_{3}\rangle$, where each
$\lambda\in\Lambda$ is an independent tolerance (weight) used to determine the
influence of each part of the joint objective function presented in equation
2. By minimizing the objective function we obtain the generated counterfactual
$\mathbf{x}^{c}$. Informally, our objective function is composed of three
parts:

  * •

$\mathcal{L}_{pr}$, for the maximization the probability towards the desired
class;

  * •

$\mathcal{L}_{\delta}$, to minimize the distance between the instance and a
generated counterfactual;

  * •

$\mathcal{L}_{\epsilon}$, to minimize the uncertainty of both the generated
paths and generated counterfactual examples.

Combining these terms with our weighting vector, we have

$\displaystyle\mathcal{G}(\mathbf{x})$
$\displaystyle=\underset{\mathbf{x}^{c}}{\arg\min}\
\Lambda\odot(\mathcal{L}_{pr}+\mathcal{L}_{\delta}+\mathcal{L}_{\epsilon}).$
(2)

Having constructed our generative objective function, we provide further
notation to illustrate the learning process. First, we consider an iterative
learning process such as gradient descent on $\mathbf{x}$ to produce a path
from $\mathbf{x}$ to a generated $\mathbf{x}^{c}$. Thus we aim to minimize the
developed function $\mathcal{G}(\mathbf{x})$ through the gradient descent
approach (variants include, e.g., SGD and ADAM). We initially let
$\mathbf{x}^{c}=\mathbf{x}$; $\mathbf{x}^{c}$ is updated via

$\displaystyle\mathbf{x}^{c}\Leftarrow\mathbf{x}_{\Delta_{i}},$
$\displaystyle\Delta_{i}=\varphi\nabla_{\mathbf{x}^{c}}(\mathcal{G}(\mathbf{x})),$
$\displaystyle\mathbf{x}_{\Delta_{i}}=\mathbf{x}-\Delta_{i}.$

Let $\mathbf{x}$ be updated on a loop over $i(0\leq i\leq n)$ iterations; when
$i=n$ we let have our $\mathbf{x}^{c}$ indicating our generated point. Here
$\varphi$ represents the “learning rate," a small positive multiplier value
$\varphi\in[0,1]:\varphi<<1$. We store each update on $\mathbf{x}^{c}$ as a
vector $\mathbf{x}^{\Delta}$ =
$\langle\mathbf{x}_{\Delta_{0}},\ldots,\mathbf{x}_{\Delta_{n}}\rangle$.

###  4.2 Finding Counterfactuals

####  4.2.1 Valid Counterfactuals

A key concept in finding counterfactual examples is ensuring that the
counterfactual is indeed _valid_ , and thus we aim to produce counterfactual
examples that belong to a counterfactual class.

######  Remark 1.

For the sake of simplicity we let $f(\mathbf{x})=\neg\tau$ throughout the
remainder of this work; thus $\tau$ is the counterfactual class for which we
aim to generate some $\mathbf{x}^{c}$ such that $f(\mathbf{x}^{c})=\tau$.

In light of remark 1, given a target counterfactual class $\tau$ and
probabilistic decision threshold $\vartheta$ we aim to find an instance
$\mathbf{x}^{c}$ that satisfies

$\displaystyle F(\tau|\mathbf{x}^{c})\geq\vartheta.$ (3)

Here $\vartheta$ is a probability threshold for the target class $\tau$. Thus,
we need a generator $\mathcal{G}$ that satisfies the condition in equation 3.
To achieve this, we can maximize the likelihood of an instance belonging to a
class $\tau$. The maximum log likelihood criterion is defined as:

$\displaystyle\mathcal{L}_{pr}=\bigg{[}log[F(\tau|f(\mathbf{x}^{c}))]\bigg{]}$

such that, we only accept counterfactuals where $F(\cdot|\cdot)\geq\vartheta$.
For optimisation we can rewrite this as a minimization problem, instead
minimizing the negative log-likelihood:

$\displaystyle\mathcal{L}_{pr}=\bigg{[}-log[F(\tau|f(\mathbf{x}^{c}))]\bigg{]}.$

Thus we require a valid counterfactual to have $F(\cdot|\cdot)<(1-\vartheta)$.
This constitutes the constrained optimisation problem with respect to some
target class $\tau$.

####  4.2.2 Proximity for Counterfactuals

Given an instance, in the production of counterfactual examples we often aim
to find a counterfactual example that is “similar" in feature space to the
instance. This is often termed _proximity_.

In this work, we use the $l_{2}$ norm as the proposed model focuses on
producing counterfactuals from continuous features, defining proximity as
follows:

######  Definition 5 (Proximity).

Given an instance $\mathbf{x}$ and its counterfactual example
$\mathbf{x}^{c}$, the proximity between the two instances is given by

$\displaystyle\mathcal{L}_{\delta}=\bigg{[}\frac{1}{2}||\mathbf{x}^{c}-\mathbf{x}||^{2}\bigg{]}.$
(4)

####  4.2.3 Minimally Uncertain Counterfactuals

To maximize the certainty of counterfactual examples, we examine their
complement—namely, the uncertainty associated with a counterfactual example.
To explore this, we establish the concept of _counterfactual uncertainty_.
Informally, we consider uncertainty to be the Evidence Lower Bound (ELBO) as
measured by a Variational Autoencoder (VAE) framework. The objective function
ELBO is comprised of two components, namely the Kullbeck–Leibler (KL)
divergence and a reconstruction loss. This is defined as the following:

$\displaystyle\operatorname{VAELoss}(\mathbf{x})=\mathbbm{E}_{q_{\theta}}[log\
q_{\theta}(\mathbf{z}|\mathbf{x})-log\
p_{\psi}(\mathbf{z})]-\mathbbm{E}_{q_{\theta}}\ log\
p_{\psi}(\mathbf{x}|\mathbf{z})$

where $p$ and $q$ are probability distributions and $\mathbf{z}$ is the latent
representation of $\mathbf{x}$. The aim is to find a $\theta$ that
successfully models the true training data distribution $\psi$, and thereby
satisfying the following minimization problem:

$\displaystyle\\{\theta^{*},\psi^{*}\\}=\underset{\theta,\psi}{\arg\min}\operatorname{VAELoss}(\mathbf{x})$
(5)

Intuitively, ELBO aims to produce a $p$ and $q$ that are as close as possible
while simultaneously optimizing the mapping of $\mathbf{z}$ back to its
original representation $\mathbf{x}$, through learning to model the
distributions with the parameters $\theta$ and $\psi$. Thus it suffices to say
that if the reconstructing loss is minimized, the decoded $\mathbf{z}$ should
successfully map back to $\mathbf{x}$. Posterior to the pre-trained VAE, we
now have a fixed representation shaping our $p$ and $q$ distributions with our
parameters $\theta^{*}$ and $\psi^{*}$ and thus we can provide counterfactual
uncertainty as:

######  Definition 6 (Counterfactual Uncertainty).

Following equation 4.2.3, we simply take the loss without minimization to be
the counterfactual uncertainty, thus given our fixed parameters $\theta^{*}$
and $\psi^{*}$, we have

$\displaystyle\mathcal{L}_{\epsilon}=$ (6)
$\displaystyle\mathbbm{E}_{q_{\theta^{*}}}[log\
{q_{\theta^{*}}}(\mathbf{z}|\mathbf{x}^{c})-log\
p_{\psi^{*}}(\mathbf{z})]-\mathbbm{E}_{q_{\theta^{*}}}\ log\
p_{\psi^{*}}(\mathbf{x}^{c}|\mathbf{z}).$

###  4.3 Uncertainty in Counterfactual Explanations

Definition 6 allows for the evaluation of new generated instances with a
measure of how “good" the fit of the new instance is with respect to the
training data distribution, and similarly how well a path fits into the data
distribution.

From definition 6 we have a quantifiable measure of uncertainty for the
generated counterfactual $\mathbf{x}^{c}$. Expressing this in vector form as a
difference, which is needed for modifying and updating our generated
$\mathbf{x}^{c}$ to adjust for underlying uncertainty, we define _Feature-wise
Counterfactual Uncertainty_ by simply looking at the reconstruction error,
defined for simplicity as follows:

######  Definition 7 (Feature-wise Counterfactual Uncertainty).

Given counterfactual uncertainty, we can rewrite a reconstruction error
equivalent at a feature level by calculating a vector of differences
$\epsilon_{\mathbf{d}}$ such that

$\displaystyle\epsilon_{\mathbf{d}}=\langle|d^{1}|,\ldots,|d^{J}|\rangle\ :\
\langle d^{1},\ldots,d^{J}\rangle=\mathbf{d};$ (7) $\displaystyle\mbox{where
}\ \mathbf{d}=\mathbf{x}^{c}-\hat{\mathbf{x}}^{c}.$ (8)

With this representation, we can then successfully update $\mathbf{x}^{c}$ by
both adding and subtracting this vector of feature-wise counterfactual
uncertainty as given by the reconstruction error, and thus we can calculate
_Counterfactual Explanation Uncertainty_.

######  Definition 8 (Counterfactual Explanation Uncertainty).

Given a CFA $\Phi_{CF}$ and the feature-wise counterfactual uncertainty
$\epsilon_{\mathbf{d}}$, the counterfactual explanation uncertainty is given
by

$\displaystyle\Phi_{CF}^{\epsilon_{\mathbf{d}}}=\Phi_{CF}(\mathbf{x}^{c}\pm\epsilon_{\mathbf{d}}|\mathbf{x}).$
(9)

###  4.4 Path Explanations

Figure 1: Here we illustrate two explanations produced on the Wisconsin Breast
Cancer Dataset given by the proposed QUCE method. We observe how each feature
influenced the change in the prediction in attempting to generate a
counterfactual example. We see the left explanation has almost no uncertainty
in generated explanation, whereas the right image demonstrates a large degree
of uncertainty in the generated counterfactual explanation.

The Path-Integrated Gradients [30, 24, 10, 17, 35] formulation is the only
approach to our knowledge within the landscape of feature attribution methods
that satisfies all the feature attribution axioms in section 3. Therefore we
adopt the path integral formulation and relax the straight-line constraint
seen in IG. To achieve this, recall the set of learned updates on
$\mathbf{x}$, namely $\mathbf{x}^{\Delta}$. It follows that we can produce
explanations over $\mathbf{x}^{\Delta}$ with respect to a classifier $f$.

Formally, let the function $F$ be a continuously differentiable function, the
QUCE explanation takes the path integral formulation such that given a smooth
function
$\psi=\langle\psi_{1},\ldots,\psi_{m}\rangle:[0,1]\rightarrow\mathbb{R}^{m}$
defining a path in $\mathbb{R}^{m}$, where $\psi(\alpha)$ is a point along a
path at $\alpha\in[0,1]$ with $\psi(0)=\mathbf{x}_{\Delta_{0}}$ and
$\psi(1)=\mathbf{x}_{\Delta_{n}}$, the single-path QUCE explainer is defined
as

$\displaystyle\Phi_{\mbox{QUCE}}(\mathbf{x}^{\Delta}):=(\mathbf{x}_{\Delta_{n}}-\mathbf{x}_{\Delta_{0}})\times\bigg{(}\int^{\mathbf{x}_{\Delta_{n}}}_{\mathbf{x}_{\Delta_{0}}}\nabla
F(\psi(\alpha))d\alpha\bigg{)}.$ (10)

It follows that explanation uncertainty with respect to a single generated
counterfactual $\mathbf{x}^{c}$ is given as

$\displaystyle\Phi^{\pm\epsilon_{\mathbf{d}}}_{\mbox{QUCE}}(\mathbf{x}^{c}):=((\mathbf{x}^{c}\pm\epsilon_{\mathbf{d}})-\mathbf{x}^{c})\times\bigg{(}\int^{\mathbf{x}^{c}\pm\epsilon_{\mathbf{d}}}_{\mathbf{x}^{c}}\nabla
F(\psi(\alpha))d\alpha\bigg{)}.$ (11)

We now show through the proof of proposition 1 that the QUCE explanations are
easily computed.

######  Proposition 1.

The QUCE explainer has a computable Riemann approximation solution for each
feature.

###### Proof.

Proof provided in the supplementary material. ∎

Both attributed values from equation 11 illustrate uncertainty in feature
attribution values given by the QUCE explainer. In proposition 2 and its
associated proof we show the implications of weighting uncertainty.

######  Definition 9 ($\lambda$-tolerance).

An increase in $\lambda$-tolerance refers to the reduction of any component
$\lambda$ of $\Lambda$. Likewise, a decrease in $\lambda$-tolerance refers to
the increase of any component $\lambda$ of $\Lambda$.

Naturally, as the weight for uncertainty decreases, we lean towards an
increased tolerance of the effects of uncertainty in our explanations. The
reasoning behind this is that we may sometimes accept a higher degree of
uncertainty, depending on the purpose of generating counterfactuals and the
volatility of the task.

######  Proposition 2.

Increasing the $\lambda$-tolerance of uncertainty provides a more flexible
search space for possible paths to a generative counterfactual example.

###### Proof.

This is trivial subject to $\lambda_{3}$ approaching zero. Proof provided in
the supplementary material. ∎

Due to the stochastic nature of our model under nondeterministic variants of
gradient descent (e.g. optimization with stochastic gradient descent [28]),
and with potentially multiple minima (e.g. we may have two points equally
“close" to the decision bound with different values), we consider a set of
generated counterfactual examples to be given as
$C=\langle\mathbf{x}^{c}_{1,1},\ldots,\mathbf{x}^{c}_{1,k}\rangle$, where $k$
is the number of generated counterfactual examples over some set $\mathbf{x}$.
Given $C$, we can accumulate attribution over many counterfactuals by avoiding
the specification of $\mathbf{x}^{c}$; our lower limit is implicitly assumed
to be our instance to explain $\mathbf{x}$, so that we have

$\displaystyle\Phi_{\mbox{exQUCE}}(\mathbf{x}):=\int_{\mathbf{x}^{c}}\bigg{(}\Phi_{\mbox{QUCE}}(\mathbf{x}^{\Delta})\bigg{)}p_{C}(\mathbf{x}^{c})d\mathbf{x}^{c}$
(12)

where we integrate over $p_{C}$ the distribution of $C$ for all
$\mathbf{x}^{c}\in C$ and we can instead rewrite the integral as an
expectation as follows:

$\displaystyle\Phi_{\mbox{exQUCE}}(\mathbf{x}):=\underset{{\mathbf{x}^{c}\sim
C,\alpha\sim\mathcal{U}(0,1)}}{\mathbbm{E}}\bigg{[}\Phi_{\mbox{QUCE}}(\mathbf{x}^{\Delta})\bigg{]}.$
(13)

Here we let $\alpha\sim\mathcal{U}(0,1)$ indicate interpolation over $\alpha$
for $m$ counterfactual steps in the generator function. Informally, we get the
expectation of the gradients over the piecewise linear path between
counterfactual steps of the generator. We take a similar approach to the
Expected Gradients [12] formulation, except that instead of sampling from a
background set of baselines, we sample from a set of generative counterfactual
examples. We make two arguments as to why we use this approach:

  * •

In explaining a counterfactual outcome, we do not know the specific path taken
and thus we can average over many paths.

  * •

We can invert the path to explain $\mathbf{x}$ and therefore we can have many
generative baselines. This relaxes the specified baseline of many existing
path-based explanation methods.

To further extend on the axiomatic guarantees of path-based explainers, we
show via proposition 3 that completeness holds when working with the many-
paths approach for expected values.

######  Proposition 3.

Given the function $\Phi_{\mbox{exQUCE}}(\mathbf{x})$, the expected difference
in prediction probabilities between generated counterfactuals in the set $C$
with respect to the prediction probability given by $F(\mathbf{x})$, the
following equality holds:

$\displaystyle\underset{{\mathbf{x}^{c}\sim
C,\alpha\sim\mathcal{U}(0,1)}}{\mathbbm{E}}\bigg{[}\Phi_{\mbox{QUCE}}(\mathbf{x}^{\Delta})\bigg{]}$
(14) $\displaystyle=\underset{{\mathbf{x}^{c}\sim
C}}{\mathbbm{E}}\bigg{[}F(\mathbf{x}^{c})-F(\mathbf{x})\bigg{]}$ (15)

###### Proof.

This is a direct consequence of the _completeness_ axiom; the proof is
provided in the supplementary material. ∎

It follows that the expected value approach is equally computable and is a
direct extension of proposition 2. We illustrate this in corollary 4’s simple
proof.

######  Corollary 4.

The expected QUCE variant has a computable Riemann approximation solution for
each feature.

###### Proof.

This follows from proposition 1; the proof is provided in the supplementary
material. ∎

Going further, we demonstrate monotonic relationships for generated
counterfactual instances that are given by the multiple paths approach. This
is a further consequence of the completeness axiom and is expressed in
corollary 5.

######  Corollary 5.

Given two sets of counterfactual examples $C^{1}$ and $C^{2}$ for an instance
$\mathbf{x}$,

$\displaystyle\mbox{if}\ \underset{{\mathbf{x}^{c}\sim
C^{1}}}{\mathbbm{E}}\bigg{[}F(\mathbf{x}^{c})-F(\mathbf{x})\bigg{]}\leq\underset{{\mathbf{x}^{c}\sim
C^{2}}}{\mathbbm{E}}\bigg{[}F(\mathbf{x}^{c})-F(\mathbf{x})\bigg{]},$
$\displaystyle\mbox{then}\ \underset{{\mathbf{x}^{c}\sim
C^{1}}}{\mathbbm{E}}\bigg{[}\Phi_{\mbox{exQUCE}}(\mathbf{x})\bigg{]}\leq\underset{{\mathbf{x}^{c}\sim
C^{2}}}{\mathbbm{E}}\bigg{[}\Phi_{\mbox{exQUCE}}(\mathbf{x})\bigg{]}.$

###### Proof.

This is a direct consequence of proposition 3 and the _completeness_ axiom. ∎

Given we can compute many-paths explanations, it follows that we can also take
the expected gradients for the explanation uncertainty computed by QUCE along
each path, such that

$\displaystyle\Phi^{\pm\epsilon_{\mathbf{d}}}_{\mbox{exQUCE}}(\mathbf{x}^{c}):=\underset{\mathbf{x}^{c}\sim
C,\alpha\sim\mathcal{U}(0,1)}{\mathbbm{E}}\bigg{[}\Phi^{\pm\epsilon_{\mathbf{d}}}_{\mbox{QUCE}}(\mathbf{x}^{c})\bigg{]}.$
(16)

##  5 Experimental Setup

Table 2: Comparison of the average path uncertainty on the generated counterfactual instances. This is experimented over 100 instances from the training and testing sets of each dataset. Here we have 1000 steps (path interpolation instances) for the Riemann approximation of every path-based approach, thus effectively 100$\times$1000 instances. Here the lower value the better. The proposed QUCE method shows superior performance on average when comparing counterfactual path-based approaches. Path $\mathcal{L}_{\epsilon}$ | Lung | Breast | Skin | Lymph | Rectal | COVID | W-BC
---|---|---|---|---|---|---|---
Train |  |  |  |  |  |  |
QUCE | 0.92$\pm$0.32 | 0.82$\pm$0.26 | 0.94$\pm$0.41 | 0.74$\pm$0.06 | 0.86$\pm$0.28 | 1.36$\pm$0.15 | 0.82$\pm$0.16
IG-QUCE | 0.95$\pm$0.32 | 0.84$\pm$0.27 | 0.96$\pm$0.41 | 0.76$\pm$0.06 | 0.86$\pm$0.29 | 1.39$\pm$0.11 | 0.84$\pm$0.20
AGI | 1.94$\pm$1.86 | 1.49$\pm$0.95 | 1.80$\pm$1.47 | 0.92$\pm$0.28 | 2.19$\pm$2.23 | 2.01$\pm$0.15 | 0.93$\pm$0.36
Test |  |  |  |  |  |  |
QUCE | 0.82$\pm$0.34 | 0.91$\pm$0.28 | 0.82$\pm$0.31 | 0.69$\pm$0.19 | 0.80$\pm$0.29 | 0.61$\pm$0.07 | 0.83$\pm$0.29
IG-QUCE | 0.83$\pm$0.33 | 0.91$\pm$0.28 | 0.82$\pm$0.31 | 0.70$\pm$0.19 | 0.79$\pm$0.29 | 0.67$\pm$0.05 | 0.85$\pm$0.33
AGI | 1.22$\pm$1.16 | 1.83$\pm$0.97 | 1.34$\pm$1.20 | 0.89$\pm$0.28 | 1.57$\pm$1.35 | 0.82$\pm$0.11 | 0.96$\pm$0.55
Table 3: Comparison of the average reconstruction error between original instances and their generated counterfactual examples. This is experimented over 100 instances on each dataset. Here we observe that the proposed QUCE method performs best across all datasets. Counterfactual $\mathcal{L}_{\epsilon}$ | Lung | Breast | Skin | Lymph | Rectal | COVID | W-BC
---|---|---|---|---|---|---|---
Train |  |  |  |  |  |  |
QUCE | 1.01 | 0.78 | 0.97 | 0.71 | 0.85 | 1.25 | 0.93
DiCE | 1.97 | 1.02 | 1.48 | 1.10 | 1.27 | 1.57 | 3.63
AGI | 2.93 | 2.07 | 2.61 | 1.01 | 3.48 | 2.40 | 1.22
Test |  |  |  |  |  |  |
QUCE | 0.93 | 0.91 | 0.86 | 0.67 | 0.83 | 0.76 | 1.08
DiCE | 1.95 | 1.09 | 1.33 | 1.18 | 1.22 | 0.92 | 3.00
AGI | 1.68 | 2.75 | 1.87 | 1.02 | 2.38 | 0.84 | 1.41

###  5.1 Datasets

####  5.1.1 The Simulacrum

The Simulacrum222https://simulacrum.healthdatainsight.org.uk/ is a synthetic
dataset used in this study, the Simulacrum is a large dataset developed by
Health Data Insight CiC and derived from anonymous cancer data provided by the
National Disease Registration Service, NHS England. We produce five subsets of
patient records based on ICD-10 codes corresponding to lung cancer, breast
cancer, skin cancer, lymphoma and rectal cancer. These datasets are organised
as survival time classification problems, where patients are predicted a
survival time of either at least 6 months or less than 6 months.

####  5.1.2 COVID Rate of Infection

The COVID rate of infection dataset contains details on control measures,
temperature, humidity and the daily rate of infection for different regions of
the UK. Details on data collection are provided in [16]. This dataset is a
binary classification task identifying an increased rate of infection against
a non-increased rate of infection.

####  5.1.3 Wisconsin Breast Cancer

The Wisconsin Breast Cancer (W-BC) [33] dataset, provided in the scikit-learn
library333https://scikit-
learn.org/stable/modules/generated/sklearn.datasets.load_breast_cancer.html,
is a binary classification dataset that classifies malginant and benign
tumours given a set of independent features from breast mass measurements.

###  5.2 Baseline Methods

For comparison, we consider a selection of methods that aim to generate
counterfactual examples and also a collection of path-based explainers.

####  5.2.1 Diverse Counterfactual Explanations

DiCE, a counterfactual generator, provides feature attribution values for an
instance with respect to its counterfactual examples. We use the DiCE method
as a comparison for generating counterfactual examples, as DiCE is not a path-
based explainer, we can only compare the generated counterfactuals.

####  5.2.2 Integrated Gradients

The IG feature attribution method produces explanations for instances in a
given dataset. To achieve this, the approach is to integrate over the
gradients of a straight-line path from an all-zero baseline vector to the
instance to be explained. We modify this in our experiments so that our
baseline becomes the instance to be explained, while the target instance is
the counterfactual generated by QUCE, so we can evaluate the straight path
solution against the QUCE-generated path.

####  5.2.3 Adversarial Gradient Integration

The Adversarial Gradient Integration (AGI) [24] approach provides an
alternative form of feature attribution that also produces generative
counterfactual examples. The AGI method is the only path-based generative
counterfactual method currently available to our knowledge and thus forms the
primary focus of our comparison. We use the Individual AGI algorithm presented
in their paper.

##  6 Quantitative Evaluation

To evaluate the implementation of generative counterfactual examples, we
propose using the VAE loss to determine how well the counterfactual examples
fit to the underlying data distribution.

###  6.1 Path-Based Uncertainty Comparison

To evaluate the QUCE method, we provide a comparison of uncertainty along a
path. To do this, we use a pre-trained VAE, feeding all generated instances
along any path into the VAE to determine the reconstruction error for all
given instances along a path. The intuition behind this is that a smaller
reconstruction error is associated with a path that better follows the data
distribution and is therefore more “realistic".

In Table 2 we evaluate the path uncertainty across 100 instances from both the
training and test data. From this, we observe that the QUCE method provides
paths that better follow the data distribution when compared against both IG
and AGI on average. Here we reiterate that IG is used with the generated QUCE
instance, which already aims to minimize VAE loss, and is therefore used to
show the minor differences in relaxing the straight-line path requirements,
although this may not always be necessary.

###  6.2 Counterfactual Uncertainty

We use counterfactual uncertainty as a measure to evaluate generative
counterfactual examples given by the QUCE method. This is a simple measure of
the average reconstruction error across the generative counterfactual examples
across each instance in a dataset. We measure this against the DiCE and AGI
methods, as both provide counterfactual examples in their generative process.

In Table 3 we present the counterfactual uncertainty over 100 instances from
both the training and test datasets over each of the datasets. Here, we
observe that QUCE provides a lower value with respect to uncertainty
measurements compared to both the DiCE and AGI methods, implying that the
instances generated by QUCE better follow the data distribution and can thus
be thought of as more likely, subject to the dataset. Further experiments on
the deletion game seen in [34, 1], reconstruction error and a theoretical
evaluation against further explainability axioms presented in [2] are provided
in the supplementary material.

##  7 Conclusion

In this paper, we provide a novel approach that combines generative
counterfactual methods and path-based explainers, minimizing uncertainty along
generated paths and for generated counterfactual examples. We provide an
analysis of the proposed QUCE method on path uncertainty, generative
counterfactual example uncertainty, and proximity. Our approach provides paths
that are less uncertain in their interpolations, so that more reliable
gradients and explanations can be extracted. Similarly, we provide a clear
explanation of uncertainty, including when and where it exists, as seen in the
example explanations provided in Figure 1.

In future work we aim to relax the assumption that all feature values are
continuous to provide more realistic and reliable generated examples.
Exploration of optimal parameters in the QUCE framework is currently a manual
process. Automating this approach would provide greater flexibility and ease
of application upon distributing the method to end users.

Appendix: Supplementary Materials for QUCE

##  Appendix A Computing QUCE Explanations

Proposition 1. _The QUCE explainer has a computable Riemann approximation
solution for each feature._

###### Proof.

Given an instance $\mathbf{x}$ that is the origin of our instance for our
counterfactual explanation, we consider the $\mathbf{x}^{\Delta}$ steps
produced by our learning process. In order to compute QUCE for feature
explanations, we deconstruct the definition to focus on a single step of QUCE.
This can be directly extracted from the additive property of integrals. Recall
the definition:

$\displaystyle\Phi_{\mbox{QUCE}}(\mathbf{x}^{\Delta})$
$\displaystyle:=(\mathbf{x}_{\Delta_{n}}-\mathbf{x}_{\Delta_{0}})\times\bigg{(}\int^{\mathbf{x}_{\Delta_{n}}}_{\mathbf{x}_{\Delta_{0}}}\nabla
F(\psi(\alpha))d\alpha\bigg{)}.$

Expanding this, we obtain

$\displaystyle\Phi_{\mbox{QUCE}}(\mathbf{x}^{\Delta}):=\bigg{(}(\mathbf{x}_{\Delta_{1}}-\mathbf{x}_{\Delta_{0}})\times\bigg{(}\int^{\mathbf{x}_{\Delta_{1}}}_{\mathbf{x}_{\Delta_{0}}}\nabla
F(\psi(\alpha))d\alpha\bigg{)}\bigg{)}$
$\displaystyle+\ldots+\bigg{(}(\mathbf{x}_{\Delta_{n}}-\mathbf{x}_{\Delta_{n-1}})\times\bigg{(}\int^{\mathbf{x}_{\Delta_{n}}}_{\mathbf{x}_{\Delta_{n-1}}}\nabla
F(\psi(\alpha))d\alpha\bigg{)}\bigg{)}.$

Rewriting in terms of partial derivatives we get,

$\displaystyle\Phi_{\mbox{QUCE}}(\mathbf{x}^{\Delta}):=\sum_{j=1}^{J}\bigg{(}\int^{\mathbf{x}_{\Delta_{1}}}_{\mathbf{x}_{\Delta_{0}}}\frac{\partial
F(\psi(\alpha))}{\partial\psi^{j}(\alpha)}\frac{\partial\psi^{j}(\alpha)}{\partial\alpha}d\alpha\bigg{)}$
$\displaystyle+\ldots+\sum_{j=1}^{J}\bigg{(}\int^{\mathbf{x}_{\Delta_{n}}}_{\mathbf{x}_{\Delta_{n-1}}}\frac{\partial
F(\psi(\alpha))}{\partial\psi^{j}(\alpha)}\frac{\partial\psi^{j}(\alpha)}{\partial\alpha}d\alpha\bigg{)}$

where $\psi^{j}$ is along a single feature dimension of a path $\psi$. With
this decomposition, we can further extrapolate by defining each step in
$\mathbf{x}^{\Delta}$ as a piecewise linear path, such that

$\displaystyle\Phi_{\mbox{QUCE}}(\mathbf{x}^{\Delta}):=$
$\displaystyle\sum_{j=1}^{J}\bigg{(}(x_{\Delta_{1}}^{j}-x_{\Delta_{0}}^{j})\times\int^{1}_{\alpha=0}\frac{\partial
F(\mathbf{x}_{\Delta_{0}}+\alpha(\mathbf{x}_{\Delta_{1}}-\mathbf{x}_{\Delta_{0}}))}{\partial
x^{j}}d\alpha\bigg{)}$ $\displaystyle+\ldots+$
$\displaystyle\bigg{(}(x_{\Delta_{n}}^{j}-x_{\Delta_{n-1}}^{j})\times\int^{1}_{\alpha=0}\frac{\partial
F(\mathbf{x}_{\Delta_{n-1}}+\alpha(\mathbf{x}_{\Delta_{n}}-\mathbf{x}_{\Delta_{n-1}}))}{\partial
x^{j}}d\alpha\bigg{)}$

and rewriting for a single feature $j$, we simply remove the summation and
define the explanation over $x^{\Delta,j}$:

$\displaystyle\Phi_{\mbox{QUCE}}(x^{\Delta,j}):=$
$\displaystyle\bigg{(}(x_{\Delta_{1}}^{j}-x_{\Delta_{0}}^{j})\times\int^{1}_{\alpha=0}\frac{\partial
F(\mathbf{x}_{\Delta_{0}}+\alpha(\mathbf{x}_{\Delta_{1}}-\mathbf{x}_{\Delta_{0}}))}{\partial
x^{j}}d\alpha\bigg{)}$ $\displaystyle+\ldots+$
$\displaystyle\bigg{(}(x_{\Delta_{n}}^{j}-x_{\Delta_{n-1}}^{j})\times\int^{1}_{\alpha=0}\frac{\partial
F(\mathbf{x}_{\Delta_{n-1}}+\alpha(\mathbf{x}_{\Delta_{n}}-\mathbf{x}_{\Delta_{n-1}}))}{\partial
x^{j}}d\alpha\bigg{)}.$

We can then rewrite this as a computable Riemann approximation for $K$ steps,
as

$\displaystyle\Phi_{\mbox{QUCE}^{\mathcal{R}}}(x^{\Delta,j}):=$
$\displaystyle\bigg{(}(x_{\Delta_{1}}^{j}-x_{\Delta_{0}}^{j})\times\frac{1}{K}\sum_{k=1}^{K}\frac{\partial
F(\mathbf{x}_{\Delta_{0}}+\frac{k}{K}(\mathbf{x}_{\Delta_{1}}-\mathbf{x}_{\Delta_{0}}))}{\partial
x^{j}}\bigg{)}$ $\displaystyle+\ldots+$
$\displaystyle\bigg{(}(x_{\Delta_{n}}^{j}-x_{\Delta_{n-1}}^{j})\times\frac{1}{K}\sum_{k=1}^{K}\frac{\partial
F(\mathbf{x}_{\Delta_{n-1}}+\frac{k}{K}(\mathbf{x}_{\Delta_{n}}-\mathbf{x}_{\Delta_{n-1}}))}{\partial
x^{j}}\bigg{)}$

which yields a computable explanation over the $j^{th}$ dimension of an
instance $\mathbf{x}$ along a piecewise linear path. This can simply be
executed across all features $j$. ∎

The Riemann approximation is the algorithm for computing explanations that we
use in our implementation carried out on each dimension $j$, which returns a
vector containing the overall attribution over piecewise linear path integral
formulation. It follows that the expected gradients variant can be easily
computed by averaging explanations over each counterfactual example in the set
$C$.

Corollary 1. _The expected QUCE variant has a Riemann approximation solution
for each feature._

###### Proof.

Given a set of $k$ counterfactual examples in the set $C$, we can simply
compute the expected attribution over the $j^{th}$ feature of each generative
counterfactual example as

$\displaystyle\Phi_{\mbox{exQUCE}^{\mathcal{R}}}(x^{\Delta,j})=\frac{1}{k}\sum_{\mathbf{x^{c}}\in
C}\Phi_{\mbox{QUCE}^{\mathcal{R}}}(x^{\Delta,j}).$

∎

Table 4: Comparison of the average sum of feature-wise reconstruction error between original instances and their generated counterfactual examples. This is experimented over 100 instances on each dataset. Here we observe that the QUCE method performs best in generating counterfactuals with minimal uncertainty across all datasets. Counterfactual Reconstruction Error | Lung | Breast | Skin | Lymph | Rectal | COVID | W-BC
---|---|---|---|---|---|---|---
Train |  |  |  |  |  |  |
QUCE | 0.95 | 0.73 | 0.90 | 0.66 | 0.78 | 1.16 | 0.73
DiCE | 1.80 | 0.95 | 1.33 | 1.01 | 1.18 | 1.38 | 1.04
AGI | 2.55 | 1.89 | 2.25 | 0.91 | 2.91 | 2.19 | 0.80
Test |  |  |  |  |  |  |
QUCE | 0.88 | 0.85 | 0.80 | 0.63 | 0.78 | 0.57 | 0.76
DiCE | 1.81 | 0.95 | 1.32 | 0.99 | 1.18 | 1.38 | 1.04
AGI | 1.53 | 2.41 | 1.66 | 0.92 | 2.08 | 0.79 | 0.81
Table 5: Comparison of the deletion scores for counterfactual generative methods that provide feature attribution values. This is experimented over 100 instances on each dataset. Here the lower the value the better. We observe that the proposed QUCE method performs best across a larger fraction of datasets. Deletion | Lung | Breast | Skin | Lymph | Rectal | COVID | W-BC
---|---|---|---|---|---|---|---
QUCE | 0.556 | 0.689 | 0.656 | 0.611 | 0.669 | 0.699 | 0.632
DiCE | 0.561 | 0.688 | 0.649 | 0.619 | 0.669 | 0.710 | 0.637
AGI | 0.559 | 0.683 | 0.659 | 0.607 | 0.670 | 0.728 | 0.648

##  Appendix B Proof of Proposition 2

Proposition 2._Increasing the $\lambda$-tolerance of uncertainty provides a
more flexible search space for possible paths to a generative counterfactual
example. _

###### Proof.

Recall the objective function in equation 2. For simplicity we let
$\lambda_{1}=0$, $0<\lambda_{2}\leq 1$ and $0\leq\lambda_{3}\leq 1$, such that

$\displaystyle\mathcal{G}(\mathbf{x})=\lambda_{2}\mathcal{L}_{\delta}+\lambda_{3}\mathcal{L}_{\epsilon}$
$\displaystyle=\frac{\lambda_{2}}{2}||\mathbf{x}^{c}-\mathbf{x}||^{2}+$
$\displaystyle\lambda_{3}(\mathbbm{E}_{q_{\theta^{*}}}[log\
q_{q_{\theta^{*}}}(\mathbf{z}|\mathbf{x}^{c})-log\
p_{\psi^{*}}(\mathbf{z})]-\mathbbm{E}_{q_{\theta^{*}}}\ log\
p_{\psi^{*}}(\mathbf{x}^{c}|\mathbf{z})).$

Trivially, as $\lambda_{3}$ decreases toward zero (we accept more
uncertainty), the freedom in the distance function increases, as it not
constrained by uncertainty, since

$\displaystyle\lim_{\lambda_{3}\rightarrow
0^{+}}(\frac{\lambda_{2}}{2}||\mathbf{x}^{c}-\mathbf{x}||^{2}+$
$\displaystyle\lambda_{3}(\mathbbm{E}_{q_{\theta^{*}}}[log\
q_{q_{\theta^{*}}}(\mathbf{z}|\mathbf{x}^{c})-log\
p_{\psi^{*}}(\mathbf{z})]-\mathbbm{E}_{q_{\theta^{*}}}\ log\
p_{\psi^{*}}(\mathbf{x}^{c}|\mathbf{z})))$
$\displaystyle=\frac{\lambda_{2}}{2}||\mathbf{x}^{c}-\mathbf{x}||^{2}.$

Here the eigenvalues of the Hessian are given by $\lambda_{2}$ and since
$\lambda_{2}>0$, the Hessian is positive definite and thus the search space is
convex, implying the global minima are conditioned on $\mathcal{L}_{\delta}$
when uncertainty $\mathcal{L}_{\epsilon}$ is relaxed. ∎

##  Appendix C Proof of Proposition 3

Proposition 3. _Given the function $\Phi_{\mbox{exQUCE}}(\mathbf{x})$, the
expected difference in prediction probabilities between generated
counterfactuals in the set $C$ with respect to the prediction probability
given by $F(\mathbf{x})$, the following equality holds: _

$\displaystyle\underset{{\mathbf{x}^{c}\sim
C,\alpha\sim\mathcal{U}(0,1)}}{\mathbbm{E}}\bigg{[}\Phi_{\mbox{QUCE}}(\mathbf{x}^{\Delta})\bigg{]}$
(17) $\displaystyle=\underset{{\mathbf{x}^{c}\sim
C}}{\mathbbm{E}}\bigg{[}F(\mathbf{x}^{c})-F(\mathbf{x})\bigg{]}.$ (18)

###### Proof.

Due to the _completeness_ axiom the following holds true:

$\displaystyle
F(\mathbf{x}_{\Delta_{n}})-F(\mathbf{x}_{\Delta_{0}})=(\mathbf{x}_{\Delta_{n}}-\mathbf{x}_{\Delta_{0}})\times\bigg{(}\int^{\mathbf{x}_{\Delta_{n}}}_{\mathbf{x}_{\Delta_{0}}}\nabla
F(\psi(\alpha))d\alpha\bigg{)}$ (19)

and extending this via $\mathbf{x}_{\Delta_{n}}=\mathbf{x}^{c}$ and
$\mathbf{x}_{\Delta_{0}}=\mathbf{x}$ respectively, and thus we have

$\displaystyle F(\mathbf{x}^{c})-F(\mathbf{x})$

while rearranging the RHS of equation 19 we have

$\displaystyle\underbrace{\mathbf{x}^{c}\times\bigg{(}\int^{\mathbf{x}^{c}}_{\mathbf{x}}\nabla
F(\psi(\alpha))d\alpha\bigg{)}}_{F(\mathbf{x}^{c})}$ (20)
$\displaystyle-\underbrace{\mathbf{x}\times\bigg{(}\int^{\mathbf{x}^{c}}_{\mathbf{x}}\nabla
F(\psi(\alpha))d\alpha\bigg{)}}_{F(\mathbf{x})}$ (21)

which follows by relaxing a strict definition of $\mathbf{x}^{c}$, where we
instead use a set of generated counterfactuals in the set $C$. Our proposed
approach takes the integral over $\mathbf{x}^{c}$ with respect to a change in
the set $C$, yielding

$\displaystyle\int_{\mathbf{x}^{c}}\bigg{(}\mathbf{x}^{c}\times\big{(}\int^{\mathbf{x}^{c}}_{\mathbf{x}}\nabla
F(\psi(\alpha))\big{)}d\alpha\bigg{)}p_{C}d\mathbf{x}^{c}$ (22)
$\displaystyle=\underset{\mathbf{x}^{c}\sim
C,\alpha\sim\mathcal{U}(0,1)}{\mathbbm{E}}\bigg{[}\mathbf{x}^{c}\times\nabla
F(\psi(\alpha))\bigg{]}$ (23) $\displaystyle=\underset{{\mathbf{x}^{c}\sim
C}}{\mathbbm{E}}\bigg{[}F(\mathbf{x}^{c})\bigg{]}$ (24)

for equation 20 and $F(\mathbf{x})$ for equation 21. and since $F(\mathbf{x})$
is a constant, namely

$\displaystyle\underset{{\mathbf{x}^{c}\sim
C}}{\mathbbm{E}}\bigg{[}F(\mathbf{x}^{c})\bigg{]}-F(\mathbf{x})=\underset{{\mathbf{x}^{c}\sim
C}}{\mathbbm{E}}\bigg{[}F(\mathbf{x}^{c})-F(\mathbf{x})\bigg{]},$ (25)

equations 17 and 18 are equivalent. ∎

##  Appendix D Counterfactual Reconstruction Error

In addition to evaluating the VAE loss, we also analyze the average
reconstruction error per instance across all 100 instances on both the
training and test datasets. This highlights the closeness of the reconstructed
sample against the ground truth counterfactual generated by different methods.
In Table 4 we observe that our proposed QUCE method provides better reproduced
counterfactuals through the VAE than either the DiCE or AGI methods.

##  Appendix E QUCE Evaluated against Further Properties of Explainability

In the work of [2] the authors present a desirable set of axiomatic
foundations for XAI methods. As a brief informal overview, we consider the
following proposed axioms:

  * •

Success: The explainer method should be able to produce explanations for any
instance.

  * •

Explainability: An explanation method should provide informative explanations.
An empty explanation here is not recommended.

  * •

Irreducability: An explanation should not contain irrelevant information.

  * •

Representativity: An explanation should be possible on unseen instances.

  * •

Relevance: Information should only be included if it impacts the prediction.

We evaluate QUCE against these properties as interpreted with respect to our
model design. We show that it is inherently straightforward to prove that our
proposed QUCE method satisfies these desirable axioms.

######  Proposition 6.

The QUCE method can always provide an explanation satisfying the success
axiom.

###### Proof.

Whilst a counterfactual may not always be _valid_ , as a direct implication of
the generative learning process QUCE will achieve an explanation. ∎

######  Corollary 7.

The QUCE method satisifies the Explainability axiom.

###### Proof.

Assuming that different instances generated by QUCE do not have the same
prediction probability with respect to the target class, it follows
immediately from proposition 3 that explainability holds. ∎

######  Corollary 8.

The QUCE method satisfies the Irreducability axiom.

###### Proof.

Here, we characterize irrelevance under our own interpretation: since a
feature that does not change does not affect the predicted outcome, it should
be assigned zero attribution. Then directly from the definition of QUCE it is
clear that irreducability holds, as the gradients are multiplied by a zero-
value scalar for the same valued features. ∎

######  Proposition 9.

The QUCE method satisfies the Representativity axiom.

###### Proof.

It is easy to see that any instance with the same dimensionality of the
instances from a training dataset can utilize the QUCE approach. ∎

######  Corollary 10.

The QUCE method satisfies the relevance axiom.

###### Proof.

Relevance holds as a direct implication of irreducability as seen in corollary
8 and the fact that gradients are traced over the change in the predictions
along paths, thereby guaranteeing model-specific relevance. ∎

##  Appendix F Experimental Setup

For the experiments presented in this paper, the details of hyper-parameters
and experimental setup are found in the notebook file at:
https://github.com/jamie-duell/QUCE. For the comparative experiments in this
paper we set up QUCE with single-path solutions using the Adam optimizer. The
empirical intuition for using Adam for a single-path approach is as follows:
as the loss is minimized and points along the path become more reliable there
will be an increase in the frequency of points as we approach the solution;
thus averaging the gradient along such paths become more reliable. Similarly
the approach is deterministic, therefore does not require multiple paths as
with alternative optimizers.

##  Appendix G Deletion Experiments

To compare counterfactual feature attribution methods we evaluate the deletion
score, a common metric used for evaluating feature attribution methods for
identifying important features. The deletion score is used in various studies
[34, 1]; here, a lower value indicates better performance. In Table 5 we
observe that the QUCE method performs better on average than both DiCE and AGI
for counterfactual feature attribution performance.

## References

  * Akhtar and Jalwana [2023] N. Akhtar and M. A. A. K. Jalwana.  Towards credible visual model interpretation with path attribution.  In _Proceedings of the 40th International Conference on Machine Learning_ , volume 202 of _Proceedings of Machine Learning Research_ , pages 439–457. PMLR, 23–29 Jul 2023.  URL https://proceedings.mlr.press/v202/akhtar23a.html.
  * Amgoud and Ben-Naim [2022] L. Amgoud and J. Ben-Naim.  Axiomatic foundations of explainability.  In _IJCAI_ , pages 636–642, 2022.
  * Angiulli et al. [2023] F. Angiulli, F. Fassetti, and L. Ferragina.  Reconstruction error-based anomaly detection with few outlying examples, 2023.
  * Antoran et al. [2021] J. Antoran, U. Bhatt, T. Adel, A. Weller, and J. M. Hernández-Lobato.  Getting a CLUE: A method for explaining uncertainty estimates.  In _International Conference on Learning Representations_ , 2021.  URL https://openreview.net/forum?id=XSLF1XFq5h.
  * Beutler et al. [2016] L. E. Beutler, K. Someah, S. Kimpara, and K. Miller.  Selecting the most appropriate treatment for each patient.  _International Journal of Clinical and Health Psychology_ , 16(1):99–108, Jan. 2016.  ISSN 1697-2600.  10.1016/j.ijchp.2015.08.001.  URL http://dx.doi.org/10.1016/j.ijchp.2015.08.001.
  * Bharati et al. [2023] S. Bharati, M. R. H. Mondal, and P. Podder.  A review on explainable artificial intelligence for healthcare: Why, how, and when?  _IEEE Transactions on Artificial Intelligence_ , pages 1–15, 2023.  10.1109/TAI.2023.3266418.
  * Bull et al. [2020] J. W. Bull, N. Strange, R. J. Smith, and A. Gordon.  Reconciling multiple counterfactuals when evaluating biodiversity conservation impact in social-ecological systems.  _Conservation Biology_ , 35(2):510–521, Sept. 2020.  ISSN 1523-1739.  10.1111/cobi.13570.  URL http://dx.doi.org/10.1111/cobi.13570.
  * Chen et al. [2023] X.-Q. Chen, C.-Q. Ma, Y.-S. Ren, Y.-T. Lei, N. Q. A. Huynh, and S. Narayan.  Explainable artificial intelligence in finance: A bibliometric review.  _Finance Research Letters_ , 56:104145, 2023.  ISSN 1544-6123.  https://doi.org/10.1016/j.frl.2023.104145.  URL https://www.sciencedirect.com/science/article/pii/S1544612323005172.
  * Dejl et al. [2023] A. Dejl, H. Ayoobi, M. Williams, and F. Toni.  Cafe: Conflict-aware feature-wise explanations, 2023.
  * Duell et al. [2023a] J. Duell, X. Fan, H. Fu, and M. Seisenberger.  Batch integrated gradients: Explanations for temporal electronic health records.  In _Artificial Intelligence in Medicine_ , pages 120–124, Cham, 2023a. Springer Nature Switzerland.  ISBN 978-3-031-34344-5.
  * Duell et al. [2023b] J. Duell, M. Seisenberger, T. Zhong, H. Fu, and X. Fan.  A formal introduction to batch-integrated gradients for temporal explanations.  In _2023 IEEE 35th International Conference on Tools with Artificial Intelligence (ICTAI)_ , pages 452–459, Los Alamitos, CA, USA, nov 2023b. IEEE Computer Society.  10.1109/ICTAI59109.2023.00072.  URL https://doi.ieeecomputersociety.org/10.1109/ICTAI59109.2023.00072.
  * Erion et al. [2021] G. Erion, J. D. Janizek, P. Sturmfels, S. M. Lundberg, and S.-I. Lee.  Improving performance of deep learning models with axiomatic attribution priors and expected gradients.  _Nature Machine Intelligence_ , 3(7):620–631, May 2021.  ISSN 2522-5839.  10.1038/s42256-021-00343-w.  URL http://dx.doi.org/10.1038/s42256-021-00343-w.
  * Guidotti [2022] R. Guidotti.  Counterfactual explanations and how to find them: literature review and benchmarking.  _Data Mining and Knowledge Discovery_ , Apr. 2022.  ISSN 1573-756X.  10.1007/s10618-022-00831-6.  URL http://dx.doi.org/10.1007/s10618-022-00831-6.
  * Guyomard et al. [2021] V. Guyomard, F. Fessant, T. Bouadi, and T. Guyet.  Post-hoc counterfactual generation with supervised autoencoder.  In _Machine Learning and Principles and Practice of Knowledge Discovery in Databases_ , pages 105–114, Cham, 2021. Springer International Publishing.  ISBN 978-3-030-93736-2.
  * Höfler [2005] M. Höfler.  Causal inference based on counterfactuals.  _BMC Medical Research Methodology_ , 5(1), Sept. 2005.  ISSN 1471-2288.  10.1186/1471-2288-5-28.  URL http://dx.doi.org/10.1186/1471-2288-5-28.
  * Kapcia et al. [2021] M. Kapcia, H. Eshkiki, J. Duell, X. Fan, S. Zhou, and B. Mora.  Exmed: An ai tool for experimenting explainable ai techniques on medical data analytics.  In _2021 IEEE 33rd International Conference on Tools with Artificial Intelligence (ICTAI)_ , pages 841–845, 2021.  10.1109/ICTAI52525.2021.00134.
  * Kapishnikov et al. [2021] A. Kapishnikov, S. Venugopalan, B. Avci, B. Wedin, M. Terry, and T. Bolukbasi.  Guided integrated gradients: an adaptive path method for removing noise.  In _2021 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 5048–5056, Los Alamitos, CA, USA, jun 2021. IEEE Computer Society.  10.1109/CVPR46437.2021.00501.  URL https://doi.ieeecomputersociety.org/10.1109/CVPR46437.2021.00501.
  * Kovvuri et al. [2022] V. R. R. Kovvuri, S. Liu, M. Seisenberger, X. Fan, B. Müller, and H. Fu.  On understanding the influence of controllable factors with a feature attribution algorithm: a medical case study.  In _INISTA_ , pages 1–8, 2022.
  * Li et al. [2023] X. Li, D. Pan, C. Li, Y. Qiang, and D. Zhu.  Negative flux aggregation to estimate feature attributions.  In E. Elkind, editor, _Proceedings of the Thirty-Second International Joint Conference on Artificial Intelligence, IJCAI-23_ , pages 446–454. International Joint Conferences on Artificial Intelligence Organization, 8 2023.  10.24963/ijcai.2023/50.  URL https://doi.org/10.24963/ijcai.2023/50.  Main Track.
  * Lundberg and Lee [2017] S. M. Lundberg and S. Lee.  A unified approach to interpreting model predictions.  In _NeurIPS_ , pages 4765–4774, 2017.
  * Ma et al. [2022] J. Ma, R. Guo, S. Mishra, A. Zhang, and J. Li.  Clear: Generative counterfactual explanations on graphs.  In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, _Advances in Neural Information Processing Systems_ , volume 35, pages 25895–25907. Curran Associates, Inc., 2022.  URL https://proceedings.neurips.cc/paper_files/paper/2022/file/a69d7f3a1340d55c720e572742439eaf-Paper-Conference.pdf.
  * Madaan et al. [2023] A. Madaan, T. Chowdhury, N. Rana, J. Allan, and T. Chakraborty.  Uncertainty in additive feature attribution methods, 2023.
  * Mothilal et al. [2020] R. K. Mothilal, A. Sharma, and C. Tan.  Explaining machine learning classifiers through diverse counterfactual explanations.  In _FAccT_ , page 607–617, 2020.
  * Pan et al. [2021] D. Pan, X. Li, and D. Zhu.  Explaining deep neural network models with adversarial gradient integration.  In Z.-H. Zhou, editor, _Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence, IJCAI-21_ , pages 2876–2883. International Joint Conferences on Artificial Intelligence Organization, 8 2021.  10.24963/ijcai.2021/396.  URL https://doi.org/10.24963/ijcai.2021/396.  Main Track.
  * Prosperi et al. [2020] M. Prosperi, Y. Guo, M. Sperrin, J. S. Koopman, J. S. Min, X. He, S. Rich, M. Wang, I. E. Buchan, and J. Bian.  Causal inference and counterfactual prediction in machine learning for actionable healthcare.  _Nature Machine Intelligence_ , 2(7):369–375, July 2020.  ISSN 2522-5839.  10.1038/s42256-020-0197-y.  URL http://dx.doi.org/10.1038/s42256-020-0197-y.
  * Ribeiro et al. [2016] M. T. Ribeiro, S. Singh, and C. Guestrin.  "Why should I trust you?": Explaining the predictions of any classifier.  In _ACM SIGKDD_ , page 1135–1144, 2016.
  * Sagar [2022] A. Sagar.  Uncertainty quantification using variational inference for biomedical image segmentation.  In _2022 IEEE/CVF Winter Conference on Applications of Computer Vision Workshops (WACVW)_ , pages 44–51, Los Alamitos, CA, USA, jan 2022. IEEE Computer Society.  10.1109/WACVW54805.2022.00010.  URL https://doi.ieeecomputersociety.org/10.1109/WACVW54805.2022.00010.
  * Sinha and Griscik [1971] N. K. Sinha and M. P. Griscik.  A stochastic approximation method.  _IEEE Transactions on Systems, Man, and Cybernetics_ , SMC-1(4):338–344, 1971.  10.1109/TSMC.1971.4308316.
  * Slack et al. [2021] D. Slack, A. Hilgard, S. Singh, and H. Lakkaraju.  Reliable post hoc explanations: Modeling uncertainty in explainability.  In M. Ranzato, A. Beygelzimer, Y. Dauphin, P. Liang, and J. W. Vaughan, editors, _Advances in Neural Information Processing Systems_ , volume 34, pages 9391–9404. Curran Associates, Inc., 2021.  URL https://proceedings.neurips.cc/paper_files/paper/2021/file/4e246a381baf2ce038b3b0f82c7d6fb4-Paper.pdf.
  * Sundararajan et al. [2017] M. Sundararajan, A. Taly, and Q. Yan.  Axiomatic attribution for deep networks.  page 3319–3328, 2017.
  * Torabi et al. [2023] H. Torabi, S. L. Mirtaheri, and S. Greco.  Practical autoencoder based anomaly detection by using vector reconstruction error.  _Cybersecurity_ , 6(1), Jan. 2023.  ISSN 2523-3246.  10.1186/s42400-022-00134-9.  URL http://dx.doi.org/10.1186/s42400-022-00134-9.
  * Wachter et al. [2017] S. Wachter, B. Mittelstadt, and C. Russell.  "counterfactual explanations without opening the black box: Automated decisions and the gdpr.  In _arXiv_ , 2017.
  * Wolberg et al. [1995] W. Wolberg, O. Mangasarian, N. Street, and W. Street.  Breast Cancer Wisconsin (Diagnostic).  UCI Machine Learning Repository, 1995.  DOI: https://doi.org/10.24432/C5DW2B.
  * Yang et al. [2023a] P. Yang, N. Akhtar, Z. Wen, and A. Mian.  Local path integration for attribution.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 37(3):3173–3180, Jun. 2023a.  10.1609/aaai.v37i3.25422.  URL https://ojs.aaai.org/index.php/AAAI/article/view/25422.
  * Yang et al. [2023b] R. Yang, B. Wang, and M. Bilgic.  IDGI: A framework to eliminate explanation noise from integrated gradients.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2023b.  URL http://www.cs.iit.edu/~ml/pdfs/yang-cvpr23.pdf.