# On Minimizing the Impact of Dataset Shifts on Actionable Explanations

Anna P. Meyer ‡ University of Wisconsin - Madison. Email<EMAIL_ADDRESS>Dan Ley Harvard University. Emails: {dley<EMAIL_ADDRESS><EMAIL_ADDRESS>contribution  Suraj Srinivas† Himabindu Lakkaraju†

###### Abstract

The Right to Explanation is an important regulatory principle that allows
individuals to request actionable explanations for algorithmic decisions.
However, several technical challenges arise when providing such actionable
explanations in practice. For instance, models are periodically retrained to
handle dataset shifts. This process may invalidate some of the previously
prescribed explanations, thus rendering them unactionable. But, it is unclear
if and when such invalidations occur, and what factors determine explanation
stability i.e., if an explanation remains unchanged amidst model retraining
due to dataset shifts. In this paper, we address the aforementioned gaps and
provide one of the first theoretical and empirical characterizations of the
factors influencing explanation stability. To this end, we conduct rigorous
theoretical analysis to demonstrate that model curvature, weight decay
parameters while training, and the magnitude of the dataset shift are key
factors that determine the extent of explanation (in)stability. Extensive
experimentation with real-world datasets not only validates our theoretical
results, but also demonstrates that the aforementioned factors dramatically
impact the stability of explanations produced by various state-of-the-art
methods.

##  1 Introduction

Machine learning (ML) models have recently witnessed increased utility in
critical real-world applications. This, in turn, led to the introduction of
several regulatory principles that aim to safeguard the practice of
algorithmic decision making in such settings [13, 32]. The Right to
Explanation is one such important regulatory principle that allows individuals
to request actionable explanations for algorithmic decisions that adversely
impact them. Post-hoc explanation methods such as LIME, SHAP, input gradients,
and Smoothgrad [27, 37, 40, 42] have commonly been employed in practice to
operationalize this principle. These methods explain complex model predictions
by assigning importance scores to input features, typically via learning local
linear approximations of the underlying functions [18].

Providing actionable explanations to end-users is often hindered by
operational challenges in practice. For instance, models are periodically
retrained to handle dataset shifts, and the original explanations may no
longer be valid under the new model. If so, the original explanations do not
remain actionable. For example, a user may be informed that their low salary
was the primary reason for the rejection of their loan application, thus
prompting them to increase their income. However, an update to the model could
result in a shift of model internals such that the user’s credit score, and
not their income, is now the primary factor for the rejection. In this case,
the user’s action of increasing their income would be less likely to result in
a positive outcome once the updated model was deployed.

The aforementioned scenario could have been avoided if the underlying model
and its explanations remained (relatively) unchanged after model retraining
due to dataset shifts. Such _explanation stability_ implies actionability
because if explanations remain (relatively) unchanged despite model
retraining, they are likely to be actionable over an extended period of time.
Therefore, it is necessary to develop a systematic understanding of
explanation stability in the face of dataset shifts, characterize the
conditions that lead to unstable explanations, and find ways to mitigate
explanation instability after retraining. Despite its significance, there is
little to no research on characterizing the factors that influence explanation
stability.

Our work addresses the aforementioned critical gaps by providing one of the
first theoretical and empirical characterizations of the factors influencing
explanation stability after dataset shifts. To this end, we conduct rigorous
theoretical analysis to demonstrate that model curvature, weight decay
parameters while training, and the magnitude of the dataset shift are key
factors that determine the extent of explanation (in)stability. Our
theoretical analysis emphasizes how seemingly inconsequential modeling
decisions that may not impact predictive accuracy can heavily influence other
critical aspects such as explanation stability. We also conduct extensive
experimentation with real-world datasets to validate our theoretical insights
with various state-of-the-art explanation methods (e.g., gradient and
perturbation based methods such as LIME, SHAP, SmoothGrad etc.) and
explanation stability metrics (e.g., $\ell_{2}$, top-k consistency). We also
empirically analyze how other training decisions (e.g., learning rate and
number of training epochs) impact explanation stability. In summary, our work
makes the following key contributions:

  1. 1.

We provide one of the first theoretical and empirical characterizations of the
factors influencing explanation stability in the face of dataset shifts.

  2. 2.

We conduct a rigorous theoretical analysis to demonstrate that model
curvature, weight decay parameters during training, and the magnitude of the
dataset shift are key factors influencing the extent of explanation
(in)stability (§3).

  3. 3.

We carry out extensive experimentation with multiple real-world datasets to
validate our theoretical results (§4.1). Our empirical findings suggest that
standard neural network training pipelines exhibit low explanation stability,
and confirm that our theoretical results are valid even when some of the
underlying assumptions do not hold in practice (§4.2).

  4. 4.

We empirically analyze the impact of other training hyperparameters such as
learning rate, batch size, and number of training epochs on explanation
stability (§4.3).

##  2 Related Work

##### Explanation Methods

A variety of post-hoc techniques have been proposed to explain complex models
[10, 23, 38]. These techniques differ in their access to the complex model
(i.e., black box vs. access to internals), scope of approximation (e.g.,
global vs. local), search technique (e.g., perturbation-based vs. gradient-
based), explanation families (e.g., linear vs. non-linear), etc. For instance,
LIME [37] and SHAP [27] are _perturbation-based_ , _local explanation_
approaches that learn a linear model locally around each prediction. Other
_local explanation_ methods capture feature importances by computing the
gradient with respect to the input [39, 40, 42, 44]. Counterfactual
explanation methods, on the other hand, capture the changes that need to be
made to a given instance in order to flip its prediction [21, 22, 26, 34, 46,
47, 48]. In this work, we focus on analyzing the stability of explanations
output by perturbation-based and gradient-based local explanation methods.

##### Explanation Stability and Robustness

Various notions of explanation stability and robustness have been suggested in
the literature.111The literature is inconsistent on the distinction between
explanation robustness and stability; for clarity we use “robustness” when the
model is fixed (e.g., robustness to input perturbations) and “stability” when
the model is changing. One line of work shows that explanations are not robust
by crafting adversarial manipulations of test inputs that change their
explanations but not predictions [8, 17]. Another set of approaches studies
explanation stability by altering the underlying model to maintain accuracy,
but change explanations to a desired target [1, 20, 41]. Multiple works
suggest using smooth or low-curvature models to improve both model and
explanation stability [8, 43]. However, all of these works focus on
adversarial manipulations of either the model or the input, while we are
interested in more realistic model shifts occurring due to natural dataset
shifts. In contrast, other works have proposed algorithmic techniques to
generate counterfactual explanations or algorithmic recourses that are stable
under model shifts [6, 12, 16, 30, 31, 36, 45]. However, these works differ
from ours in three ways: first, we focus on feature attribution-based
explanations that highlight influential features, while the existing work
focuses on counterfactual explanations or algorithmic recourse. Second, we
adapt the model to yield more stable explanations after retraining due to
dataset shifts, whereas most of the existing work aims to adapt the
explanation method. And finally, we theoretically characterize the factors
that impact explanation stability.

More generally, neural networks are known to be non-robust to small training
modifications, e.g., different random initializations or alternate model
selection due to underspecification [2, 7, 25, 28]. However, relatively few
works have studied how _explanations_ change in the same setting. These works
only show empirically that explanations are not stable with respect to
underspecification [5] or random model initializations [3]. These works are
orthogonal to ours, as we focus on theoretically and empirically
characterizing specific properties of models and datasets that impact
explanation stability.

##### Algorithmic Stability

Algorithmic stability is a classical learning-theoretic framework that
characterizes the consistency of outputs of learning algorithms when trained
on different but similarly distributed data [4]. In contrast, our study does
not make the assumption that the shifted data is similarly distributed. [19]
studied the hypothesis stability of stochastic gradient descent (SGD) and
showed that practical choices such as learning rate and number of epochs of
training can impact the stability of model parameters. Recent work has applied
the algorithmic stability perspective to post-hoc explanations [14] by
requiring that feature attributions be stable to a re-sampling of the
underlying dataset. However, none of these prior works study realistic model
shifts occurring due to natural dataset shifts or characterize the factors
influencing explanation stability, which is the main focus of our work.

##  3 Theoretical Analysis

In this section, we provide an analytical characterization of explanation
shifts $e_{1}\rightarrow e_{2}$ and model shifts $f_{1}\rightarrow f_{2}$ that
happen due to underlying dataset shifts
$\mathcal{D}_{1}\rightarrow\mathcal{D}_{2}$. Concretely, we address the
following question:  How much do the explanations for a model $f_{1}$ trained
on dataset $\mathcal{D}_{1}$ change when fine-tuning222 Here, fine-tuning
refers to initializing model parameters of $f_{2}$ with those of $f_{1}$, and
then training $f_{2}$ using $D_{2}$ until its loss converges. We provide
practical details of this in §4.1.on a slightly shifted dataset
$\mathcal{D}_{2}$, resulting in a new model $f_{2}$?

Let the dataset $\mathcal{D}_{1}$ be composed of input-output pairs
($\mathbf{x},y$), where $\mathbf{x}\in\mathbb{R}^{d}$, and $y$ is either a
real value or a one-hot vector depending on the task (regression or
classification). For simplicity of analysis, in this section we shall focus on
gradient explanations $e=\nabla_{\mathbf{x}}f\in\mathbb{R}^{d}$, where $f$ is
a scalar-valued neural network function
$f(\cdot,\theta):\mathbb{R}^{d}\rightarrow\mathbb{R}$ mapping $D$-dimensional
inputs to scalar outputs. Let
$\ell(\cdot,\theta):\mathbb{R}^{d}\rightarrow\mathbb{R}$ denote $f$ composed
with a loss function such as cross-entropy loss. Here, $\theta$ represents the
parameters of the function being minimized, e.g., a deep neural network.

Our goal here is to relate the explanation shift to the dataset shift (i.e.,
$d(\mathcal{D}_{1},\mathcal{D}_{2})$, for some distance measure $d$) and the
learning algorithm. To this end, our overall strategy is:

  1. 1.

In §3.1, we bound the _parameter shift_ $\|\theta_{2}-\theta_{1}\|_{2}$ in
terms of the _dataset shift_ $d(\mathcal{D}_{1},\mathcal{D}_{2})$.

  2. 2.

In §3.2, we bound the _explanation shift_
$\|\nabla_{\mathbf{x}}f(\mathbf{x},\theta_{2})-\nabla_{\mathbf{x}}f(\mathbf{x},\theta_{1})\|_{2}$
in terms of the _parameter shift_ $\|\theta_{2}-\theta_{1}\|_{2}$ derived
above.

While the existence of these shifts are qualitatively evident, we make novel
contributions in that (a) we quantify the parameter and gradient changes and
(b) we show that these changes are affected by factors in the modelling
process.

###  3.1 Bounding the Parameter Shift

We begin by defining the distance between datasets
$d(\mathcal{D}_{1},\mathcal{D}_{2})$ using Hungarian distance, as defined in
Defn. 1. The idea behind the Hungarian distance is simple: we find the
”alignment” between data points of two datasets that minimizes the average
$\ell_{2}$ distance between a point in $\mathcal{D}_{1}$ and its counterpart
in $\mathcal{D}_{2}$.

######  Definition 1.

Given datasets $\mathcal{D}_{1}=\\{(x_{i},y_{i})_{i=1,\ldots,N}\\}\ $ and
$\mathcal{D}_{2}=\\{(x^{\prime}_{i},y^{\prime}_{i})_{i=1,\ldots,N}\\}$, their
_Hungarian distance_ is

$\displaystyle
d(\mathcal{D}_{1},\mathcal{D}_{2})=\min_{P(\mathcal{D}_{2})}\sum_{i=1}^{N}\|x_{i}-x^{\prime}_{i}\|_{2}$

where the minimum is taken over $P(\mathcal{D}_{2})$, all permutations on the
data points of $\mathcal{D}_{2}$.

Using this distance metric, we derive an expression for the parameter change
$\|\theta_{2}-\theta_{1}\|_{2}$ in terms of the dataset distance
$d(\mathcal{D}_{1},\mathcal{D}_{2})$. We use the following assumptions to
derive this result: (1) we minimize a regularized loss $\ell_{reg}$ that
involves weight decay, i.e,
$\ell_{reg}(\theta)=\ell(\theta)+\gamma\|\theta\|^{2}_{2}$; (2) $\ell$ is
locally quadratic; (3) the learning algorithm returns a unique minimum
$\theta$ given a dataset $\mathcal{D}$.

######  Theorem 1.

Given the assumptions stated above, and that $\mathcal{L}_{x}(\theta_{1})$ is
the Lipschitz constant of the model with parameters $\theta_{1}$, we have

$\displaystyle\|\theta_{2}-\theta_{1}\|_{2}$
$\displaystyle\leq\sqrt{\frac{\mathcal{L}_{x}(\theta_{1})d(\mathcal{D}_{1},\mathcal{D}_{2})}{\gamma}}+C$

where $\gamma$ is the weight decay regularization constant, and $C$ is a small
problem-dependent constant.

###### Proof Idea.

The idea of the proof is to first estimate
$\ell_{\mathcal{D}_{2}}(\theta_{1})$, i.e., the initial loss of the model
$\theta_{1}$, before fine-tuning on $\mathcal{D}_{2}$. Assuming, for the sake
of explanation, that the optimal loss value obtained by $\theta_{2}$ on
$\mathcal{D}_{2}$ is zero, we have obtained the change in loss value from
$\theta_{1}\rightarrow\theta_{2}$. To derive the change in parameters from the
change in loss derived above, we use a second order Taylor series expansion
(with the assumption that the loss is locally quadratic). In particular, we
lower bound the Hessian with its lowest eigenvalue using the weight decay term
$\gamma$. The “problem-dependent constant” $C$ arises from deviating from the
assumption that the optimal loss value is zero. The complete proof is provided
in the supplementary material. ∎

The above theorem achieves our first goal of relating parameter change to
dataset shift. Intuitively, this tells us that parameter shift depends
directly on dataset shift (as expected), directly on the Lipschitz constant
(which depends on the model’s robustness), and inversely on the weight decay
parameter.

###  3.2 Bounding the Explanation Shift

We now consider the problem of bounding the explanation shift
$\|\nabla_{\mathbf{x}}f(\cdot,\theta_{2})-\nabla_{\mathbf{x}}f(\cdot,\theta_{1})\|_{2}$
given an estimate of the parameter shift $\|\theta_{2}-\theta_{1}\|_{2}$. To
this end, we first define a quantity called the gradient-parameter Lipschitz
constant, defined in Defn. 2. To the best of our knowledge, this quantity has
not been considered in previous works.

######  Definition 2.

We define the “gradient-parameter Lipschitz constant”
$\mathcal{L}_{\Theta,\mathcal{D}}$ w.r.t. an input distribution $\mathcal{D}$,
and parameter set $\Theta$ as follows:

$\displaystyle\mathcal{L}_{\Theta,\mathcal{D}}=\mathbb{E}_{\theta\in\Theta}\mathbb{E}_{\mathbf{x}\sim\mathcal{D}}\|\nabla_{\theta}\nabla_{x}f(x;\theta)\|_{2}$

Intuitively, this quantity captures the sensitivity of gradient explanations
to small changes in the parameter $\theta$. In contrast to usual definitions
of Lipschitz constants, this is defined (1) locally with respect a particular
distribution and parameter set, and (2) using the mean instead of the $\sup$.
This allows us to derive the following relationship between average gradient
difference and parameter shift.

######  Lemma 1.

The gradient-parameter Lipschitz has the following property:

$\displaystyle\mathbb{E}_{\mathbf{x}\sim\mathcal{D}}\|\nabla_{x}f(x;\theta_{1})-\nabla_{x}f(x;\theta_{2})\|_{2}\leq\leavevmode\nobreak\
\mathcal{L}_{\Theta,\mathcal{D}}\times\|\theta_{2}-\theta_{1}\|_{2}$

where
$\Theta=\\{\lambda\theta_{1}+(1-\lambda)\theta_{2}\mid\lambda\in[0,1]\\}$.

This relation follows immediately from the fundamental theorem of integral
calculus, and the full proof is given in the supplementary material. This
result reveals an intuitive fact that larger parameter shifts lead to larger
gradient shifts, but mediated by the gradient-parameter Lipschitz. This
quantity is difficult to analyse for general neural networks, and we provide
below an analysis for the case of 1-hidden layer neural networks, with a
specific simplifying assumption on the data distribution.

######  Theorem 2.

Assume that we have a 1-hidden layer neural network with weights $\theta$, and
random inputs $\mathbf{x}\sim\mathcal{N}(0,I)$333Covariance of I is chosen for
notational brevity. Further assume that we use an activation function $\sigma$
with well-defined second derivatives (e.g: softplus). For this case, the
gradient-parameter Lipschitz constant is

$\displaystyle\mathbb{E}_{\mathbf{x},\theta}\|\nabla_{\theta}\nabla_{\mathbf{x}}f(\mathbf{x},\theta)\|_{2}\leq\left(\mathbb{E}_{\theta}\|\theta\|_{2}\right)+\beta\left(\mathbb{E}_{\theta}\phi(\theta)\right)$

where $\beta$ is the maximum curvature of activation function $\sigma$, and
$\phi(\theta)$ is the path-norm [29] of the model.

###### Proof Idea.

The proof involves computing the gradient-parameter derivatives for the
1-hidden layer neural network case. Our assumption regarding smoothness of
activation function ensures that $\sigma^{\prime\prime}\leq\beta$, i.e., that
the second derivative of the activation function is upper bounded, which is
true for smooth non-linearities like softplus. Another simplifying assumption
regarding the distribution of the inputs ($\mathbf{x}\sim\mathcal{N}(0,I)$)
helps us dramatically simplify the expression for the expected value of the
gradient norm terms.

The path-norm for a 1-hidden layer neural network (for weights $W_{2},W_{1}$)
is given by $\phi(\theta)=\sqrt{\sum_{i,j}(W_{2}^{j}W_{1}^{i,j})^{2}}$ and has
been linked to model generalization [29]. The complete proof is provided in
the supplementary material. ∎

These results achieve our goal of relating the change in gradients to
parameter change. Taken together with the results in §3.1, this achieves our
overall goal of relating gradient change to dataset shift. Intuitively, this
tells us the following regarding gradient shift:

  1. 1.

It depends directly on the _dataset distance_
$d(\mathcal{D}_{2},\mathcal{D}_{1})$, which is as expected: the larger the
dataset shift, the larger can be the parameter shift, and explanation shift.

  2. 2.

It depends inversely on the _weight decay_ parameter $\gamma$, and directly on
the norm of weights $\|\theta\|_{2}$. Intuitively, large weight decay shifts
the optima to be closer to zero, reducing the norm of weights, which makes all
optima closer to each other.

  3. 3.

It depends directly on the _Lipschitz constant_ of the first model, which is
connected to the model’s robustness. Intuitively, if the model is already
robust to small input changes, we don’t expect it to shift too much upon fine-
tuning on slightly shifted data.

  4. 4.

It depends directly on the _smoothness constant_ $\beta$ of the activation
function used in the model. The smoother the model, the smaller are its second
derivatives, and thus its gradient-input Lipschitz.

These insights motivate the core thesis of the paper, i.e., that explanation
stability depends on properties beyond predictive performance. Specifically,
we see that specific algorithm choices that do not necessarily improve
accuracy, such as weight decay, smoothness constant and robustness-inducing
losses, affect explanation and model stability to data shifts.

We conclude by making two remarks regarding the theory. First, it only
highlights _sufficient conditions_ for explanation stability, not necessary
conditions. Meaning, there can be ways to achieve gradient stability without
the parameters $\theta_{1},\theta_{2}$ being close, or without neural nets
being Lipschitz, or when ReLU or weight decay are not used, but these lie
outside the scope of this theory. Second, the focus of these bounds is not to
produce numerically tight estimates for explanation stability, rather, it is
intended as a tool to make qualitative predictions regarding specific
modelling interventions that can improve explanation stability in practice.
Section 4 explores this in more detail.

##  4 Experimental Evaluation

In this section, we describe our general evaluation setup. In §4.1, we
experimentally validate that the modeling interventions identified in §3
improve explanation stability. Next, in §4.2, we show that these relationships
extend to realistic training setups and more-sophisticated explanation
techniques (i.e., top-K feature attributions of techniques like SmoothGrad,
LIME, and SHAP). Finally, in §4.3, we show that other hyperparameters not
discussed in our theory – e.g., learning rate, learning rate decay, batch
size, and number of training epochs – also impact gradient stability.

##### Datasets

We evaluate explanation stability on three binary datasets with continuous
features. We use the WHO life expectancy dataset ($n$=2928, number of features
$d$=18) [35], which evaluates whether a country’s life expectancy is above the
global median based on health and economic factors; the HELOC dataset
($n$=9871, $d$=23) [15], which evaluates loan acceptance given the applicant’s
financial information; and the Adult Income dataset ($n$=32561, $d$=6) [11],
which evaluates whether a person’s income is above $50,000. We model dataset
shifts in two ways:

  * •

Synthetic Noise. We add Gaussian $\mathcal{N}(0,\sigma^{2})$ noise to all
training samples, with varying levels of noise $\sigma$.

  * •

Temporal Shift. We create a temporal shift dataset from the WHO data by using
pre-2012 data as the “original” dataset and the full dataset as the “shifted”
dataset.

##### Models

We evaluate the stability of medium-sized neural networks containing 5 layers
of 50 nodes. By default, we use ReLU activation; however, to see the effects
of model curvature on explanation stability, we use softplus (SP) for some
experiments. For the experiments using SP, we vary the parameter $\beta$ from
2 to 10 depending on the dataset and the experiment’s goals. Smaller values of
$\beta$ correspond to lower-curvature models, while SP with $\beta=\infty$ is
equivalent to ReLU. Unless otherwise specified, all models we compare within
each section achieve similar predictive accuracy.

##### Experimental procedure

For each experiment, we average over 10 trials to reduce the impact that
random seeds have on the outcomes. Each trial in the synthetic noise
experiments consists of training one “base” model plus 10 “noisy” models, each
corresponding to a different random perturbation of the data. In the temporal
shift experiments, each trial compares the base model (trained on the partial
dataset) to a new model (trained on the full dataset). We compute explanation
stability for all test samples, except for LIME and SHAP results, where we
limit computations to 100 randomly-chosen test samples. In the results, we
display the mean of all trials along with the middle 50% of values.

###  4.1 Validating Our Theoretical Results

##### Experiment setup

In this section, we test our theoretical results from §3 by evaluating the
impact of _weight decay_ , _model curvature_ , and _size of the data shift_ on
explanation stability. To ensure that our analysis is sound, we adhere to the
theoretical assumptions as closely as possible. In particular, we focus on
_model parameters_ ($\theta$) and _gradients_ ($g$) rather than top-K feature
attributions, since top-K explanations are not differentiable. We train each
base model for a large number of epochs to ensure convergence of the training
loss to a minimum (i.e., we ignore overfitting on the test set). To fine-tune,
we copy the parameters of the base model and then train on the shifted data.
Also, we start fine-tuning with a lower learning rate and train for fewer
epochs, given that the initial loss is already close to a minimum.

Figure 1: Effect of weight decay value on the mean gradient and parameter changes for HELOC after fine-tuning on a shifted dataset. All models use ReLU activation. The x-axis is the size of the data shift, represented by the standard deviation of noise, $\sigma$. Table 1: Effect of weight decay ($\gamma$) on WHO (with fine-tuning). | $\gamma=0$ | $0.001$ | $0.01$
---|---|---|---
$\|\theta_{1}-\theta_{2}\|_{2}$ | 1.69±0.09 | 1.53±0.03 | 0.95±0.02
$\|g_{1}-g_{2}\|_{2}$ | 0.201±0.065 | 0.122±0.034 | 0.010±0.002

For brevity, we omit similar results from the Adult dataset, with complete
results in the appendix. We also provide the training hyperparameters for each
dataset in the appendix.

##### Weight decay results

Figure 1 and Table 1 show the effect of weight decay for the synthetically and
naturally-shifted datasets, respectively. We see that larger weight decay
values correspond to orders of magnitude smaller mean gradient and parameter
changes, which is in line with the theoretical relationship between weight
decay and stability. As the amount of synthetic noise grows, we observe that
both parameters and gradients diverge more between the base model and fine-
tuned model. This pattern supports our hypotheses that (a) the loss landscape
changes more with larger data shift and (b) a larger change in model
parameters is correlated with a larger change in test sample gradients.

##### Curvature results

Figure 2 and Table 2 show that using low-curvature training techniques
increases gradient stability. In particular, SP does better than ReLU,
particularly with smaller values of $\beta$. For both the real-world and
synthetic shifts, the gradient changes are less distinct (i.e., the confidence
intervals overlap more) than the parameter changes, but still show a trend in
the predicted direction.

Figure 2: Effect of model curvature on the mean gradient and parameter changes for the HELOC dataset after fine-tuning on a shifted dataset. The x-axis is the size of the data shift, represented by the standard deviation of noise, $\sigma$. Table 2: Effect of curvature on WHO (with fine-tuning). | ReLU |  SP ($\beta=10$) |  SP ($\beta=5$)
---|---|---|---
$\|\theta_{1}-\theta_{2}\|_{2}$ | 1.53±0.03 | 1.37±0.05 | 1.23±0.02
$\|g_{1}-g_{2}\|_{2}$ | 0.122±0.034 | 0.096±0.026 | 0.070±0.024

###  4.2 Going Beyond Theoretical Results

In this section, we evaluate whether our theoretical predictions still hold if
we deviate from the theoretical assumptions, namely: fine-tuning, exactly
reaching a local minimum, and using non-differentiable top-K feature
attribution metrics rather than gradients $\ell_{2}$ distance. The top-K
features for a model’s predictions are the $k$ features with attribution
values (e.g., gradients or functions of gradients for gradient-based methods;
weights of linear models for LIME and SHAP) of highest magnitude. For example,
if an instance $\mathbf{x}$ comprises of four features indexed by $[1\cdots
4]$, and the gradient (feature attribution) of a given model w.r.t.
$\mathbf{x}$ is $[-0.1,-0.4,0.2,0.3]$, then the top-2 features are features 2
and 4 since their corresponding feature attribution magnitudes -0.4 and 0.3
are the highest.

##### Experiment setup

We maintain our experimental setup from §4.1, except that we train each model
for a smaller number of epochs (until accuracy stabilizes), before retraining
a model from scratch on the shifted data, rather than fine-tuning
(hyperparameters are located in the appendix). We use the same random network
initialization for both training and retraining, as prior work has shown that
feature attributions are not stable to changes in random seed [3].

Figure 3: Effect of weight decay on parameter stability when retraining on HELOC and Adult datasets. The x-axis is the size of the data shift, represented by the standard deviation of noise, $\sigma$. Table 3: Effect of weight decay ($\gamma$) on WHO (with retraining). | $\gamma=0$ | $0.001$ | $0.01$
---|---|---|---
$\|\theta_{1}-\theta_{2}\|_{2}$ | 3.97±0.13 | 3.50±0.19 | 1.79±0.17
$\|g_{1}-g_{2}\|_{2}$ | 0.295±0.075 | 0.150±0.034 | 0.003±0.001
Table 4: Effect of curvature on WHO (with retraining). | ReLU |  SP ($\beta=10$) |  SP ($\beta=5$)
---|---|---|---
$\|\theta_{1}-\theta_{2}\|_{2}$ | 3.50±0.19 | 2.60±0.05 | 1.97±0.12
$\|g_{1}-g_{2}\|_{2}$ | 0.150±0.034 | 0.102±0.040 | 0.076±0.024
Figure 4: Effect of model curvature on gradient stability when retraining on
HELOC and Adult datasets. The x-axis is the size of the data shift,
represented by the standard deviation of noise, $\sigma$. Figure 5:  Top-5 SA
for the HELOC dataset under model retraining. The left two graphs show the
effect of weight decay and the right two graphs show the effect of curvature.
Graphs 1 and 3 show Top-5 SA for saliency, while graphs 2 and 4 show Top-5 SA
for SmoothGrad. Confidence intervals represent the middle 50% of values. See
the appendix for other datasets and explanation techniques.

To complete the analysis, we consider four feature attribution techniques:
input gradients, input gradients with SmoothGrad, LIME, and kernel SHAP [27,
37, 42, 40]. We use three top-K stability metrics adapted from Brunet et al.
[5]. They each take in a pair of top-K feature sets corresponding to two
models. Sign Agreement (SA) returns the fraction of top-K features that appear
in both models’ top-K features and have the same signed value. Consistent
Direction of Contribution (CDC) is binary (per-sample) and measures whether
all features in the top-K (for either model) have the same signed value in the
other model. Signed-Set Agreement (SSA) is also binary (per-sample) and is 1
if the two model’s top-K feature sets contain the same features and the
features have the same signed value (the order of the top-K features does not
matter).

##### Relaxing Training Assumptions

From Table 3, we see that increasing the algorithm’s weight decay
significantly decreases both parameters and gradients. Figure 3 confirms that
a larger weight decay during training also increases gradient stability for
models with synthetic noise (analogous results for parameter stability are
included in the appendix). Figure 4 and Table 4 show that the curvature trends
discovered in §4.1 still hold, i.e., lower model curvature induces higher
gradient stability.

Although curvature and weight decay increase explanation stability upon
retraining, we note parameter stability scores are much lower for retraining
than for fine-tuning given fixed curvatures and weight decays. For example,
the mean change in parameters when using ReLU with a weight decay of 0.001 is
1.53±0.03 when fine-tuning, versus 3.50±0.19 when retraining. The difference
in gradient change between the two techniques is less stark, however,
suggesting that even though the fine-tuned models are more similar to each
other overall, the retrained models are sufficiently similar on the parts of
the input domain represented in the test set.

##### Using Diverse Explanation Techniques

Figure 5 shows Top-5 SA stability for two popular explanation techniques,
namely, Saliency (input gradients) and SmoothGrad [24]. As in the previous
experiments, models with low curvature and large weight decay are more
effective at preserving explanation stability. In case of saliency, the
explanation stability starts to decrease sharply for shift magnitudes larger
than $\mathcal{N}(0,0.1^{2})$. For SmoothGrad, the noise has less of an
effect, which we hypothesize is because SmoothGrad is more stable overall as
it computes an averaged gradient over a small local neighborhood of a given
instance (see Table 5). Lowering the model’s curvature and increasing the
weight decay similarly improve explanation stability under data shifts for the
other metrics and explanation techniques (see the appendix for a full set of
graphs).

Table 5 gives a more detailed breakdown of top-K feature attribution stability
across all datasets and explanation techniques. The table corresponds to a
single modeling choice – ReLU with weight decay 0.001. We notice trends based
on the _explanation technique_ , _explanation metric_ , and _dataset_.

The different explanation techniques exhibit vastly different stabilities
under dataset shifts. Notably, SmoothGrad outperforms all other explanation
techniques, obtaining scores of over 90% for both Top-3 and Top-5 explanation
stability for all of the metrics we used. That is, over 90% of test samples’
attributions had the same signs for all of each model’s top-K features (CDC),
over 90% of test samples’ attributions had the same signed features in the
top-K (SSA), and on average, over 90% of the top-K features for each test
sample appeared in the other model’s top-K features (SA).

To compare the explanation metrics, we look at the relative difference between
SSA and each of SA and CDC (since SSA is strictly stronger than either of the
other metrics). We see that neither SA nor CDC is consistently more stable,
suggesting that the failure to have identical top-K feature sets stems both
from different features ranking as the most influential (i.e., what SA
captures) and also whether the top features have the same sign in the other
model (i.e., CDC).

We also see that different datasets exhibit different levels of explanation
stability. The top-3 and top-5 scores for the Adult datasets are very high;
however, this dataset only has six features, so there are fewer possible
feature combinations. We see that SSA scores for Adult are much lower than the
SA scores, indicating that even with only 6 features, the learnt models often
disagree on different features’ impacts.

Table 5: Explanation stability scores for all datasets. The synthetically-shifted datasets (HELOC and Adult) are modified by $\mathcal{N}(0,0.1^{2})$ noise. Models use ReLU activation and weight decay $\gamma=0.001$. All values averaged across all test samples and then across all trials. Error bounds are such that the lower and upper bound values represent the 25${}^{\text{th}}$ and 75${}^{\text{th}}$ percentile respectively. | Explanation | Top-3 | Top-5
---|---|---|---
Dataset | technique | SA | CDC | SSA | SA | CDC | SSA
WHO | Saliency | 0.63±0.01 | 0.83±0.03 | 0.18±0.03 | 0.61±0.03 | 0.61±0.10 | 0.03±0.01
SmoothGrad | 0.94±0.00 | 0.94±0.01 | 0.91±0.00 | 0.94±0.00 | 0.92±0.00 | 0.90±0.00
LIME | 0.69±0.09 | 0.35±0.24 | 0.27±0.20 | 0.59±0.57 | 0.08±0.05 | 0.06±0.05
K.SHAP | 0.58±0.10 | 0.60±0.39 | 0.13±0.08 | 0.61±0.10 | 0.41±0.31 | 0.05±0.05
HELOC | Saliency | 0.51±0.02 | 0.88±0.03 | 0.07±0.01 | 0.55±0.02 | 0.71±0.05 | 0.01±0.00
SmoothGrad | 0.94±0.00 | 0.96±0.00 | 0.90±0.00 | 0.94±0.00 | 0.94±0.00 | 0.90±0.00
LIME | 0.63±0.02 | 0.45±0.04 | 0.18±0.03 | 0.66±0.02 | 0.16±0.03 | 0.07±0.02
K.SHAP | 0.64±0.02 | 0.91±0.03 | 0.19±0.03 | 0.70±0.02 | 0.79±0.05 | 0.09±0.02
Adult | Saliency | 0.87±0.01 | 0.92±0.01 | 0.62±0.04 | 0.85±0.01 | 0.57±0.04 | 0.38±0.03
SmoothGrad | 0.98±0.00 | 0.98±0.00 | 0.96±0.00 | 0.98±0.00 | 0.96±0.00 | 0.94±0.00
LIME | 0.95±0.00 | 0.86±0.02 | 0.86±0.02 | 0.97±0.00 | 0.84±0.02 | 0.84±0.02
K.SHAP | 0.88±0.01 | 0.88±0.04 | 0.65±0.04 | 0.84±0.02 | 0.47±0.10 | 0.32±0.08

###  4.3 Sensitivity Analysis

In this section, we analyze the impact of other training hyperparameters on
explanation stability, epoch-by-epoch throughout the entire training process.

Given our observation that parameter stability promotes explanation stability,
we expect a lower number of epochs and a lower learning rate to both increase
explanation stability under dataset shift by decreasing the divergence between
optimization paths in parameter space. Based on the geometry of two divergent
optimization paths, we anticipate that if the base model is trained for a long
period of time, the retrained model’s parameters that are most similar to the
base model will occur relatively sooner in the retraining process. This
emphasizes the challenge of retraining models that both perform well and
remain close to the base model.

We track explanations during each epoch of retraining, measuring similarity
using gradient $\ell_{2}$ distance and top-K SA consistency. Each trial
consists of varying a single hyperparameter while keeping all others fixed.
For each value of the varied hyperparameter, we compare one base model with 10
models retrained on synthetically shifted HELOC data
$\sim\mathcal{N}(0,0.1^{2})$. Recall that [19] studied the hypothesis
stability of stochastic gradient descent (SGD) and showed that practical
choices such as learning rate and number of epochs of training can impact the
stability of model parameters; the following experiments confirm these
findings in the context of explanation stability.

##### Learning Rate

We verify in Figure 6a that increasing the learning rate of all models
effectively worsens gradient stability. We observe that the maximal gradient
stability occurs sooner when the learning rate is higher. We hypothesize that
this is because a high learning rate stretches the optimization path of the
base model, taking it further away from the initialization in parameter space.
Since the retrained models take divergent paths from the base model (due to
the dataset shifts), the likelihood that they achieve similar parameters
decreases as learning rate increases.

##### Learning Rate Decay

Building upon the understanding of the effect of learning rate on gradient
stability, we delve into its counterpart, learning rate decay, and its
influence, as depicted in Figure 6b. The results are as expected: decreasing
the decay value, i.e. increasing the amount of decay, shows similar behaviour
to decreasing the learning rate. Namely, the optimization trajectories have
less chance to diverge, and the resultant base model and retrained models are
thus more likely to have similar parameters. While a higher amount of decay
(lower decay value) is preferable, this can inhibit the training process
(notice how a decay value of 0.9 causes training to stop after around 20
epochs).

##### Batch Size

In Figure 6c, we see that decreasing batch size worsens gradient stability. We
hypothesize that this behaviour is because there are a greater number of steps
in each epoch, inducing a higher level of randomness in the optimization
trajectory. Each step adds an opportunity for the retrained model’s path to
deviate from the initial base model’s path, due to the dataset shifts. As
such, the likelihood of achieving similar parameters between the retrained
models and the base models decreases with smaller batch sizes. Additionally,
the points at which explanation stability is maximized occur after fewer
epochs for smaller batch sizes, since those correspond to both longer and more
divergent optimization trajectories.

(a)

(b)

(c)

(d)

Figure 6: Raw gradient similarity (y-axis) of retrained models are shown per
epoch of training (x-axis). For gradient norm, lower values are better, while
for top-5 SA, higher values (closer to 1) are better. Default values are:
learning rate of 0.2; learning rate decay of 1.0 (no decay); batch size of
128; base models are trained for 30 epochs (besides the bottom right); and
synthetic noise, $\sigma=0.1$.

##### Base Model Epochs

Figure 6d verifies that as the number of base model epochs is increased,
explanation stability suffers, with significantly larger gradient changes and
lower top-5 consistency scores. For a sufficiently small base model epoch of
20, we observe that similarity also peaks at around 20 training epochs.
However, if we were to increase the base model epoch to 60, as shown in green,
we observe that the peak still occurs at around the same value. This means
that, after just 20 epochs, the optimization path of the retrained model must
have diverged so much from the base model that any further training leads only
to a decrease in stability.

##  5 Conclusions

In this paper, we theoretically and empirically characterize how modelling
choices impact model explanations’ stability to retraining under data shifts.
Our formal results suggest, and our experimental results confirm, that model
explanations can be made more stable to retraining after data shifts by using
low-curvature activation functions and large weight decay values. We also took
a first step towards empirically characterizing how different explanation
techniques yield variable explanation stabilities under dataset shifts, noting
that SmoothGrad outperforms other methods.

The message for ML practitioners is clear: in order to preserve explanations
upon model retraining due to data shifts at little cost to accuracy, use
softplus with the lowest possible $\beta$, and the largest weight decay that
does not affect predictive accuracy. If possible, practitioners should fine-
tune the model (§4.1) rather than retrain from scratch to avoid instability
due to the multiplicity of optimal solutions. However, as shown in §4.3,
training hyperparameters can also play a large role in explanation stability,
so care must be taken to choose hyperparameters in a way that balances
accuracy and explanation stability.

There are several directions for future work. First, our theoretical results
are only applicable to differentiable models. In the appendix, we empirically
evaluate explanation stability when retraining decision trees, random forests,
and XGBoost models, but additional work is needed to formalize explanation
stability to retraining for these and other classes of models. Second, while
we do perform experiments with a temporal dataset shift, we do not evaluate
other types of naturally-occurring data shifts like geographic shifts or data
corrections (e.g., as with the German Credit Statlog dataset [11]). Future
work should focus on understanding explanation stability for various types of
real-world data shifts. Third, we limit our evaluation to tabular datasets due
to the commonality of these datasets where laws like a Right to Explanation
are relevant. However, other sources of data may become more relevant to user
explanations in the future and validating our theoretical results on non-
tabular data is of independent interest. Finally, rather than sticking to
existing explanation techniques, future work can explore how to develop novel
explanation techniques that preserve explanation stability despite small
changes in the underlying models and data.

#### Acknowledgements

The authors would like to thank the anonymous reviewers for their helpful
feedback and all the funding agencies listed below for supporting this work.
This work is supported in part by the NSF awards #IIS-2008461, #IIS-2040989,
and #IIS-2238714, and research awards from Google, JP Morgan, Amazon, Harvard
Data Science Initiative, and D3̂ Institute at Harvard. The views expressed
here are those of the authors and do not reflect the official policy or
position of the funding agencies.

## References

  * [1] C. J. Anders, P. Pasliev, A.-K. Dombrowski, K.-R. Müller, and P. Kessel.  Fairwashing Explanations with Off-Manifold Detergent.  In Proceedings of the 37th International Conference on Machine Learning, ICML’20. JMLR.org, 2020.

  * [2] E. Black and M. Fredrikson.  Leave-One-out Unfairness.  In Proceedings of the 2021 ACM Conference on Fairness, Accountability, and Transparency, FAccT ’21, page 285–295. Association for Computing Machinery, 2021.
  * [3] E. Black, K. Leino, and M. Fredrikson.  Selective Ensembles for Consistent Predictions.  In International Conference on Learning Representations, 2022.
  * [4] O. Bousquet and A. Elisseeff.  Stability and Generalization.  The Journal of Machine Learning Research, 2:499–526, 2002.
  * [5] M.-E. Brunet, A. Anderson, and R. Zemel.  Implications of Model Indeterminacy for Explanations of Automated Decisions.  In A. H. Oh, A. Agarwal, D. Belgrave, and K. Cho, editors, Advances in Neural Information Processing Systems, 2022.
  * [6] N. Bui, D. Nguyen, and V. A. Nguyen.  Counterfactual Plans under Distributional Ambiguity.  In International Conference on Learning Representations, 2022.
  * [7] A. D’Amour, K. Heller, D. Moldovan, B. Adlam, B. Alipanahi, A. Beutel, C. Chen, J. Deaton, J. Eisenstein, M. D. Hoffman, F. Hormozdiari, N. Houlsby, S. Hou, G. Jerfel, A. Karthikesalingam, M. Lucic, Y. Ma, C. McLean, D. Mincu, A. Mitani, A. Montanari, Z. Nado, V. Natarajan, C. Nielson, T. F. Osborne, R. Raman, K. Ramasamy, R. Sayres, J. Schrouff, M. Seneviratne, S. Sequeira, H. Suresh, V. Veitch, M. Vladymyrov, X. Wang, K. Webster, S. Yadlowsky, T. Yun, X. Zhai, and D. Sculley.  Underspecification Presents Challenges for Credibility in Modern Machine Learning.  Journal of Machine Learning Research, 23(226):1–61, 2022.
  * [8] A.-K. Dombrowski, M. Alber, C. Anders, M. Ackermann, K.-R. Müller, and P. Kessel.  Explanations can be manipulated and geometry is to blame.  Advances in Neural Information Processing Systems, 32, 2019.
  * [9] R. Dominguez-Olmedo, A. H. Karimi, and B. Schölkopf.  On the Adversarial Robustness of Causal Algorithmic Recourse.  In K. Chaudhuri, S. Jegelka, L. Song, C. Szepesvari, G. Niu, and S. Sabato, editors, Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 5324–5342. PMLR, 17–23 Jul 2022.
  * [10] F. Doshi-Velez and B. Kim.  Towards A Rigorous Science of Interpretable Machine Learning.  arXiv: Machine Learning, 2017.
  * [11] D. Dua and C. Graff.  UCI Machine Learning Repository, 2017.
  * [12] S. Dutta, J. Long, S. Mishra, C. Tilli, and D. Magazzeni.  Robust Counterfactual Explanations for Tree-Based Ensembles.  In International Conference on Machine Learning, pages 5742–5756. PMLR, 2022.
  * [13] European Commission.  General Data Protection Regulation (GDPR).
  * [14] T. Fel, D. Vigouroux, R. Cadène, and T. Serre.  How Good Is Your Explanation? Algorithmic Stability Measures To Assess the Quality of Explanations for Deep Neural Networks.  In Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision (WACV), pages 720–730, January 2022.
  * [15] FICO.  Explainable Machine Learning Challenge, 2018.
  * [16] A. Forel, A. Parmentier, and T. Vidal.  Robust Counterfactual Explanations for Random Forests.  arXiv preprint arXiv:2205.14116, 2022.
  * [17] A. Ghorbani, A. Abid, and J. Zou.  Interpretation of Neural Networks is Fragile.  In Proceedings of the Thirty-Third AAAI Conference on Artificial Intelligence and Thirty-First Innovative Applications of Artificial Intelligence Conference and Ninth AAAI Symposium on Educational Advances in Artificial Intelligence, AAAI’19/IAAI’19/EAAI’19. AAAI Press, 2019.
  * [18] T. Han, S. Srinivas, and H. Lakkaraju.  Which Explanation Should I Choose? A Function Approximation Perspective to Characterizing Post Hoc Explanations.  In A. H. Oh, A. Agarwal, D. Belgrave, and K. Cho, editors, Advances in Neural Information Processing Systems, 2022.
  * [19] M. Hardt, B. Recht, and Y. Singer.  Train faster, generalize better: Stability of stochastic gradient descent.  In International conference on machine learning, pages 1225–1234. PMLR, 2016.
  * [20] J. Heo, S. Joo, and T. Moon.  Fooling Neural Network Interpretations via Adversarial Model Manipulation.  Advances in Neural Information Processing Systems, 32, 2019.
  * [21] A.-H. Karimi, B. Schölkopf, and I. Valera.  Algorithmic Recourse: From Counterfactual Explanations to Interventions.  In Proceedings of the 2021 ACM Conference on Fairness, Accountability, and Transparency, FAccT ’21, page 353–362, New York, NY, USA, 2021. Association for Computing Machinery.
  * [22] A.-H. Karimi, J. von Kügelgen, B. Schölkopf, and I. Valera.  Algorithmic recourse under imperfect causal knowledge: a probabilistic approach.  In H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 265–277. Curran Associates, Inc., 2020.
  * [23] P. W. Koh and P. Liang.  Understanding Black-Box Predictions via Influence Functions.  In Proceedings of the 34th International Conference on Machine Learning - Volume 70, ICML’17, page 1885–1894. JMLR.org, 2017.
  * [24] N. Kokhlikyan, V. Miglani, M. Martin, E. Wang, B. Alsallakh, J. Reynolds, A. Melnikov, N. Kliushkina, C. Araya, S. Yan, and O. Reblitz-Richardson.  Captum: A unified and generic model interpretability library for PyTorch, 2020.
  * [25] J. F. Kolen and J. B. Pollack.  Back Propagation is Sensitive to Initial Conditions.  In Proceedings of the 1990 Conference on Advances in Neural Information Processing Systems 3, NIPS-3, page 860–867, San Francisco, CA, USA, 1990. Morgan Kaufmann Publishers Inc.
  * [26] A. V. Looveren and J. Klaise.  Interpretable Counterfactual Explanations Guided by Prototypes.  ArXiv, abs/1907.02584, 2019.
  * [27] S. M. Lundberg and S.-I. Lee.  A Unified Approach to Interpreting Model Predictions.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
  * [28] J. Mehrer, C. J. Spoerer, N. Kriegeskorte, and T. C. Kietzmann.  Individual differences among deep neural network models.  Nature communications, 11(1):5725, November 2020.
  * [29] B. Neyshabur, R. Tomioka, and N. Srebro.  Norm-based capacity control in neural networks.  In Conference on learning theory, pages 1376–1401. PMLR, 2015.
  * [30] D. Nguyen, N. Bui, and V. A. Nguyen.  Distributionally Robust Recourse Action, 2022.
  * [31] T.-D. H. Nguyen, N. Bui, D. Nguyen, M.-C. Yue, and V. A. Nguyen.  Robust Bayesian Recourse.  In Uncertainty in Artificial Intelligence, pages 1498–1508. PMLR, 2022.
  * [32] Office of Science and Technology Policy.  Blueprint for an AI Bill of Rights.
  * [33] M. Pawelczyk, C. Agarwal, S. Joshi, S. Upadhyay, and H. Lakkaraju.  Exploring Counterfactual Explanations Through the Lens of Adversarial Examples: A Theoretical and Empirical Analysis.  In G. Camps-Valls, F. J. R. Ruiz, and I. Valera, editors, Proceedings of The 25th International Conference on Artificial Intelligence and Statistics, volume 151 of Proceedings of Machine Learning Research, pages 4574–4594. PMLR, 28–30 Mar 2022.
  * [34] R. Poyiadzi, K. Sokol, R. Santos-Rodriguez, T. De Bie, and P. Flach.  FACE: Feasible and Actionable Counterfactual Explanations.  In Proceedings of the AAAI/ACM Conference on AI, Ethics, and Society, AIES ’20, page 344–350, New York, NY, USA, 2020. Association for Computing Machinery.
  * [35] K. Rajarshi.  Life Expectancy (WHO), 2017.
  * [36] K. Rawal, E. Kamar, and H. Lakkaraju.  Algorithmic Recourse in the Wild: Understanding the Impact of Data and Model Shifts, 2020.
  * [37] M. T. Ribeiro, S. Singh, and C. Guestrin.  ”Why Should I Trust You?”: Explaining the Predictions of Any Classifier.  In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’16, page 1135–1144, New York, NY, USA, 2016. Association for Computing Machinery.
  * [38] M. T. Ribeiro, S. Singh, and C. Guestrin.  Anchors: High-Precision Model-Agnostic Explanations.  In Proceedings of the Thirty-Second AAAI Conference on Artificial Intelligence and Thirtieth Innovative Applications of Artificial Intelligence Conference and Eighth AAAI Symposium on Educational Advances in Artificial Intelligence, AAAI’18/IAAI’18/EAAI’18. AAAI Press, 2018.
  * [39] R. R. Selvaraju, M. Cogswell, A. Das, R. Vedantam, D. Parikh, and D. Batra.  Grad-CAM: Visual Explanations from Deep Networks via Gradient-Based Localization.  In 2017 IEEE International Conference on Computer Vision (ICCV), pages 618–626, 2017.
  * [40] K. Simonyan, A. Vedaldi, and A. Zisserman.  Deep Inside Convolutional Networks: Visualising Image Classification Models and Saliency Maps.  CoRR, abs/1312.6034, 2013.
  * [41] D. Slack, S. Hilgard, E. Jia, S. Singh, and H. Lakkaraju.  Fooling LIME and SHAP: Adversarial Attacks on Post hoc Explanation Methods.  In Proceedings of the AAAI/ACM Conference on AI, Ethics, and Society, pages 180–186, 2020.
  * [42] D. Smilkov, N. Thorat, B. Kim, F. Viegas, and M. M. Wattenberg.  SmoothGrad: removing noise by adding noise.  In ICML Workshop on Visualization for Deep Learning, Sydney, Australia, 2017.
  * [43] S. Srinivas, K. Matoba, H. Lakkaraju, and F. Fleuret.  Efficient Training of Low-Curvature Neural Networks.  In A. H. Oh, A. Agarwal, D. Belgrave, and K. Cho, editors, Advances in Neural Information Processing Systems, 2022.
  * [44] M. Sundararajan, A. Taly, and Q. Yan.  Axiomatic Attribution for Deep Networks.  In D. Precup and Y. W. Teh, editors, Proceedings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research, pages 3319–3328. PMLR, 06–11 Aug 2017.
  * [45] S. Upadhyay, S. Joshi, and H. Lakkaraju.  Towards Robust and Reliable Algorithmic Recourse.  In M. Ranzato, A. Beygelzimer, Y. Dauphin, P. Liang, and J. W. Vaughan, editors, Advances in Neural Information Processing Systems, volume 34, pages 16926–16937. Curran Associates, Inc., 2021.
  * [46] B. Ustun, A. Spanghler, and Y. Liu.  Actionable Recourse in Linear Classification.  In ACM Conference on Fairness, Accountability, and Transparency (FAT ’19), 2019.
  * [47] S. Wachter, B. Mittelstadt, and C. Russell.  Counterfactual Explanations Without Opening the Black Box: Automated Decisions and the GDPR.  Harvard Journal of Law & Technology, 31:841–887, 04 2018.
  * [48] W. Yang, J. Li, C. Xiong, and S. C. H. Hoi.  MACE: An Efficient Model-Agnostic Framework for Counterfactual Explanation.  ArXiv, abs/2205.15540, 2022.

##  Appendix A Proofs

In this section, we provide the proofs for the theory presented in the main
paper. We present the notations here again for convenience in Table 6. As with
the main paper, we first bound parameter shift in §A.1 and then explanation
shift in §A.2.

Table 6: Notations $z_{i}$ | $(x_{i},y_{i});\leavevmode\nobreak\ x_{i}\in\mathbb{R}^{d};\leavevmode\nobreak\ y\in\mathbb{R}$ | data point
---|---|---
$\mathcal{D}_{1}$ | $\\{z_{1},..z_{n}\\}$ | original training set
$\mathcal{D}_{2}$ | $\\{z^{\prime}_{1},...z^{\prime}_{n}\\}$ | new training set
$\ell_{reg}(\theta,z_{i})$ |  $\ell(\theta,z_{i})+\gamma\|\theta\|^{2}_{2}$; $\ell(\theta,z_{i})\in\mathbb{R}^{+}$ | regularized loss function
$\ell_{\mathcal{D}_{1}}(\theta)$ | $\frac{1}{n}\sum_{i=1}^{n}\ell_{reg}(\theta,z_{i})$ |  loss incurred by $\theta$ on dataset $D_{1}$
$\theta_{1}$ | $\arg\min_{\theta}\ell_{\mathcal{D}_{1}}(\theta)$ |  optimal weights on dataset $D_{1}$
$\epsilon_{1}+\gamma\|\theta_{1}\|^{2}$ | $\ell_{\mathcal{D}_{1}}(\theta_{1})=\min_{\theta}\ell_{\mathcal{D}_{1}}(\theta)$ |  optimal loss value on dataset $\mathcal{D}_{1}$

###  A.1 Bounding the parameter shift

In order to bound the parameter shift, we first prove an intermediate result
connecting the loss on $\mathcal{D}_{2}$ of the optimum obtained by minimizing
the loss on $\mathcal{D}_{1}$.

######  Lemma 2.

The loss $\ell_{\mathcal{D}_{2}}$ on $\theta_{1}$ is given by

$\displaystyle\ell_{\mathcal{D}_{2}}(\theta_{1})\leq\epsilon_{1}+\gamma\|\theta_{1}\|^{2}+\mathcal{L}_{x}(\theta_{1})d(\mathcal{D}_{1},\mathcal{D}_{2})$

###### Proof.

$\displaystyle\ell_{\mathcal{D}_{2}}(\theta_{1})$
$\displaystyle=\ell_{\mathcal{D}_{1}}(\theta_{1})+\left(\ell_{\mathcal{D}_{2}}(\theta_{1})-\ell_{\mathcal{D}_{1}}(\theta_{1})\right)$
$\displaystyle\leq\ell_{\mathcal{D}_{1}}(\theta_{1})+|\ell_{\mathcal{D}_{2}}(\theta_{1})-\ell_{\mathcal{D}_{1}}(\theta_{1})|$
$\displaystyle\leq\epsilon_{1}+\gamma\|\theta_{1}\|^{2}+\frac{1}{n}\sum_{i=1}^{n}|\ell(\theta_{1},z_{i})-\ell(\theta_{1},z^{\prime}_{i})|\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\text{(Convexity of abs fn)}$
$\displaystyle\leq\epsilon_{1}+\gamma\|\theta_{1}\|^{2}+\frac{1}{n}\sum_{i=1}^{n}\mathcal{L}_{x}(\theta_{1})\|z_{i}-z^{\prime}_{i}\|_{2}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \text{(Defn of
Lipschitz)}$
$\displaystyle\leq\epsilon_{1}+\gamma\|\theta_{1}\|^{2}+\mathcal{L}_{x}(\theta_{1})\underbrace{\frac{1}{n}\sum_{i=1}^{n}\|z_{i}-z^{\prime}_{i}\|}_{d(\mathcal{D}_{1},\mathcal{D}_{2})}$

∎

In the result above, $\epsilon_{1}+\gamma\|\theta_{1}\|^{2}$ is the optimal
value of the regularized loss, and the term following that denotes the shift
from the optimal value incurred. As mentioned in the paper, we use the
following additional assumptions to derive the main result: (1) we minimize a
regularized loss $\ell_{reg}$ that involves weight decay, i.e,
$\ell_{reg}(\theta)=\ell(\theta)+\gamma\|\theta\|^{2}_{2}$; (2) $\ell$ is
locally quadratic; (3) the learning algorithm returns a unique minimum
$\theta$ given a dataset $\mathcal{D}$.

######  Theorem 3.

Given the assumptions stated above, and that $\mathcal{L}_{x}(\theta_{1})$ is
the Lipschitz constant of the model with parameters $\theta_{1}$, we have

$\displaystyle\|\theta_{2}-\theta_{1}\|_{2}$
$\displaystyle\leq\sqrt{\frac{\mathcal{L}_{x}(\theta_{1})d(\mathcal{D}_{1},\mathcal{D}_{2})}{\gamma}}+C$

where $\gamma$ is the weight decay regularization constant, and $C$ is a small
problem-dependent constant.

###### Proof.

Assuming the loss $\ell_{D_{2}}$ is locally quadratic around the optimal
solution $\theta_{2}$, we can employ a second order Taylor series expansion as
follows:

$\displaystyle\ell_{D_{2}}(\theta_{1})$
$\displaystyle=\ell_{D_{2}}(\theta_{2})+\nabla_{\theta}\ell_{D_{2}}(\theta_{2})^{\top}(\theta_{1}-\theta_{2})+\frac{1}{2}(\theta_{1}-\theta_{2})^{\top}H_{\theta_{2}}(\theta_{1}-\theta_{2})$
(1)
$\displaystyle=\epsilon_{2}+\gamma\|\theta_{2}\|_{2}^{2}+\frac{1}{2}(\theta_{1}-\theta_{2})^{\top}H_{\theta_{2}}(\theta_{1}-\theta_{2})\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ (\text{global
optimality of $\theta_{2}$ on $S_{2}$})$ (2)

In the statement above we have used the fact that
$\nabla_{\theta}\ell_{\mathcal{D}_{2}}(\theta_{2})=0$ because of first order
optimality criteria, and that
$\epsilon_{2}+\gamma\|\theta_{2}\|^{2}=\ell_{\mathcal{D}_{2}}(\theta_{2})$, is
the optimal value. Now we use Lemma 2 to substitute
$\ell_{\mathcal{D}_{2}}(\theta_{1})$, and we have:

$\displaystyle\epsilon_{2}+\gamma\|\theta_{2}\|_{2}^{2}+\frac{1}{2}(\theta_{1}-\theta_{2})^{\top}H_{\theta_{2}}(\theta_{1}-\theta_{2})$
$\displaystyle\leq\epsilon_{1}+\gamma\|\theta_{1}\|_{2}^{2}+\mathcal{L}_{x}(\theta_{1})d(\mathcal{D}_{1},\mathcal{D}_{2})$
(3)
$\displaystyle\frac{1}{2}\lambda_{\min}(H_{\theta_{2}})\|\theta_{1}-\theta_{2}\|^{2}_{2}$
$\displaystyle\leq(\epsilon_{1}-\epsilon_{2})+\gamma(\|\theta_{1}\|_{2}^{2}-\|\theta_{2}\|_{2}^{2})+\mathcal{L}_{x}(\theta_{1})d(\mathcal{D}_{1},\mathcal{D}_{2})$
(4)

To derive the second step we have used the following fact that a lower bound
on the quadratic term is given by the lowest eigenvalue i.e.,
$\arg\min_{\theta;\|\theta\|=1}\theta^{\top}H\theta=\lambda_{\min}\implies\theta^{\top}H\theta\geq\lambda_{\min}\|\theta\|^{2}$.
We now observe the following relationship for $\lambda_{min}$:

$\displaystyle\lambda_{\min}(\nabla_{\theta}^{2}\ell_{reg}(\theta_{2}))=\lambda_{\min}(\nabla_{\theta}^{2}\ell(\theta_{2}))+2\gamma\geq
2\gamma$ (5)

as $\nabla_{\theta}^{2}\ell(\theta_{2})$ is positive definite due to global
optimality. Putting eqn. 5 in eqn. 4, we have the final result. Here, the
“problem dependent constant”
$C=\sqrt{\frac{\ell_{\mathcal{D}_{1}}(\theta_{1})-\ell_{\mathcal{D}_{2}}(\theta_{2})}{\gamma}}$,
which we can expect to be quite small in practice if the optimal values
$\ell_{\mathcal{D}_{1}}(\theta_{1})\approx\ell_{\mathcal{D}_{2}}(\theta_{2})$.
Thus in practice we can assume $C$ to be very small.

∎

###  A.2 Bounding the explanation shift

######  Lemma 3.

The parameter-input Lipschitz has the following property:

$\displaystyle\mathbb{E}_{\mathbf{x}\sim\mathcal{D}}\|\nabla_{x}f(x;\theta_{1})-\nabla_{x}f(x;\theta_{2})\|_{2}\leq\leavevmode\nobreak\
\mathbb{E}_{\theta\in\Theta}\mathbb{E}_{\mathbf{x}\in\mathcal{D}}\|\nabla_{\theta}\nabla_{\mathbf{x}}f(x,\theta)\|_{2}\times\|\theta_{2}-\theta_{1}\|_{2}$

where
$\Theta=\\{\theta_{\lambda}=\lambda\theta_{2}+(1-\lambda)\theta_{1}\mid\lambda\in[0,1]\\}$

###### Proof.

The proof follows from the fundamental theorem of integral calculus. Let
$g(\theta,\mathbf{x})=\nabla_{\mathbf{x}}f(\mathbf{x},\theta)$ for
convenience.

$\displaystyle g(\theta_{2},\mathbf{x})-g(\theta_{1},\mathbf{x})$
$\displaystyle=\left(\int_{\lambda=0}^{1}\nabla_{\theta}g(\theta_{\lambda},\mathbf{x})d\lambda\right)^{\top}(\theta_{2}-\theta_{1})$
$\displaystyle\|g(\theta_{2},\mathbf{x})-g(\theta_{1},\mathbf{x})\|_{2}$
$\displaystyle\leq\|\int_{\lambda=0}^{1}\nabla_{\theta}g(\theta_{\lambda},\mathbf{x})d\lambda\|_{2}\|\theta_{2}-\theta_{1}\|_{2}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ (\text{Cauchy Schwartz})$
$\displaystyle\leq\left(\int_{\lambda=0}^{1}\|\nabla_{\theta}g(\theta_{\lambda},\mathbf{x})\|_{2}d\lambda\right)\|\theta_{2}-\theta_{1}\|_{2}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ (\text{Jensen's
inequality})$
$\displaystyle\mathbb{E}_{\mathbf{x}\in\mathcal{D}}\|g(\theta_{2},\mathbf{x})-g(\theta_{1},\mathbf{x})\|_{2}$
$\displaystyle\leq\left(\int_{\lambda=0}^{1}\mathbb{E}_{\mathcal{D}}\|\nabla_{\theta}g(\theta_{\lambda},\mathbf{x})\|_{2}d\lambda\right)\|\theta_{2}-\theta_{1}\|_{2}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ (\text{Swapping
expectation and integral})$
$\displaystyle\mathbb{E}_{\mathbf{x}\in\mathcal{D}}\|g(\theta_{2},\mathbf{x})-g(\theta_{1},\mathbf{x})\|_{2}$
$\displaystyle\leq\mathbb{E}_{\Theta}\mathbb{E}_{\mathcal{D}}\|\nabla_{\theta}g(\theta,\mathbf{x})\|_{2}\|\theta_{2}-\theta_{1}\|_{2}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ (\text{defn of
expectation})$

∎

######  Theorem 4.

Assume that a 1-hidden layer neural network with weights $\theta$, and random
inputs $\mathbf{x}\sim\mathcal{N}(0,I)$444Covariance of I is chosen for
notational brevity. Further assume that we have use an activation function
$\sigma$ with well-defined second derivatives (e.g: softplus). For this case,
the parameter-input derivatives have the following form:

$\displaystyle
E_{\mathbf{x}}\|\nabla_{\theta}\nabla_{\mathbf{x}}\ell(\mathbf{x},\theta)\|_{2}\leq
E_{\mathbf{x}}\|\nabla_{\theta}\nabla_{\mathbf{x}}\ell(\mathbf{x},\theta)\|_{F}\leq\|\theta\|_{2}+\beta\phi(\theta)$

where $\beta$ is an maximum curvature, and $\phi(\theta)$ is the path-norm
[29] of the model.

###### Proof.

We derive this expression for the case of a scalar valued neural network
$f(\mathbf{x})=W_{2}\sigma(W_{1}\mathbf{x})$, where
$W_{2}\in\mathbb{R}^{h\times d},W_{1}\in\mathbb{R}^{1\times h}$, where $h$ is
the number of hidden layers, and $d$ is the input dimensionality. Here,
$\sigma:\mathbb{R}\rightarrow\mathbb{R}$ is a point-wise smooth non-linearity
with a well-defined curvature, such as softplus, but not ReLU. For such a non-
linearity, let
$\sigma^{\prime\prime}(x)\vcentcolon=\frac{\partial^{2}\sigma(x)}{\partial
x^{2}}\leq M$. For softplus, this $M=\beta$ hyper-parameter is implicit in its
definition [43].

By straightforward calculus, its gradient and the “parameter-input
derivatives” are given by

$\displaystyle\frac{\partial f(\mathbf{x})}{\partial\mathbf{x}_{i}}$
$\displaystyle=\sum_{k=1}^{h}W_{2}^{k}\sigma^{\prime}(W_{1}\mathbf{x})^{k}W_{1}^{i,k}$
$\displaystyle\frac{\partial^{2}f(\mathbf{x})}{\partial\mathbf{x}_{i}\partial
W_{2}^{j}}$ $\displaystyle=\sigma^{\prime}(W_{1}\mathbf{x})^{j}W_{1}^{i,j}$
$\displaystyle\frac{\partial^{2}f(\mathbf{x})}{\partial\mathbf{x}_{i}\partial
W_{1}^{i,j}}$
$\displaystyle=W^{j}_{2}\sigma^{\prime}(W\mathbf{x})^{j}+W^{j}_{2}\sigma^{\prime\prime}(W_{1}\mathbf{x})^{j}W_{1}^{i,j}\mathbf{x}_{i}$
$\displaystyle\frac{\partial^{2}f(\mathbf{x})}{\partial\mathbf{x}_{i}\partial
W_{1}^{m,j}}$
$\displaystyle=W^{j}_{2}\sigma^{\prime\prime}(W_{1}\mathbf{x})^{j}W_{1}^{i,j}\mathbf{x}_{m}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ (\forall m\neq i)$

Now, assume that $\mathbf{x}\sim\mathcal{N}(0,I)$, i.e., the input is an
independent normal distribution. Let us now compute the squared terms for the
parameter-input derivatives to eventually be able to compute its Frobenius
norm.

$\displaystyle\left(\frac{\partial^{2}f(\mathbf{x})}{\partial\mathbf{x}_{i}\partial
W_{2}^{j}}\right)^{2}$
$\displaystyle\leq(\sigma^{\prime}(W_{1}\mathbf{x})^{j}W_{1}^{i,j})^{2}\leq(W_{1}^{i,j})^{2}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
(\sigma^{\prime}\leq 1\leavevmode\nobreak\ \text{for softplus})$
$\displaystyle\mathbb{E}_{\mathbf{x}}\left(\frac{\partial^{2}f(\mathbf{x})}{\partial\mathbf{x}_{i}\partial
W_{1}^{i,j}}\right)^{2}$
$\displaystyle=\mathbb{E}_{\mathbf{x}}(W^{j}_{2}\sigma^{\prime}(W\mathbf{x})^{j}+W^{j}_{2}\sigma^{\prime\prime}(W_{1}\mathbf{x})^{j}W_{1}^{i,j}\mathbf{x}_{i})^{2}$
$\displaystyle\leq(W^{j}_{2})^{2}+(W^{j}_{2}\sigma^{\prime\prime}(W_{1}\mathbf{x})^{j}W_{1}^{i,j})^{2}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
(\text{Using}\leavevmode\nobreak\ \mathbf{x}\sim\mathcal{N}(0,I))$
$\displaystyle\leq(W^{j}_{2})^{2}+\beta^{2}(W^{j}_{2}W_{1}^{i,j})^{2}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
(\sigma^{\prime\prime}\leq\beta)$
$\displaystyle\mathbb{E}_{\mathbf{x}}\left(\frac{\partial^{2}f(\mathbf{x})}{\partial\mathbf{x}_{i}\partial
W_{1}^{m,j}}\right)^{2}$
$\displaystyle\leq\beta^{2}(W^{j}_{2}W_{1}^{i,j})^{2}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ (\forall m\neq i)$

Computing the Frobenius norm, we have

$\displaystyle\mathbb{E}_{\mathbf{x}}\|\nabla_{\theta}\nabla_{\mathbf{x}}f(\mathbf{x})\|_{F}^{2}$
$\displaystyle=\sum_{i,j,k}\mathbb{E}_{X}\left(\frac{\partial^{2}f(\mathbf{x})}{\partial\mathbf{x}_{i}\partial
W_{1}^{j,k}}\right)+\sum_{i,j}\mathbb{E}_{X}\left(\frac{\partial^{2}f(\mathbf{x})}{\partial\mathbf{x}_{i}\partial
W_{2}^{j}}\right)$
$\displaystyle\leq\left(\sum_{j,k}(W_{2}^{k})^{2}+\sum_{j=i}\beta^{2}(W_{2}^{k}W_{1}^{j,k})^{2}+\sum_{j\neq
i}\beta^{2}(W_{2}^{k}W_{1}^{j,k})^{2}\right)+\sum_{i,j}(W_{1}^{i,j})^{2}$
$\displaystyle\leq\left(\sum_{k}(W_{2}^{k})^{2}+\sum_{i,j}(W_{1}^{i,j})^{2}\right)+\beta^{2}\sum_{j,k}(W_{2}^{k}W_{1}^{j,k})^{2}$
$\displaystyle\leq\|\theta\|^{2}+\beta^{2}\phi^{2}(\theta)$

Here, we use the fact that
$\phi(\theta)=\sqrt{\sum_{j,k}(W_{2}^{k}W_{1}^{j,k})^{2}}$ is the 2-path-norm
[29] for the purpose of characterizing generalization in neural networks.
Finally, we re-write the above using square roots:

$\displaystyle\mathbb{E}_{\mathbf{x}}\|\nabla_{\theta}\nabla_{\mathbf{x}}f(\mathbf{x})\|_{F}$
$\displaystyle\leq\sqrt{\mathbb{E}_{\mathbf{x}}\|\nabla_{\theta}\nabla_{\mathbf{x}}f(\mathbf{x})\|_{F}^{2}}$
$\displaystyle\leq\sqrt{\|\theta\|_{2}^{2}+\beta^{2}\phi^{2}(\theta)}$
$\displaystyle\leq\|\theta\|_{2}+\beta\phi(\theta)$

∎

##  Appendix B Additional experiments

##### Additional data from §4.1

The specific hyperparameters that we use to a) train base models, and b) fine-
tune these models, are shown in Table 7. These hyperparameters are chosen such
that a minimum can be attained by the optimizer (SGD), as is required to
verify the theory. The learning rates shown in the fine-tuning process for
HELOC and Adult (indicated by a star in Table 7) are multiplied by a scalar
that is dependent on the amount of noise added to the dataset. Intuitively,
adding more noise causes a larger increase in loss, so we increase the
learning rate used during fine-tuning in order to converge to the new minimum
within 50 epochs. The alternative approach, used in the WHO experiments, is to
use more epochs during fine-tuning, and instead start with a constant learning
rate that is higher but decays more, such that we again see convergence.

Table 7: Details of (dataset specific) hyperparameters used in §4.1 fine-
tuning. Base models are trained for a large enough number of epochs to ensure
convergence to a minimum. Decay value indicates the amount of decay applied to
the learning rate, and step size indicates the number of epochs between
applying decay.

| WHO | HELOC | Adult
---|---|---|---
| Base | Fine-Tune | Base | Fine-Tune | Base | Fine-Tune
Epochs | 4000 | 1000 | 500 | 50 | 500 | 50
Learning Rate | 0.5 | 0.1 | 0.5 | 0.01* | 0.5 | 0.01*
Decay Value | 0.9 | 0.9 | 0.95 | 0.95 | 0.95 | 0.95
Step Size | 40 | 40 | 100 | 100 | 100 | 100

* these learning rates vary as described in this appendix.

Accuracy data for WHO in the fine-tuning experiments is shown in Table 8.
Larger weight decay values and lower model curvatures do reduce train
accuracy, as these modifications make it more difficult to overfit the model
in reaching an exact minimum. However, test accuracy remains very similar
across techniques. For the experiments that add synthetic noise, the final
accuracy after fine-tuning is consistent regardless of how much noise we add,
up to a standard deviation of 0.1. The accuracy is in the mid-70% range for
HELOC and the mid-80% range for Adult.

Figure 7 shows the effect of weight decay and curvature on the Adult dataset
for the fine-tuning experiments (omitted in the main text). In the case of
curvature, the trend towards higher gradient stability given less curvature is
not very pronounced, which we attribute to the difficulty associated with
training a model to an exact minimum. In the presence of larger noise, the
closest minimum in the loss landscape during fine-tuning may also shift to a
different minimum.

Table 8: Accuracy for fine-tuning in the WHO dataset. Confidence intervals are the middle 50% of values. For softplus, $\gamma=0.001$. ReLU / $\gamma$=0 |  ReLU / $\gamma$=0.001 |  ReLU / $\gamma$=0.01 |  SP / $\beta=10$ |  SP / $\beta=5$
---|---|---|---|---
Train | Test | Train | Test | Train | Test | Train | Test | Train | Test
99.1 ±0.1 | 92.8±0.5 | 98.9 ±0.2 | 92.9 ±0.4 | 96.5±0.2 | 92.7 ±0.3 | 97.2 ±0.2 | 93.2 ±0.2 | 95.3±0.1 | 91.8±0.3

Figure 7: Effect of weight decay value (left) and model curvature (right) on
gradient and parameter stability for the Adult dataset. All models use ReLU
activation. Shifted models were trained via fine-tuning the original model.
The x-axis is the size of the data shift, represented by the standard
deviation $\sigma$ of Gaussian $\mathcal{N}(0,\sigma^{2})$ noise.

##### Additional data from §4.2

For a given dataset, the hyperparameters for the retraining experiments are
fixed between base models and retrained models in this section. In the HELOC
and Adult datasets, we train any particular model for 30 epochs with a
learning rate of 0.2, and no decay. For the WHO dataset, we use 80 epochs, an
initial learning rate 0.8, with decay value 0.8 and step size 10 (i.e. the
learning rate is multiplied by 0.8 every 10 epochs).

Accuracy data for the retraining experiments is in Table 9 for WHO, and Figure
8 for HELOC and Adult. In the latter, the base model accuracy is essentially
identical to the lower extreme on the x-axis (when very little/no noise has
been added to the training data). In all datasets, while there exists some
disparity between training and test data, this is much smaller than in the
fine-tuning experiments. Importantly, test accuracies remain approximately
constant across model types (e.g. lower or higher weight decay or curvature).
We see a slight drop in accuracy for the Adult dataset using the softplus
model with $\beta=2$. Despite how similar the test accuracies of all shifted
models are, we observe in the main text and in the following figures that
gradient stability across test inputs can still vary greatly, independent of
test accuracy.

The left and right sections of Figure 9 show how parameters change given
different weight decays and curvatures, respectively, as the data shift grows.
This illustrates the same trends that we see in the main text, namely, that
increasing weight decay or reducing curvature results in a smaller change in
parameters, while adding more noise to the data results in a larger change.

Figure 8: Accuracy for retraining (HELOC and Adult datasets). Confidence intervals are the middle 50% of values. Table 9: Accuracy for of base models (upper row) and shifted models (lower row) in retraining experiments. Confidence intervals are the middle 50% of values. For softplus, $\gamma=0.001$. |  ReLU / $\gamma$=0 |  ReLU / $\gamma$=0.001 |  ReLU / $\gamma$=0.01 |  SP / $\beta=10$ large∗ |  SP / $\beta=5$
---|---|---|---|---|---
| Train | Test | Train | Test | Train | Test | Train | Test | Train | Test
Base | 98.3 ±0.2 | 93.3 ±0.2 | 97.7±0.2 | 93.9±0.1 | 94.7±0.0 | 93.4±0.4 | 95.8 ±0.0 | 92.9 ±0.4 | 94.5±0.0 | 91.5±0.2
Shifted | 98.0 ±0.2 | 93.7 ±0.4 | 97.4±0.2 | 93.7±0.4 | 94.7±0.2 | 92.5±0.3 | 95.7 ±0.0 | 93.2±0.2 | 94.7 ±0.1 | 92.6 ±0.2

Figure 9: Effect of weight decay (left) and model curvature (right) on
parameter stability when retraining on the HELOC and Adult datasets. The
x-axis is the size of the data shift, represented by the standard deviation
$\sigma$ of $\mathcal{N}(0,1)$ noise.

Figures 10 and 11 show the Top-5 consistency scores of HELOC and Adult,
respectively, each for Saliency and SmoothGrad techniques. We show the same
for LIME and K.SHAP in Figure 12 for HELOC and Figure 13 for Adult. Overall,
SmoothGrad demonstrates consistently strong performance with respect to
explanation stability, while LIME and K.SHAP show poorer performance with
higher variability. For all explanation techniques, we see the same consistent
trends with regards to weight decay and curvature, though the effect of each
of these is occasionally less distinct in the cases of LIME and K.SHAP (i.e.
when the explanation technique is inherently less stable, which we attribute
mostly to the number of samples used in LIME, the countermeasures we propose
are sometimes less effective).

Figure 10: Top-5 consistency. From the top, HELOC SA, HELOC CDC, HELOC SSA
(the graph for HELOC SA appears in the main text). Each row shows the effects
of weight decay and curvature as data shift grows for salience and SmoothGrad
top-5 metrics. Confidence intervals represent the middle 50% of values.

Figure 11: Top-5 consistency. From the top, Adult SA, Adult CDC, Adult SSA.
Each row shows the effects of weight decay and curvature as data shift grows
for salience and SmoothGrad top-5 metrics. Confidence intervals represent the
middle 50% of values.

Figure 12:  Top-5 consistency. From the top, HELOC SA, HELOC CDC, HELOC SSA.
Each row shows the effects of weight decay and curvature as data shift grows
for LIME and K.SHAP top-5 metrics. Confidence intervals represent the middle
50% of values.

Figure 13:  Top-5 consistency. From the top, Adult SA, Adult CDC, Adult SSA.
Each row shows the effects of weight decay and curvature as data shift grows
for LIME and K.SHAP top-5 metrics. Confidence intervals represent the middle
50% of values.

Table 10 shows the explanation stability scores for various metrics and
explanation techniques for the WHO dataset. We make the following
observations: first, for salience and SHAP (and to a lesser extent, LIME),
using a weight decay of 0.01 is much more stable than using smaller weight
decays. SmoothGrad is already very stable, so increasing the weight decay has
no effect. Second, lowering the model curvature greatly increases the
stability of CDC and SSA across explanation techniques, but has less effect on
SA. This trend suggests that increasing the model curvature has limited impact
on what features are in the top-K (i.e., what SA measures), but has a large
impact on keeping the same sign for influential features (even though they
might not be in the top K), which is what CDC measures.

Table 10: Explanation Stability for WHO |  |  ReLU ($\gamma=0$) |  ReLU ($\gamma$=0.001) |  ReLU ($\gamma$=0.01) |  SP ($\beta=10$) |  SP ($\beta=5$)
---|---|---|---|---|---|---
SA | Salience | 0.63±0.01 | 0.64±0.01 | 0.73±0.01 | 0.72±0.02 | 0.77±0.02
| SmoothGrad | 0.94±0.00 | 0.94±0.00 | 0.95±0.00 | 0.95±0.00 | 0.95±0.00
| LIME | 0.59±0.02 | 0.59±0.03 | 0.70±0.03 | 0.65±0.03 | 0.69±0.02
| SHAP | 0.52±0.01 | 0.56±0.02 | 0.67±0.03 | 0.61±0.02 | 0.65±0.02
CDC | Salience | 0.80±0.04 | 0.83±0.04 | 0.95±0.01 | 0.94±0.02 | 0.97±0.01
| SmoothGrad | 0.94±0.00 | 0.94±0.01 | 0.95±0.00 | 0.95±0.00 | 0.96±0.00
| LIME | 0.52±0.03 | 0.52±0.05 | 0.57±0.05 | 0.58±0.05 | 0.65±0.04
| SHAP | 0.62±0.02 | 0.70±0.04 | 0.87±0.03 | 0.77±0.03 | 0.86±0.04
SSA | Salience | 0.19±0.04 | 0.18±0.03 | 0.30±0.02 | 0.29±0.01 | 0.39±0.03
| SmoothGrad | 0.91±0.00 | 0.91±0.00 | 0.91±0.00 | 0.91±0.00 | 0.92±0.00
| LIME | 0.27±0.02 | 0.25±0.02 | 0.32±0.03 | 0.30±0.03 | 0.35±0.03
| SHAP | 0.19±0.04 | 0.21±0.01 | 0.29±0.04 | 0.23±0.03 | 0.28±0.03
Figure 14: Various decision tree depths for the HELOC dataset. Similarity of
LIME explanations between base models (original dataset) and retrained models
(averaged across 10 random noise seeds for each value of $\sigma$). Left to
right show Top-5 SA, CDC, and SSA metrics.

##### Additional experiments with tree-based models

Our theoretical results primarily focus on differentiable models, which
encompass popular classes like neural networks, though the trends we have
identified may not be universally applicable to all model classes. That said,
it is worth noting that prior works analyzing the behavior and stability of ML
models and explanation methods also typically rely on linear or differentiable
models [9, 33]. We show that differentiable models with softplus activation
functions and higher weight decay values, both of which can be viewed as a way
to limit the functional complexity of the model, tend to promote explanation
stability. Thus, we posit that explanations of tree-based models could
similarly be made more robust by imposing analogous limits on their
complexity.

We empirically evaluate explanation stability when retraining decision trees,
random forests, and XGBoost models, controlling for tree depth in the first
two cases and the $\lambda$ regularization term in XGBoost. Higher $\lambda$
values indicate greater regularization and lower complexity. We evaluated the
stability of post-hoc explanations on the HELOC dataset with synthetic noise
$\sim\mathcal{N}(0,\sigma^{2})$, conducting empirical evaluations using LIME
and SHAP. Please refer to §4 for detailed information on the setup of
synthetic noise experiments. Interestingly, our findings show mixed results.
While larger values of $\lambda$ (corresponding to higher regularization
constants) or smaller tree depths tend to yield more robust explanations in
some cases, the outcomes vary for different explanation methods and dataset
shifts.

Figures 14 and 15 demonstrate the impact of decision tree depth on LIME and
SHAP, respectively, for Top-5 SA, CDC, and SSA metrics. Decision trees with
lower depth exhibit improved explanation stability for smaller noise values,
though the effects can break down as the dataset shifts by larger amounts.
While these results support our findings for neural networks to a moderate
degree, further research is required to investigate the behaviour that occurs
for larger noise values (recall that our theory applies to small dataset
shifts).

Figures 16 and 17 demonstrate the impact of random forest tree depth on LIME
and SHAP, respectively, for Top-5 SA, CDC, and SSA metrics. Observe how the
explanation stability of random forests is in general much higher than for
decision trees– this would support the idea that the ensemble approach of
random forests, which reduces functional diversity and complexity through
averaging across multiple decision trees, indeed promotes explanation
stability. However, the relative effects of depth are mixed. For LIME, at
depth 5, explanations are relatively less stable, though for higher depths,
explanations retained strong similarity until sufficiently large synthetic
shifts were added ($\sigma>0.1$). For SHAP, we observe the trends that we
expect (higher depth reduces stability), though in the case of CDC this was
reversed, albeit at very high stability values. Again, further research is
required to investigate the behaviour that occurs for larger shifts.

Figures 18 and 19 demonstrate the impact of the XGBoost $\ell_{2}$
regularization term $\lambda$ on LIME and SHAP, respectively, for Top-5 SA,
CDC, and SSA metrics. For LIME, the outcomes appear to be quite volatile,
showing a relatively unpredictable range of optimal $\lambda$ values for
explanation stability. Interestingly, the added synthetic noise didn’t
significantly destabilize the explanations, and in some instances, it
surprisingly seemed to bolster stability. This counter-intuitive finding
underscores the need for further investigation to understand the underlying
phenomena driving these trends (though LIME itself may be the significant
driving factor). In contrast, for SHAP, the results are much clearer,
revealing that higher regularization indeed promotes more stable explanations.
Nonetheless, as anticipated, the stability decreases as the dataset undergoes
more significant shifts. Notably, the performance of SHAP appeared to be
superior to LIME.

We want to emphasize that these results are preliminary - more work is needed
to generalize these results to larger experiments, to extend these results to
real-world (not synthetic) data shifts, and to develop theory for when
explanations are be more stable for tree-based models. Notably, we did not
perform experiments with the WHO dataset (our real-world shift example)
because the results of performing a single experiment were too noisy (with
synthetic shifts, the random initialization of the noise can be used to
average over multiple trials).

Figure 15: Various decision tree depths for the HELOC dataset. Similarity of
SHAP explanations between base models (original dataset) and retrained models
(averaged across 10 random noise seeds for each value of $\sigma$). Left to
right show Top-5 SA, CDC, and SSA metrics. Figure 16: Various random forest
depths for the HELOC dataset. Similarity of LIME explanations between base
models (original dataset) and retrained models (averaged across 10 random
noise seeds for each value of $\sigma$). Left to right show Top-5 SA, CDC, and
SSA metrics. Figure 17: Various random forest depths for the HELOC dataset.
Similarity of SHAP explanations between base models (original dataset) and
retrained models (averaged across 10 random noise seeds for each value of
$\sigma$). Left to right show Top-5 SA, CDC, and SSA metrics. Figure 18:
Various XGBoost $\lambda$ values for the HELOC dataset. Similarity of LIME
explanations between base models (original dataset) and retrained models
(averaged across 10 random noise seeds for each value of $\sigma$). Left to
right show Top-5 SA, CDC, and SSA metrics. Figure 19: Various XGBoost
$\lambda$ values for the HELOC dataset. Similarity of SHAP explanations
between base models (original dataset) and retrained models (averaged across
10 random noise seeds for each value of $\sigma$). Left to right show Top-5
SA, CDC, and SSA metrics.