# On the Robustness of Removal-Based Feature Attributions

Chris Lin
University of Washington
<EMAIL_ADDRESS>&Ian Covert∗†
Stanford University
<EMAIL_ADDRESS>&Su-In Lee
University of Washington
<EMAIL_ADDRESS>Equal contribution. †Work done while at the
University of Washington.

###### Abstract

To explain predictions made by complex machine learning models, many feature
attribution methods have been developed that assign importance scores to input
features. Some recent work challenges the robustness of these methods by
showing that they are sensitive to input and model perturbations, while other
work addresses this issue by proposing robust attribution methods. However,
previous work on attribution robustness has focused primarily on gradient-
based feature attributions, whereas the robustness of removal-based
attribution methods is not currently well understood. To bridge this gap, we
theoretically characterize the robustness properties of removal-based feature
attributions. Specifically, we provide a unified analysis of such methods and
derive upper bounds for the difference between intact and perturbed
attributions, under settings of both input and model perturbations. Our
empirical results on synthetic and real-world data validate our theoretical
results and demonstrate their practical implications, including the ability to
increase attribution robustness by improving the model’s Lipschitz regularity.

##  1 Introduction

In recent years, machine learning has shown great promise in a variety of
real-world applications. An obstacle to its widespread deployment is the lack
of transparency, an issue that has prompted a wave of research on
interpretable or explainable machine learning [56, 67, 53, 60, 45, 51]. One
popular way of explaining a machine learning model is feature attribution,
which assigns importance scores to input features of the model [56, 53, 60,
45, 51]. Many feature attribution methods can be categorized as either
gradient-based methods that compute gradients of model predictions with
respect to input features [6], or removal-based methods that remove features
to quantify each feature’s influence [17]. While substantial progress has been
made and current tools are widely used for model debugging and scientific
discovery [10, 18, 37, 50], some important concerns remain. Among them is
unclear robustness properties: feature attributions are vulnerable to
adversarial attacks, and even without an adversary appear unstable under small
changes to the input or model.

As a motivating example, consider the work of Ghorbani et al. [28], which
shows that minor perturbations to an image can lead to substantially different
feature attributions while preserving the original prediction. Similar to
adversarial examples designed to alter model predictions [61, 29, 43], this
phenomenon violates the intuition that explanations should be invariant to
imperceptible changes. A natural question, and the focus of this work, is
therefore: When can we guarantee that feature attributions are robust to small
changes in the input or small changes in the model? In answering this
question, we find that previous theoretical work on this topic has primarily
focused on gradient-based methods. Hence, in this work we provide a unified
analysis to characterize the robustness of removal-based feature attributions,
under both input and model perturbations, and considering a range of existing
methods in combination with different approaches to implement feature removal.

Related work. Recent work has demonstrated that gradient-based attributions
for neural networks are susceptible to small input perturbations that can lead
to markedly different attributions [21, 28, 41]. Theoretical analyses have
established that this issue relates to model smoothness, and approaches like
stochastic smoothing of gradients, weight decay regularization, smoothing
activation functions, and Hessian minimization have been shown to generate
more robust gradient-based attributions [21, 22, 65, 4]. As for the robustness
of removal-based attributions under input perturbations, Alvarez-Melis and
Jaakkola [5] empirically assess the robustness of LIME, SHAP, and Occlusion
with the notion of Lipschitz continuity in local neighborhoods. Khan et al.
[40] provide theoretical guarantees on the robustness of SHAP, RISE, and
leave-one-out attributions using Lipschitz continuity, in the specific setting
where held-out features are replaced with zeros. Agarwal et al. [3] study the
robustness of discrete attributions for graph neural networks and establish
Lipschitz upper bounds for the robustness of GraphLIME and GraphMASK. Our
analysis of input perturbations is most similar to Khan et al. [40], but we
generalize this work by considering alternative techniques for both feature
removal and importance summarization (e.g., LIME’s weighted least squares
approach [53]).

Other works have considered model perturbations and manipulated neural
networks to produce targeted, arbitrary gradient-based attributions [7, 34].
Anders et al. [7] explain this phenomenon through the insight that neural
network gradients are underdetermined with many degrees of freedom to exploit,
and they demonstrate better robustness when gradient-based attributions are
projected onto the tangent space of the data manifold. As for removal-based
attributions, it has been demonstrated that models can be modified to hide
their reliance on sensitive features in LIME and SHAP [58, 20]; however, Frye
et al. [27] show that such hidden features can be discovered if features are
removed using their conditional distribution, suggesting a key role for the
feature removal approach in determining sensitivity to such attacks. Our work
formalizes these results through the lens of robustness under model
perturbations, and provides new theoretical guarantees.

Contribution. The main contribution of this work is to develop a comprehensive
characterization of the robustness properties of removal-based feature
attributions, focusing on two notions of robustness that have received
attention in the literature: stability under input changes and stability under
model changes. Our specific contributions are the following: (1) We provide
theoretical guarantees for the robustness of model predictions with the
removal of arbitrary feature subsets, under both input and model
perturbations. (2) We analyze the robustness of techniques for summarizing
each feature’s influence (e.g., Shapley values) and derive best- and worst-
case robustness within several classes of such techniques. (3) We combine the
analyses above to prove unified attribution robustness properties to both
input and model perturbations, where changes in removal-based attributions are
controlled by the scale of perturbations in either input space or function
space. (4) We validate our theoretical results and demonstrate practical
implications using synthetic and real-world datasets.

##  2 Background

Here, we introduce the notation used in the paper and review removal-based
feature attributions.

###  2.1 Notation

Let $f:\mathbb{R}^{d}\mapsto\mathbb{R}$ be a model whose predictions we seek
to explain. We consider that the input variable ${\bm{x}}$ consists of $d$
separate features, or ${\bm{x}}=({\bm{x}}_{1},\ldots,{\bm{x}}_{d})$. Given an
index set $S\subseteq[d]\equiv\\{1,\ldots,d\\}$, we define the corresponding
feature subset as ${\bm{x}}_{S}\equiv\\{{\bm{x}}_{i}:i\in S\\}$. The power set
of $[d]$ is denoted by ${\mathcal{P}}(d)$, and we let
$\bar{S}\equiv[d]\setminus S$ denote the set complement. We consider the
model’s output space to be one-dimensional, which is not restrictive because
attributions are typically calculated for scalar predictions (e.g., the
probability for a given class). We use the bold symbols
${\bm{x}},{\bm{x}}_{S}$ to denote random variables, the symbols $x,x_{S}$ to
denote specific values, and $p({\bm{x}})$ to represent the data distribution
with support on ${\mathcal{X}}\subseteq\mathbb{R}^{d}$. We assume that all
explanations are generated for inputs such that $x\in{\mathcal{X}}$.

###  2.2 Removal-based feature attributions

Our work focuses on a class of methods known as removal-based explanations
[17]. Intuitively, these are algorithms that remove subsets of inputs and
summarize how each feature affects the model. This framework describes a large
number of methods that are distinguished by two main implementation choices:
(1) how feature information is removed from the model, and (2) how the
algorithm summarizes each feature’s influence.111The framework also considers
the choice to explain different model behaviors (e.g., the loss calculated
over one or more examples) [17], but we focus on methods that explain
individual predictions. This perspective shows a close connection between
methods like leave-one-out [67], LIME [53] and Shapley values [45], and
enables a unified analysis of their robustness properties. Below, we describe
the two main implementation choices in detail.

Feature removal. Most machine learning models require all feature values to
generate predictions, so we must specify a convention for depriving the model
of feature information. We denote the prediction given partial inputs by
$f(x_{S})$. Many implementations have been discussed in the literature [17],
but we focus here on three common choices. Given a set of observed values
$x_{S}$, these techniques can all be viewed as averaging the prediction over a
distribution $q({\bm{x}}_{\bar{S}})$ for the held-out feature values:

$f(x_{S}):=\mathbb{E}_{q({\bm{x}}_{\bar{S}})}\left[f(x_{S},{\bm{x}}_{\bar{S}})\right]=\int
f(x_{S},x_{\bar{S}})q(x_{\bar{S}})dx_{\bar{S}}.$ (1)

Specifically, the three choices we consider are:

  * •

(Baseline values) Given a baseline input $b\in\mathbb{R}^{d}$, we can set the
held-out features to their corresponding values $b_{\bar{S}}$ [59]. This is
equivalent to letting $q({\bm{x}}_{\bar{S}})$ be a Dirac delta centered at
$b_{\bar{S}}$.

  * •

(Marginal distribution) Rather than using a single replacement value, we can
average across values sampled from the input’s marginal distribution, or let
$q({\bm{x}}_{\bar{S}})=p({\bm{x}}_{\bar{S}})$ [45].

  * •

(Conditional distribution) Finally, we can average across replacement values
while conditioning on the available features, which is equivalent to setting
$q({\bm{x}}_{\bar{S}})=p({\bm{x}}_{\bar{S}}\mid x_{S})$ [27].

The first two options are commonly implemented in practice because they are
simple to estimate, but the third choice is viewed by some work as more
informative to users [1, 27, 17]. The conditional distribution approach
requires more complex and error-prone estimation procedures (see [12] for a
discussion), but we assume here that all three versions of $f(x_{S})$ can be
calculated exactly.

Summary technique. Given a feature removal technique that allows us to query
the model with arbitrary feature sets, we must define a convention for
summarizing each feature’s influence. This is challenging due to the
exponential number of feature sets, or because $|{\mathcal{P}}(d)|=2^{d}$.
Again, many techniques have been discussed in the literature [17], with the
simplest option comparing the prediction with all features included and with a
single feature missing [67]. We take a broad perspective here, considering a
range of approaches that define attributions as a linear combination of
predictions with different feature sets. These methods include leave-one-out
[67], RISE [51], LIME [53], Shapley values [45], and Banzhaf values [13],
among other possible options.

All of these methods yield per-feature attribution scores
$\phi_{i}(f,x)\in\mathbb{R}$ for $i\in[d]$. For example, the Shapley value
calculates feature attribution scores as follows [45]:

$\phi_{i}(f,x)=\frac{1}{d}\sum_{S\subseteq[d]\setminus\\{i\\}}\binom{d-1}{|S|}^{-1}\left(f(x_{S\cup\\{i\\}})-f(x_{S})\right).$
(2)

The specific linear combinations for the other approaches we consider are
shown in Table 3. The remainder of the paper uses the notation $f(x_{S})$ to
refer to predictions with partial information, and
$\phi(f,x)=[\phi_{1}(f,x),\ldots,\phi_{d}(f,x)]\in\mathbb{R}^{d}$ to refer to
feature attributions. We do not introduce separate notation to distinguish
between implementation choices, opting instead to make the relevant choices
clear in each result.

###  2.3 Problem formulation

The problem formulation in this work is straightforward: our goal is to
understand the stability of feature attributions under input perturbations and
model perturbations. Formally, we aim to study

  1. 1.

whether $\lVert\phi(f,x)-\phi(f,x^{\prime})\rVert$ is controlled by $\lVert
x-x^{\prime}\rVert$ (input perturbation), and

  2. 2.

whether $\lVert\phi(f,x)-\phi(f^{\prime},x)\rVert$ is controlled by $\lVert
f-f^{\prime}\rVert$ (model perturbation).

The following sections show how this can be guaranteed using certain distance
metrics, and under certain assumptions about the model and/or the data
distribution.

##  3 Preliminary results

As an intermediate step towards understanding explanation robustness, we first
provide results for sub-components of the algorithms. We begin by addressing
the robustness of predictions with partial information to input and model
perturbations, and we then discuss robustness properties of the summary
technique. Proofs for all results are in the Appendix.

Before proceeding, we introduce two assumptions that are necessary for our
analysis.

######  Assumption 1.

We assume that the model $f$ is globally $L$-Lipschitz continuous, or that we
have $|f(x)-f(x^{\prime})|\leq L\cdot\lVert x-x^{\prime}\rVert_{2}$ for all
$x,x^{\prime}\in\mathbb{R}^{d}$.

######  Assumption 2.

We assume that the model $f$ has bounded predictions, or that there exists a
constant $B$ such that $|f(x)|\leq B$ for all $x\in\mathbb{R}^{d}$.

The first assumption holds for most deep learning architectures, because they
typically compose a series of Lipschitz continuous layers [63]. The second
assumption holds with $B=1$ for classification models. These are therefore
mild assumptions, but they are discussed further in Section 6.

We also define the notion of a functional norm, which is useful for several of
our results.

######  Definition 1.

The $L^{p}$ norm for a function $g:\mathbb{R}^{d}\mapsto\mathbb{R}$ is defined
as $\lVert g\rVert_{p}\equiv(\int|g(x)|^{p}dx)^{1/p}$, where the integral is
taken over $\mathbb{R}^{d}$. $\lVert g\rVert_{p,{\mathcal{X}}^{\prime}}$
denotes the same integral taken over the domain
${\mathcal{X}}^{\prime}\subseteq\mathbb{R}^{d}$.

Our results adopt this to define a notion of distance between functions
$\lVert g-g^{\prime}\rVert_{p}$, focusing on the cases where $p=1$ for the
$L^{1}$ distance and $p=\infty$ for the Chebyshev distance.

###  3.1 Prediction robustness to input perturbations

Our goal here is to understand how the prediction function $f(x_{S})$ for a
feature set $x_{S}$ behaves under small input perturbations. We first consider
the case where features are removed using the baseline or marginal approach,
and we find that Lipschitz continuity is inherited from the original model.

######  Lemma 1.

When removing features using either the baseline or marginal approaches, the
prediction function $f(x_{S})$ for any feature set $x_{S}$ is $L$-Lipschitz
continuous:

$|f(x_{S})-f(x_{S}^{\prime})|\leq L\cdot\lVert
x_{S}-x_{S}^{\prime}\rVert_{2}\quad\forall\;x_{S},x_{S}^{\prime}\in\mathbb{R}^{|S|}.$

Next, we consider the case where features are removed using the conditional
distribution approach. We show that the continuity of $f(x_{S})$ depends not
only on the original model, but also on the similarity of the conditional
distribution for the held-out features.

######  Lemma 2.

When removing features using the conditional approach, the prediction function
$f(x_{S})$ for a feature set $x_{S}$ satisfies

$|f(x_{S})-f(x_{S}^{\prime})|\leq L\cdot\lVert
x_{S}-x_{S}^{\prime}\rVert_{2}+2B\cdot d_{TV}\Big{(}p({\bm{x}}_{\bar{S}}\mid
x_{S}),p({\bm{x}}_{\bar{S}}\mid x_{S}^{\prime})\Big{)},$

where the total variation distance is defined via the $L^{1}$ functional
distance as

$d_{TV}\Big{(}p({\bm{x}}_{\bar{S}}\mid x_{S}),p({\bm{x}}_{\bar{S}}\mid
x_{S}^{\prime})\Big{)}\equiv\frac{1}{2}\Big{\lVert}p({\bm{x}}_{\bar{S}}\mid
x_{S})-p({\bm{x}}_{\bar{S}}\mid x_{S}^{\prime})\Big{\rVert}_{1}.$

Lemma 2 does not immediately imply Lipschitz continuity for $f(x_{S})$,
because we have not bounded the total variation distance in terms of $\lVert
x_{S}-x_{S}^{\prime}\rVert_{2}$. To address this, we require the following
Lipschitz-like continuity property in the total variation distance.

######  Assumption 3.

We assume that there exists a constant $M$ such that for all $S\subseteq[d]$,
we have

$d_{TV}\Big{(}p({\bm{x}}_{\bar{S}}\mid x_{S}),p({\bm{x}}_{\bar{S}}\mid
x_{S}^{\prime})\Big{)}\leq M\cdot\lVert
x_{S}-x_{S}^{\prime}\rVert_{2}\quad\forall
x_{S},x_{S}^{\prime}\in\mathbb{R}^{|S|}.$

Intuitively, this says that the conditional distribution
$p({\bm{x}}_{\bar{S}}\mid x_{S})$ cannot change too quickly as a function of
the observed features. This property does not hold for all data distributions
$p({\bm{x}})$, and it may hold in some cases only for large $M$; in such
scenarios, the function $f(x_{S})$ is not guaranteed to change slowly.
However, there are cases where it holds. For example, we have $M=0$ if the
features are independent, and the following example shows a case where it
holds with dependent features.

######  Example 1.

For a Gaussian random variable ${\bm{x}}\sim{\mathcal{N}}(\mu,\Sigma)$ with
mean $\mu\in\mathbb{R}^{d}$ and covariance $\Sigma\in\mathbb{R}^{d\times d}$,
Assumption 3 holds with
$M=\frac{1}{2}\sqrt{\lambda_{\max}(\Sigma^{-1})-\lambda_{\min}(\Sigma^{-1})}$.
If ${\bm{x}}$ is assumed to be standardized, this captures the case where
independent features yield $M=0$. Intuitively, it also means that if one
dimension is roughly a linear combination of the others, or if
$\lambda_{\max}(\Sigma^{-1})$ is large, the conditional distribution can
change quickly as a function of the conditioning variables.

Now, assuming that this property holds for $p({\bm{x}})$, we show that we can
improve upon Lemma 2.

######  Lemma 3.

Under Assumption 3, the prediction function $f(x_{S})$ defined using the
conditional approach is Lipschitz continuous with constant $L+2BM$ for any
feature set $x_{S}$.

Between Lemmas 1 and 3, we have established that the function $f(x_{S})$
remains Lipschitz continuous for any feature set $x_{S}$, although in some
cases with a larger constant that depends on the data distribution. These
results are summarized in Table 2. In Appendix B we show that our analysis can
be extended to account for sampling in eq. 1 when the expectation is not
calculated exactly.

Table 1: Lipschitz constants induced by each feature removal technique. Baseline | Marginal | Conditional
---|---|---
$L$ | $L$ | $L+2BM$

Table 2: Relevant functional distances for each feature removal technique. Baseline | Marginal | Conditional
---|---|---
$\lVert f-f^{\prime}\rVert_{\infty}$ | $\lVert f-f^{\prime}\rVert_{\infty}$ | $\lVert f-f^{\prime}\rVert_{\infty,{\mathcal{X}}}$

###  3.2 Prediction robustness to model perturbations

Our next goal is to understand how the function $f(x_{S})$ for a feature set
$x_{S}$ behaves under small changes to the model. The intuition is that if two
models make very similar predictions with all features, they should continue
to make similar predictions with a subset of features. We first derive a
general result involving the proximity between two models $f$ and $f^{\prime}$
within a subdomain.

######  Lemma 4.

For two models $f,f^{\prime}:\mathbb{R}^{d}\mapsto\mathbb{R}$ and a subdomain
${\mathcal{X}}^{\prime}\subseteq\mathbb{R}^{d}$, the prediction functions
$f(x_{S}),f^{\prime}(x_{S})$ for any feature set $x_{S}$ satisfy

$|f(x_{S})-f^{\prime}(x_{S})|\leq\lVert
f-f^{\prime}\rVert_{\infty,{\mathcal{X}}^{\prime}}\cdot
Q_{x_{S}}({\mathcal{X}}^{\prime})+2B\cdot\big{(}1-Q_{x_{S}}({\mathcal{X}}^{\prime})\big{)},$

where $Q_{x_{S}}({\mathcal{X}}^{\prime})$ is the probability of imputed
samples lying in ${\mathcal{X}}^{\prime}$ based on the distribution
$q({\bm{x}}_{\bar{S}})$:

$Q_{x_{S}}({\mathcal{X}}^{\prime})\equiv\mathbb{E}_{q({\bm{x}}_{\bar{S}})}\left[{\mathbb{I}}\left\\{(x_{S},{\bm{x}}_{\bar{S}})\in{\mathcal{X}}^{\prime}\right\\}\right].$

The difference in predictions therefore depends on the distance between the
models only within the subdomain ${\mathcal{X}}^{\prime}$, as well as the
likelihood of imputed sampled lying in this subdomain. This result illustrates
a point shown by Slack et al. [58]: that two models which are equivalent on a
small subdomain, or even on the entire data manifold ${\mathcal{X}}$, can lead
to different attributions if we use a removal technique that yields a low
value for $Q_{x_{S}}({\mathcal{X}}^{\prime})$. In fact, when
$Q_{x_{S}}({\mathcal{X}}^{\prime})\to 0$, Lemma 4 reduces to a trivial bound
$|f(x_{S})-f^{\prime}(x_{S})|\leq 2B$ that follows directly from Assumption 2.

The general result in Lemma 4 allows us to show two simpler ones. In both
cases, we choose the subdomain ${\mathcal{X}}^{\prime}$ and removal technique
$q({\bm{x}}_{\bar{S}})$ so that $Q_{x_{S}}({\mathcal{X}}^{\prime})=1$.

######  Lemma 5.

When removing features using the conditional approach, the prediction
functions $f(x_{S}),f^{\prime}(x_{S})$ for two models $f,f^{\prime}$ and any
feature set $x_{S}$ satisfy

$|f(x_{S})-f^{\prime}(x_{S})|\leq\lVert
f-f^{\prime}\rVert_{\infty,{\mathcal{X}}}.$

The next result is similar but involves a potentially larger upper bound
$\lVert f-f^{\prime}\rVert_{\infty}\geq\lVert
f-f^{\prime}\rVert_{\infty,{\mathcal{X}}}$.

######  Lemma 6.

When removing features using the baseline or marginal approach, the prediction
functions $f(x_{S}),f^{\prime}(x_{S})$ for two models $f,f^{\prime}$ and any
feature set $x_{S}$ satisfy

$|f(x_{S})-f^{\prime}(x_{S})|\leq\lVert f-f^{\prime}\rVert_{\infty}.$

The relevant functional distances for Lemma 5 and Lemma 6 are summarized in
Table 2.

######  Remark 1.

Lemma 6 implies that if two models $f,f^{\prime}$ are functionally equivalent,
or $f(x)=f(x^{\prime})$ for all $x\in\mathbb{R}^{d}$, they remain equivalent
with any feature subset $x_{S}$. In Section 4, we will see that this implies
equal attributions for the two models—a natural property that, perhaps
surprisingly, is not satisfied by all feature attribution methods [55, 47],
and that has been described in the literature as implementation invariance
[60]. This property holds automatically for all removal-based methods.

######  Remark 2.

In contrast, Lemma 5 implies the same for models that are equivalent only on
the data manifold ${\mathcal{X}}\subseteq\mathbb{R}^{d}$. This is a less
stringent requirement to guarantee equal attributions, and it is perhaps
reasonable that models with equal predictions for all realistic inputs receive
equal attributions. To emphasize this, we propose distinguishing between a
notion of weak implementation invariance and strong implementation invariance,
depending on whether equal attributions are guaranteed when we have
$f(x)=f^{\prime}(x)$ everywhere ($x\in\mathbb{R}^{d}$) or almost everywhere
($x\in{\mathcal{X}}$).

###  3.3 Summary technique robustness

We previously focused on the effects of the feature removal choice, so our
goal is now to understand the role of the summary technique. Given a removal
approach that lets us query the model with any feature set, the summary
generates attribution scores $\phi(f,x)\in\mathbb{R}^{d}$, often using a
linear combination of the outputs $f(x_{S})$ for each $S\subseteq[d]$. We
formalize this in the following proposition.

Table 3: Summary technique linear operators used by various removal-based explanations. Summary | Method |  $A_{iS}$ ($i\in S$) |  $A_{iS}$ ($i\notin S$) | $\lVert A\rVert_{1,\infty}$ | $\lVert A\rVert_{2}$
---|---|---|---|---|---
Leave-one-out |  Occlusion [67] | ${\mathbb{I}}\left\\{|S|=d\right\\}$ | $-{\mathbb{I}}\left\\{|S|=d-1\right\\}$ | $2\sqrt{d}$ | $\sqrt{d+1}$
Shapley value |  SHAP [45] | $\frac{(|S|-1)!(d-|S|)!}{d!}$ | $-\frac{|S|!(d-|S|-1)!}{d!}$ | $2\sqrt{d}$ | $\sqrt{2/d}$
Banzhaf value |  Banzhaf [13] | $1/2^{d-1}$ | $-1/2^{d-1}$ | $2\sqrt{d}$ | $1/2^{d/2-1}$
Mean when included |  RISE [51] | $1/2^{d-1}$ | $0$ | $\sqrt{d}$ | $\sqrt{(d+1)/2^{d}}$
Weighted least squares |  LIME [53] |  Depends on implementation choices (see Appendix D)

######  Proposition 1.

The attributions for each method can be calculated by applying a linear
operator $A\in\mathbb{R}^{d\times 2^{d}}$ to a vector $v\in\mathbb{R}^{2^{d}}$
representing the predictions with each feature set, or

$\phi(f,x)=Av,$

where the linear operator $A$ for each method is listed in Table 3, and $v$ is
defined as $v_{S}=f(x_{S})$ for each $S\subseteq[d]$ based on the chosen
feature removal technique.

In this notation, the entries of $v$ are indexed as $v_{S}$ for
$S\subseteq[d]$, and the entries of $A$ are indexed as $A_{iS}$ for $i\in[d]$
and $S\subseteq[d]$. The operation $Av$ can be understood as summing across
all subsets, or $(Av)_{i}=\sum_{S\subseteq[d]}A_{iS}\cdot v_{S}$. Representing
the attributions this way is convenient for our next results.

The entries for the linear operator in Table 3 are straightforward for the
first four methods, but LIME depends on several implementation choices: these
include the choice of weighting kernel, regularization, and intercept term.
The first three methods are in fact special cases of LIME [17]. The class of
weighting kernel used in practice is difficult to characterize analytically,
but its default parameters for tabular, image and text data approach limiting
cases where LIME reduces to other methods (Appendix H). For simplicity, the
remainder of this section focuses on the other methods.

Perturbing either the input or model induces a change in
$v\in\mathbb{R}^{2^{d}}$, and we must consider how the attributions differ
between the original vector $v$ and a perturbed version $v^{\prime}$.
Proposition 1 suggests that we can bound the attribtuion difference via the
change in model outputs $\lVert v-v^{\prime}\rVert$ and properties of the
matrix $A$. We can even do this under different distance metrics, as we show
in the next result.

######  Lemma 7.

The difference in attributions given the same summary technique $A$ and
different model outputs $v,v^{\prime}$ can be bounded as

$\lVert Av-Av^{\prime}\rVert_{2}\leq\lVert A\rVert_{2}\cdot\lVert
v-v^{\prime}\rVert_{2}\quad\quad\text{or}\quad\quad\lVert Av-
Av^{\prime}\rVert_{2}\leq\lVert A\rVert_{1,\infty}\cdot\lVert
v-v^{\prime}\rVert_{\infty},$

where $\lVert A\rVert_{2}$ is the spectral norm, and the operator norm $\lVert
A\rVert_{1,\infty}$ is the square root of the sum of squared row 1-norms, with
values for each $A$ given in Table 3.

For both inequalities, smaller norms $\lVert A\rVert$ represent stronger
robustness to perturbations. Neither bound is necessarily tighter, and the
choice of which to use depends on how $v-v^{\prime}$ is bounded. The first
bound using the Euclidean distance $\lVert v-v^{\prime}\rVert_{2}$ is perhaps
more intuitive, but the second bound is more useful here because the results
in Sections 3.1 and 3.2 effectively relate to $\lVert
v-v^{\prime}\rVert_{\infty}$.

This view of the summary technique’s robustness raises a natural question,
which is whether the approaches used in practice have relatively good or bad
robustness (Table 3). The answer is not immediately clear, because within the
class of linear operators we can control the robustness arbitrarily: we can
achieve $\lVert A\rVert_{2}=\lVert A\rVert_{1,\infty}=0$ by setting $A=0$
(strong robustness), and we can likewise get $\lVert A\rVert_{2},\lVert
A\rVert_{1,\infty}\to\infty$ by setting entries of $A$ to a large value (weak
robustness). These results are not useful, however, because neither limiting
case results in meaningful attributions.

Rather than considering robustness within the class of all linear operators
$A\in\mathbb{R}^{d\times 2^{d}}$, we therefore restrict our attention to
attributions that are in some sense meaningful. There are multiple ways to
define this, but we consider solutions that satisfy one or more of the
following properties, which are motivated by axioms from the literature on
game-theoretic credit allocation [54, 46].

######  Definition 2.

We define the following properties for linear operators
$A\in\mathbb{R}^{d\times 2^{d}}$ that are used for removal-based explanations:

  * •

(Boundedness) For all $v\in\mathbb{R}^{2^{d}}$, the attributions are bounded
by each feature’s smallest and largest contributions, or $\min_{S\not\ni
i}(v_{S\cup\\{i\\}}-v_{S})\leq(Av)_{i}\leq\max_{S\not\ni
i}(v_{S\cup\\{i\\}}-v_{S})$ for all $i\in[d]$.

  * •

(Symmetry) For all $v\in\mathbb{R}^{2^{d}}$, two features that make equal
marginal contributions to all feature sets must have equal attributions, or
$(Av)_{i}=(Av)_{j}$ if $v_{S\cup\\{i\\}}=v_{S\cup\\{j\\}}$ for all
$S\subseteq[d]$.

  * •

(Efficiency) For all $v\in\mathbb{R}^{2^{d}}$, the attributions sum to the
difference in predictions for the complete and empty sets, or
${\bm{1}}^{\top}Av=v_{[d]}-v_{\\{\\}}$.

Among these properties, we prioritize boundedness because this leads to a
class of solutions that are well-known in game theory, and which are called
probabilistic values [46]. We now show that constraining the summary technique
to satisfy different combinations of these properties yields clear bounds for
the robustness. We first address the operator norm $\lVert
A\rVert_{1,\infty}$, which is the simpler case.

######  Lemma 8.

When the linear operator $A$ satisfies the boundedness property, we have
$\lVert A\rVert_{1,\infty}=2\sqrt{d}$.

The above result applies to several summary techniques shown in Table 3,
including leave-one-out [67], Shapley values [45] and Banzhaf values [13]. We
now address the case of the spectral norm $\lVert A\rVert_{2}$, beginning with
the case where boundedness and symmetry are satisfied.

######  Lemma 9.

When the linear operator $A$ satisfies the boundedness and symmetry
properties, the spectral norm is bounded as follows,

$\frac{1}{2^{d/2-1}}\leq\lVert A\rVert_{2}\leq\sqrt{d+1},$

with $1/2^{d/2-1}$ achieved by the Banzhaf value and $\sqrt{d+1}$ achieved by
leave-one-out.

This shows that the Banzhaf value is the most robust summary within a class of
linear operators, which are known as semivalues [46]. Its robustness is even
better than the Shapley value, a result that was recently discussed by Wang
and Jia [64]. However, the Banzhaf value is criticized for failing to satisfy
the efficiency property [25], so we now address the case where this property
holds.

######  Lemma 10.

When the linear operator $A$ satisfies the boundedness and efficiency
properties, the spectral norm is bounded as follows,

$\sqrt{\frac{2}{d}}\leq\lVert
A\rVert_{2}\leq\sqrt{2+2\cdot\cos\left(\frac{\pi}{d+1}\right)},$

with $\sqrt{2/d}$ achieved by the Shapley value.

We therefore observe that the Shapley value is the most robust choice within a
different class of linear operators, which are known as random-order values
[46]. In comparison, the worst-case robustness is achieved by a method that
has not been discussed in the literature, but which is a special case of a
previously proposed method [26]. The upper bound in Lemma 10 is approximately
equal to $2$, whereas the upper bound in Lemma 9 grows with the input
dimensionality, suggesting that the efficiency property imposes a minimum
level of robustness but limits robustness in the best case.

Finally, we can trivially say that $\lVert A\rVert_{2}=\sqrt{2/d}$ when
boundedness, symmetry and efficiency all hold, because the Shapley value is
the only linear operator to satisfy all three properties [54, 46].

######  Lemma 11.

When the linear operator $A$ satisfies the boundedness, symmetry and
efficiency properties, the spectral norm is $\lVert A\rVert_{2}=\sqrt{2/d}$.

##  4 Feature attribution robustness: main results

We now shift our attention to the robustness properties of entire feature
attribution algorithms. These results build on those presented in Section 3,
combining them to provide general results that apply to existing methods, and
to new methods that combine their implementation choices arbitrarily.

Our first main result relates to the robustness to input perturbations.

######  Theorem 1.

The robustness of removal-based explanations to input perturbations is given
by the following meta-formula,

$\lVert\phi(f,x)-\phi(f,x^{\prime})\rVert_{2}\leq g(\textnormal{removal})\cdot
h(\textnormal{summary})\cdot\lVert x-x^{\prime}\rVert_{2},$

where the factors for each method are defined as follows:

$\displaystyle g(\textnormal{removal})$ $\displaystyle=\begin{cases}L&\text{if
{removal} $=$ {baseline} or {marginal}}\\\ L+2BM&\text{if {removal} $=$
{conditional}},\end{cases}$ $\displaystyle h(\textnormal{summary})$
$\displaystyle=\begin{cases}2\sqrt{d}&\text{if {summary} $=$ {Shapley},
{Banzhaf}, or {leave-one-out}}\\\ \sqrt{d}&\text{if {summary} $=$ {mean when
included}}.\end{cases}$

We therefore observe that the explanation functions themselves are Lipschitz
continuous, but with a constant that combines properties of the original model
with implementation choices of the attribution method. The model’s inherent
robustness interacts with the removal choice to yield the $g(\text{removal})$
factor, and the summary technique is accounted for by $h(\text{summary})$.

Our second result takes a similar form, but relates to the robustness to model
perturbations.

######  Theorem 2.

The robustness of removal-based explanations to model perturbations is given
by the following meta-formula,

$\lVert\phi(f,x)-\phi(f^{\prime},x)\rVert_{2}\leq
h(\textnormal{summary})\cdot\lVert f-f^{\prime}\rVert,$

where the functional distance and factor associated with the summary technique
are specified as follows:

$\displaystyle\lVert f-f^{\prime}\rVert$ $\displaystyle=\begin{cases}\lVert
f-f^{\prime}\rVert_{\infty}&\text{if {removal} $=$ {baseline} or
{marginal}}\\\ \lVert f-f^{\prime}\rVert_{\infty,{\mathcal{X}}}&\text{if
{removal} $=$ {conditional}},\end{cases}$ $\displaystyle
h(\textnormal{summary})$ $\displaystyle=\begin{cases}2\sqrt{d}&\text{if
{summary} $=$ {Shapley}, {Banzhaf}, or {leave-one-out}}\\\ \sqrt{d}&\text{if
{summary} $=$ {mean when included}}.\end{cases}$

Similarly, we see here that the attribution difference is determined by the
strength of the model perturbation, as measured by a functional distance
metric that depends on the feature removal technique. This represents a form
of Lipschitz continuity with respect to the model $f$. These results address
the case with either an input perturbation or model perturbation, but we can
also account for simultaneous input and model perturbations, as shown in
Corollary 1 in Appendix E. The remainder of the paper empirically verifies
these results and discusses their practical implications.

##  5 Summary of experiments

Figure 1: Shapley attribution difference for networks trained with increasing
weight decay, under input perturbations with varying perturbation norms. The
results include (a) the wine quality dataset with FCNs and baseline, marginal,
and conditional feature removal; and (b) MNIST with CNNs and baseline feature
removal with either training set means or zeros. Error bars show the mean and
$95\%$ confidence intervals across explicand-perturbation pairs. Figure 2:
Sanity checks for attributions using cascading randomization for the FCN
trained on the wine quality dataset. Attribution similarity is measured by
Pearson correlation and Spearman rank correlation. We show the mean and $95\%$
confidence intervals across $10$ random seeds.

Due to space constraints, we defer our experiments to Appendix A but briefly
describe the results here. First, we empirically validate our results using
synthetic data and a logistic regression classifier, where the model’s
Lipschitz constant and the upper bound for the total variation constant $M$
can be readily computed. This allows us to verify the robustness bound under
input perturbations (Theorem 1). Using the same setup, we construct logistic
regression classifiers where the functional distance on the data manifold is
small, and we show that feature removal with the conditional distribution is
indeed more robust to model perturbations than baseline or marginal removal,
as implied by Theorem 2.

Next, we demonstrate two practical implications of our findings with the UCI
wine quality [14, 24], MNIST [44], CIFAR-10 [42] and Imagenette datasets [19,
36]. (1) Removal-based attributions are more robust under input perturbations
for networks trained with weight decay, which encourages a smaller Lipschitz
constant. Formally, this is because a network’s Lipschitz constant is upper
bounded by the product of its layers’ spectral norms [61], which are
respectively upper bounded by each layer’s Frobenius norm. These results are
shown for Shapley value attributions in Figure 1 for the wine quality and
MNIST datasets, where we see a consistent effect of weight decay across
feature removal approaches. We speculate that a similar effect may be observed
with other techniques that encourage Lipschitz regularity [63, 30, 8]. (2)
Viewing the sanity check with parameter randomization proposed by Adebayo et
al. [2] as a form of model perturbation, we show that removal-based
attributions empirically pass the sanity check, in contrast to certain
gradient-based methods [55, 60]. Figure 2 shows results for the wine quality
dataset, where the similarity for removal-based methods decays more rapidly
than for IntGrad and Grad $\times$ Input. Results for the remaining datasets
are in Appendix A.

##  6 Discussion

Our analysis focused on the robustness properties of removal-based feature
attributions to perturbations in both the model $f$ and input ${\bm{x}}$. For
input perturbations, Theorem 1 shows that the attribution change is
proportional to the input perturbation strength. The Lipchitz constant is
determined by the attribution’s removal and summary technique, as well as the
model’s inherent robustness. We find that both implementation choices play an
important role, but our analysis in Section 3 reveals that the associated
factors in our main result can only be so small for meaningful attributions,
implying that attribution robustness in some sense reduces to model
robustness. Fortunately, there are many existing works on improving a neural
network’s Lipschitz regularity [63, 30, 8].

For model perturbations, Theorem 2 shows that the attribution difference is
proportional to functional distance between the perturbed and original models,
as measured by the Chebyshev distance (or infinity norm). Depending on the
removal approach, the distance may depend on only a subdomain of the input
space. One practical implication is that attributions should remain similar
under minimal model changes, e.g., accurate model distillation or
quantization. Another implication is that bringing the model closer to its
correct form (e.g., the Bayes classifier) results in attributions closer to
those of the correct model; several recent works have demonstrated this by
exploring ensembled attributions [23, 37, 48, 31], which are in some cases
mathematically equivalent to attributions of an ensembled model.

Our robustness results provide a new lens to compare different removal-based
attribution methods. In terms of the summary technique, our main proof
technique relies on the associated operator norm, which does not distinguish
between several approaches (Shapley, leave-one-out, Banzhaf); however, our
results for the spectral norm suggest that the Banzhaf value is in some sense
most robust. This result echoes recent work [64] and has implications for
other uses of game-theoretic credit allocation, but we also find that the
Shapley value is most robust within a specific class of summary techniques.
Our findings regarding the removal approach are more complex: the conditional
distribution approach has been praised by some works for being more
informative than baseline or marginal [27, 17], but this comes at a cost of
worse robustness to input perturbations. On the other hand, it leads to
improved robustness to model perturbations, which has implications for fooling
attributions with imperceptible off-manifold model changes [58].

Next, we make connections between current understanding of gradient-based
methods and our findings for removal-based methods. Overall, the notion of
Lipschitz continuity is important for robustness guarantees under input
perturbations. Our results for removal-based attributions rely on Lipschitz
continuity of the model itself, whereas gradient-based attributions rely on
Lipschitz continuity of the model’s gradient (i.e., Lipschitz smoothness)
[22]. Practically, both are computationally intensive to determine for
realistic networks [63, 30]. Empirically, weight decay is an effective
approach for generating robust attributions, for both removal-based methods
(our results) and gradient-based methods (see Dombrowski et al. [22]), as the
Frobenius norms of a neural network’s weights upper bound both the network’s
Lipschitz constant and Hessian norm. Under model perturbations, removal-based
attributions are more robust when features are removed such that
${\bm{x}}_{S}$ stay on the data manifold, whereas gradient-based attributions
are more robust when gradients are projected on the tangent space of the data
manifold [21].

Finally, we consider limitations of our analysis. Future work may consider
bounding attribution differences via different distance measures (e.g.,
$\ell_{p}$ norms for $p\neq 2$). The main limitation to our results is that
they are conservative. By relying on global properties of the model, i.e., the
Lipschitz constant and global prediction bounds, we arrive at worst-case
bounds that are in many cases overly conservative, as reflected in our
experiments. An alternative approach can focus on localized versions of our
bounds. Analogous to recent work on certified robustness [39, 62, 52, 57, 38],
tighter bounds may be achievable with specific perturbation strengths, unlike
our approach that provides bounds simultaneously for any perturbation strength
via global continuity properties.

## Acknowledgements

We thank members of the Lee Lab and Jonathan Hayase for helpful discussions.
This work was funded by the National Science Foundation [DBI-1552309 and
DBI-1759487]; and the National Institutes of Health [R35 GM 128638 and R01 NIA
AG 061132].

## References

  * Aas et al. [2021] Kjersti Aas, Martin Jullum, and Anders Løland.  Explaining individual predictions when features are dependent: More accurate approximations to Shapley values.  _Artificial Intelligence_ , 298:103502, 2021.
  * Adebayo et al. [2018] Julius Adebayo, Justin Gilmer, Michael Muelly, Ian Goodfellow, Moritz Hardt, and Been Kim.  Sanity checks for saliency maps.  _Advances in Neural Information Processing Systems_ , 31, 2018.
  * Agarwal et al. [2022] Chirag Agarwal, Marinka Zitnik, and Himabindu Lakkaraju.  Probing gnn explainers: A rigorous theoretical and empirical analysis of gnn explanation methods.  In _International Conference on Artificial Intelligence and Statistics_ , pages 8969–8996. PMLR, 2022.
  * Agarwal et al. [2021] Sushant Agarwal, Shahin Jabbari, Chirag Agarwal, Sohini Upadhyay, Steven Wu, and Himabindu Lakkaraju.  Towards the unification and robustness of perturbation and gradient based explanations.  In _International Conference on Machine Learning_ , pages 110–119. PMLR, 2021.
  * Alvarez-Melis and Jaakkola [2018] David Alvarez-Melis and Tommi S Jaakkola.  On the robustness of interpretability methods.  _arXiv preprint arXiv:1806.08049_ , 2018.
  * Ancona et al. [2017] Marco Ancona, Enea Ceolini, Cengiz Öztireli, and Markus Gross.  Towards better understanding of gradient-based attribution methods for deep neural networks.  _arXiv preprint arXiv:1711.06104_ , 2017.
  * Anders et al. [2020] Christopher Anders, Plamen Pasliev, Ann-Kathrin Dombrowski, Klaus-Robert Müller, and Pan Kessel.  Fairwashing explanations with off-manifold detergent.  In _International Conference on Machine Learning_ , pages 314–323. PMLR, 2020.
  * Araujo et al. [2023] Alexandre Araujo, Aaron Havens, Blaise Delattre, Alexandre Allauzen, and Bin Hu.  A unified algebraic perspective on Lipschitz neural networks.  _arXiv preprint arXiv:2303.03169_ , 2023.
  * Banzhaf III [1964] John F Banzhaf III.  Weighted voting doesn’t work: A mathematical analysis.  _Rutgers Law Review_ , 19:317, 1964.
  * Bhatt et al. [2020] Umang Bhatt, Alice Xiang, Shubham Sharma, Adrian Weller, Ankur Taly, Yunhan Jia, Joydeep Ghosh, Ruchir Puri, José MF Moura, and Peter Eckersley.  Explainable machine learning in deployment.  In _Proceedings of the 2020 Conference on Fairness, Accountability, and Transparency_ , pages 648–657, 2020.
  * Boyd et al. [2004] Stephen Boyd, Stephen P Boyd, and Lieven Vandenberghe.  _Convex Optimization_.  Cambridge University Press, 2004.
  * Chen et al. [2022] Hugh Chen, Ian C Covert, Scott M Lundberg, and Su-In Lee.  Algorithms to estimate Shapley value feature attributions.  _arXiv preprint arXiv:2207.07605_ , 2022.
  * Chen and Jordan [2020] Jianbo Chen and Michael Jordan.  LS-Tree: Model interpretation when the data are linguistic.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 34, pages 3454–3461, 2020.
  * Cortez et al. [2009] Paulo Cortez, António Cerdeira, Fernando Almeida, Telmo Matos, and José Reis.  Modeling wine preferences by data mining from physicochemical properties.  _Decision Support Systems_ , 47(4):547–553, 2009.
  * Covert et al. [2020] Ian Covert, Scott M Lundberg, and Su-In Lee.  Understanding global feature contributions with additive importance measures.  _Advances in Neural Information Processing Systems_ , 33:17212–17223, 2020.
  * Covert et al. [2022] Ian Covert, Chanwoo Kim, and Su-In Lee.  Learning to estimate Shapley values with vision transformers.  _arXiv preprint arXiv:2206.05282_ , 2022.
  * Covert et al. [2021] Ian C Covert, Scott Lundberg, and Su-In Lee.  Explaining by removing: A unified framework for model explanation.  _The Journal of Machine Learning Research_ , 22(1):9477–9566, 2021.
  * DeGrave et al. [2021] Alex J DeGrave, Joseph D Janizek, and Su-In Lee.  Ai for radiographic covid-19 detection selects shortcuts over signal.  _Nature Machine Intelligence_ , 3(7):610–619, 2021.
  * Deng et al. [2009] Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei.  Imagenet: A large-scale hierarchical image database.  In _2009 IEEE conference on computer vision and pattern recognition_ , pages 248–255. Ieee, 2009.
  * Dimanov et al. [2020] Botty Dimanov, Umang Bhatt, Mateja Jamnik, and Adrian Weller.  You shouldn’t trust me: Learning models which conceal unfairness from multiple explanation methods.  2020\.
  * Dombrowski et al. [2019] Ann-Kathrin Dombrowski, Maximillian Alber, Christopher Anders, Marcel Ackermann, Klaus-Robert Müller, and Pan Kessel.  Explanations can be manipulated and geometry is to blame.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Dombrowski et al. [2022] Ann-Kathrin Dombrowski, Christopher J Anders, Klaus-Robert Müller, and Pan Kessel.  Towards robust explanations for deep neural networks.  _Pattern Recognition_ , 121:108194, 2022.
  * Dong and Rudin [2020] Jiayun Dong and Cynthia Rudin.  Exploring the cloud of variable importance for the set of all good models.  _Nature Machine Intelligence_ , 2(12):810–824, 2020.
  * Dua and Graff [2017] Dheeru Dua and Casey Graff.  UCI machine learning repository, 2017.  URL http://archive.ics.uci.edu/ml.
  * Dubey and Shapley [1979] Pradeep Dubey and Lloyd S Shapley.  Mathematical properties of the Banzhaf power index.  _Mathematics of Operations Research_ , 4(2):99–131, 1979.
  * Frye et al. [2020] Christopher Frye, Colin Rowat, and Ilya Feige.  Asymmetric Shapley values: incorporating causal knowledge into model-agnostic explainability.  _Advances in Neural Information Processing Systems_ , 33:1229–1239, 2020.
  * Frye et al. [2021] Christopher Frye, Damien de Mijolla, Tom Begley, Laurence Cowton, Megan Stanley, and Ilya Feige.  Shapley explainability on the data manifold.  In _International Conference on Learning Representations_ , 2021.
  * Ghorbani et al. [2019] Amirata Ghorbani, Abubakar Abid, and James Zou.  Interpretation of neural networks is fragile.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 33, pages 3681–3688, 2019.
  * Goodfellow et al. [2014] Ian J Goodfellow, Jonathon Shlens, and Christian Szegedy.  Explaining and harnessing adversarial examples.  _arXiv preprint arXiv:1412.6572_ , 2014.
  * Gouk et al. [2021] Henry Gouk, Eibe Frank, Bernhard Pfahringer, and Michael J Cree.  Regularisation of neural networks by enforcing Lipschitz continuity.  _Machine Learning_ , 110:393–416, 2021.
  * Gyawali et al. [2022] Prashnna K Gyawali, Xiaoxia Liu, James Zou, and Zihuai He.  Ensembling improves stability and power of feature selection for deep learning models.  In _Machine Learning in Computational Biology_ , pages 33–45. PMLR, 2022.
  * Hammer and Holzman [1992] Peter L Hammer and Ron Holzman.  Approximations of pseudo-boolean functions; applications to game theory.  _Zeitschrift für Operations Research_ , 36(1):3–21, 1992.
  * He et al. [2016] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 770–778, 2016.
  * Heo et al. [2019] Juyeon Heo, Sunghwan Joo, and Taesup Moon.  Fooling neural network interpretations via adversarial model manipulation.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Hoeffding [1994] Wassily Hoeffding.  Probability inequalities for sums of bounded random variables.  _The Collected Works of Wassily Hoeffding_ , pages 409–426, 1994.
  * Howard and Gugger [2020] Jeremy Howard and Sylvain Gugger.  FastAI: A layered API for deep learning.  _Information_ , 11(2):108, 2020.
  * Janizek et al. [2021] Joseph D Janizek, Ayse B Dincer, Safiye Celik, Hugh Chen, William Chen, Kamila Naxerova, and Su-In Lee.  Uncovering expression signatures of synergistic drug response using an ensemble of explainable AI models.  _bioRxiv_ , pages 2021–10, 2021.
  * Jordan et al. [2022] Matt Jordan, Jonathan Hayase, Alex Dimakis, and Sewoong Oh.  Zonotope domains for lagrangian neural network verification.  _Advances in Neural Information Processing Systems_ , 35:8400–8413, 2022.
  * Katz et al. [2017] Guy Katz, Clark Barrett, David L Dill, Kyle Julian, and Mykel J Kochenderfer.  Reluplex: An efficient SMT solver for verifying deep neural networks.  In _Computer Aided Verification: 29th International Conference, CAV 2017, Heidelberg, Germany, July 24-28, 2017, Proceedings, Part I 30_ , pages 97–117. Springer, 2017.
  * Khan et al. [2022] Zulqarnain Khan, Aria Masoomi, Davin Hill, and Jennifer Dy.  Analyzing the effects of classifier lipschitzness on explainers.  _arXiv preprint arXiv:2206.12481_ , 2022.
  * Kindermans et al. [2019] Pieter-Jan Kindermans, Sara Hooker, Julius Adebayo, Maximilian Alber, Kristof T Schütt, Sven Dähne, Dumitru Erhan, and Been Kim.  The (un) reliability of saliency methods.  _Explainable AI: Interpreting, explaining and visualizing deep learning_ , pages 267–280, 2019.
  * Krizhevsky et al. [2009] Alex Krizhevsky, Geoffrey Hinton, et al.  Learning multiple layers of features from tiny images.  2009\.
  * Kurakin et al. [2018] Alexey Kurakin, Ian J Goodfellow, and Samy Bengio.  Adversarial examples in the physical world.  In _Artificial Intelligence Safety and Security_ , pages 99–112. Chapman and Hall/CRC, 2018.
  * LeCun et al. [1998] Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner.  Gradient-based learning applied to document recognition.  _Proceedings of the IEEE_ , 86(11):2278–2324, 1998.
  * Lundberg and Lee [2017] Scott M Lundberg and Su-In Lee.  A unified approach to interpreting model predictions.  _Advances in Neural Information Processing Systems_ , 30, 2017.
  * Monderer and Samet [2002] Dov Monderer and Dov Samet.  Variations on the Shapley value.  _Handbook of Game Theory with Economic Applications_ , 3:2055–2076, 2002.
  * Montavon et al. [2019] Grégoire Montavon, Alexander Binder, Sebastian Lapuschkin, Wojciech Samek, and Klaus-Robert Müller.  Layer-wise relevance propagation: an overview.  _Explainable AI: Interpreting, Explaining and Visualizing Deep Learning_ , pages 193–209, 2019.
  * Ning et al. [2022] Yilin Ning, Mingxuan Liu, and Nan Liu.  Shapley variable importance cloud for machine learning models.  _arXiv preprint arXiv:2212.08370_ , 2022.
  * Noschese et al. [2013] Silvia Noschese, Lionello Pasquini, and Lothar Reichel.  Tridiagonal toeplitz matrices: properties and novel applications.  _Numerical Linear Algebra With Applications_ , 20(2):302–326, 2013.
  * Novakovsky et al. [2022] Gherman Novakovsky, Nick Dexter, Maxwell W Libbrecht, Wyeth W Wasserman, and Sara Mostafavi.  Obtaining genetics insights from deep learning via explainable artificial intelligence.  _Nature Reviews Genetics_ , pages 1–13, 2022.
  * Petsiuk et al. [2018] Vitali Petsiuk, Abir Das, and Kate Saenko.  RISE: Randomized input sampling for explanation of black-box models.  _arXiv preprint arXiv:1806.07421_ , 2018.
  * Raghunathan et al. [2018] Aditi Raghunathan, Jacob Steinhardt, and Percy S Liang.  Semidefinite relaxations for certifying robustness to adversarial examples.  _Advances in Neural Information Processing Systems_ , 31, 2018.
  * Ribeiro et al. [2016] Marco Tulio Ribeiro, Sameer Singh, and Carlos 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_ , pages 1135–1144, 2016.
  * Shapley [1953] Lloyd S Shapley.  A value for n-person games.  _Contributions to the Theory of Games_ , 2(28):307–317, 1953.
  * Shrikumar et al. [2017] Avanti Shrikumar, Peyton Greenside, and Anshul Kundaje.  Learning important features through propagating activation differences.  In _International Conference on Machine Learning_ , pages 3145–3153. PMLR, 2017.
  * Simonyan et al. [2013] Karen Simonyan, Andrea Vedaldi, and Andrew Zisserman.  Deep inside convolutional networks: Visualising image classification models and saliency maps.  _arXiv preprint arXiv:1312.6034_ , 2013.
  * Singh et al. [2019] Gagandeep Singh, Timon Gehr, Markus Püschel, and Martin Vechev.  An abstract domain for certifying neural networks.  _Proceedings of the ACM on Programming Languages_ , 3(POPL):1–30, 2019.
  * Slack et al. [2020] Dylan Slack, Sophie Hilgard, Emily Jia, Sameer Singh, and Himabindu 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.
  * Sundararajan and Najmi [2020] Mukund Sundararajan and Amir Najmi.  The many Shapley values for model explanation.  In _International Conference on Machine Learning_ , pages 9269–9278. PMLR, 2020.
  * Sundararajan et al. [2017] Mukund Sundararajan, Ankur Taly, and Qiqi Yan.  Axiomatic attribution for deep networks.  In _International Conference on Machine Learning_ , pages 3319–3328. PMLR, 2017.
  * Szegedy et al. [2013] Christian Szegedy, Wojciech Zaremba, Ilya Sutskever, Joan Bruna, Dumitru Erhan, Ian Goodfellow, and Rob Fergus.  Intriguing properties of neural networks.  _arXiv preprint arXiv:1312.6199_ , 2013.
  * Tjeng et al. [2017] Vincent Tjeng, Kai Xiao, and Russ Tedrake.  Evaluating robustness of neural networks with mixed integer programming.  _arXiv preprint arXiv:1711.07356_ , 2017.
  * Virmaux and Scaman [2018] Aladin Virmaux and Kevin Scaman.  Lipschitz regularity of deep neural networks: analysis and efficient estimation.  _Advances in Neural Information Processing Systems_ , 31, 2018.
  * Wang and Jia [2022] Tianhao Wang and Ruoxi Jia.  Data Banzhaf: A data valuation framework with maximal robustness to learning stochasticity.  _arXiv preprint arXiv:2205.15466_ , 2022.
  * Wang et al. [2020] Zifan Wang, Haofan Wang, Shakul Ramkumar, Piotr Mardziel, Matt Fredrikson, and Anupam Datta.  Smoothed geometry for robust attribution.  _Advances in Neural Information Processing Systems_ , 33:13623–13634, 2020.
  * Yoshida and Miyato [2017] Yuichi Yoshida and Takeru Miyato.  Spectral norm regularization for improving the generalizability of deep learning.  _arXiv preprint arXiv:1705.10941_ , 2017.
  * Zeiler and Fergus [2014] Matthew D Zeiler and Rob Fergus.  Visualizing and understanding convolutional networks.  In _Computer Vision–ECCV 2014: 13th European Conference, Zurich, Switzerland, September 6-12, 2014, Proceedings, Part I 13_ , pages 818–833. Springer, 2014.

##  Appendix A Experiments

This section conducts a range of experiments to provide empirical validation
for our theoretical results, and to examine their practical implications for
common machine learning models. Code is available at
https://github.com/suinleelab/removal-robustness.

###  A.1 Synthetic data with input perturbation

First, we empirically validate our results from Section 4 using synthetic data
and a logistic regression classifier, where the model’s Lipschitz constant $L$
and the data distribution’s total variation constant $M$ can be readily
computed for the bound in Theorem 1.

Setup. For each input $x\in\mathbb{R}^{4}$ in the synthetic dataset, features
are generated from a multivariate Gaussian distribution with zero mean,
${\bm{x}}_{3},{\bm{x}}_{4}$ independent from each other and from
${\bm{x}}_{1},{\bm{x}}_{2}$, and with a correlation $\rho$ between
${\bm{x}}_{1},{\bm{x}}_{2}$. That is, the covariance matrix
$\Sigma\in\mathbb{R}^{4\times 4}$ has the following form:

$\Sigma=\begin{pmatrix}1&\rho&0&0\\\ \rho&1&0&0\\\ 0&0&1&0\\\ 0&0&0&1\\\
\end{pmatrix}.$

The binary label ${\bm{y}}\in\\{0,1\\}$ for each input is sampled from a
Bernoulli distribution with probability $\sigma(\beta^{\top}{\bm{x}})$, where
$\sigma(\cdot)$ is the sigmoid function and $\beta^{\top}=[5,0,3,1]$. Unless
stated otherwise, we consider $\rho=1$, zeros are used for baseline feature
removal, and the exact Gaussian distributions are used for marginal and
conditional feature removal with $20{,}000$ samples.

A logistic regression classifier was trained with $500$ samples, and removal-
based feature attributions were computed for $50$ test samples (i.e.,
explicands) with and without input perturbation. In our experiment,
perturbations are randomly sampled from standard Gaussian distributions with
$\ell_{2}$ norms normalized to values ranging between $0$ to $2$ with a total
of $50$ different values. For each perturbation norm, $50$ perturbations are
generated and applied to each explicand, resulting in $2{,}500$ total
explicand-perturbation pairs. As illustrated in Figure 3, a line is drawn to
map the maximum $\ell_{2}$ difference between original and perturbed
attributions at each perturbation norm. The slope of such a line estimates the
dependency of the upper bound of attribution differences with respect to the
perturbation norm. Theorem 1 guarantees a specific slope, and this experiment
aims to verify this property.

Furthermore, note that the sigmoid function $\sigma(\cdot)$ is $1/4$-Lipschitz
because its largest derivative is $1/4$. A logistic regression classifier is a
composition of a linear function and sigmoid function $\sigma(\cdot)$, so the
model is $L$-Lipschitz with $L=\lVert\hat{\beta}\rVert_{2}/4$, where
$\hat{\beta}$ is the trained model coefficients. Finally, because the input
features follow a Gaussian distribution, the total variation constant $M$
derived in Example 1 applies here. Hence, we have specific values for all
terms in Theorem 1 bound, which allows a direct comparison between theoretical
and empirical results.

Figure 3: Shapley attribution differences under input perturbation with
various perturbation norms in our synthetic data experiment. Lines bounding
the maximum attribution difference at each perturbation norm are shown as
empirical upper bounds. Theoretical upper bounds are included for baseline and
marginal feature removal and omitted for conditional feature removal because
it is infinite for $\rho=1$.

Results. As shown in Figure 3, for all three feature removal approaches, the
theoretical upper bounds for the Shapley value are indeed greater than the
corresponding empirical upper bounds. We also observe that the empirical upper
bound with conditional feature removal exceeds those with baseline and
marginal feature removal, as expected due to our results involving the
conditional distribution’s total variation distance (Lemma 2).

Moreover, we vary the Lipschitz constant $L$ of the trained logistic
regression classifier by a factor of $\alpha=\\{0.05,0.1,0.5,1,2\\}$, which is
done by multiplying the trained coefficients $\hat{\beta}$ by $\alpha$. Figure
4 shows that the perturbation norm-dependent slope of the empirical bound
increases linearly with $\alpha$, aligning with our theoretical results. Also,
Example 1 suggests that increasing $\rho$ would lead to a larger upper bound
for conditional feature removal, which is observed in the empirical bound with
$\rho=\\{0,0.25,0.5,0.75,1\\}$ (see Figure 4 right).

Overall, these empirical results validate our theoretical bound from Theorem
1. We note that gaps between our theoretical and empirical bounds exist,
likely because we use the global Lipschitz constant to provide guarantee under
all input perturbation strengths, instead of local constants for specific
input regions and perturbation strengths, and also due to looseness in the
total variation constant $M$. For example, our bound for multivariate Gaussian
distributions in Example 1 suggests that $M=\infty$ when $\rho=1$, which is
far more conservative than our empirical results (see Figure 3 right); this is
due to looseness in Pinsker’s inequality (Appendix B). Similar results are
observed for Occlusion, Banzhaf, RISE, and LIME attributions, and these are
included in Appendix I.

Figure 4: Slopes of empirical and theoretical upper bounds with respect to the
input perturbation norm, for Shapley value attribution difference in the
synthetic experiment. The parameter $\alpha$ varies the logistic regression
Lipschitz constant, and $\rho$ varies the correlation between the two
synthetic features ${\bm{x}}_{1},{\bm{x}}_{2}$.

###  A.2 Synthetic data with model perturbation

Here, we construct a logistic regression classifier as well as a perturbed
version such that their functional distance on the data manifold is small. We
use this construction to empirically verify that removing features with the
conditional distribution leads to attributions that are more robust against
model perturbation, as suggested by Theorem 2.

Setup. The synthetic data generation follows that of Section A.1. A logistic
regression classifier $f$ is constructed with coefficients $[5,0,3,1]$, and
its perturbed counterpart $f^{\prime}$ has coefficients $[0,5,3,1]$. Because
the synthetic features ${\bm{x}}_{1},{\bm{x}}_{2}$ are identical with the
default value $\rho=1$, we have the functional distance $\lVert
f-f^{\prime}\rVert_{\infty,{\mathcal{X}}}=0$, so our theoretical bound
suggests that $\lVert\phi(f,x)-\phi(f^{\prime},x)\rVert_{2}=0$ for conditional
feature removal. When we consider smaller values for $\rho$, we still expect
better robustness for conditional rather than marginal removal due to Lemma 4,
because significant probability mass lies in a subdomain where $f$ and
$f^{\prime}$ are functionally close (near the hyperplane where
${\bm{x}}_{1}={\bm{x}}_{2}$). Feature attributions are computed for $100$
explicands with the same feature removal settings as in Section A.1.

Figure 5: Shapley value attribution differences between the logistic
regression classifiers $f,f^{\prime}$ as constructed in Section A.2 with
baseline, marginal, and conditional feature removal, and varying correlation
$\rho$.

Results. As suggested by our theoretical bound for model perturbation when
$\rho=1$, the Shapley attribution difference between $f,f^{\prime}$ is near
zero when features are removed using their conditional distribution, whereas
the attribution difference can be large with baseline and marginal feature
removal (Figure 5). When we decrease $\rho$ in the synthetic data generation,
the models $f,f^{\prime}$ become more likely to generate different outputs,
and we observe that the robustness gap between conditional vs. baseline and
marginal feature removal reduces. These empirical results corroborate Lemma 4.
Similar results are observed for Occlusion, Banzhaf, RISE, and LIME
attributions and are included in Appendix I. It is worth noting that the
arbitrarily large bound for baseline and marginal feature removal underpins
the adversarial attacks on LIME and SHAP demonstrated by Slack et al. [58].
Our theory and synthetic experiments for model perturbation are similar to
those in Frye et al. [27], but we provide a unified analysis for removal-based
attribution methods instead of focusing on the Shapley value, and we consider
cases where models disagree even on the data manifold (Lemma 4).

###  A.3 Implication I: robust attributions for regularized neural networks

Our theoretical upper bound for input perturbation decreases with the model’s
Lipschitz constant (Theorem 1), suggesting that training a model with better
Lipschitz regularity can help obtain more robust removal-based feature
attributions. Weight decay is one way to achieve this, as it encourages lower
Frobenius norm for the individual layers in a neural network, and the product
of these Frobenius norms upper bounds the network’s Lipschitz constant [66].
We therefore examine whether neural networks trained with weight decay produce
attributions that are more robust against input perturbation.

Datasets. The UCI white wine quality dataset [14, 24] consists of $4{,}898$
samples of white wine, each with 11 input features for predicting the wine
quality (whether the rating $>6$). Two subsets of $500$ samples are randomly
chosen as the validation and test set. For datasets with larger input
dimension, we use MNIST [44], CIFAR-10 [42], and ten classes of ImageNet [19]
(i.e., Imagenette222https://github.com/fastai/imagenette). From the official
training set of MNIST and CIFAR-10, $10{,}000$ images are randomly chosen as
the validation set.

Setup. Fully connected networks (FCNs) and convolutional networks (CNNs) are
trained with a range of weight decay values for the wine quality dataset and
MNIST, respectively, where we use the values
$\\{0,0.001,0.0025,0.005,0.01\\}$. ResNet-18 networks [33] are trained for
CIFAR-10, with weight decay values $\\{0.0075,0.01,0.025,0.05,0.075\\}$.
ResNet-50 networks pre-trained on ImageNet are fine-tuned for Imagenette, with
weight decay values $\\{0.005,0.01,0.025,0.05\\}$. The test set accuracy does
not differ much with these weight decay values for the wine quality dataset,
MNIST, and CIFAR-10 (Figure 20 and Figure 21), but large weight decay values
can worsen the performance of ResNet-50 on Imagenette (Figure 21). Details on
model architectures and other hyperparameters are in Appendix J. Removal-based
attributions are computed for the $500$ test samples and all features in the
wine quality dataset, and for $100$ test images with non-overlapping patches
of size $4\times 4$ as features in MNIST and CIFAR-10, and of size $16\times
16$ as features in Imagenette (see Appendix G for a discussion of how our
theory extends to the use of feature groups). For the wine quality dataset,
removal-based attributions are computed exactly with all $2^{11}$ feature
subsets. For MNIST, certain attributions must be estimated rather than
calculated exactly: Occlusion attributions are computed exactly, whereas
Shapley, Banzhaf, LIME, and RISE are estimated with $10{,}000$ random feature
sets, and Shapley and Banzhaf use a least squares regression similar to
KernelSHAP [45]. For CIFAR-10 and Imagenette, we focus on Occlusion, RISE, and
Shapley values, computed with the same settings as for MNIST, except RISE and
Shapley values are estimated with $20{,}000$ random feature sets for
Imagenette due to its large number of features.

Figure 6: Identical to Figure 1 and shown here again for readability. Shapley
attribution difference for networks trained with increasing weight decay,
under input perturbations with perturbation norms $\\{0.1,0.5,1,2\\}$. The
results include (a) the wine quality dataset with FCNs and baseline (replacing
with training set minimums), marginal, and conditional feature removal; and
(b) MNIST with CNNs and baseline feature removal with either training set
means or zeros. Error bars show the mean and $95\%$ confidence intervals
across explicand-perturbation pairs. Figure 7: Attribution difference for
ResNet-18 networks trained on CIFAR-10 with increasing weight decay, under
input perturbations with perturbation norms $\\{0.1,0.5,1,2\\}$. Error bars
show the means and $95\%$ confidence intervals across explicand-perturbation
pairs.

When using baseline feature removal, features are replaced with the training
set minimums for the wine quality dataset, with the training set means or
zeros for MNIST and CIFAR-10, and with zeros (identical to the training set
means due to normalization) for Imagenette. The wine quality training set is
used for removing features with their marginal distributions, while a
conditional variational autoencoder (CVAE) is trained following the procedure
in Frye et al. [27] to approximate feature removal with the conditional
distribution (details in Appendix J). For both marginal and conditional
feature removal in the wine quality dataset, $250$ samples are used. Marginal
and conditional feature removal are computationally expensive and less common
for image data, so we use only baseline feature removal as is often done in
practice. Finally, input perturbations are generated as in Section A.1, with
$10$ random perturbations applied per explicand and with perturbation norms
$\\{0.1,0.5,1,2\\}$ for the wine quality dataset, MNIST, and CIFAR-10. For
Imagenette, $5$ random perturbations with norms $\\{50,100,200,400\\}$ are
applied per explicand to obtain similar perturbation strength per input.

Results. Figure 6 shows that Shapley value attribution differences decrease
with increasing weight decay, across different feature removal techniques, for
both the wine quality dataset and MNIST, and for all perturbation norms in
$\\{0.1,0.5,1,2\\}$. This demonstrates the practical implication that removal-
based attributions can be made more robust against input perturbation by
training with weight decay. Interestingly, weight decay has also been shown to
improve the robustness of gradient-based attributions [22]. In Appendix K,
similar results are observed for Occlusion, Banzhaf, RISE, and LIME
attributions, and also when we measure explanation similarity using Pearson or
Spearman rank correlation rather than $\ell_{2}$ distance. For ResNet-18
networks trained on CIFAR-10 and ResNet-50 networks trained on Imagenette,
similar trends are observed for Occlusion, RISE, and Shapley attributions
(Figure 7 and Figure 8).

Figure 8: Attribution difference for ResNet-50 networks trained on Imagenette
with increasing weight decay, under input perturbations with perturbation
norms $\\{50,100,200,400\\}$. Error bars show the means and $95\%$ confidence
intervals across explicand-perturbation pairs.

###  A.4 Implication II: sanity checks for removal-based attributions

The robustness bound for model perturbation in Theorem 2 suggests that the
attribution difference between an intact and perturbed model should be small
if the model perturbation is small. Similarly, we expect that if the
perturbation intensifies, the attribution differences may become larger. One
way to perturb a neural network with increasing intensity is to sequentially
randomize its parameters from the output to input layer, a procedure known as
cascading randomization introduced by Adebayo et al. [2] for sanity checking
feature attributions. Here, we study whether removal-based attributions indeed
pass the parameter randomization check, as suggested by our theoretical bound.

Figure 9: Identical to Figure 2 and shown here again for readability. Sanity
checks for attributions using cascading randomization for the FCN trained on
the wine quality dataset. Attribution similarity is measured by Pearson
correlation and Spearman rank correlation. We show the mean and $95\%$
confidence intervals across $10$ random seeds.  Figure 10: Sanity checks for
attributions using cascading randomization for the CNN trained on MNIST.
Attribution similarity is measured by Pearson correlation and Spearman rank
correlation. We show the mean and $95\%$ confidence intervals across $10$
random seeds.

Datasets. We use the same datasets and splits as in Section A.3.

Setup. The best performing networks from Section A.3 are used here. The
settings for removal-based attributions are the same as in Section A.3, except
$50{,}000$ random feature sets are used for more stable Shapley estimations in
Imagenette. Parameters for each network are randomized in a cascading fashion
from the output to input layer as in Adebayo et al. [2], with attributions
computed after randomizing each layer and compared to the original
attributions using Pearson correlation and Spearman rank correlation.
Randomization is performed using the layer’s default initialization. Feature
attributions with the original model are computed twice and compared to
capture variations due to sampling in feature removal and in attribution
estimation using feature sets. This experiment is run $10$ times with
different random seeds for the wine quality dataset, MNIST, and CIFAR-10, and
$3$ times for Imagenette. Gradient-based attribution methods such as Vanilla
Gradients [56], Integrated Gradients [60], and Grad $\times$ Input [55] are
included as references, as they were the focus of the original work on sanity
checks [2].

Results. As suggested by our robustness bound for model perturbation, the
original and perturbed removal-based attributions become dissimilar as more
parameters are randomized, leading to near-zero correlations when the entire
network is re-initialized (Figure 9). This shows that removal-based
attributions reliably pass the model randomization sanity check. In contrast,
we observe that for Integrated Gradients and Grad $\times$ Input, the
correlation between intact and perturbed attributions remains around $0.5$
even when the entire network is randomized. Similarly, for the CNN trained on
MNIST, the original and perturbed removal-based attributions have lower
correlations as more parameters are randomized, whereas the correlations
decrease only slightly for Integrated Gradients and Grad $\times$ Input
(Figure 10). Similar trends are observed for the ResNet-18 trained on CIFAR-10
and ResNet-50 trained on Imagenette (Figure 11, Figure 12, and Figure 29).

Figure 11: Sanity checks for attributions using cascading randomization for
the ResNet-18 trained on CIFAR-10. Features are removed by replacing them with
training set means. Attribution similarity is measured by Pearson correlation
and Spearman rank correlation. We show the means and $95\%$ confidence
intervals across $10$ random seeds. Missing points correspond to undefined
correlations due to constant attributions under model randomization.  Figure
12: Sanity checks for attributions using cascading randomization for the
ResNet-50 trained on Imagenette. Features are removed by replacing them with
zeros (identical to training set means due to normalization). Attribution
similarity is measured by Pearson correlation and Spearman rank correlation.
We show the means and $95\%$ confidence intervals across $3$ random seeds.
Missing points correspond to undefined correlations due to constant
attributions under model randomization.

##  Appendix B Proofs: prediction robustness to input perturbations

We now re-state and prove our results from Section 3.1.

######  Lemma 1.

When removing features using either the baseline or marginal approaches, the
prediction function $f(x_{S})$ for any feature set $x_{S}$ is $L$-Lipschitz
continuous:

$|f(x_{S})-f(x_{S}^{\prime})|\leq L\cdot\lVert
x_{S}-x_{S}^{\prime}\rVert_{2}\quad\forall\;x_{S},x_{S}^{\prime}\in\mathbb{R}^{|S|}.$

###### Proof.

The key insight for this result is that removing features using either the
baseline or marginal approach is equivalent to using a distribution
$q({\bm{x}}_{\bar{S}})$ that does not depend on the observed features. Because
of that, we have:

$\displaystyle|f(x_{S})-f(x_{S}^{\prime})|$ $\displaystyle=\Big{|}\int
f(x_{S},x_{\bar{S}})q(x_{\bar{S}})dx_{\bar{S}}-\int
f(x_{S}^{\prime},x_{\bar{S}})q(x_{\bar{S}})dx_{\bar{S}}\Big{|}$
$\displaystyle\leq\int\big{|}f(x_{S},x_{\bar{S}})-f(x_{S}^{\prime},x_{\bar{S}})\big{|}\cdot
q(x_{\bar{S}})dx_{\bar{S}}$ $\displaystyle\leq L\cdot\lVert
x_{S}-x_{S}^{\prime}\rVert_{2}\cdot\int q(x_{\bar{S}})dx_{\bar{S}}$
$\displaystyle=L\cdot\lVert x_{S}-x_{S}^{\prime}\rVert_{2}.$

∎

######  Lemma 2.

When removing features using the conditional approach, the prediction function
$f(x_{S})$ for a feature set $x_{S}$ satisfies

$|f(x_{S})-f(x_{S}^{\prime})|\leq L\cdot\lVert
x_{S}-x_{S}^{\prime}\rVert_{2}+2B\cdot d_{TV}\Big{(}p({\bm{x}}_{\bar{S}}\mid
x_{S}),p({\bm{x}}_{\bar{S}}\mid x_{S}^{\prime})\Big{)},$

where the total variation distance is defined via the $L^{1}$ functional
distance as

$d_{TV}\Big{(}p({\bm{x}}_{\bar{S}}\mid x_{S}),p({\bm{x}}_{\bar{S}}\mid
x_{S}^{\prime})\Big{)}\equiv\frac{1}{2}\Big{\lVert}p({\bm{x}}_{\bar{S}}\mid
x_{S})-p({\bm{x}}_{\bar{S}}\mid x_{S}^{\prime})\Big{\rVert}_{1}.$

###### Proof.

In this case, we must be careful to account for the different distributions
used when marginalizing out the removed features.

$\displaystyle|f(x_{S})-f(x_{S}^{\prime})|$ $\displaystyle=\Big{|}\int
f(x_{S},x_{\bar{S}})p(x_{\bar{S}}\mid x_{S})dx_{\bar{S}}-\int
f(x_{S}^{\prime},x_{\bar{S}})p(x_{\bar{S}}\mid
x_{S}^{\prime})dx_{\bar{S}}\Big{|}$
$\displaystyle\leq\int\big{|}f(x_{S},x_{\bar{S}})p(x_{\bar{S}}\mid
x_{S})-f(x_{S}^{\prime},x_{\bar{S}})p(x_{\bar{S}}\mid
x_{S}^{\prime})\big{|}dx_{\bar{S}}$
$\displaystyle\leq\int\big{|}f(x_{S},x_{\bar{S}})p(x_{\bar{S}}\mid
x_{S})-f(x_{S}^{\prime},x_{\bar{S}})p(x_{\bar{S}}\mid
x_{S})\big{|}dx_{\bar{S}}$
$\displaystyle\quad+\int\big{|}f(x_{S}^{\prime},x_{\bar{S}})p(x_{\bar{S}}\mid
x_{S})-f(x_{S}^{\prime},x_{\bar{S}})p(x_{\bar{S}}\mid
x_{S}^{\prime})\big{|}dx_{\bar{S}}$
$\displaystyle\leq\int\big{|}f(x_{S},x_{\bar{S}})-f(x_{S}^{\prime},x_{\bar{S}})\big{|}\cdot
p(x_{\bar{S}}\mid x_{S})dx_{\bar{S}}$
$\displaystyle\quad+\int\big{|}f(x_{S}^{\prime},x_{\bar{S}})\big{|}\cdot\big{|}p(x_{\bar{S}}\mid
x_{S})-p(x_{\bar{S}}\mid x_{S}^{\prime})\big{|}dx_{\bar{S}}$
$\displaystyle\leq L\cdot\lVert x_{S}-x_{S}^{\prime}\rVert_{2}\cdot\int
p(x_{\bar{S}}\mid x_{S})dx_{\bar{S}}+B\cdot\int\big{|}p(x_{\bar{S}}\mid
x_{S})-p(x_{\bar{S}}\mid x_{S}^{\prime})\big{|}dx_{\bar{S}}$
$\displaystyle=L\cdot\lVert x_{S}-x_{S}^{\prime}\rVert_{2}+2B\cdot
d_{TV}\Big{(}p(x_{\bar{S}}\mid x_{S}),p(x_{\bar{S}}\mid
x_{S}^{\prime})\Big{)}.$

∎

######  Example 1.

For a Gaussian random variable ${\bm{x}}\sim{\mathcal{N}}(\mu,\Sigma)$ with
mean $\mu\in\mathbb{R}^{d}$ and covariance $\Sigma\in\mathbb{R}^{d\times d}$,
Assumption 3 holds with
$M=\frac{1}{2}\sqrt{\lambda_{\max}(\Sigma^{-1})-\lambda_{\min}(\Sigma^{-1})}$.
If ${\bm{x}}$ is assumed to be standardized, this captures the case where
independent features yield $M=0$. Intuitively, it also means that if one
dimension is roughly a linear combination of the others, or if
$\lambda_{\max}(\Sigma^{-1})$ is large, the conditional distribution can
change quickly as a function of the conditioning variables.

###### Proof.

Consider a Gaussian random variable ${\bm{x}}$ with mean
$\mu\in\mathbb{R}^{d}$ and positive definite covariance
$\Sigma\in\mathbb{R}^{d\times d}$. We can partition the variable into two
parts, denoted ${\bm{x}}_{1}$ and ${\bm{x}}_{2}$, and then partition the
variable and parameters as follows:

$\displaystyle{\bm{x}}=\begin{pmatrix}{\bm{x}}_{1}\\\
{\bm{x}}_{2}\end{pmatrix}\quad\quad\quad\mu$
$\displaystyle=\begin{pmatrix}\mu_{1}\\\
\mu_{2}\end{pmatrix}\quad\quad\quad\Sigma=\begin{pmatrix}\Sigma_{11}&\Sigma_{12}\\\
\Sigma_{21}&\Sigma_{22}\end{pmatrix}.$

We can then consider the distribution for ${\bm{x}}_{1}$ conditioned on a
specific $x_{2}$ value. The conditional distribution is ${\bm{x}}_{1}\mid
x_{2}\sim{\mathcal{N}}(\bar{\mu},\bar{\Sigma})$ with the parameters defined as

$\displaystyle\bar{\mu}$
$\displaystyle=\mu_{1}+\Sigma_{12}\Sigma_{22}^{-1}(x_{2}-\mu_{2})\quad\quad\quad\bar{\Sigma}=\Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}.$

Our goal is to understand the total variation distance when conditioning on
different values $x_{2},x_{2}^{\prime}$. Pinsker’s inequality gives us the
following upper bound based on the KL divergence:

$\displaystyle d_{TV}\Big{(}p({\bm{x}}_{1}\mid x_{2}),p({\bm{x}}_{1}\mid
x_{2}^{\prime})\Big{)}\leq\sqrt{\frac{1}{2}D_{\mathrm{KL}}\Big{(}p({\bm{x}}_{1}\mid
x_{2})\mid\mid p({\bm{x}}_{1}\mid x_{2}^{\prime})\Big{)}}.$

The KL divergence between Gaussian distributions has a closed-form solution,
and in the case with equal covariance matrices reduces to the following:

$\displaystyle D_{\mathrm{KL}}\Big{(}p({\bm{x}}_{1}\mid x_{2})\mid\mid
p({\bm{x}}_{1}\mid
x_{2}^{\prime})\Big{)}=\frac{1}{2}(x_{2}-x_{2}^{\prime})^{\top}\Sigma_{22}^{-1}\Sigma_{21}\left(\Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}\right)^{-1}\Sigma_{12}\Sigma_{22}^{-1}(x_{2}-x_{2}^{\prime})$

To provide an upper bound on this term based on $\lVert
x_{2}-x_{2}^{\prime}\rVert_{2}$, it suffices to find the maximum eigenvalue of
the above matrix. To characterize the matrix, we observe that it is present in
the formula for the block inversion of the joint covariance matrix:

$\displaystyle\Sigma^{-1}=\begin{pmatrix}\left(\Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}\right)^{-1}&-\left(\Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}\right)^{-1}\Sigma_{12}\Sigma_{22}^{-1}\\\
-\Sigma_{22}^{-1}\Sigma_{21}\left(\Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}\right)^{-1}&\Sigma_{22}^{-1}+\Sigma_{22}^{-1}\Sigma_{21}\left(\Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}\right)^{-1}\Sigma_{12}\Sigma_{22}^{-1}\end{pmatrix}$

We therefore have the following, because the matrix we’re interested in is a
principal minor:

$\displaystyle\lambda_{\max}\Big{(}\Sigma_{22}^{-1}+\Sigma_{22}^{-1}\Sigma_{21}\left(\Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}\right)^{-1}\Sigma_{12}\Sigma_{22}^{-1}\Big{)}\leq\lambda_{\max}(\Sigma^{-1}).$

Weyl’s inequality also gives us:

$\displaystyle\lambda_{\min}(\Sigma_{22}^{-1})+\lambda_{\max}\Big{(}\Sigma_{22}^{-1}\Sigma_{21}\left(\Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}\right)^{-1}\Sigma_{12}\Sigma_{22}^{-1}\Big{)}$
$\displaystyle\leq\;$
$\displaystyle\lambda_{\max}\Big{(}\Sigma_{22}^{-1}+\Sigma_{22}^{-1}\Sigma_{21}\left(\Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}\right)^{-1}\Sigma_{12}\Sigma_{22}^{-1}\Big{)}$

Combining this with the previous inequality, we arrive at the following:

$\displaystyle\lambda_{\max}\Big{(}\Sigma_{22}^{-1}\Sigma_{21}\left(\Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}\right)^{-1}\Sigma_{12}\Sigma_{22}^{-1}\Big{)}\leq\lambda_{\max}(\Sigma^{-1})-\lambda_{\min}(\Sigma_{22}^{-1})$

For the subtractive term, we can write:

$\displaystyle\lambda_{\min}(\Sigma_{22}^{-1})=\frac{1}{\lambda_{\max}(\Sigma_{22})}\geq\frac{1}{\lambda_{\max}(\Sigma)}=\lambda_{\min}(\Sigma^{-1}).$

We therefore have the following upper bound on the KL divergence:

$\displaystyle D_{\mathrm{KL}}\Big{(}p({\bm{x}}_{1}\mid x_{2})\mid\mid
p({\bm{x}}_{1}\mid x_{2}^{\prime})\Big{)}\leq\frac{1}{2}\cdot\lVert
x_{2}-x_{2}^{\prime}\rVert_{2}^{2}\cdot\Big{(}\lambda_{\max}(\Sigma^{-1})-\lambda_{\min}(\Sigma^{-1})\Big{)}.$

Combining this with Pinsker’s inequality, we arrive at the final result:

$\displaystyle d_{TV}\big{(}p({\bm{x}}_{1}\mid x_{2}),p({\bm{x}}_{1}\mid
x_{2}^{\prime})\big{)}\big{)}\leq\frac{1}{2}\cdot\lVert
x_{2}-x_{2}^{\prime}\rVert_{2}\cdot\sqrt{\lambda_{\max}(\Sigma^{-1})-\lambda_{\min}(\Sigma^{-1})}.$

This means that we have
$M=\frac{1}{2}\sqrt{\lambda_{\max}(\Sigma^{-1})-\lambda_{\min}(\Sigma^{-1})}$.

Notice that this does not imply $M=0$ when the features ${\bm{x}}$ are
independent: it only captures this for the case of isotropic variance.
However, if we assume that ${\bm{x}}$ is standardized, so that each dimension
has variance equal to 1, $M=0$ is implied for independent features because we
have $\Sigma=I$.

∎

######  Lemma 3.

Under Assumption 3, the prediction function $f(x_{S})$ defined using the
conditional approach is Lipschitz continuous with constant $L+2BM$ for any
feature set $x_{S}$.

###### Proof.

The result follows directly from Lemma 2 and Assumption 3:

$\displaystyle|f(x_{S})-f(x_{S}^{\prime})|$ $\displaystyle\leq L\cdot\lVert
x_{S}-x_{S}^{\prime}\rVert_{2}+2B\cdot d_{TV}\Big{(}p(x_{\bar{S}}\mid
x_{S}),p(x_{\bar{S}}\mid x_{S}^{\prime})\Big{)}$ $\displaystyle\leq
L\cdot\lVert x_{S}-x_{S}^{\prime}\rVert_{2}+2BM\cdot\lVert
x_{S}-x_{S}^{\prime}\rVert_{2}$ $\displaystyle=(L+2BM)\cdot\lVert
x_{S}-x_{S}^{\prime}\rVert_{2}.$

∎

##  Appendix C Proofs: prediction robustness to model perturbations

We now re-state and prove our results from Section 3.2.

######  Lemma 4.

For two models $f,f^{\prime}:\mathbb{R}^{d}\mapsto\mathbb{R}$ and a subdomain
${\mathcal{X}}^{\prime}\subseteq\mathbb{R}^{d}$, the prediction functions
$f(x_{S}),f^{\prime}(x_{S})$ for any feature set $x_{S}$ satisfy

$|f(x_{S})-f^{\prime}(x_{S})|\leq\lVert
f-f^{\prime}\rVert_{\infty,{\mathcal{X}}^{\prime}}\cdot
Q_{x_{S}}({\mathcal{X}}^{\prime})+2B\cdot\big{(}1-Q_{x_{S}}({\mathcal{X}}^{\prime})\big{)},$

where $Q_{x_{S}}({\mathcal{X}}^{\prime})$ is the probability of imputed
samples lying in ${\mathcal{X}}^{\prime}$ based on the distribution
$q({\bm{x}}_{\bar{S}})$:

$Q_{x_{S}}({\mathcal{X}}^{\prime})\equiv\mathbb{E}_{q({\bm{x}}_{\bar{S}})}\left[{\mathbb{I}}\left\\{(x_{S},{\bm{x}}_{\bar{S}})\in{\mathcal{X}}^{\prime}\right\\}\right].$

###### Proof.

The main idea here is that when integrating over the difference in function
outputs, we can separate the integral over $(x_{S},{\bm{x}}_{\bar{S}})$ values
falling in the domains ${\mathcal{X}}^{\prime}$ and
$\mathbb{R}^{d}\setminus{\mathcal{X}}^{\prime}$. Because we have a single
value for the observed variables $x_{S}$, the distribution
$q({\bm{x}}_{\bar{S}})$ is identical for both terms regardless of the removal
technique.

$\displaystyle|f(x_{S})-f^{\prime}(x_{S})|$ $\displaystyle=\Big{|}\int
f(x_{S},x_{\bar{S}})q(x_{\bar{S}})dx_{\bar{S}}-\int
f^{\prime}(x_{S},x_{\bar{S}})q(x_{\bar{S}})dx_{\bar{S}}\Big{|}$
$\displaystyle\leq\int\big{|}f(x_{S},x_{\bar{S}})-f^{\prime}(x_{S},x_{\bar{S}})\big{|}\cdot
q(x_{\bar{S}})dx_{\bar{S}}$
$\displaystyle=\int\big{|}f(x_{S},x_{\bar{S}})-f^{\prime}(x_{S},x_{\bar{S}})\big{|}\cdot{\mathbb{I}}\left\\{(x_{S},{\bm{x}}_{\bar{S}})\in{\mathcal{X}}^{\prime}\right\\}\cdot
q(x_{\bar{S}})dx_{\bar{S}}$
$\displaystyle\quad+\int\big{|}f(x_{S},x_{\bar{S}})-f^{\prime}(x_{S},x_{\bar{S}})\big{|}\cdot{\mathbb{I}}\left\\{(x_{S},{\bm{x}}_{\bar{S}})\notin{\mathcal{X}}^{\prime}\right\\}\cdot
q(x_{\bar{S}})dx_{\bar{S}}$ $\displaystyle\leq\lVert
f-f^{\prime}\rVert_{\infty,{\mathcal{X}}^{\prime}}\int{\mathbb{I}}\left\\{(x_{S},{\bm{x}}_{\bar{S}})\in{\mathcal{X}}^{\prime}\right\\}\cdot
q(x_{\bar{S}})dx_{\bar{S}}$
$\displaystyle\quad+2B\cdot\int{\mathbb{I}}\left\\{(x_{S},{\bm{x}}_{\bar{S}})\notin{\mathcal{X}}^{\prime}\right\\}\cdot
q(x_{\bar{S}})dx_{\bar{S}}$ $\displaystyle=\lVert
f-f^{\prime}\rVert_{\infty,{\mathcal{X}}^{\prime}}\cdot
Q_{x_{S}}({\mathcal{X}}^{\prime})+2B\cdot\big{(}1-Q_{x_{S}}({\mathcal{X}}^{\prime})\big{)}$

∎

######  Lemma 5.

When removing features using the conditional approach, the prediction
functions $f(x_{S}),f^{\prime}(x_{S})$ for two models $f,f^{\prime}$ and any
feature set $x_{S}$ satisfy

$|f(x_{S})-f^{\prime}(x_{S})|\leq\lVert
f-f^{\prime}\rVert_{\infty,{\mathcal{X}}}.$

###### Proof.

The result follows directly from Lemma 4 when we set
${\mathcal{X}}^{\prime}={\mathcal{X}}$. Note that when we remove features
using the conditional distribution, we have $Q_{x_{S}}({\mathcal{X}})=1$ for
all $x_{S}$ lying on the data manifold, or where there exists $x_{\bar{S}}$
such that $p(x_{S},x_{\bar{S}})>0$. ∎

######  Lemma 6.

When removing features using the baseline or marginal approach, the prediction
functions $f(x_{S}),f^{\prime}(x_{S})$ for two models $f,f^{\prime}$ and any
feature set $x_{S}$ satisfy

$|f(x_{S})-f^{\prime}(x_{S})|\leq\lVert f-f^{\prime}\rVert_{\infty}.$

###### Proof.

The result follows directly from Lemma 4 when we set
${\mathcal{X}}^{\prime}=\mathbb{R}^{d}$, because we have
$Q_{x_{S}}(\mathbb{R}^{d})=1$ for all $x_{S}\in\mathbb{R}^{|S|}$. ∎

##  Appendix D Proofs: summary technique robustness

We now re-state and prove our results from Section 3.3.

######  Proposition 1.

The attributions for each method can be calculated by applying a linear
operator $A\in\mathbb{R}^{d\times 2^{d}}$ to a vector $v\in\mathbb{R}^{2^{d}}$
representing the predictions with each feature set, or

$\phi(f,x)=Av,$

where the linear operator $A$ for each method is listed in Table 3, and $v$ is
defined as $v_{S}=f(x_{S})$ for each $S\subseteq[d]$ based on the chosen
feature removal technique.

###### Proof.

This result follows directly from the definition of each method. Following the
approach of Covert et al. [17], we disentangle between each method’s feature
removal implementation and how it generates attribution scores. Assuming a
fixed input $x$ where $f(x_{S})$ denotes the prediction with a feature set,
each method’s attributions are defined as follows.

  * •

Occlusion [67]: this method simply compares the prediction given all features
and the prediction with a single missing feature. The attribution score
$\phi_{i}(f,x)$ is defined as

$\phi_{i}(f,x)=f(x)-f(x_{[d]\setminus\\{i\\}}).$

We refer to this summary technique as “leave-one-out” in Table 3.

  * •

Shapley value [54, 45]: this method calculates the impact of removing a single
feature, and then averages this over all possible preceding subsets. The
attribution score $\phi_{i}(f,x)$ is defined as

$\phi_{i}(f,x)=\frac{1}{d}\sum_{S\subseteq[d]\setminus\\{i\\}}\binom{d-1}{|S|}^{-1}\Big{(}f(x_{S\cup\\{i\\}})-f(x_{S})\Big{)}.$

  * •

Banzhaf value [9, 13]: this method is similar to the Shapley value, but it
applies a uniform weighting when averaging over all preceding subsets. The
attribution score $\phi_{i}(f,x)$ is defined as

$\phi_{i}(f,x)=\frac{1}{2^{d-1}}\sum_{S\subseteq[d]\setminus\\{i\\}}\Big{(}f(x_{S\cup\\{i\\}})-f(x_{S})\Big{)}.$

  * •

RISE [51]: this method is similar to the Banzhaf value, but it does not
subtract the value for subsets where a feature is not included. The
attribution score $\phi_{i}(f,x)$ is defined as

$\phi_{i}(f,x)=\frac{1}{2^{d-1}}\sum_{S\subseteq[d]\setminus\\{i\\}}f(x_{S\cup\\{i\\}}).$

We refer to this summary technique as “mean when included” in Table 3

  * •

LIME [53]: this method is more flexible than the others and generally involves
solving a regression problem that treats the features as inputs and the
predictions as labels. LIME can be viewed as a removal-based explanation when
the features are represented in a binary fashion, indicating whether they are
provided to the model or not [17]. There is also a choice of weighting kernel
$w(S)\geq 0$, of a regularization term, and of whether to fit an intercept
term.

We assume that the regularization is omitted and that we fit an intercept
term. Denoting the model’s coefficients as $(\beta_{0},\ldots,\beta_{d})$, the
problem is the following:

$\min_{\beta_{0},\ldots,\beta_{d}}\;\sum_{S\subseteq[d]}w(S)\Big{(}\beta_{0}+\sum_{i\in
S}\beta_{i}-f(x_{S})\Big{)}^{2}.$

We can rewrite this problem in a simpler fashion to derive its solution.
First, we let $W\in[0,1]^{2^{d}\times 2^{d}}$ be a diagonal matrix containing
the weights $w(S)$ for each subset. Next, we let the predictions be
represented by a vector $v\in\mathbb{R}^{2^{d}}$ where $v_{S}=f(x_{S})$, and
we let the matrix $Z\in\\{0,1\\}^{2^{d}\times d}$ enumerate all subsets in a
binary fashion. Finally, to accommodate the intercept term, we let
$\beta=(\beta_{0},\ldots,\beta_{d})\in\mathbb{R}^{d+1}$ represent the
parameters, and we prepend a columns of ones to $Z$, resulting in
$Z^{\prime}\in\\{0,1\\}^{2^{d}\times(d+1)}$. The problem then becomes:

$\min_{\beta}\;(Z^{\prime}\beta-v)^{\top}W(Z^{\prime}\beta-v^{\prime})=\big{\lVert}Z^{\prime}\beta-v\big{\rVert}_{W}^{2}.$

Assuming that $Z^{\prime\top}WZ^{\prime}$ is invertible, which is guaranteed
when $w(S)>0$ for all $S\subseteq[d]$, this problem has a closed-form
solution:

$\beta^{*}=(Z^{\prime\top}WZ^{\prime})^{-1}Z^{\prime\top}Wv.$

LIME’s attribution scores are given by $\phi_{i}(f,x)=\beta^{*}_{i}$, for
$i=1,...,d$, where we discard the intercept term. We can therefore conclude
that each attribution is a linear function of $v$, with coefficients given by
all but the first row of $(Z^{\prime\top}WZ^{\prime})^{-1}Z^{\prime\top}W$. If
we denote all but the first row of $(Z^{\prime\top}WZ^{\prime})^{-1}$ as
$(Z^{\prime\top}WZ^{\prime})^{-1}_{[1:]}$, then we have
$A=(Z^{\prime\top}WZ^{\prime})^{-1}_{[1:]}Z^{\prime\top}W$.

The formula for the attributions remains roughly the same if we do not fit an
intercept term: if the intercept is zero, we have
$A=(Z^{\top}WZ)^{-1}Z^{\top}W$. We can similarly allow for ridge
regularization: given a penalty parameter $\lambda>0$, we then have
$A=(Z^{\top}WZ+\lambda I_{2^{d}})^{-1}Z^{\top}W$.

These coefficients depend on the specific choice of weighting kernel $w(S)$,
and previous work has shown that leave-one-out, Shapley and Banzhaf values all
correspond to specific weighting kernels [17]. Our analysis of those methods
therefore pertains to LIME, at least in certain special cases. However, LIME
often uses a heuristic exponential kernel in practice, which is more difficult
to characterize analytically. We therefore omit the specific entries and norms
for LIME in Table 3, but show norms for the exponential kernel in our computed
results.

We refer to this summary technique as “weighted least squares” in Table 3.

We can therefore see that in each method, $\phi_{i}(f,x)$ is a linear function
of the predictions for each $S\subseteq[d]$. Across all methods, we can let
the predictions be represented by a vector $v\in\mathbb{R}^{2^{d}}$ and the
coefficients by a matrix $A\in\mathbb{R}^{d\times 2^{d}}$. ∎

######  Lemma 7.

The difference in attributions given the same summary technique $A$ and
different model outputs $v,v^{\prime}$ can be bounded as

$\lVert Av-Av^{\prime}\rVert_{2}\leq\lVert A\rVert_{2}\cdot\lVert
v-v^{\prime}\rVert_{2}\quad\quad\text{or}\quad\quad\lVert Av-
Av^{\prime}\rVert_{2}\leq\lVert A\rVert_{1,\infty}\cdot\lVert
v-v^{\prime}\rVert_{\infty},$

where $\lVert A\rVert_{2}$ is the spectral norm, and the operator norm $\lVert
A\rVert_{1,\infty}$ is the square root of the sum of squared row 1-norms, with
values for each $A$ given in Table 3.

###### Proof.

We first derive the two matrix inequalities. The first inequality follows from
the definition of the spectral norm, which is defined as the maximum singular
value:

$\displaystyle\lVert A\rVert_{2}$ $\displaystyle\equiv\sup_{u\neq
0}\frac{\lVert Au\rVert_{2}}{\lVert
u\rVert_{2}}=\sigma_{\max}(A)=\sqrt{\lambda_{\max}(AA^{\top})}.$

We prove the second bound using Holder’s inequality as follows, where $A_{i}$
denotes the $i$th row:

$\displaystyle\lVert Au\rVert_{2}$
$\displaystyle=\sqrt{\sum_{i=1}^{d}(A_{i}^{\top}u)^{2}}\leq\sqrt{\sum_{i=1}^{d}\lVert
A_{i}\rVert_{1}^{2}\cdot\lVert u\rVert_{\infty}^{2}}=\lVert
u\rVert_{\infty}\cdot\sqrt{\sum_{i=1}^{d}\lVert A_{i}\rVert_{1}^{2}}=\lVert
u\rVert_{\infty}\cdot\lVert A\rVert_{1,\infty}.$

In the special case where all of $A$’s row 1-norms are equal, or where we have
$\lVert A_{i}\rVert_{1}=\lVert A_{j}\rVert_{1}$ for all $i,j\in[d]$, we can
use the simplified expression $\lVert A\rVert_{1,\infty}=\lVert
A_{1}\rVert_{1}\cdot\sqrt{d}$.

Next, we derive the specific norm values for each linear operator in Table 3.

Beginning with the operator norm $\lVert A\rVert_{1,\infty}$, we remark that
several methods have equal 1-norms across all rows: following their
interpretation as semivalues [46], leave-one-out, Shapley and Banzhaf all have
$\lVert A_{i}\rVert_{1}=2$ for all $i\in[d]$. Similarly, RISE has $\lVert
A_{i}\rVert_{1}=1$ for all $i\in[d]$. This yields the norms $\lVert
A\rVert_{1,\infty}=2\sqrt{d}$ for leave-one-out, Banzhaf and Shapley, and
$\lVert A\rVert_{1,\infty}=\sqrt{d}$ for RISE.

For the spectral norm, we consider each method separately. We find that for
all methods, it is simplest to calculate $AA^{\top}$ and then
$\lambda_{\max}(AA^{\top})$ to find the maximum singular value
$\sigma_{\max}(A)$. Note that the entries of the matrix are determined by the
inner product $(AA^{\top})_{ij}=A_{i}^{\top}A_{j}$.

For leave-one-out, we derive $AA^{\top}$ as follows:

$\displaystyle(AA^{\top})_{ij}=\begin{cases}2&\text{if}\;i=j\\\
1&\text{otherwise}.\end{cases}$

The matrix $AA^{\top}$ is therefore the sum of a rank-one component and an
identity term, and we have $AA^{\top}={\bm{1}}_{d}{\bm{1}}_{d}+I_{d}$.
Recognizing that ${\bm{1}}_{d}$ is the leading eigenvector, we get the
following result:

$\displaystyle\sigma_{\max}(A)=\sqrt{\lambda_{\max}(AA^{\top})}=\sqrt{\frac{({\bm{1}}_{d}{\bm{1}}_{d}+I_{d}){\bm{1}}_{d}}{\lVert{\bm{1}}_{d}\rVert_{2}}}=\sqrt{\frac{{\bm{1}}_{d}d+{\bm{1}}_{d}}{\lVert{\bm{1}}_{d}\rVert_{2}}}=\sqrt{d+1}.$

For the Banzhaf value, we derive $AA^{\top}$ as follows:

$\displaystyle(AA^{\top})_{ij}=\begin{cases}1/2^{2d-2}&\text{if}\;i=j\\\
0&\text{otherwise}.\end{cases}$

The max eigenvalue is therefore $\lambda_{\max}(AA^{\top})=1/2^{d-2}$, which
yields $\lVert A\rVert_{2}=1/2^{d/2-1}$.

For the mean when included technique used by RISE, we derive $AA^{\top}$ as
follows:

$\displaystyle(AA^{\top})_{ij}=\begin{cases}1/2^{d-1}&\text{if}\;i=j\\\
1/2^{d}&\text{otherwise}.\end{cases}$

The matrix $AA^{\top}$ is therefore the sum of a rank-one component and an
identity term, or we have
$AA^{\top}={\bm{1}}_{d}{\bm{1}}_{d}/2^{d}+I_{d}/2^{d}$. Recognizing that
${\bm{1}}_{d}$ is the leading eigenvector, we get the following result:

$\displaystyle\sigma_{\max}(A)=\sqrt{\lambda_{\max}(AA^{\top})}=\sqrt{\frac{({\bm{1}}_{d}{\bm{1}}_{d}/2^{d}+I_{d}/2^{d}){\bm{1}}_{d}}{\lVert{\bm{1}}_{d}\rVert_{2}}}=\sqrt{\frac{{\bm{1}}_{d}d/2^{d}+{\bm{1}}_{d}/2^{d}}{\lVert{\bm{1}}_{d}\rVert_{2}}}=\sqrt{\frac{d+1}{2^{d}}}.$

For the Shapley value, we derive $AA^{\top}$ by first considering the diagonal
entries:

$\displaystyle(AA^{\top})_{ii}=A_{i}^{\top}A_{i}$ $\displaystyle=\sum_{S:i\in
S}\left(\frac{(|S|-1)!(d-|S|)!}{d!}\right)^{2}+\sum_{S:i\notin
S}\left(\frac{|S|!(d-|S|-1)!}{d!}\right)^{2}$
$\displaystyle=\sum_{k=0}^{d}\left[\binom{d-1}{k-1}\left(\frac{(k-1)!(d-k)!}{d!}\right)^{2}+\binom{d-1}{k}\left(\frac{k!(d-k-1)!}{d!}\right)^{2}\right]$
$\displaystyle=\frac{1}{d}\sum_{k=1}^{d}\frac{(k-1)!(d-k)!}{d!}+\frac{1}{d}\sum_{k=0}^{d-1}\frac{k!(d-k-1)!}{d!}$
$\displaystyle=\sum_{k=1}^{d-1}\frac{(k-1)!(d-k-1)!}{d!}+\underbrace{\frac{2}{d^{2}}}_{k=0,k=d}$
$\displaystyle=\frac{1}{d(d-1)}\sum_{k=1}^{d-1}\binom{d-2}{k-1}^{-1}+\frac{2}{d^{2}}$

We next consider the off-diagonal entries:

$\displaystyle(AA^{\top})_{ij}=A_{i}^{\top}A_{j}$ $\displaystyle=\sum_{S:i\in
S,j\in S}\frac{(|S|-1)!(d-|S|)!}{d!}\frac{(|S|-1)!(d-|S|)!}{d!}$
$\displaystyle\quad-\sum_{S:i\in S,j\notin
S}\frac{(|S|-1)!(d-|S|)!}{d!}\frac{|S|!(d-|S|-1)!}{d!}$
$\displaystyle\quad-\sum_{S:i\notin S,j\in
S}\frac{|S|!(d-|S|-1)!}{d!}\frac{(|S|-1)!(d-|S|)!}{d!}$
$\displaystyle\quad+\sum_{S:i\notin S,j\notin
S}\frac{|S|!(d-|S|-1)!}{d!}\frac{|S|!(d-|S|-1)!}{d!}$
$\displaystyle=\sum_{k=0}^{d}\binom{d-2}{k-2}\frac{(k-1)!(d-k)!}{d!}\frac{(k-1)!(d-k)!}{d!}$
$\displaystyle\quad-2\sum_{k=0}^{d}\binom{d-2}{k-1}\frac{(k-1)!(d-k)!}{d!}\frac{k!(d-k-1)!}{d!}$
$\displaystyle\quad+\sum_{k=0}^{d}\binom{d-2}{k}\frac{k!(d-k-1)!}{d!}\frac{k!(d-k-1)!}{d!}$
$\displaystyle=\sum_{k=2}^{d}\frac{(k-1)}{d(d-1)}\frac{(k-1)!(d-k)!}{d!}$
$\displaystyle\quad-2\sum_{k=1}^{d-1}\frac{1}{d(d-1)}\frac{k!(d-k)!}{d!}$
$\displaystyle\quad+\sum_{k=0}^{d-2}\frac{(d-k-1)}{d(d-1)}\frac{k!(d-k-1)!}{d!}$

We first focus on the case where $d>2$, and we later return to the case where
$d=2$. We can group terms that share the same $k$ values to simplify the
above:

$\displaystyle A_{i}^{\top}A_{j}$
$\displaystyle=\frac{1}{d(d-1)}\left[-\sum_{k=2}^{d-2}\frac{(k-1)!(d-k-1)!}{(d-1)!}+\underbrace{\frac{2(d-1)}{d}}_{k=0,k=d}-\underbrace{\frac{2}{d-1}}_{k=1,k=d-1}\right]$
$\displaystyle=-\frac{1}{d(d-1)}\sum_{k=2}^{d-2}\frac{(k-1)!(d-k-1)!}{(d-1)!}+\frac{2}{d^{2}}-\frac{2}{d(d-1)^{2}}$
$\displaystyle=-\frac{1}{d(d-1)^{2}}\sum_{k=2}^{d-2}\binom{d-2}{k-1}^{-1}+\frac{2}{d^{2}}-\frac{2}{d(d-1)^{2}}.$

By convention, we let the summation above evaluate to zero for $d=3$. We
therefore have the following result for $AA^{\top}$:

$\displaystyle(AA^{\top})_{ij}=\begin{cases}\frac{1}{d(d-1)}\sum_{k=1}^{d-1}\binom{d-2}{k-1}^{-1}+\frac{2}{d^{2}}&\text{if}\;i=j\\\
-\frac{1}{d(d-1)^{2}}\sum_{k=2}^{d-2}\binom{d-2}{k-1}^{-1}+\frac{2}{d^{2}}-\frac{2}{d(d-1)^{2}}&\text{otherwise}.\end{cases}$

From this, we can see that $AA^{\top}$ is the sum of a rank-one component and
an identity term, or that we have
$AA^{\top}={\bm{1}}_{d}{\bm{1}}_{d}b+I_{d}(c-b)$ with values for
$b=(AA^{\top})_{ij}$ and $c=(AA^{\top})_{ii}$ given above. This matrix has two
eigenvalues, $c-b$ and $db+(c-b)$, and we must consider which of the two is
larger. Under our assumption of $d>2$, we can observe that $b>0$:

$\displaystyle b=(AA^{\top})_{ij}$
$\displaystyle=\frac{2}{d^{2}}-\frac{1}{d(d-1)^{2}}\sum_{k=2}^{d-2}\binom{d-2}{k-1}^{-1}-\frac{2}{d(d-1)^{2}}$
$\displaystyle\geq\frac{2}{d^{2}}-\frac{1}{d(d-1)^{2}}\sum_{k=2}^{d-2}\binom{d-2}{1}^{-1}-\frac{2}{d(d-1)^{2}}$
$\displaystyle=\frac{2}{d^{2}}-\frac{d-3}{d(d-1)^{2}(d-2)}-\frac{2}{d(d-1)^{2}}$
$\displaystyle=\frac{2(d-1)^{2}(d-2)-d(d-3)-2d(d-2)}{d^{2}(d-1)^{2}(d-2)}$

It can be verified that the numerator of the above fraction is positive for
$d>2$. This implies that $db+(c-b)$ is the larger eigenvalue. We can therefore
calculate $\lambda_{\max}(AA^{\top})$ as follows:

$\displaystyle\lambda_{\max}(AA^{\top})$ $\displaystyle=db+(c-b)=c+(d-1)b$
$\displaystyle=\frac{2}{d^{2}}+\underbrace{\frac{2}{d(d-1)}}_{k=1,k=d-1}+\frac{2(d-1)}{d^{2}}-\frac{2}{d(d-1)}$
$\displaystyle=\frac{2(d-1)+2(d-1)^{2}}{d^{2}(d-1)}$
$\displaystyle=\frac{2}{d}$

On the other hand, when we have $d=2$ we have the following entries for
$AA^{\top}$:

$\displaystyle(AA^{\top})_{ij}=\begin{cases}1&\text{if}\;i=j\\\
0&\text{otherwise}.\end{cases}$

This yields $\lambda_{\max}(AA^{\top})=1$, which coincides with the formula
$\lambda_{\max}(AA^{\top})=2/d$ derived above. Finally, this let us conclude
that the spectral norm for the Shapley value is the following:

$\displaystyle\lVert
A\rVert_{2}=\sqrt{\lambda_{\max}(AA^{\top})}=\sqrt{\frac{2}{d}}.$

∎

######  Lemma 8.

When the linear operator $A$ satisfies the boundedness property, we have
$\lVert A\rVert_{1,\infty}=2\sqrt{d}$.

###### Proof.

Following Theorem 1 in Monderer and Samet [46], solution concepts satisfying
linearity and the boundedness property (also known as the Milnor axiom) are
probabilistic values. This means that each attribution is the average of a
feature’s marginal contributions, or that for each feature $i\in[d]$ there
exists a distribution $p_{i}(S)$ with support on the power set of
$[d]\setminus\\{i\\}$, such that the attributions are given by

$\displaystyle\phi_{i}(f,x)=\mathbb{E}_{p_{i}(S)}\big{[}f(x_{S\cup\\{i\\}})-f(x_{S})\big{]}=\sum_{S\subseteq[d]\setminus\\{i\\}}p_{i}(S)\big{(}f(x_{S\cup\\{i\\}})-f(x_{S})\big{)}.$
(D.0.1)

The corresponding matrix $A\in\mathbb{R}^{d\times 2^{d}}$ therefore has rows
$A_{i}$ with entries $A_{iS}\geq 0$ when $i\in S$ and $A_{iS}\leq 0$ when
$i\notin S$, where the positive and negative entries sum to 1 and -1
respectively. We can therefore conclude that each row satisfies $\lVert
A_{i}\rVert_{1}=2$, and that $\lVert A\rVert_{1,\infty}=2\sqrt{d}$.

∎

######  Lemma 9.

When the linear operator $A$ satisfies the boundedness and symmetry
properties, the spectral norm is bounded as follows,

$\frac{1}{2^{d/2-1}}\leq\lVert A\rVert_{2}\leq\sqrt{d+1},$

with $1/2^{d/2-1}$ achieved by the Banzhaf value and $\sqrt{d+1}$ achieved by
leave-one-out.

###### Proof.

Following Monderer and Samet [46], solution concepts satisfying linearity,
boundedness and symmetry are semivalues. This means that they are
probabilistic values for which the probability distribution $p_{i}(S)$ depends
only on the cardinality $|S|$ (see the definition of probabilistic values
above). As a result, we can parameterize the summary technique by how much
weight it places on each cardinality.

Let $\alpha\in\mathbb{R}^{d}$ denote a probability distribution over the
cardinalities $k=1,\ldots,d$. We can write the attributions generated by any
semivalue as follows:

$\displaystyle\phi_{i}(f,x)=\sum_{k=1}^{d}\alpha_{k}\sum_{\begin{subarray}{c}S\subseteq[d]\setminus\\{i\\}\\\
|S|=k-1\end{subarray}}\binom{d-1}{k-1}^{-1}\Big{(}f(x_{S\cup\\{i\\}})-f(x_{S})\Big{)}.$

Note that this resembles eq. D.0.1, only with coefficients that are shared
between all subsets with the same cardinality. Due to each attribution being a
linear combination of the predictions, we can associate a matrix
$A\in\mathbb{R}^{d\times 2^{d}}$ with each semivalue, similar to Proposition
1. Furthermore, we can let $A$ be a linear combination of basis elements
corresponding to each cardinality. We define the basis matrix
$A^{(k)}\in\mathbb{R}^{d\times 2^{d}}$ for each cardinality $k=1,\ldots,d$ as
follows:

$\displaystyle
A^{(k)}_{iS}=\begin{cases}\binom{d-1}{|S|-1}^{-1}&\text{if}\;i\in
S\;\text{and}\;|S|=k\\\ -\binom{d-1}{|S|}^{-1}&\text{if}\;i\notin
S\;\text{and}\;|S|=k-1\\\ 0&\text{otherwise}.\end{cases}$

With this, we can write the matrix associated with any semivalue as a function
of the cardinality weights $\alpha$:

$\displaystyle A(\alpha)=\sum_{k=1}^{d}\alpha_{k}A^{(k)}.$

An important insight is that by composing a linear function with the convex
spectral norm function [11], we can see that $\lVert A(\alpha)\rVert_{2}$ is a
convex function of $\alpha$. This lets us find both the upper and lower bound
of $\lVert A(\alpha)\rVert_{2}$ subject to $\alpha$ being a valid probability
distribution.

For the lower bound, we refer readers to Theorem C.3 from Wang and Jia [64],
who prove that the Banzhaf value achieves the minimum with $\alpha^{*}$ such
that $\alpha^{*}_{k}=\binom{d-1}{k-1}/2^{d-1}$ and spectral norm $\lVert
A(\alpha^{*})\rVert_{2}=1/2^{d/2-1}$. Notice that this is a valid probability
distribution because we have $\sum_{k=1}^{d}\binom{d-1}{k-1}=2^{d-1}$. Also
note that the proof in [64] parameterizes $\alpha$ slightly differently.

For the upper bound, we can leverage convexity to make the following argument.
The fact that spectral norm is convex with $\alpha$ representing probabilities
lets us write the following for any $\alpha$ value:

$\displaystyle\sigma_{\max}\Big{(}A(\alpha)\Big{)}=\sigma_{\max}\left(\sum_{k=1}^{d}\alpha_{k}\cdot
A^{(k)}\right)\leq\sum_{k=1}^{d}\alpha_{k}\cdot\sigma_{\max}\big{(}A^{(k)}\big{)}\leq\max_{k}\Big{\\{}\sigma_{\max}\big{(}A^{(k)}\big{)}\Big{\\}}.$

Finding the upper bound therefore reduces to finding the basis element with
the largest spectral norm. For each basis element $A^{(k)}$, we can derive the
spectral norm by first finding the maximum eigenvalue of
$A^{(k)}{A^{(k)}}^{\top}\in\mathbb{R}^{d\times d}$. In deriving this matrix,
we first consider the diagonal entries:

$\displaystyle\big{(}A^{(k)}{A^{(k)}}^{\top}\big{)}_{ii}$
$\displaystyle=\sum_{S:i\in S,|S|=k}\binom{d-1}{|S|-1}^{-2}+\sum_{S:i\notin
S,|S|=k-1}\binom{d-1}{|S|}^{-2}$
$\displaystyle=\binom{d-1}{k-1}\binom{d-1}{k-1}^{-2}+\binom{d-1}{k-1}\binom{d-1}{k-1}^{-2}$
$\displaystyle=2\binom{d-1}{k-1}^{-1}.$

Next, we consider the off-diagonal entries:

$\displaystyle\big{(}A^{(k)}{A^{(k)}}^{\top}\big{)}_{ij}$
$\displaystyle=\sum_{S:i\in S,j\in
S,|S|=k}\binom{d-1}{|S|-1}^{-2}+\sum_{S:i\notin S,j\notin
S,|S|=k-1}\binom{d-1}{|S|}^{-2}$
$\displaystyle=\binom{d-2}{k-2}\binom{d-1}{k-1}^{-2}+\binom{d-2}{k-1}\binom{d-1}{k-1}^{-2}$
$\displaystyle=\frac{k-1}{d-1}\binom{d-1}{k-1}^{-1}+\frac{d-k}{d-1}\binom{d-1}{k-1}^{-1}$
$\displaystyle=\binom{d-1}{k-1}^{-1}.$

To summarize, the entries of this matrix are the following:

$\displaystyle\big{(}A^{(k)}{A^{(k)}}^{\top}\big{)}_{ij}=\begin{cases}2\binom{d-1}{k-1}^{-1}&\text{if}\;i=j\\\
\binom{d-1}{k-1}^{-1}&\text{otherwise.}\end{cases}$

We therefore see that $A^{(k)}{A^{(k)}}^{\top}$ is equal to a rank-one matrix
plus an identity matrix, or
$A^{(k)}{A^{(k)}}^{\top}=b{\bm{1}}_{d}{\bm{1}}_{d}^{\top}+(c-b)I_{d}$, with
values for $b=\big{(}A^{(k)}{A^{(k)}}^{\top}\big{)}_{ij}$ and
$c=\big{(}A^{(k)}{A^{(k)}}^{\top}\big{)}_{ii}$ given above. Recognizing that
$b>0$ and therefore that ${\bm{1}}_{d}$ is the leading eigenvector, we have
the following expression for the maximum eigenvalue:

$\displaystyle\lambda_{\max}\big{(}A^{(k)}{A^{(k)}}^{\top}\big{)}$
$\displaystyle=c+(d-1)b$
$\displaystyle=\big{(}A^{(k)}{A^{(k)}}^{\top}\big{)}_{ii}+(d-1)\cdot\big{(}A^{(k)}{A^{(k)}}^{\top}\big{)}_{ij}$
$\displaystyle=(d+1)\binom{d-1}{k-1}^{-1}.$

Finally, we must consider how to maximize this as a function of $k$. Setting
$k=1$ or $k=d$ achieves the same value of $d+1$, which yields the following
spectral norm:

$\displaystyle\lVert
A^{(d)}\rVert_{2}=\sqrt{\lambda_{\max}\big{(}A^{(d)}{A^{(d)}}^{\top}\big{)}}=\sqrt{d+1}.$

The case with $k=d$ corresponds to leave-one-out (see Lemma 7), and $k=1$
corresponds to an analogous leave-one-in approach that has been discussed in
prior work [15, 17]. ∎

######  Lemma 10.

When the linear operator $A$ satisfies the boundedness and efficiency
properties, the spectral norm is bounded as follows,

$\sqrt{\frac{2}{d}}\leq\lVert
A\rVert_{2}\leq\sqrt{2+2\cdot\cos\left(\frac{\pi}{d+1}\right)},$

with $\sqrt{2/d}$ achieved by the Shapley value.

###### Proof.

Following Monderer and Samet [46], solution concepts satisfying linearity,
boundedness and efficiency are random-order values (or quasivalues). This
means that the attributions are the average marginal contributions across a
distribution of feature orderings. Reasoning about such orderings requires new
notation, so we use $\pi$ to denote an ordering over $[d]$, and $\pi(i)$ to
denote the position of index $i$ within the ordering. We also let
$\text{Pre}(\pi,i)$ denote the set of all elements preceding $i$ in the
ordering, or $\text{Pre}(\pi,i)\equiv\\{j:\pi(j)<\pi(i)\\}$. Similarly, we let
$\text{Pre}(\pi,i,j)$ be an indicator of whether $j$ directly precedes $i$ in
the ordering, or $\text{Pre}(\pi,i,j)\equiv{\mathbb{I}}\\{\pi(j)=\pi(i)-1\\}$.
Finally, we use $\Pi$ to denote the set of all orderings, where $|\Pi|=d!$.

Now, given a distribution over orderings $p(\pi)$, we can write the
attributions for any random-order value as follows:

$\displaystyle\phi_{i}(f,x)$
$\displaystyle=\sum_{\pi\in\Pi}p(\pi)\Big{(}f(x_{\text{Pre}(\pi,i)\cup\\{i\\}})-f(x_{\text{Pre}(\pi,i)})\Big{)}$
$\displaystyle=\sum_{S\subseteq[d]\setminus\\{i\\}}\Big{(}\sum_{\pi:\text{Pre}(\pi,i)=S}p(\pi)\Big{)}\Big{(}f(x_{S\cup\\{i\\}})-f(x_{S})\Big{)}.$

Based on this formulation, we can associate a matrix $A\in\mathbb{R}^{d\times
2^{d}}$ with each random-order value based on the coefficients applied to each
prediction. The entries are defined as follows:

$\displaystyle
A_{iS}=\begin{cases}\sum_{\pi:\text{Pre}(\pi,i)=S\setminus\\{i\\}}p(\pi)&\text{if}\;i\in
S\\\ -\sum_{\pi:\text{Pre}(\pi,i)=S}p(\pi)&\text{otherwise}.\end{cases}$

Furthermore, we can view $A$ as a linear combination of basis elements
$A^{(\pi)}\in\mathbb{R}^{d\times 2^{d}}$ corresponding to each individual
ordering. We define the basis matrix $A^{(\pi)}$ associated with each
$\pi\in\Pi$ as follows:

$\displaystyle
A^{(\pi)}_{iS}=\begin{cases}1&\text{if}\;\text{Pre}(\pi,i)=S\setminus\\{i\\}\\\
-1&\text{if}\;\text{Pre}(\pi,i)=S\\\ 0&\text{otherwise}.\end{cases}$

With this, we can write the linear operator $A$ corresponding to a
distribution $p(\pi)$ as the following linear combination:

$\displaystyle A(p)=\sum_{\pi\in\Pi}p(\pi)\cdot A^{(\pi)}.$

By composing a linear operation with the convex spectral norm function [11],
we can see that $\lVert A(p)\rVert_{2}$ is a convex function of the
probabilities $p(\pi)$ assigned to each ordering. This will allow us to derive
a form for both the upper and lower bound on $\lVert A(p)\rVert_{2}$.

Beginning with the upper bound, we can use convexity to write the following:

$\displaystyle\sigma_{\max}\big{(}A(p)\big{)}=\sigma_{\max}\Big{(}\sum_{\pi\in\Pi}p(\pi)\cdot
A^{(\pi)}\Big{)}\leq\sum_{\pi\in\Pi}p(\pi)\cdot\sigma_{\max}\big{(}A^{(\pi)}\big{)}\leq\max_{\pi\in\Pi}\Big{\\{}\sigma_{\max}\big{(}A^{(\pi)}\big{)}\Big{\\}}.$

Finding the upper bound therefore reduces to finding the basis element with
the largest spectral norm. For an arbitrary ordering $\pi$, we can derive the
spectral norm via the matrix $A^{(\pi)}{A^{(\pi)}}^{\top}$, whose entries are
given by:

$\displaystyle\Big{(}A^{(\pi)}{A^{(\pi)}}^{\top}\Big{)}_{ij}=\begin{cases}2&\text{if}\;i=j\\\
-1&\text{if}\;i\neq j\;\text{and}\;\text{Pre}(\pi,i,j)\\\ -1&\text{if}\;i\neq
j\;\text{and}\;\text{Pre}(\pi,j,i)\\\ 0&\text{otherwise}.\end{cases}$

Within these matrices, the diagonal entries are always equal to $2$, most rows
contain two $-1$’s, and two rows (corresponding to the first and last element
in the ordering $\pi$) contain a single $-1$. We can see that the matrices
$A^{(\pi)}{A^{(\pi)}}^{\top}$ are identical up to a permutation of the rows
and columns, or that they are similar matrices. They therefore share the same
eigenvalues.

Without loss of generality, we choose to focus on the ordering $\pi^{\prime}$
that places the indices in their original order, or $\pi^{\prime}(i)=i$ for
$i\in[d]$. The corresponding matrix
$A^{(\pi^{\prime})}{A^{(\pi^{\prime})}}^{\top}$ has the following form:

$\displaystyle
A^{(\pi^{\prime})}{A^{(\pi^{\prime})}}^{\top}=\begin{pmatrix}2&-1&0&0&\cdots&0\\\
-1&2&-1&0&\cdots&0\\\ 0&-1&2&-1&\cdots&0\\\ 0&0&\ddots&\ddots&\ddots&\vdots\\\
\vdots&\vdots&\ddots&-1&2&-1\\\ 0&\cdots&\cdots&0&-1&2\end{pmatrix}$

This is a tridiagonal Toeplitz matrix, a type of matrix whose eigenvalues have
a known closed-form solution [49]. The maximum eigenvalue is the following:

$\displaystyle\lambda_{\max}\Big{(}A^{(\pi^{\prime})}{A^{(\pi^{\prime})}}^{\top}\Big{)}$
$\displaystyle=2+2\sqrt{(-1)\cdot(-1)}\cos\left(\frac{\pi}{d+1}\right)$
$\displaystyle=2+2\cdot\cos\left(\frac{\pi}{d+1}\right).$

We therefore have the following upper bound on the spectral norm for a random-
order value:

$\displaystyle\lVert
A\rVert_{2}\leq\sqrt{2+2\cdot\cos\left(\frac{\pi}{d+1}\right)}.$

This value tends asymptotically towards 2, and it is achieved by any random-
order value that places all its probability mass on a single ordering. Such a
single-order value approach has not been discussed in the literature, but it
is a special case of a proposed variant of the Shapley value [26].

Next, we consider the lower bound. We now exploit the convexity of $\lVert
A(p)\rVert_{2}$ to argue that a specific distribution achieves a lower
spectral norm than any other distribution.

Consider an arbitrary distribution $p(\pi)$ over orderings, with the
associated linear operator $A(p)$ and spectral norm $\lVert A(p)\rVert_{2}$.
Our proof technique is to find other distributions that result in different
linear operators $A$ that share the same spectral norm. For any ordering
$\pi^{\prime}\in\Pi$, we can imagine permuting the distribution $p(\pi)$ in
the following sense: we generate a new distribution $p_{\pi^{\prime}}(\pi)$
such that $p_{\pi^{\prime}}(\pi)=p(\pi(\pi^{\prime}))$, where
$\pi(\pi^{\prime})$ is a modified ordering with
$\pi(\pi^{\prime})(i)=\pi(\pi^{\prime}(i))$. That is, the probability mass is
reassigned based on a permutation of the indices.

Crucially, the linear operator $A(p_{\pi^{\prime}})$ is different from $A(p)$
but shares the same spectral norm. This is because
$A(p_{\pi^{\prime}})A(p_{\pi^{\prime}})^{\top}$ is similar to
$A(p)A(p)^{\top}$, or they are equal up to a permutation of the rows and
columns. We therefore have $\lVert A(p_{\pi^{\prime}})\rVert_{2}=\lVert
A(p)\rVert_{2}$, and due to convexity we have

$\displaystyle\big{\lVert}A\big{(}p/2+p_{\pi^{\prime}}/2\big{)}\big{\rVert}_{2}\leq\frac{1}{2}\lVert
A(p_{\pi^{\prime}})\rVert_{2}+\frac{1}{2}\lVert A(p)\rVert_{2}=\lVert
A(p)\rVert_{2}.$

This does not imply that $p/2+p_{\pi^{\prime}}/2$ is a minimizer, but we can
now generalize this technique to identify a lower bound. Rather than
generating a single permuted distribution, we propose generating all possible
permuted distributions. That is, we consider the set of all $d!$ permutations
$\pi^{\prime}\in\Pi$ and generate the corresponding distributions
$\\{p_{\pi^{\prime}}:\pi^{\prime}\in\Pi\\}$. Notice that when we take the mean
of these distributions, the result is identical regardless of the original
distributions $p(\pi)$, and that it assigns uniform probability mass to all
orderings, because we have:

$\displaystyle\frac{1}{d!}\sum_{\pi^{\prime}\in\Pi}p_{\pi^{\prime}}(\pi)=\frac{1}{d!}\sum_{\pi^{\prime}\in\Pi}p(\pi^{\prime})=\frac{1}{d!}.$

We can then exploit convexity to write the following inequality:

$\displaystyle\Big{\lVert}A\big{(}\sum_{\pi^{\prime}\in\Pi}p_{\pi^{\prime}}/d!\big{)}\Big{\rVert}_{2}\leq\lVert
A(p)\rVert_{2}.$

We can make this argument for any original distribution $p$, so this provides
a global lower bound. The random-order value with uniform weighting is the
Shapley value, whose spectral norm was derived in Lemma 7. We therefore
conclude with the following lower bound for random-order values, which is
achieved by the Shapley value:

$\displaystyle\lVert A\rVert_{2}\geq\sqrt{\frac{2}{d}}.$

∎

######  Lemma 11.

When the linear operator $A$ satisfies the boundedness, symmetry and
efficiency properties, the spectral norm is $\lVert A\rVert_{2}=\sqrt{2/d}$.

###### Proof.

The result follows from the spectral norm result in Lemma 7, and the fact that
the Shapley value is the only solution concept to satisfy linearity,
boundedness, symmetry and efficiency [46]. ∎

##  Appendix E Proofs: main results

We now re-state and prove our results from Section 4.

######  Theorem 1.

The robustness of removal-based explanations to input perturbations is given
by the following meta-formula,

$\lVert\phi(f,x)-\phi(f,x^{\prime})\rVert_{2}\leq g(\textnormal{removal})\cdot
h(\textnormal{summary})\cdot\lVert x-x^{\prime}\rVert_{2},$

where the factors for each method are defined as follows:

$\displaystyle g(\textnormal{removal})$ $\displaystyle=\begin{cases}L&\text{if
{removal} $=$ {baseline} or {marginal}}\\\ L+2BM&\text{if {removal} $=$
{conditional}},\end{cases}$ $\displaystyle h(\textnormal{summary})$
$\displaystyle=\begin{cases}2\sqrt{d}&\text{if {summary} $=$ {Shapley},
{Banzhaf}, or {leave-one-out}}\\\ \sqrt{d}&\text{if {summary} $=$ {mean when
included}}.\end{cases}$

###### Proof.

Given two inputs $x$ and $x^{\prime}$, the results in Lemma 1 and Lemma 3
imply the following bound for any feature set ${\bm{x}}_{S}$,

$\displaystyle|f(x_{S}^{\prime})-f(x_{S}^{\prime})|\leq L^{\prime}\cdot\lVert
x_{S}-x_{S}^{\prime}\rVert_{2},$

where we have $L^{\prime}=L$ when using the baseline or marginal approach and
$L^{\prime}=L+2BM$ when using the conditional approach. This bound depends on
the specific subset, we also have a bound across all subsets because $\lVert
x_{S}-x_{S}^{\prime}\rVert_{2}\leq\lVert x-x^{\prime}\rVert_{2}$ for all
$S\subseteq[d]$. In terms of the corresponding vectors
$v,v^{\prime}\in\mathbb{R}^{2^{d}}$, this implies the following bound:

$\displaystyle\lVert
v-v^{\prime}\rVert_{\infty}=\max_{S\subseteq[d]}|v_{S}-v^{\prime}_{S}|\leq
L^{\prime}\cdot\lVert x-x^{\prime}\rVert_{2}.$

Next, Lemma 7 shows that we have the following bound on the difference in
attributions:

$\displaystyle\lVert\phi(f,x)-\phi(f,x^{\prime})\rVert_{2}=\lVert Av-
Av^{\prime}\rVert_{2}\leq\lVert A\rVert_{1,\infty}\cdot\lVert
v-v^{\prime}\rVert_{\infty},$

where the norm $\lVert A\rVert_{1,\infty}$ depends on the chosen summary
technique. Combining these bounds, we arrive at the following result:

$\displaystyle\lVert\phi(f,x)-\phi(f,x^{\prime})\rVert_{2}\leq\lVert
A\rVert_{1,\infty}\cdot L^{\prime}\cdot\lVert x-x^{\prime}\rVert_{2}.$

Substituting in the specific values for $L^{\prime}$ and $\lVert
A\rVert_{1,\infty}$ completes the proof.

∎

######  Theorem 2.

The robustness of removal-based explanations to model perturbations is given
by the following meta-formula,

$\lVert\phi(f,x)-\phi(f^{\prime},x)\rVert_{2}\leq
h(\textnormal{summary})\cdot\lVert f-f^{\prime}\rVert,$

where the functional distance and factor associated with the summary technique
are specified as follows:

$\displaystyle\lVert f-f^{\prime}\rVert$ $\displaystyle=\begin{cases}\lVert
f-f^{\prime}\rVert_{\infty}&\text{if {removal} $=$ {baseline} or
{marginal}}\\\ \lVert f-f^{\prime}\rVert_{\infty,{\mathcal{X}}}&\text{if
{removal} $=$ {conditional}},\end{cases}$ $\displaystyle
h(\textnormal{summary})$ $\displaystyle=\begin{cases}2\sqrt{d}&\text{if
{summary} $=$ {Shapley}, {Banzhaf}, or {leave-one-out}}\\\ \sqrt{d}&\text{if
{summary} $=$ {mean when included}}.\end{cases}$

###### Proof.

Given two models $f$ and $f^{\prime}$, the results in Lemmas 5 and 6 imply the
following bound for any feature set ${\bm{x}}_{S}$,

$\displaystyle|f(x_{S})-f^{\prime}(x_{S})|\leq\lVert f-f^{\prime}\rVert,$

where the specific norm is $\lVert f-f^{\prime}\rVert_{\infty}$ when using the
baseline or marginal approach and $\lVert
f-f^{\prime}\rVert_{\infty,{\mathcal{X}}}$ when using the conditional
approach. In terms of the corresponding vectors
$v,v^{\prime}\in\mathbb{R}^{2^{d}}$, this implies the following bound:

$\displaystyle\lVert
v-v^{\prime}\rVert_{\infty}=\max_{S\subseteq[d]}|v_{S}-v^{\prime}_{S}|\leq\lVert
f-f^{\prime}\rVert.$

Next, Lemma 7 showed that we have the following bound on the difference in
attributions:

$\displaystyle\lVert\phi(f,x)-\phi(f^{\prime},x)\rVert_{2}=\lVert Av-
Av^{\prime}\rVert_{2}\leq\lVert A\rVert_{1,\infty}\cdot\lVert
v-v^{\prime}\rVert_{\infty},$

where the norm $\lVert A\rVert_{1,\infty}$ depends on the chosen summary
technique. Combining these bounds, we arrive at the following result:

$\displaystyle\lVert\phi(f,x)-\phi(f^{\prime},x)\rVert_{2}\leq\lVert
A\rVert_{1,\infty}\cdot\lVert f-f^{\prime}\rVert.$

Substituting in the specific values for $\lVert f-f\rVert$ and $\lVert
A\rVert_{1,\infty}$ completes the proof.

∎

Next, we present a corollary that considers the case of simultaneous input and
model perturbations. While this case has received less attention in the
literature, we include this result for completeness because it is easy to
analyze given our previous analysis.

######  Corollary 1.

The robustness of removal-based explanations to simultaneous input and model
perturbations is given by the following meta-formula,

$\lVert\phi(f,x)-\phi(f^{\prime},x^{\prime})\rVert_{2}\leq
g(\textnormal{removal})\cdot h(\textnormal{summary})\cdot\lVert
x-x^{\prime}\rVert_{2}+h(\textnormal{summary})\cdot\lVert f-f^{\prime}\rVert,$

where the functional distance and factors associated with the removal and
summary technique are given in Theorems 1 and 2.

###### Proof.

This result follows from triangle inequality:

$\displaystyle\lVert\phi(f,x)-\phi(f^{\prime},x^{\prime})\rVert_{2}\leq\lVert\phi(f,x)-\phi(f,x^{\prime})\rVert_{2}+\lVert\phi(f,x^{\prime})-\phi(f^{\prime},x^{\prime})\rVert_{2}.$

The bounds for the first and second terms are given by Theorem 1 and Theorem
2, respectively. ∎

##  Appendix F The role of sampling when removing features

In this section, we present results analyzing the impact of sampling when
removing features. Sampling is often required in practice to integrate across
the distribution $q({\bm{x}}_{\bar{S}})$ for removed features (see eq. 1), and
this can create further discrepancy between attributions for similar inputs.

###  F.1 Theory

Here, we provide a probabilistic bound for the difference between the model
prediction $f(x_{S})$ with a subset of features $x_{S}$ and its estimator
$\hat{f}(x_{S})$ with $m$ independent samples from $q({\bm{x}}_{\bar{S}})$.

######  Lemma 12.

Let $\hat{f}(x_{S})\equiv\frac{1}{m}\sum_{i=1}^{m}f(x_{S},x_{\bar{S}}^{(i)})$
be the empirical estimator of $f(x_{S})$ in eq. 1, where
$x_{\bar{S}}^{(i)}\overset{\text{i.i.d.}}{\sim}q({\bm{x}}_{\bar{S}})$. For any
$\delta\in(0,1)$ and $\varepsilon>0$, if the sample size $m$ satisfies
$m\geq\frac{2B^{2}\log(2/\delta)}{\varepsilon^{2}}$, then
$\mathbb{P}\left(|\hat{f}(x_{S})-f(x_{S})|\leq\varepsilon\right)\geq
1-\delta$.

###### Proof.

With Assumption 2, we have
$-\frac{B}{m}\leq\frac{f(x_{S},x_{\bar{S}}^{(i)})}{m}\leq\frac{B}{m}$. Also,
$\mathbb{E}[\hat{f}(x_{\bar{S}})]=f(x_{\bar{S}})$ because $x_{\bar{S}}^{(i)}$
are sampled from $q({\bm{x}}_{\bar{S}})$. Applying Hoeffding’s inequality
[35], we obtain the probabilistic bound

$\displaystyle\mathbb{P}(|\hat{f}(x_{S})-f(x_{S})|\geq\varepsilon)$
$\displaystyle\leq
2\exp\left(-\frac{2\varepsilon^{2}}{\sum_{i=1}^{m}\left(\frac{2B}{m}\right)^{2}}\right)$
$\displaystyle=2\exp\left(-\frac{m\varepsilon^{2}}{2B^{2}}\right)\leq
2\exp\left(-\frac{\left(\frac{2B^{2}\log(2/\delta)}{\varepsilon^{2}}\right)\varepsilon^{2}}{2B^{2}}\right)=\delta,$

and hence for the complement event we have

$\mathbb{P}(|\hat{f}(x_{S})-f(x_{S})|\leq\varepsilon)\geq 1-\delta.$

∎

Lemma 12 can be applied to bound the difference between the original and
perturbed attributions when the integral in eq. 1 is estimated via sampling.
For example, under input perturbations, the difference
$|\hat{f}(x_{S})-\hat{f}(x_{S}^{\prime})|$ decomposes into
$|\hat{f}(x_{S})-f(x_{S})|+|f(x_{S})-f(x_{S}^{\prime})|+|f(x_{S}^{\prime})-\hat{f}(x_{S}^{\prime})|$
by the triangle inequality. The first and last terms can be bounded with high
probability using Lemma 12 and the union bound, and the middle term is bounded
using results from Section 3.1 (e.g., Lemma 1 or Lemma 3). Intuitively, by
combining probabilistic bounds over all the feature subsets, it is possible to
derive a final probabilistic upper bound with the following meta-formula:

$\lVert\phi(f,x)-\phi(f,x^{\prime})\rVert_{2}\leq g(\text{removal})\cdot
h(\text{summary})\cdot\lVert
x-x^{\prime}\rVert_{2}+g_{\text{sampling}}(\text{removal}),$ (F.1.1)

where $g_{\text{sampling}}(\text{removal})$ is a function of
$\delta,\varepsilon,m$ when removal $=$ marginal or conditional, and
$g_{\text{sampling}}(\textbf{baseline})=0$ because no sampling is needed.
Similarly, for our bound under model perturbations, the term
$g_{\text{sampling}}(\text{removal})$ can be included to account for sampling