# Data-centric Prediction Explanation via Kernelized Stein Discrepancy

Firstname1 Lastname1  Firstname2 Lastname2  Firstname3 Lastname3  Firstname4
Lastname4  Firstname5 Lastname5  Firstname6 Lastname6  Firstname7 Lastname7
Firstname8 Lastname8  Firstname8 Lastname8

###### Abstract

Existing example-based prediction explanation methods often bridge test and
training data points through the model’s parameters or latent representations.
While these methods offer clues to the causes of model predictions, they often
exhibit innate shortcomings, such as incurring significant computational
overhead or producing coarse-grained explanations. This paper presents a
Highly-precise and Data-centric Explanation (HD-Explain), a straightforward
prediction explanation method exploiting properties of Kernelized Stein
Discrepancy (KSD). Specifically, the KSD uniquely defines a parameterized
kernel function for a trained model that encodes model-dependent data
correlation. By leveraging the kernel function, one can identify training
samples that provide the best predictive support to a test point efficiently.
We conducted thorough analyses and experiments across multiple classification
domains, where we show that HD-Explain outperforms existing methods from
various aspects, including 1) preciseness (fine-grained explanation), 2)
consistency, and 3) computation efficiency, leading to a surprisingly simple,
effective, and robust prediction explanation solution. Furthermore, we
demonstrated that the proposed method can reveal underlying distribution
mismatches between training samples, making it an effective tool for data
debugging.

Machine Learning, ICML

##  1 Introduction

As one of the decisive factors affecting the performance of a Machine Learning
(ML) model, training data points are of great value in promoting the model’s
transparency and trustworthiness, including explaining prediction results,
tracing sources of errors, or summarizing the characteristics of the model
(cai2019effects; anik2021data; nam2022human). The challenges of example-based
prediction explanation mainly come from retrieving relevant data points from a
vast pool of training samples or justifying the rationale of such explanations
(lim2019these; zhou2021evaluating).

Table 1: Summary of existing Post-hoc Example-based Prediction Explanation
Methods that work with deep neural networks. Practicality of the whole model
explanation is measured by the feasibility of explaining the prediction of
ResNet-18 trained on CIFAR-10 with a single A100 GPU machine. CIFAR-10 is a
small benchmark data with 50000 training samples.

Method |  Explanation of |  Need optimization as sub-routine | Whole model explanation |  Inference computation complexity bounded by |  Memory/cache (of each training sample) bounded by
---|---|---|---|---|---
Theoretical | Practical
Influence Function | Original Model |  Yes (Iterative HVP approximation) | Yes | No |  1.$H_{\theta}^{-1}\nabla_{\theta}L(\mathbf{x}_{t},\theta)$ approximation 2\. $<\nabla_{\theta}L(\mathbf{x},\theta),H_{\theta}^{-1}\nabla_{\theta}L(\mathbf{x}_{t},\theta)>$ | Size of model parameters
RPS | Fine-tuned Model |  Yes (L2 regularized last layer retrain) | No | No |  1.last layer representation $\mathbf{f}_{t}$ 2.$<\alpha_{i}\mathbf{f}_{i},\mathbf{f}_{t}>$ |  Size of model parameters of the last layer
TracIn* | Original Model | No | Yes | No |  1.$\nabla_{\theta}L(\mathbf{x}_{t},\theta)$ approximation 2\. $<\nabla_{\theta}L(\mathbf{x},\theta),\nabla_{\theta}L(\mathbf{x}_{t},\theta)>$ | Size of model parameters
HD-Explain | Original Model | No | Yes | Yes |  1.$\nabla_{\mathbf{x}_{t}}f(\mathbf{x}_{t},\theta)_{y_{t}}$ 2\. Closed-form $k_{\theta}(\mathbf{x},\mathbf{x}_{t})$ defined by KSD | Size of data dimension

* TracIn typically requires to access training process. Here, TracIn* refers to a special case that only use the last training checkpoint.

Modern example-based prediction explanation methods commonly approach the
above challenges by constructing an influence chain between training and test
data points (li2020survey; nam2022human; tsai2023sample). The influence chain
could be either data points’ co-influence on model parameters or their
similarity in terms of latent representations. In particular, Influence
Function (koh2017understanding), one of the representative, model-aware
explanation methods, looks for the shift of the model parameters (due to up-
weighting each training sample) as the sample’s influence score. Since
computing the inverse Hessian matrix is challenging, the approach adapts
Conjugate Gradients Stochastic Estimation and the Perlmutter trick to reduce
its computation cost. Representer Point Selection (RPS) (yeh2018representer),
as another example, reproduces the representer theorem by refining the trained
neural network model with $L2$ regularization, such that the influence score
of each training sample can be represented as the gradient of the predictive
layer. While computationally efficient, RPS is criticized for producing
coarse-grained explanations that are more class-level rather than instance-
level explanations (sui2021representer) (In this paper we use the instance
level explanation and example-based explanation interchangeably.). Multiple
later variants (pruthi2020estimating; sui2021representer) attempted to
mitigate the drawbacks above, but their improvements were often limited by the
cause of their shared theoretical scalability bounds.

This paper presents Highly-precise and Data-centric Explanation (HD-Explain),
a post-hoc, model-aware, example-based explanation solution for neural
classifiers. Instead of relying on data co-influence on model parameters or
feature representation similarity, HD-Explain retains the influence chain
between training and test data points by exploiting the underrated properties
of Kernelized Stein Discrepancy (KSD) between the trained predictive model and
its training dataset. Specifically, we note that the Stein operator augmented
kernel uniquely defines a pairwise data correlation (in the context of a
trained model) whose expectation on the training dataset results in the
minimum KSD (as a discrete approximation) compared to that of the dataset
sampled from different distributions. By exploiting this property, we can 1)
reveal a subset of training data points that provides the best predictive
support to the test point and 2) identify the potential distribution mismatch
among training data points. Jointly leveraging these advantages, HD-Explain
can produce explanations faithful to the original trained model with a
corresponding confidence level.

We conducted a wide range of qualitative and quantitative experiments and
analyses to demonstrate the effectiveness and efficiency of HD-Explain. Our
empirical results show that the HD-Explain offers fine-grained, instance-level
explanations with remarkable computational efficiency on multiple
applications, compared to well-known example-based prediction explanation
solutions. Given its algorithmic simplicity and scalability, we believe HD-
Explain is an excellent addition to the literature on example-based prediction
explanation research.

##  2 Preliminary and Related Work

###  2.1 Post-hoc Classifier Explanation by Examples

Post-hoc Classifier Explanation by Examples (a.k.a prototypes, representers)
refers to a category of classifier explanation approaches that pick a subset
of training data points as prediction explanations without accessing the model
training process. Its research history spans from the model-intrinsic approach
to the recent impact-based approach.

Model-inherent approaches are essentially a feature of self-explanatory
machine learning models such as $k$-nearest neighbor (peterson2009k) or
decision tree; For a given test data point, similar data points on the raw
feature space can be efficiently selected as explanations through the inherent
decision making mechanism of the self-explanatory machine learning models. In
fact, attracted by their inherent explanatory power, multiple well-known works
attempted to compile complex black-box models into self-explanatory models for
enabling prediction explanation (frosst2017distilling), while computationally
inefficient.

To unlock the general explanatory power applicable to black-box models,
multiple later studies suggest to fall back to statistics-based solutions,
looking for prototype samples that could represent common data points or play
critical roles in data distribution. MMD-critic (kim2016examples) and
Normative and Comprarative explanations (cai2019effects) are the well-known
examples in this category. Unfortunately, those approaches often overlook the
characteristics of trained models, making their prediction explanations
general to the training dataset rather than a trained model instance.

Recently, influence-based methods have emerged as the prevailing technique in
model explanation (li2020survey; nam2022human; bae2022if; park2023trak).
Influence function (koh2017understanding), as one of the earliest influence-
based solutions, bridges the outcome of a prediction task to training data
points by, first, evaluating training data’s influence on the model parameters
and, then, estimating how model parameter changes affect prediction.
Similarly, Representative Point Selection (RPS) (yeh2018representer) builds
such an influence chain by fitting the representation theorem, where the
weighted product between the representations of test and training samples
comes into play. Concerning the computational overhead of previous work, the
later solution TracIn (pruthi2020estimating) proposed a simple approximation
of the influence function via a first-order Taylor expansion (essentially
Neural Tangent Kernel (jacot2018neural)), successfully discarding the inverse
Hessian matrix from the influence chain formulation. BoostIn
(brophy2023adapting) further extends TracIn and is dedicated to interpreting
the predictions of gradient-boosted decision trees. RPS-LJE
(sui2021representer), on the other hand, alleviated the inconsistent
explanation problem of RPS through Local Jacobian Expansion. In the latest
publication (tsai2023sample), all of the above methods described in this
paragraph are identified as special cases of Generalized Representers but with
different chosen kernels.

One limitation of the current influence-based methods is that they attribute
the influence of each training data point to the parameters of the trained
model as an essential intermediate step. Indeed, as the nature of stochastic
gradient descent (the dominating training strategy of neural networks),
isolating such contribution is barely possible without 1) relying on
approximations or 2) accessing the training process. Unfortunately, either
solution would result in performance degradation or heavy computational
overhead (schioppa2022scaling). See Table 1. Hence, this work delves into the
exploration of an alternative influence connection between training and test
data points without exploiting the perturbation of model parameters.

###  2.2 Kernelized Stein Discrepancy

The idea of Kernelized Stein Discrepancy (KSD) (liu2016kernelized) can be
traced back to a theorem called Stein Identity (kattumannil2009stein) that
states, if a smooth distribution $p(x)$ and a function $\phi(x)$ satisfy
$\lim_{||x||\rightarrow\infty}p(x)\phi(x)=0$, we have

$\mathbb{E}_{x\sim p}[\phi(x)\nabla_{x}\log
p(x)+\nabla_{x}\phi(x)]=0,\quad\forall f.$

The above expression could be further abstracted to use function operator
$\mathcal{A}_{p}$ (a.k.a Stein operator) such that

$\mathcal{A}_{p}\phi(x)=\phi(x)\nabla_{x}\log p(x)+\nabla_{x}\phi(x),$

where the operator encodes distribution $p(x)$ in the form of derivative to
input (a.k.a score function).

Stein’s identity offers a mechanism to measure the gap of two distributions by
assuming the variable $x$ is sampled from a different distribution $q\neq p$
such that

$\sqrt{\mathbb{S}(q,p)}=\max_{\phi\in\mathcal{F}}\mathbb{E}_{x\sim
q}[\mathcal{A}_{p}\phi(x)],$

where the expression takes the most discriminant $\phi$ that maximizes the
violation of Stein’s identity to quantify the distribution discrepancy. This
discrepancy is, accordingly, referred as Stein Discrepancy.

The challenge of computing Stein Discrepancy comes from the selection of
function set $\mathcal{F}$, which motivates the later innovation of KSD that
takes $\mathcal{F}$ to be the unit ball of a reproducing kernel Hilbert space
(RKHS). By leveraging the reproducing property of RKHS, the KSD could be
eventually transformed into

$\mathbb{S}(q,p)=\mathbb{E}_{x,x^{\prime}\sim q}[\kappa_{p}(x,x^{\prime})]$

where
$\kappa_{p}(x,x^{\prime})=\mathcal{A}_{p}^{x}\mathcal{A}_{p}^{x^{\prime}}k(x,x^{\prime})$
that can work with arbitrary kernel function $k(x,x^{\prime})$. See Appendix B

In the literature, KSD has been adopted for tackling three types of
application tasks – 1) parameter inference (barp2019minimum), 2) Goodness-of-
fit tests (chwialkowski2016kernel; liu2016kernelized; yang2018goodness), and
3) particle filtering (sampling) (gorham2020stochastic; korba2021kernel).
However, to the best of our knowledge, its innate property that uniquely
defines model-dependent data correlation has never been exploited, which, we
note, is valuable to interpret model behaviour from various aspects, including
instance-level prediction explanation and global prototypical explanations.

Figure 1: Varying of Kernelized Stein Discrepancy given the shift of training
data distribution on Two Moon dataset.

##  3 Example-based Prediction Explanation via Kernelized Stein Discrepancy
(HD-Explain)

Consider a trained classifier $f_{\theta}$ as the outcome of a training
process with Maximum Likelihood Estimation (MLE)

$\operatorname*{argmax}_{\theta}\mathbb{E}_{(\mathbf{x},y)\sim P_{D}}[\log
P_{\theta}(y|\mathbf{x})].$

Theoratically, maximizing observation likelihood is equivalent to minimizing a
KL divergence between data distribution $P_{D}$ and the parameterized
distribution $P_{\theta}$ such that

$\displaystyle\mathbb{D}_{\text{KL}}(P_{D},P_{\theta})=\mathbb{E}_{(\mathbf{x},y)\sim
P_{D}}\left[\log\frac{P_{D}(\mathbf{x},y)}{P_{\theta}(\mathbf{x},y)}\right]$
$\displaystyle=-\underbrace{\mathbb{E}_{(\mathbf{x},y)\sim P_{D}}[\log
P_{\theta}(y|\mathbf{x})]}_{\text{Likelihood}}+\underbrace{\mathbb{E}_{(\mathbf{x},y)\sim
P_{D}}[\log P_{D}(y|\mathbf{x})]}_{\text{constant}},$

which, in turns, is proven recently to align with minimizing Stein Discrepancy
in the form of gradient descent (liu2016stein)

$\nabla_{\theta}\mathbb{D}_{\text{KL}}(P_{D},P_{\theta})=-\mathbb{S}(P_{D},P_{\theta}).$

The chain of reasoning above shows that a well-trained classifier $f_{\theta}$
through gradient-descent should lead to minimum Stein Discrepancy between the
training dataset distribution and the model encoded distribution
$\mathbb{S}(P_{D},P_{\theta})$ 111Since $P_{D}$ is discrete distribution while
$P_{\theta}$ is continuous, the Stein Discrepancy between the two
distributions will not recap Stein Identity ($=0$) with limited number of
training data points.. We can empirically verify the connection through simple
examples as shown in Figure 1, where the changes of training data distribution
would result in larger KSD compared to that of the original training data
distribution. Intuitively, the connection shows that there is a tie between a
model and its training data points, encoded in the form of a stein kernel
function $k_{\theta}(\cdot,\cdot)$ defined on each pair of data points. As the
kernel function is conditioned on model $f_{\theta}$, we note it is an
encoding of data correlation under the context of a trained model, which paves
the foundation of the example-based prediction explanation as shown next.

###  3.1 KSD between Model and Training Data

Recall that KSD, $\mathbb{S}(P_{D},P_{\theta})$, defines the correlation
between pairs of training samples through model $\theta$ dependent kernel
function with closed-form decomposition

$\displaystyle\kappa_{\theta}((\mathbf{x}_{a},y_{a}),(\mathbf{x}_{b},y_{b}))=\mathcal{A}_{\theta}^{a}\mathcal{A}_{\theta}^{b}k(a,b)$
(1) $\displaystyle=\nabla_{a}\nabla_{b}k(a,b)+k(a,b)\nabla_{a}\log
P_{\theta}(a)\nabla_{b}\log P_{\theta}(b)$
$\displaystyle+\nabla_{a}k(a,b)\nabla_{b}\log
P_{\theta}(b)+\nabla_{b}k(a,b)\nabla_{a}\log P_{\theta}(a),$

where we denote data point $(\mathbf{x}_{a},y_{a})$ with $a$ for clean
notation. The only model-dependent factor in the above decomposition is a
derivative $\nabla_{\mathbf{x},y}\log P_{\theta}(\mathbf{x},y)$ (for both data
$a$ or $b$). Hence, to utilize KSD for prediction explanation, we first need
to unravel the derivative of the trained model with respect to a whole data
point (including its label). Indeed, while the gradient of a classifier
$\nabla_{\mathbf{x}}\log P_{\theta}(y|\mathbf{x})$ on its input feature
$\mathbf{x}$ is straightforward, the gradient to the entire data points is
less evident due to the lack of definition of $P_{\theta}(\mathbf{x})$.

To obtain the joint distribution of $P_{\theta}(\mathbf{x},y)$ for
facilitating KSD computation, we propose to set $P_{\theta}(\mathbf{x})\equiv
P_{D}(\mathbf{x})$ as uniform distribution. Although this setting appears
hasty, we note $P_{D}$ represents the uniform data distribution in dataset
$D$, which, while reflects, but is not the complex unknown distribution from
which the data is sampled. In particular, for a data point $(\mathbf{x},y)$ in
dataset D with i.i.d assumption, $P_{D}(\mathbf{x})$ is indeed uniformly
distributed, whereas $P_{D}(y|\mathbf{x})$ is a delta distribution with
probability $1$.

With the above setting, the score function $\nabla_{\mathbf{x},y}\log
P_{\theta}(\mathbf{x},y)$ in the Stein operator $\mathcal{A}_{\theta}$ could
be derived as a concatenation of the gradient of model
$f_{\theta}(\mathbf{x})_{y}$ to its input $\mathbf{x}$ and its probabilistic
prediction $f_{\theta}(\mathbf{x})$, since

$\displaystyle\nabla_{\mathbf{x},\mathbf{y}}\log
P_{\theta}(\mathbf{x},\mathbf{y})=\nabla_{\mathbf{x},\mathbf{y}}[\log
P_{\theta}(\mathbf{y}|\mathbf{x})+\log P_{d}(\mathbf{x})]$ (2)
$\displaystyle=\nabla_{\mathbf{x},\mathbf{y}}\log
P_{\theta}(\mathbf{y}|\mathbf{x})+[\nabla_{\mathbf{x}}\log
P_{d}(\mathbf{x})||\nabla_{\mathbf{y}}P_{d}(\mathbf{x})]$
$\displaystyle=[\nabla_{\mathbf{x}}\log\mathbf{y}^{\top}f_{\theta}(\mathbf{x})||\nabla_{\mathbf{y}}\log\mathbf{y}^{\top}f_{\theta}(\mathbf{x})]+[\mathbf{0}||0]$
$\displaystyle=[\nabla_{\mathbf{x}}\log
f_{\theta}(\mathbf{x})_{y}||f_{\theta}(\mathbf{x})],$

where $[\cdot||\cdot]$ denotes concatenation operation, and we use one-hot
encoding $\mathbf{y}$ of $y$ to clarify the gradient computation.

Combining Equation 1 and 2, we can estimate correlation of any pairs of
training data points conditioned on the trained machine learning model.
Computationally, since a score function $\nabla_{\mathbf{x},y}\log
P_{\theta}(\mathbf{x},y)$ depends on single data point, its outputs of the
training set could be pre-computed and cached for accelerate the kernel
computation. In particular, the output dimension of the score function is
simply $m+k$ for data with $m$ dimensional features and $k$ class labels.
Compared to the existing solutions, whose training data cache (or influence)
are bounded by the dimension of model parameters (such as Influence function,
TracIn, RPS, RPS-JLE), the explanation method built on KSD would come with a
significant advantage in terms of scalability. We note this statement is
generally true for neural network based classifiers, whenever the size of
model parameters is far larger than the data dimension.

Figure 2: Demonstration of HD-Explain on 2D Rectangular synthetic dataset.
Top-left shows the training dataset with three classes. Top-right shows the
explanation support of training data points to a given test point (as black
cross), where green shows higher KSD kernel value. The bottom shows the
distribution of KSD kernel value to the test point, where only a small number
of training data points provide strong support to this prediction.

###  3.2 Prediction Explanation via KSD

Now, we present how to extend the previously described kernel function to
support the prediction explanation. As stated in Equation 1, computing the
kernel function requires to access features and labels of a data point. While
the ground-truth label information is provided in the training set, it is
inaccessible during inference time given a test data point. Hence, to utilize
the KSD kernel function, we place predicted class $\hat{y}_{t}$ of the test
point $\mathbf{x}_{t}$ as label to construct a complete data point
$(\mathbf{x}_{t},\hat{y}_{t})$.

At this stage, searching supportive training examples for a given test point
$\mathbf{x}_{t}$ could be reduced to look for top-k training data points that
maximize the KSD defined kernel. For top-1 explanation sample, in particular,
we look for

$\operatorname*{argmax}_{i}\kappa_{\theta}([\mathbf{x}_{t}||\hat{\mathbf{y}}_{t}],[\mathbf{x}_{i}||\mathbf{y}_{i}]),$
(3)

where, again, we use one-hot encoding for data labels.

Figure 2 demonstrates a simple example on 2d synthetic dataset. By iterating
over all training data points and evaluating
$\kappa_{\theta}(\mathbf{d},\cdot)$ , we highlight that only a small number of
training data points provide strong influence on a particular prediction (see
Figure 2 (bottom)).

(a) Highly-Confident Correct Prediction

(b) Low-Confident Correct Prediction

(c) Low-Confident Incorrect Prediction

Figure 3: Qualitative explanation comparison among candidate example-based
explanation methods. We show three scenarios where the target model makes 1) a
highly-confident prediction that matches ground truth label, 2) a low-
confident prediction that matches ground truth label, 3) low-confident
prediction that does not match ground truth label (which is a bird). For each
sub plot, we show top-3 influential training data points picked by the
explanation methods for the test example.

(a) Overian Cancer Histopathology

(b) Brain Tumor MRI

Figure 4: Qualitative explanation comparison among candidate example-based
explanation methods on Overian Cancer histopathology and Brain Tumor MRI
datasets. We show two test data points that are predicted to belong to the
same class in each dataset. Red triangle in the top right corner of an image
shows the duplicate explanations across test samples.

##  4 Evaluations and Analysis

In this section, we conduct several qualitative and quantitative experiments
to demonstrate various properties of the proposed HD-Explain solution and
analyze it with existing example-based solutions to understand where and how
the approach can be used in practice.

Datasets: We are particularly interested in classification tasks that are
highly sensitive to trustworthiness and with a clear demand for prediction
explanations. To this end, we designed our experiments to include multiple
disease classification datasets, where diagnosis explanation is highly
desired. We also introduced synthetic and benchmark classification datasets to
deliver the main idea without the need for medical background knowledge.
Concretely, the primary datasets used in our experiments include CIFAR-10
($32\times 32\times 3$), Brain Tumor (Magnetic Resonance Imaging, $128\times
128\times 3$), and Ovarian Cancer (Histopathology Images, $128\times 128\times
3$). More details of the datasets are listed in the Appendix.

Baselines: The baseline explainers used in our experiments include Influence
Function, Representer Point Selection, and TracIn. While other variants of
these baseline explainers exist (barshan2020relatif; sui2021representer;
chen2021hydra), we note they don’t offer fundamental performance improvements
over the classic ones. In addition, as Influence Function and TracIn face
scalability issues, we limit the influence of parameters to the last layer of
the model so that they can work with models that contain a large number of
parameters. Our experiments use ResNet-18 as the backbone model architecture
(with around $11$ million trainable parameters) for all image datasets.
Finally, we also introduce an HD-Explain variant (HD-Explain*) to match the
last layer setting of other baseline models, even though HD-Explain can scale
up to the whole model without computation pressure (see Appendix F).

Metrics: In existing example-based explanation works, the experimental results
are often demonstrated qualitatively, as visualized explanation instances,
without quantitative evaluation, which, we believe, is subjective and hardly
convincing. To obtain a comprehensive understanding of the effectiveness of
explanation methods, we proposed several quantitative metrics in this paper.

  * •

Hit Rate: Hit rate shows how likely an explanation sample hits the desired
example cases. Here, the desired examples are guaranteed to be undisputed.
Specifically, we modify a training data point with minor augmentations (adding
noise or flipping horizontally) and use it as a test data point, such that the
best explanation for the generated test data point should be the original data
point in the training set.

  * •

Coverage: Given n test data points, we show how many unique explanation
samples an explanation method can produce when configuring to return top-k
training samples. Formally,

$\text{Coverage}=\frac{|\cup_{i=1}^{n}e_{i}|}{n\times k}$

where $e_{i}$ is the set of top-k explanations for test data point $i$. Larger
coverage is preferred as it suggests the higher granularity (per test point)
of the explanation.

  * •

Run Time: Measuring the run time of explanation methods. We measure the wall
clock time.

###  4.1 Qualitative Evaluation

Following the convention of previous contributions in the literature, we show
the prediction explanation instances for all candidate approaches side by
side. Figure 3 shows three test cases of the CIFAR10 classification task that
cover different classification outcomes, including high-confident correct
prediction, low-confident correct prediction, and low-confident incorrect
prediction. For both correct prediction cases, we are confident that HD-
Explain provides a better explanation than others in terms of visually
matching test data points. E.g. brown frogs in Figure 3 (a) and deer on the
grass in Figure 3 (b). In contrast, for the misclassified prediction case
(Figure 3 (c)), we note the HD-Explain produces examples that do not even
belong to the same class as the predicted one. It is a good outcome as the
method is not confident in providing consistent explanations. RPS also shows
such inconsistency, which aligns with its claim (yeh2018representer). The
other two baseline methods do not offer such properties and still produce
explanations that match the predicted label well. It is hard to justify how
those samples support such prediction visually (since no clear shared pattern
is obvious to us). In addition, it is interesting to see that Influence
Function and TracIn produce near identical explanations, reflecting their
similarity in leveraging the perturbation of model parameters.

Figure 4 provides additional insights into Ovarian Cancer histopathology and
Brain Tumor MRI datasets. HD-Explain again shows a better explanation for
producing training samples that appear similar to the test samples (note, for
the semantic similarity should be referred to a medical practitioner). E.g.
respecting the scanning orientation of test points in MRI as shown in Figure 4
(b). We note all baseline approaches tend to produce similar explanations to
test samples belonging to the same classes. Rather than providing individual
prediction explanations, those approaches act closer to per-class interpreters
that look for class prototypes. To verify this observation further, we
conducted a quantitative evaluation, as described in the next section.

(a) Hit Rate

(b) Coverage

(c) Execution Time

Figure 5: Quantitative explanation comparison among candidate example-based
explanation methods. Data augmentation strategy used is Noise Injection. Error
bar shows 95% confidence interval.

(a) Hit Rate

(b) Coverage

(c) Execution Time

Figure 6: Quantitative explanation comparison among candidate example-based
explanation methods. Data augmentation strategy used is Horizontal Flip. Error
bar shows 95% confidence interval. We reuse the legend of Figure 5.

###  4.2 Quantitative Evaluation

Now, we aim to evaluate the performance of the candidate explanation method
through quantitative analysis. As explainability is often subjective, we limit
our experiments to training datasets where ground-truth explanation samples
are available to obtain fair evaluation results. Specifically, given a
training data sample $(\mathbf{x}_{i},y_{i})$, we generate a test point
$\mathbf{x}_{t}$ by adopting two image data augmentation methods:

  * •

Noise Injection:
$\mathbf{x}_{t}=\mathbf{x}_{i}+\epsilon\quad\text{s.t.}\quad\epsilon\sim\mathcal{N}(0,0.01\sigma)$,
where $\sigma$ is the element-wise standard deviation of features in the
entire training dataset.

  * •

Horizontal Flip: $\mathbf{x}_{t}=\text{flip}(\mathbf{x}_{i})$, where we flip
images horizontally that do not compromise the semantic meaning of images for
the three datasets.

Since the data augmentation is guaranteed to maintain prediction consistency,
the ideally best explanation for the generated test point is the original data
point $\mathbf{x}_{i}$ itself, considering their similarity. Hence, the
quantitative evaluation could be a sample retrieval evaluation where Hit Rate
measures the probability of successful retrieval. To provide reliable
statistics, we created $30$ augmented test points for each training data point
($>10,000$ data points) in each dataset, resulting in more than $300,000$
independent runs.

Figure 5(a) shows the hit rate comparison among candidate methods on the three
image classification datasets under Noise Injection data augmentation. It was
surprising to observe that existing methods face significant difficulty in
retrieving the idea explanatory sample ($\leq 10\%$), even with such a simple
problem setup; only HD-Explain (and its variant) produces a reasonable
successful rate ($>80\%$). This observation raised our curiosity to know the
underlying reason. Figure 5(b) shows the following experiments where we
compute Coverage of unique explanations. The metric, coverage, shows the ratio
of explanation samples that are unique over many test points. The purpose is
to understand how likely a user may receive instance-specific explanations if
they have many test points to interpret. It turns out existing solutions
produce only 10% - 50% coverage – many test points receive the same set of
explanations, disregarding their unique characteristic. This result reflects
our previous observation in Figure 4, where we saw multiple duplicated
explanations (as shown with red corner tags). In fact, the explanations are
often dominated by the class labels; data points predicted as the same class
would receive a similar set of explanations. In contrast, HD-Explain shows
substantially higher coverage, providing a significant advantage for instance-
level prediction explanation.

Regarding computation efficiency, while we have summarized the scalability
limitation of the candidate methods in Table 1, there was no computational
efficiency evaluation conducted in previous works. To complement this lack, we
recorded the wall clock execution time for previous quantitative experiments
as shown in Figure 5(c). As expected, the Influence Function takes longer to
return its explanation than other methods. HD-Explain and TracIn* are
comparable. However, it is worth knowing that HD-Explain works on the whole
model prediction, whereas TracIn* is configured only to compute the gradient
of the predictive layer (due to the scalability issue of memory usage).
Indeed, if we look at HD-Explain*, we note its computation time is much
shorter than TracIn*. RPS shows a significant advantage in computation
efficiency as it is 1) limited to the last layer of the neural network and 2)
no auto-differentiation is needed for computing the training data influence.

We observe a similar trend in the other data augmentation scenario, Horizontal
Flip, where computation time and coverage are roughly the same, as shown in
Figure 6. However, we do notice that, as the outcome of image flipping, the
raw feature (pixel) level similarity between $\mathbf{x}_{t}$ and
$\mathbf{x}_{i}$ is destroyed. As an outcome, the HD-Explain that works on raw
features suffers from performance deduction while other methods, including HD-
Explain*, are less impacted. This observation suggests that choosing the layer
of explanation might be considered in the practical usage of this approach.

(a) Hit Rate (Noise Injection)

(b) Coverage (both)

(c) Hit Rate (Horizontal Flip)

(d) Execution Time (both)

Figure 7: Quantitative explanation comparison among HD-Explainers with
different kernel functions on all image classification datasets. Error bar
shows 95% confidence interval.

###  4.3 Kernel Options

Our previous experiments used the Radial Basis Function (RBF) kernel by
default. One may suggest another choice over kernels that better fit a
particular application domain. In this experiment, we compare three well-known
kernels on the three image classification datasets, including Linear, RBF, and
Inverse Multi-Quadric (IMQ) kernels.

Figure 7 shows our experiment results under both data augmentation scenarios.
We note that HD-Explain equipped with an IMQ kernel generally performs better
than the RBF kernel regarding explanation quality (Hit Rate). The advantage is
significant when the data augmentation scenario is Horizontal Flip (Figure 7
(c)), which appears more challenging than Noise Injection. In addition, the
Coverage of explanation (Figure 7 (b)) aligns with the hit rate, where IMQ
leads to a higher value. However, HD-Explain equipped with linear kernel runs
much faster than the other options, as shown in Figure 7 (d), providing it
with higher efficiency on large datasets. Indeed, if we cross-check the
results with Figure 5, we note that even a Linear kernel is sufficient for HD-
Explain to stand out from other baseline models considering its performance
and efficiency. RBF and IMQ kernels often need similar computation time.

(a) Recall@K

(b) Precision@K

Figure 8: Data debugging comparison among candidate methods on CIFAR-10
dataset. Results collected from 30 independent runs. Error bar shows 95%
confidence interval.

###  4.4 Data Debugging

Now, we describe our data debugging experiments to highlight the self-
explanatory ability of candidate methods on the training data. In particular,
we generalized previous research’s binary classification-based data debugging
experiment into a multi-classification scenario, where we randomly flip labels
of 100 training data points at each run. We adopt standard information
retrieval metrics, Precision and Recall, that measure how likely the candidate
methods can retrieve the mislabeled training data points.

Figure 8 shows our experimental results. While HD-Explain on the entire model
has little data debugging ability, its variant on the last layer offers
outstanding performance compared to the other last-layer explanation methods.
Note, the data debugging functionality is a side effect/benefit of HD-Explain,
which is not our main proposal

##  5 Conclusion

This paper presents HD-Explain, a Kernel Stein Discrepancy-driven neural
classifier explanation method. The method is straightforward but powerful in
terms of prediction explanation and data debugging tasks. The computation is
efficient and can work with large machine learning models (e.g. ResNet) as a
whole model explainer. With empirical evidence demonstrated in this paper, we
are confident that the method is practically effective for improving the
transparency of Machine Learning predictive systems.

## Impact Statements

This paper presents a prediction explanation method for neural network-based
classifiers. While the ultimate goal of prediction explanation (and this work)
is to improve the model transparency and trustworthiness, we note it could be
misused to cause intentional prediction bias (e.g. removing high-influential
training data points for specific purposes).

##  Appendix A Problem Definition Recap

We considers the task of explaining the prediction of a differentiable
classifier $f:{\rm I\\!R}^{d}\rightarrow{\rm I\\!R}^{l}$, given inputs test
sample $\mathbf{x}_{t}\in{\rm I\\!R}^{d}$, where $d$ denotes the input
dimension and $l$ denotes the number of classes. Specifically, we are
interested in explaining a prediction of a model $f(\cdot)$ by returning a
subset of its training samples $D=\\{(\mathbf{x}_{i},y_{i})\\}^{n}_{i=1}$ that
has strong predictive support to the prediction of test point
$\mathbf{x}_{t}$. While not explicitly stated in the main paper, we treat
example-based prediction explanation as a function
$\psi({f,D,\mathbf{x}_{t}}):\mathcal{F}\times\mathcal{D}\times{\rm
I\\!R}^{d}\rightarrow\\{{\rm I\\!R}^{d},{\rm I\\!R}^{l}\\}^{k}$ such that it
takes a trained model $f$, a training dataset $D$, and an arbitrary test point
$\mathbf{x}_{t}$ as inputs and output top-$k$ training samples as
explanations.

##  Appendix B Additional Derivation of Kernelized Stein Discrepancy

While Stein’s Identity has been well described in many previous works
(liu2016kernelized; liu2016stein; chwialkowski2016kernel), we briefly recap
some key derivations in this paper to seek for self-contained.

As mentioned in the main paper, Stein’s Identity states that, if a smooth
distribution $p(x)$ and a function $\phi(x)$ satisfy
$\lim_{||x||\rightarrow\infty}p(x)\phi(x)=0$, we have

$\mathbb{E}_{x\sim p}[\phi(x)\nabla_{x}\log
p(x)+\nabla_{x}\phi(x)]=\mathbb{E}_{x\sim
p}[\mathcal{A}_{p}\phi(x)]=0,\quad\forall f.$

Intuitively, by using integration by part rules, we can reveal the original
assumption from the derived expression such that

$\int_{x}\phi(x)\nabla_{x}\log
p(x)+\nabla_{x}\phi(x)dx=p(x)\phi(x)\Bigr{|}^{+\infty}_{-\infty}$

Stein Discrepancy measures the difference between two distributions $q$ and
$p$ by replacing the expectation of distribution $p$ term in Stein’s Identity
expression with distribution $q$, which reveals the difference between two
distributions by projecting their score functions (gradients) with the
function $\phi(x)$

$\max_{\phi\in\mathcal{F}}\mathbb{E}_{x\sim{\color[rgb]{1,0.6,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.6,0.2}q}}[\mathcal{A}_{p}\phi(x)]=\max_{\phi\in\mathcal{F}}\mathbb{E}_{x\sim{\color[rgb]{1,0.6,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.6,0.2}q}}[\mathcal{A}_{p}\phi(x)]-\underbrace{\mathbb{E}_{x\sim{\color[rgb]{1,0.6,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.6,0.2}q}}[\mathcal{A}_{\color[rgb]{1,0.6,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.6,0.2}q}\phi(x)]}_{=0}=\max_{\phi\in\mathcal{F}}\mathbb{E}_{x\sim{\color[rgb]{1,0.6,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.6,0.2}q}}[\underbrace{\phi(x)}_{\text{projection
coefficients}}\underbrace{(\nabla_{x}\log
p(x)-\nabla_{x}\log{\color[rgb]{1,0.6,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.6,0.2}q}(x))}_{\text{score
function difference}}]$

Clearly, the choice of projection coefficients (function $\phi(x)$) term is
critical to measure the distribution difference.

Kernelized Stein Discrepancy (KSD) addresses the task of searching function
$\phi$ by treating the above challenge as an optimization task where it
decompose the target function $\phi$ with linear decomposition such that

$\max_{\phi\in\mathcal{F}}\mathbb{E}_{x\sim{\color[rgb]{1,0.6,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.6,0.2}q}}[\mathcal{A}_{p}\phi(x)]=\max_{\phi\in\mathcal{F}}\mathbb{E}_{x\sim{\color[rgb]{1,0.6,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.6,0.2}q}}[\mathcal{A}_{p}\sum_{i}w_{i}\phi_{i}(x)]=\max_{\phi\in\mathcal{F}}\sum_{i}w_{i}\mathbb{E}_{x\sim{\color[rgb]{1,0.6,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.6,0.2}q}}[\mathcal{A}_{p}\phi_{i}(x)],$

with linear property of Stein operator $\mathcal{A}_{p}$. The linear
decomposition path the way to reduce the optimization task into looking for a
finite number of the base functions $\phi_{i}\in\mathcal{F}$ whose coefficient
norm is constraint to $1$ ($||\mathbf{w}||_{\mathcal{H}}\leq 1$). KSD takes
$\mathcal{F}$ to be the unit ball of a reproducing kernel Hilbert space (RKHS)
and leverages its reproducing property such that
$\phi(x)=\langle\phi(\cdot),k(x,\cdot)\rangle$, which in turn transforms the
maximization objective of the Stein Discrepancy into

$\max_{\phi}\langle\phi(\cdot),\mathbb{E}_{x\sim{\color[rgb]{1,0.6,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.6,0.2}q}}[\mathcal{A}_{p}k(\cdot,x)]\rangle_{\mathcal{H}},\quad\textit{s.t.}||\phi||_{\mathcal{H}}\leq
1.$

The optimal $\phi$ is therefore a normalized version of
$\mathbb{E}_{x\sim{\color[rgb]{1,0.6,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.6,0.2}q}}[\mathcal{A}_{p}k(\cdot,x)]$.
Hence, KSD is defined as the optimal between the distribution $p$ and $q$ with
the optimal solution of $\phi$

$\mathbb{S}(q,p)=\mathbb{E}_{x,x^{\prime}\sim{\color[rgb]{1,0.6,0.2}\definecolor[named]{pgfstrokecolor}{rgb}{1,0.6,0.2}q}}[\kappa_{p}(x,x^{\prime})],\quad\text{where}\quad\kappa_{p}(x,x^{\prime})=\mathcal{A}_{p}^{x}\mathcal{A}_{p}^{x^{\prime}}k(x,x^{\prime}).$

##  Appendix C Discussion: Intuition on why HD-Explain works

In HD-Explain, the key metric on measuring the predictive supports of a test
point $\mathbf{x}_{t}$ given a training data $(\mathbf{x}_{i},\mathbf{y}_{i})$
is the KSD defined kernel
$\kappa_{\theta}([\mathbf{x}_{t}||\hat{\mathbf{y}}_{t}],[\mathbf{x}_{i}||\mathbf{y}_{i}])$,
where $\hat{\mathbf{y}}_{t}$ denotes the predicted class label by model
$f_{\theta}$ in one-hot encoding. By definition, the kernel
$\kappa_{\theta}((\mathbf{x}_{a},y_{a}),(\mathbf{x}_{b},y_{b}))=k_{\theta}(a,b)$
between two data points can be decomposed into four terms

$\underbrace{\text{trace}(\nabla_{a}\nabla_{b}k(a,b))}_{\text{\char
172}}+\underbrace{k(a,b)\nabla_{a}\log P_{\theta}(a)^{\top}\nabla_{b}\log
P_{\theta}(b)}_{\text{\char
173}}+\underbrace{\nabla_{a}k(a,b)^{\top}\nabla_{b}\log
P_{\theta}(b)}_{\text{\char
174}}+\underbrace{\nabla_{b}k(a,b)^{\top}\nabla_{a}\log
P_{\theta}(a)}_{\text{\char 175}}.$

We exam the effect of each term as follows:

  * •

①: We note that the first term is often a similarity bias of raw data points
given a specified kernel function. In particular, for the RBF kernel
$k(a,b)=\exp(-\gamma||a-b||^{2})$, the first term is simply
$\sum_{i}^{d+l}2\gamma k(a,b)$, where $d+l$ refers to the sum of input and
output (in one-hot) dimensions of a data point. Intuitively, the term shows
how similar the two data point is given the RBF kernel. For linear kernel
$k(a,b)=a^{\top}b$, on another hand, the first term is simply $d+l$ as a
constant bias term, which does not deliver any similarity information between
the two data points.

  * •

②: The second term reflects the similarity between two data points in the
context of trained model. In particular, considering the sub-term
$\nabla_{a}\log P_{\theta}(a)^{\top}\nabla_{b}\log P_{\theta}(b)$, based our
derivation in Equation 2 (in the main paper), we note it is equivalent to

$[\nabla_{\mathbf{x}_{a}}\log
f_{\theta}(\mathbf{x}_{a})_{y_{a}}||f_{\theta}(\mathbf{x}_{a})]^{\top}[\nabla_{\mathbf{x}_{b}}\log
f_{\theta}(\mathbf{x}_{b})_{y_{b}}||f_{\theta}(\mathbf{x}_{b})]=\underbrace{\nabla_{\mathbf{x}_{a}}\log
f_{\theta}(\mathbf{x}_{a})_{y_{a}}^{\top}\nabla_{\mathbf{x}_{b}}\log
f_{\theta}(\mathbf{x}_{b})_{y_{b}}}_{\text{similarity of scores (input
gradients)}}+\underbrace{f_{\theta}(\mathbf{x}_{a})^{\top}f_{\theta}(\mathbf{x}_{b})}_{\text{similarity
of predictions}},$

where both terms could be viewed as similarity between data points in the
context of trained model.

  * •

③-④: Both of the last two terms exam the alignment between score of one data
point and the kernel derivative of another data point. Our conjecture is that
this alignment reflects how a test prediction would change if there is a
training data point presents closer to it than before.

##  Appendix D Dataset Details

Table 2: Summary of datasets used in the paper.

Dataset | Application | Type | # Size | # Feature Dimension | # Number of Classes | Duplicated Samples | Public Dataset
---|---|---|---|---|---|---|---
Two Moons | Synthetic | 2D Numeric | 500 | 2 | 2 | No | Shared with code
Rectangulars | Synthetic | 2D Numeric | 500 | 2 | 3 | No | Shared with code
CIFAR-10 | Classification Benchmark | Image | 60,000 | $32\times 32\times 3$ | 10 | No | Yes
Overian Cancer | Histopathology (Private) | Image | 20,000 | $128\times 128\times 3$ | 5 | Yes | No
Brain Tumor | MRI Benchmark | Image | 7,023 | $128\times 128\times 3$ | 4 | Yes | Yes

In this paper, we conducted our experiments on five datasets – two synthetic
and three benchmark image classification datasets. As the work concerns the
trustworthiness of the machine learning model in high-stakes applications, we
also introduced medical diagnosis datasets to provide more insight into the
potential benefit the proposed work introduced. To train the target machine
learning models, we conducted data augmentations to increase the number of
training data samples, including random cropping, rotation, shifting,
horizontal flipping, and noise injection. Table 2 summarizes more details
about the datasets.

##  Appendix E Experiment Setup for Data Debugging

Before describing the data debugging setting of this paper, we want to recap
that the data debugging functionality is a side effect/benefit of HD-Explain,
which is not our main proposal. Indeed, using the prediction explanation
method as a data debugging tool is still under investigation since it might be
over-claimed due to the over-regularized setting in previous works (e.g.
Binary classification tasks). While we relaxed some settings, we don’t claim
it practical for real-world applications.

The data debugging task in this paper is a data sample retrieval task where we
retrieve samples that intentionally flipped their classification labels.
Higher Precision and Recall of the retrieval reflects higher performance of
data debugging.

For the HD-Explain (and its variant HD-Explain*), the retrivial order is
determined by the values of the diagonals of the KSD-defined kernel matrix,
$\kappa_{\theta}(a,a)$ for all $a\in D$. This setting is very similar to how
the Influence function does the data debugging with the self-influence of a
data sample. Indeed, $\kappa_{\theta}(a,a)$ could be treated as a self-
influence that does not rely on model parameters.

##  Appendix F HD-Explain*

In the main paper, we briefly mentioned HD-Explain* as an HD-Explain variant
to match the last layer setting of other baseline models (e.g. RPS). The HD-
Explain* is a simple change of HD-Explain in terms of using data
representations (the output of the last non-predictive layer of the neural
classifier) rather than the raw features. Specifically, we assume a neural
network model $f_{\theta}$ could be decomposed into two components
$f_{\theta_{2}}\cdot f_{\theta_{1}}$, where $f_{\theta_{1}}$ is a
representation encoder and $f_{\theta_{2}}$ is a linear model for prediction.
With this decomposition, we define the KSD kernel function for HD-Explain* as

$\displaystyle\kappa_{\theta}((f_{\theta_{1}}(\mathbf{x}_{a}),y_{a}),(f_{\theta_{1}}(\mathbf{x}_{b}),y_{b}))=\mathcal{A}_{\theta_{2}}^{a}\mathcal{A}_{\theta_{2}}^{b}k(a,b)$
$\displaystyle=\nabla_{a}\nabla_{b}k(a,b)+k(a,b)\nabla_{a}\log
P_{\theta_{2}}(a)\nabla_{b}\log
P_{\theta_{2}}(b)+\nabla_{a}k(a,b)\nabla_{b}\log
P_{\theta_{2}}(b)+\nabla_{b}k(a,b)\nabla_{a}\log P_{\theta_{2}}(a),$

where we define $a=(f_{\theta_{1}}(\mathbf{x}_{a}),y_{a})$ and
$b=(f_{\theta_{1}}(\mathbf{x}_{b}),y_{b})$ for short. This setting reduces the
prediction explanation to the last layer of the neural network in a similar
fashion to RPS.

##  Appendix G Additional Experiment Results

(a) Influence Function

(b) Representer Point Selection

(c) TracIn

(d) HD-Explain

(e) HD-Explain*

(f) Influence Function

(g) Representer Point Selection

(h) TracIn

(i) HD-Explain

(j) HD-Explain*

Figure 10: Two additional explanation examples on CIFAR-10 dataset. The first
image on the left is the test data point and the remaining three images are
top-3 explanations. Similar to the samples shown in the main paper, we note
HD-Explain’s output is visually more similar to the task data point.