On the Robustness of Removal-Based Feature Attributions
Authors: Chris Lin
Ian Connick Covert
Su-In Lee
Venue: NeurIPS 2023 poster (2023)
Keywords: explainable artificial intelligence
interpretable machine learning
feature attributions
removal-based feature attributions
robustness
TL;DR: Theoretical characterization of the robustness of removal-based feature attributions with empirical implications.

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.

On the Robustness of Removal-Based Feature Attributions
Chris Lin" University of Washington clin256cs, washington. edu
lan Covert Stanford University
icovert@stanford edu
Su-In Lee
University of Washington suinleeêcs washington. edu
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 these methods by showing that they 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 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 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.
Introduction
In recent years, machine learning has shown great promise in 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 interpretable explainable machine learning [45] 51
60 67). One popular way of explaining machine learning model is feature attribution, which assigns importance scores to 56, 60). Many feature attribution methods can be categorized as either gradient- -baiiedt-baeddd methods that compute gradients of model predictions with respect to input features (6), or eemoval-based methods that remove features to quantify each feature's '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 unclear robustness properties: feature attributions are vulnerable to adversarial attacks, and even without an adversary appear unstable under small changes to the input model.
53 56
input features of the model [45] 51
[33
As a motivating example, consider the work of Ghorbani et al. 28 which shows that minor perturbations an image can lead substantially different feature attributions while preserving the original prediction. Similar to adversarial examples designed to alter model predictions (29 43 61], 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
Equal contribution. +Work done while at the University of Washington.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

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- 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 22 [8] 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- gradient-based attributions 4 [1
22 65]. 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. [0). but we generalize this work by considering alternative techniques for both feature removal and importance summarization (e.g., LIME's '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-base 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 220 58]; however, Frye et al. [27] show that such hidden features can be discovered if features are removed using their conditional distribution, suggesting 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 to develop a comprehensive characterization of the robustness properties of removal-based feature attributions, focusing on two notions of robustness that have received attention 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 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 -asedd 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 1R. + IR be model whose predictions we seek to explain. We consider that the input variable
x consists of separate features, or x = (II,.... Zz). Given an index set S C [d] = {I,... ......dd) d), we detine the corresponding feature subset as æs {æ, i E 5). The power set of [d] is denoted by P(d), and we let S = [d] \ S denote the set complement. We consider the model's output space to be dinenimnnsional, 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 æ, Dg to denote random variables, the symbols x, gg to denote specific values, and p(x) to represent the data distribution with support on X' C Rd. We assume that all explanations are generated for inputs such that X.

---

2.2 Removal-based feature attributions
Our work focuses on a class of methods known as removal-based explanations 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
Intuitively,
@ This perspective shows a close connection between methods like leave-one-out [67]. LIME [53] and Shapley values [45], and enables 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(zs). Many implementations have been discussed in the literature [17]. but we focus here on three common choices. Given set of observed values z,, these techniques can all be viewed as averaging the prediction over a distribution ((zz) for the held-out feature values:
f(xs) := (g(ss) (f(ss æ5)] - f ((ss (()(3)d)
Specifically, the three choices we consider are:
(Baseline values) Given baseline input b Rd. we set the held-out features to their corresponding values bg [59] This is equivalent to letting ((#5) be a Dirac delta centered at bg.
(Marginal distribution) Rather than using a single replacement value, we can average across values sampled from the input's marginal distribution, let ((xg) = p(=s) [45)
(Conditional distribution) Finally, on the available features, which equivalent to setting g( = p(æg xs) 22
average across replacement values while conditioning
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 The conditional distribution approach requires more complex and error-prone estimation procedures (see 12] for discussion), but assume here that all three versions of f(xs) can be calculated exactly.
Summary technique. Given feature removal technique that allows us to query the model with arbitrary feature sets, we must define convention for summarizing each feature's 's influence. This challenging due the exponential number of feature sets, or because [P(d)]
24. Again, many techniques have been discussed in the literature [17]. with the simplest option comparing the prediction with all features included and with single feature missing 67.. We take a broad linear combination perspective here, considering range of approaches that define attributions of predictions with different feature sets. These methods include leave-one-out [67], RISE (51).
LIME [53] Shapley values 45), and Banzhaf values 113, among other possible options.
All of these methods yield per-feature attribution scores o. x Shapley value calculates feature attribution scores as follows [55]
for i [d]. For example, the
Σ d
(1)))
-1 (f(#su())- -
),(f,))=
SC(dN()
The specific linear combinations for the other approaches we consider are shown in Table The remainder of the paper uses the notation f(as) to refer to predictions with partial information, and o(f,z) x) = |d (f) ......... da(f, x]] E Rd 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
'The framework also considers the choice to explain different model behaviors (e.g.. the loss calculated over examples) [17]. but we focus on methods that explain individual predictions.
3

---

1. whether Il)(f.x) o(f.2'))l is controlled by || z'|l (input perturbation), and
2. whether [l)(f.)) (( x>] is controlled by |I - f'll (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.
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 input and model perturbations, and we then discuss robustness properties of the summary technique. Proofs for all results are in the Appendix.
Before proceeding, introduce two assumptions that necessary for our analysis.
Assumption 1. We assume that the model I is globally L-Lipschitz continuous, or that we have [fx)) [(<)) I
s L 1|x-- x' 12 for all x. x' E R.
Assumption 2. We assume that the model f has bounded predictions, or that there exists a constant B such that [x)] & B for all x R
The first assumption holds for most deep learning architectures, because they typically compose 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 norm for function R2 + IR is defined as llgl|> = (J [ 1/p, where the integral is taken over Rd_ |[ denotes the same integral taken over the domain X' _ TR..
Our results adopt this to define notion of distance between functions ||g g' llp- focusing the cases where for the L distance and oc for the Chebyshev distance.
3.1 Prediction robustness to input perturbations
Our goal here is to understand how the prediction function f(xs) for feature xy behaves under small input perturbations. We first consider the case where features are removed using the baseline or marginal approach, and 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 (xx)) for any feature set xs I.-Lipsciiz continuous:
(f(xs) f(-'s) L ||rs 2'sll2
V Vag, E R!
Next, we consider the case where features are removed using the conditional distribution approach. We show that the continuity of f(xs) 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(s) for a feature set Øg satisfies
(f(s) f(4's)] < L ||.s d'sll2 + 2B drv (p(xs xs), plas w''s).
where the total variation distance is defined via the L functional distance as
drv (p(ats xs). a's)) = p(as xs) - p(as m's))l;
Lemma does not immediately imply Lipschitz continuity for f(xs), because we have not bounded the total variation distance in terms of ||اş - 22 To address this, we require the following Lipschitz-like continuity property in the total variation distance.
A Assumption 3. We assume that there exists constant M such that for all S C [d], we have
drv (p(os) vss),(ar, v's)) S M ||æg x'sll2
Varga's E S1

---

Intuitively, this says that the conditional distribution p(23 æg) cannot change too quickly as function of the observed features. This property does not hold for all data distributions p(z). and may hold in some cases only for large M; in such scenarios, the function f(xs) is not guaranteed to change slowly. However, there are cases where it holds. For example, we have M 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 a ~ (N,,E with mean μa E Fd and covariance E Raxd, Assumption holds with M = (in/uuxx--)) Amin (Σ (E-1). If x: is assumed to be standardized, this captures the case where independent features yield M
means that if one dimension is roughly a linear combination of the others, or if Inluitively, it also is large, the conditional distribution can change quickly as a function of the conditioning variables.
Now, assuming that this property holds for p(z), we show that we can improve upon Lemmal 2
Lemma 3. Under Assumption the prediction function (xg) defined using the conditional ap- proach is Lipschitz continuous with constant L + M for any feature set zş.-
Between Lemmas I and 3 we have established that the function f(rs) remains Lipschitz continuous for any feature set Ty, although in some cases with larger constant that depends on the data distribution. These results are summarized in Table 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 I:
Lipschitz constants induced by each feature removal technique.
Table 2:
Relevant functional distances for each feature removal technique.
Baseline
Marginal
Conditional
Baseline
Marginal
Conditional
L
L+ 2BM
|| foo II - I'loo
|| - [/--f'llaa f'lbo,xx
3.2 Prediction robustness to model perturbations
Our next goal is to understand how the function (=rs) for a feature set æ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 make similar predictions with a subset of features. We first derive general result involving the proximity between two models f and f' within subdomain.
Lemma 4. For two models f. RRd + I and a subdomain X' C Rdd, the finntiooo f(xs), (as) for any feature set as satisfy
(f(#s) - (x(x)) S || - I' Qza(X') 2B (1-Qa((X'))
where Qz5 (X') is the probability of imputed samples lying in X" based on the distribution q(æ5):
Qay (X') = (By((s) [T {(05, gg)) x²))
The difference in predictions therefore depends the distance between the models only within the subdomain X", as well as the likelihood of imputed sampled lying in this subdomain. This result illustrates point shown by Slack al. [58): that two models which equivalent subdomain, or even on the entire data manifold X, can lead to different attributions if we use a removal technique that yields low value for (2'). In fact, when (2') 0, Lemma 4 reduces to a trivial bound [[(s) - s'(xs) < 2B that follows directly from Assumption 2
small
The general result in Lemma 4 allows subdomain X' and removal technique g(xg) so that Qzs (X') 1.
show two simpler ones. In both cases, choose the
Lemma 5. When removing features using the conditional approach. the prediction functions f(*s), I' (xs) for two models J and any feature set as satisfy
If(xs) -
(fae))f((xss) f((x)) < - flloo.
The next result is similar but involves a potentially larger upper bound |lf f' llf - ooooo. Lemma 6. When removing features using the baseline or marginal approach, the frndictionfuuctiinn ff fæ (xs) for two models f. and any feature set x satisfy
1f(xs) - fl (zs) s |l f
5

---

Table 3: Summary technique linear operators used by various removal- based explanations.
Summary
Method
A,s (i E S)
Ass (i S)
||
||A||2
Leave-one-out
Occlusion
67
I [[|]] = d)
-I -[[[ = d -d-1
2
vd +
Shapley value
SHAP 45
(1S) (9(-1)1 (d--1s
[8 1
2v
/2/d 1/24/2-1
Banzhaf value
Banzhaf 13
1/2"-1 1/2d-1
- -1/2"
2vd
Mean when included
RISE [1
vd
V(d+10 + 1)/2d +1/2d
Weighted least squares
LIME 53
Depends on implementation choices (see Appendix D)
The relevant functional distances for Lemma 5 and Lemma 6 are summarized in Table
1mpliess implies that if two models f" are functionally equivalent, or f(z) = f(x') for all x E IRd, they remain equivalent with any feature subset T.s- In Section we will see that this implies equal attributions for the two models - natural property that, perhaps surprisingly, is not satisfied by all feature attribution methods [77 (5), and that has been described in the literature as implementation invariance
Remark 1. Lemma 6
(60) This property holds automatically for all removal-based methods.
In contrast, Lemma implies the same for models that are equivalent only on the data manifold X C Ie 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 notion of weak implementation invariance and strong implementation invariance, depending on whether equal attributions are guaranteed when we have f( = (x) everywhere (x I2 or almost everywhere (z E X)).
Remark 2.
3.3 Summary technique robustness
We previously focused on the effects the feature removal choice, so our goal is now to understand the role of the summary technique. Given removal approach that lets query the model with any feature set, the summary generates attribution scores off, æ) E Rd, often using a linear combination of the outputs (*x) for each S [d]. We formalize this in the following proposition.
Proposition 1. The attributions for each method can be calculated by applying linear operator A E Rdx2* a vector I representing the predictions with each feature set,
o(fx)) Av,
where the linear operator for each method is listed in Table 3 and v is defined as us f(xs)) for each S [d] based on the chosen feature removal technique.
In this notation, the entries of v are indexed as us for S [d]. and the entries of A are indexed as Ais for E [d) and S C [d]. The operation Au can be understood as summing across all subsets, or (Av); = Esc(a Ais ug. 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 difficult to characterize analytically, but its default parameters for tabular, image and text data approach limiting cases where LIME reduces 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 E IR"", and we must consider how the attributions differ between the original vector and perturbed version v'. Proposition| suggests that we can bound the attribtuion difference via the change in model outputs || -v'll 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 altributions given the same summary technique A and different model outputs v, v' can be bounded as
||Ae Au'll2 S ||A||2 1|u0 o' or ||A Au'll2 S
|| 11,00 1|u - of Ilac
where ||A||2 is the spectral norm, and the operator norm squared row I I-noms,, with values for each A given in Table 3
||l 1,50 the square root of the sum of

---

For both inequalities, smaller norms ||.A]| represent stronger robustness to perturbations. Neither bound is necessarily tighter, and the choice of which to use depends on how v' bounded. The first bound using the Euclidean distance ||v v'l 2 is perhaps more intuitive, but the second bound is more useful here because the results in Sections 3.1 and and 3.2 effectively relate ||0-~'||
This view of the summary technique's robustness raises a natural question, which is whether the approaches used in practice have relatively good bad robustness (Table B) The answer is not immediately clear, because within the class of linear operators we can control the robustness arbitrarily: we can achieve || = >
by setting A 0 (strong robustness), and we can likewise get ||A||2 || - oc by setting entries of A to 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 E Rdx2", 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 the following properties, which are motivated by axioms from the literature on game-theoretic credit allocation 46 54
Definition 2. We define the following properties for linear operators removal-basee based explanations:
Rxx2 that are used for
(Boundedness) For e E IR=", the attributions are bounded by each feature's 's smallest and largest contributions, or minszi(Vsu((s us) S (Av); S maxszi(Vsutts us) for all E [d].
(Symmetry) For all TR?",, two features that make equal marginal contributions all feature sets must have equal attributions, or (Av); = (Av); if Usu(() = Usu(j) for all S C [d].
(Efficiency) For all E B²²", the attributions sum to the difference in predictions for the complete and empty sets, or 1 Av [[4] VO
Among these properties, prioritize boundedness because this leads to class of solutions that are well-known in game theory, and which are called probubilistic values [46]. 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 || which is the simpler case.
Lemma 8. When the linear operator A satisfies the boundedness property, we have ||A||1, 2vd.
The above result applies to several summary techniques shown in Table 3 including leave-one- [67]. Shapley values [45] and Banzhaf values [13. We now address the case of the spectral norm ||A||e. beginning with the case where boundedness and symmetry are satisfied.
Lemma 9.
When the linear operator satisfies the boundedness and symmetry properties, the spectral norm is bounded follows,
24/2-1 S1A Vd+1,
with 1/2d/2-1 achieved by the Banzhaf value and /d + achieved by leave-one-out.
This shows that the Banzhaf value is the most robust summary within 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 criticized for failing to satisfy the efficiency property [25), so we now address the case where this property holds.
When the linear operator satisfies the boundedness and efficiency properties, the spectral norm is bounded follows.
Lemma 10.
V2 S 1All < 2+ 2-+2.00 cos
d+1
with 2 achieved by the Shapley value.
We therefore observe that the Shapley value the most robust choice within a different class of linear operators, which are known random-orde values [46). In comparison, the worst-case robustness
7

---

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 T0 approximately equal to 2, whereas the upper bound in Lemma 9 grows with the input dimensionality. suggesting that the efficiency property imposes minimum level of robustness but limits robustness in the best case.
Finally, we can trivially say that ||Al|2 v 2/d when boundedness, symmetry and efficiency all hold, because the Shapley value is the only linear operator to satisfy all three properties [46 [4]
Lemma 11. When the linear operator satisfies the boundedness, symmetry and efficiency proper- ties, the spectral norm is ||Al||2 /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 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 the robustness input perturbations.
Theorem 1. The robustness of removal- based explanations to input perturbations is given by the following fortaaforuuaa.
[ "ff x1)l2 Σ g(removal) h(summary) || - x' al||2,
where the factors for each method are defined follows:
[L
if removal = baseline or marginal
g(removal)
I. 2BM if removal = conditional.
(2Vd
if summary Shapley, Banzhaf, or leave-one-out if summary mean when included.
h(summary)
Vd
We therefore observe that the explanation functions themselves Lipschitz continuous, but with constant that combines properties of the original model with implementation choices of the attribution method. The model's 's' inherent robustness interacts with the removal choice to yield the grremoval) factor, and the summary technique is accounted for by l (summaryy.
Our second result takes a similar form, but relates to the robustness to model perturbations.
The robustness of removal-bate based explanations to model perturbations is given by the following meta-formula.
Theorem 2.
1) -
S h(summary) ||ʃ f'll:
where the functional distance and factor associated with the summary technique are specified as follows:
if removal baseline marginal 1ll f'loo removal conditional,
[||5 - 1ox
Iff f'l
(2 if summaar Shapley, Banzhaf, or leave-one-out Vd summmay mean when included.
h(summary) =
Similarly, we see here that the attribution difference is determined by the strength of the model perturbation, as measured by 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 model perturbation, but we can also account for simultaneous input and model perturbations, as shown in Corollary ii Appendix F The remainder of the paper empirically verifies these results and discusses their practical implications.
5 Summary of experiments
Due to space constraints, we defer our experiments to Appendix but briefly describe the results here. First, empirically validate our results using synthetic data and a logistic regression classifier, where

---

||x-x||a 0.1
||x=x||a 0.5
||x-x||a
||- x 2.0
[|00,x0 stf.x2lls
0.920
0.05
9.150
0.25
Margial
Beseline
Shapley
0.125
0.20
Condit onal
0.015
0.06
0.303
0.55
0.04
0,815
0.000
0.02
5.000
0.10
dea nio
o sea n:10
šoz
a sos o
loz ana
Weight Decay
Weight Decay
Weight Decay
Weight Decay
(b)
|
|K- ==
|x-x|]== =10
| 20
[MCCx) elf xtllss
Baseline Baseline iZerol
0.10
Shapley
0.08
0.08-
0.000 6 507 6 DOE 0.06 0:000 9-515 0:000 6 507 6 DOI 0016 0:000 9:516 0.000 6 300 6 FOT 0.066 0.000 0:510 0.000 O 300 S S01 0.066 0.000 0:510 Weight Decay
Weight Decay
Weight Decay
Wright Decay
Figure I: 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% pairs.
confidence intervals across
Baseline
Marginal
Conditional
1.0
10
1.0
Peatson Correlation
0.6
0.6-
0.6
0.4
0.4
0.2
0.2
0.0
0.0
0.0
1.0
1.0
1.0
Rank Carrelation
0.8
0.8
2.
0.6
0.6
0.4
0.4
02
0.
0.2
0.0
0.0
Origiall
Unsar
Unear
Unsarz
Uneer
Original
Limerr
Linee
Gmeer
Unnar
Originl
Unew
Uneer?
Dcclusion
Shapley
Banzhaf
RSE
LIME
Vanilla Grad
IntGrad
Grad x Input
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.
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 [].
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 marginal removal, implied by Theorem 2
Next, we demonstrate two practical implications of our findings with the UCI wine quality 14 24
MNIST 44) CIFAR-10 10 42 and Imagenette datasets [99 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 's Lipschitz constant is upper bounded by the product of its layers' spectral norms [6], which are respectively upper bounded by each layer's Frobenius norm. These results are shown for Shapley value attributions in Figure [] for the wine quality and MNIST datasets, where we see a consistent effect of weight decay across feature removal approaches. We speculate that similar effect may be observed with other techniques that encourage

---

Lipschitz regularity 8 30 63 (2) Viewing the sanity check with parameter randomization proposed by Adebayo et al. 2 as form of model perturbation, we show that removal-based attributions empirically pass the sanity check, in contrast to certain gradient-based methods 155 60
Figure 2
shows results for the wine quality dataset, where the similarity for eemoval-based methods decays more rapidly than for IntGrad and Grad x Input. Results for the remaining datasets are in Appendix
6
Discussion
Our analysis focused on the robustness properties of removal-based feature attributions to perturba- tions in both the model f and input æ. For input perturbations, Theorem [ shows that the attribution change is proportional to the imput perturbation strength. The Lipchitz constant determined by the attribution's '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 B 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 's Lipschitz regularity [ 30 63
For model perturbations, Theorem shows that the attribution difference is proportional 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 B1 37 48), which are in some cases mathematically equivalent to attributions of an ensembled model.
Our robustness results provide a new lens to compare different removal-basee -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, lcave-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 nume-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 [7 27], 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-baas -ased 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) [2] Practically, both are computationally intensive to determine for realistic networks 33 63) 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 neural network's '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 æ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., Lp norms for * 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 may achievable with specific perturbation strengths, unlike our approach that provides bounds simultancously for any perturbation strength via global continuity properties.
[8 39
52
57 62 tighter bounds
10

---

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-1759487J; and the National Institutes of Health [R35 GM 128638 and RO1 NIA AG 061132].
References
[1] 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.
[2] 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.
[3] 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.
[4] Sushant Agarwal, Shahin Jabbari, Chirag Agarwal, Sohini Upadhyay, Steven Wu, and Himabindu Lakkaraju. Towards the unification and robustness of perturbation and gradi- based explanations. In International Conference Machine Learning, pages 110-119. PMLR, 2021.
[5] David Mlvarez-Melis and Tommi S Jaakkola. On the robustness of interpretability methods. arXiv preprini arXiv: 1806 886.08049, 2018.
[6] Marco Ancona, Enea Ceolini, Cengiz Öztireli, and Markus Gross. Towards better under- standing of gradient-base based attribution methods for deep neural networks. arXiv preprint arXiv:1711.06114, 2017.
[7] Christopher Anders, Plamen Pasliev, Ann-Kathrin Dombrowski, Klaus-Robert Müller, and Pan Kessel. Fairwashing explanations with off-manifold detergent, In International Conference Machine Learning, pages 314-323. PMLR, 2020.
[8] Alexandre Araujo, Aaron Havens, Blaise Delattre, Alexandre Allauzen, and Bin Hu. unified algebraic perspective on Lipschitz neural networks. arXiv preprint arXiv:2303.03169, 2023.
[9] John F Banzhaf III. Weighted voting doesn't work: A mathematical analysis. Rutgers Law Review, 19:317, 1964.
[10] 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.
[11] Stephen Boyd, Stephen Boyd, and Lieven Vandenberghe. Convex Optimization. Cambridge University Press, 2004.
[12] Hugh Chen, Ian C Covert, Scott M Lundberg, and Su-In Lee. Algorithms estimate Shapley value feature attributions. arXiv preprint arXiv:2207.07605, 2022.
[13] Jianbo Chen and Michael Jordan. LS- -S-Tree: Model interpretation when the data are linguistic. In Proceedings of the AAAI Conference Artificial Intelligence, volume 34, pages 3454-3461, 2020.
[14] 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-555: 2009.
[15] Ian Covert, Scott M Lundberg, and Su-In Lee. Understanding global feature contributions with additive importance measures. Advances Neural Information Processing Systems, 33: 17212-17223, 2020.

---

[16] Ian Covert, Chanwoo Kim, and Su-In Lee. Learning to estimate Shapley values with vision transformers. arXiv preprint arXiv:2206.0522 2022.
[17] Ian C Covert, Scott Lundberg, and Su-In Lee. Explaining by removing: A unified framework for model explanation. The Journal of Machine Learning Research, 2277799477-956 2021.
[18] Alex 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.
[19] 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 computer vision and pattern recognition, pages 248-255 Ieee, 2009.
[20] Botty Dimanov, Umang Bhatt, Mateja Jamnik, and Adrian Weller. You shouldn's 't trust me: Learning models which conceal unfairness from multiple explanation methods. 2020.
[21] Ann-Kathrin Dombrowski, Maximillinn Alber, Christopher Anders, Marcel Ackermann, Klaus- Robert Müller, and Pan Kessel. Explanations can be manipulated and geometry to blame. Advances in Neural Information Processing Systems, 32, 2019.
[22] 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.
[23] Jiayun Dong and Cynthia Rudin. Exploring the cloud of variable importance for the set of all good models. Nature Machine Intelligence, (11))880-8222 2020.
[24] Dheeru Dua and Casey Graff. UCI machine learning repository, 2017. Utthtpp://arccevv ics.uci. edu/ml
[25] Pradeep Dubey and Lloyd S Shapley. Mathematical properties of the Banzhaf power index. Mathematics of Operations Research, 4(2):99-131, 1979.
[26] Christopher Frye, Colin Rowat, and Ilya Feige. Asymmetric Shapley values: incorporat- ing causal knowledge into model-agnostic explainability. Advances in Neural Information Processing Systems, 33:1229-1239, 2020.
[27] 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.
[28] 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.
[29] Ian J Goodfellow, Jonathon Shlens, and Christian Szegedy. Explaining and harnessing adversar- ial examples. arXiv preprint 1rXiv:1412.6572, 2014.
[30] Henry Gouk, Fibe Frank, Bernhard Pfahringer, and Michael J Cree. Regularisation of neural networks by enforcing Lipschitz continuity. Machine Learning, 1 10:393 2021.
[31] 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 Computational Biology, pages 33-45. PMLR, 2022.
[32] Peter I. Hammer and Ron Holzman. Approximations of nseudo-boolean functions; applications to game theory. Zeitschrift für Operations Research, 36(1):3-21, 1992.
[33] 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.
[34] Juyeon Heo, Sunghwan Joo, and Taesup Moon. Fooling neural network interpretations via adversarial model manipulation. Advances in Neural Information Processing Systems, 32, 2019.
12

---

[35] Wassily Hoeffding. Probability inequalities for sums of bounded random variables. The Collected Works of Wassily Hoeffding, pages 409-426, 1994.
[36] Jeremy Howard and Sylvain Gugger. FastAl: A layered API for deep learning, Information, 11 (2):108, 2020.
[37] 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.
[38] Matt Jordan, Jonathan Hayase, Alex Dimakis, and Sewoong Oh. Zonotope domains for lagrangian neural network verification. Advances in Neural Information Processing Systems, 35:8400-84 1884000-813,,0 2022.
[39] 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 / 30, pages 97-117. Springer, 2017.
[40] Zulqarnain Khan, Aria Masoomi, Davin Hill, and Jennifer Dy. Analyzing the effects of classifier lipschitzness on explainers. arXiv preprint arXiv:2206. 12481, 2022.
[41] 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 Al: Interpreting, explaining and visualizing deep learning, pages 267-280, 2019.
[42] Alex Krizhevsky, Geoffrey Hinton, et al. Learning multiple layers of features from tiny images. 2009.
[43] Alexey Kurakin, Ian J Goodfellow, and Samy Bengio. Adversarial examples in the physical world. In Artificial Intelligence Safety and Security, pages 99-191--111111 112. Chapman and Hall/CRC. 2018.
[44] Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner. Gradient-based learning applied to document recognition. Proceedings of the IEEE, 86(1 86(11):227--2324, 1998.
[45] Scott M Lundberg and Su-In Lee. A unified approach to interpreting model predictions. Advances in Neural Information Processing Systems, 30, 2017.
[46] Dov Monderer and Dov Samet. Variations on the Shapley value. Handbook of Game Theory with Economic Applications, 3:2055-2011 2002.
[47] Grégoire Montavon, Alexander Binder, Sebastian Lapuschkin, Wojciech Samek, and Klaus- Robert Müller. Layer-wise relevance propagation: overview. Explainable Al: Interpreting. Explaining and Visualizing Deep Learning, pages 193-209, 2019.
[48] Yilin Ning, Mingxuan Liu, and Nan Liu. Shapley variable importance cloud for machine learning models. arXiv preprint arXiv:2212.083770 2022.
[49] Silvia Noschese, Lionello Pasquini, and Lothar Reichel. Tridiagonal toeplitz matrices: properties and novel applications. Numerical Linear Algebra With Applications, 20(2):302 326, 2013.
[50] 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.
[51] Vitali Petsiuk, Abir Das, and Kate Saenko. RISE: Randomized input sampling for explanation of black-box models. arXiv preprint 1 2018
[52] Aditi Raghunathan, Jacob Steinhardt, and Percy S Liang. Semidefinite relaxations for certifying robustness to adversarial examples. Advances in Neural Information Processing Systems, 31, 2018.
13

---

[53] Marco Tulio Ribeiro, Sameer Singh, and Carlos Guestrin. "Why should trust you?" Explaining the predictions of any classifier. In Proceedings of the 22nd ACM SIGKDD International Conference Knowledge Discovery and Data Mining, pages 1135-1144, 2016.
[54] Lloyd S Shapley. A value for n-person games. Contributions to the Theory of Games, 2(28):
307-317, 1953.
[55] Avanti Shrikumar, Peyton Greenside, and Anshul Kundaje. Learning important features through propagating activation differences. In International Conference Machine Learning, pages 3145-3153, PMLR, 2017.
[56] Karen Simonyan, Andrea Vedaldi, and Andrew Zisserman. Deep inside convolutional networks: Visualising image classification models and saliency maps. arXiv preprint arXiv:1312.6034, 2013.
[57] 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.
[58] 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.
[59] Mukund Sundararajan and Amir Najmi. The many Shapley values for model explanation. In International Conference on Machine Learning, pages 9269-9278. PMLR, 2020.
[60] Mukund Sundararajan, Ankur Taly, and Qiqi Yan. Axiomatic attribution for deep networks. In International Conference on Machine Learning. pages 3319 -3288. PMLR, 2017.
[61] Christian Szegedy, Wojciech Zaremba, Ilya Sutskever, Joan Bruna, Dumitru Erhan, Ian Goodfel- low, and Rob Fergus. Intriguing properties of neural networks. arXiv preprint 2013.
[62] Vincent Tjeng, Kai Xiao, and Russ Tedrake. Evaluating robustness of neural networks with mixed integer programming. arXiv preprint arXiv: 2017.
[63] Aladin Virmaux and Kevin Scaman. Lipschitz regularity of deep neural networks: analysis and efficient estimation. Advances in Neural Information Processing Systems, 31, 2018.
[64] Tianhao Wang and Ruoxi Jia. Data Banzhaf: A data valuation framework with maximal robustness to learning stochasticity. arXiv preprint arXiv:2205.15466. 2022.
[65] 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- -3331333333 2020.
[66] Yuichi Yoshida and Takeru Miyato. Spectral norm regularization for improving the generaliz- ability of deep learning. arXiv preprint arXiv: 2017.
[67] 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 13, pages 818-833. Springer, 2014.
14

---

A Experiments
This section conducts a range of experiments provide empirical validation for our theoretical results, and to examine their practical implications for common machine learning models. Code is available at httppp///ggguuu..co///uuulllllla
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 !
Setup. For each input R+ in the synthetic dataset, features generated from multivariate Gaussian distribution with zero mean, I3, 24 independent from each other and from ø1, 22, and with correlation p between 2 22. That is, the covariance matrix E E Rx44 has the following form:
3
0 0
0
0
The binary label y E {0, 1} for each input is sampled from a Bernoulli distribution with probability a (B æ), where (-) is the sigmoid function and 3² [5. 0, 3, 1]. Unless stated otherwise, we consider p = 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.
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 &g norms normalized values ranging between 0 to 2 with 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 E difference between original and perturbed attributions at each perturbation norm. The slope of such line estimates the dependency of the upper bound of attribution differences with respect the perturbation norm. Theorem guarantees specific slope, and this experiment aims to verify this property.
Furthermore, note that the sigmoid function (-) (-) 1 /4-Lippchitz because its largest derivative 1/4. logistic regression classifier composition a linear function and sigmoid function a--). so the model is L-Lipschitz with L || /4, where 3 is the trained model coefficients. Finally, because the input features follow a Gaussian distribution, the total variation constant M derived in Example applies here. Hence, we have specific values for all terms in Theorem bound, which allows direct comparison between theoretical and empirical results.
Baseline
Marginal
Conditional
Thaoediaa Upper Foand finaiies Upper Sound
Thapwotcll Unperibonn
Smolice upper sound
0.14
Shapley
0.10
0.08
0.06
0.04
0.0
0.
3.5
1.0
1.
25
1.0
1.5
2.0
3.0
0.5
25
2.5
20
||x
x/122
||x x|Iaa
x'Ila
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 p 1.
15

---

Results. As shown in Figure B 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
Moreover, we vary the Lipschitz constant L of the trained logistic regression classifier by factor of {0.05, 0.1, 0.5, 1, 2]. which is done by multiplying the trained coefficients 3 by a. Figure 4 shows that the perturbation norm-dependent slope of the empirical bound increases linearly with a, aligning with our theoretical results. Also, Example [ suggests that increasing P would lead to a larger upper bound for conditional feature removal, which is observed in the empirical bound with p = (0, 0.25, 0.5, 0.75, }} (see (
Overall, these empirical results validate our theoretical bound from Theorem I 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 looseness in the total variation constant M. For example, our bound for multivariate Gaussian distributions in Example Saample[[sugeessss suggests that MM oo when p which far more conservative than our empirical results (see Figure right); this is due to looseness in Pinsker's 's inequality (Appendix [B). Similar results are observed for Occlusion, Banzhaf, RISE, and LIME attributions, and these are included in Appendix
1.0-
Slope of Upper Bound
0.4
Marginal
Baseline
Baseline
Marginal
Conditional
Con ditional
0.3
Theoretical
Shapley
Marginal)
0.6
Theoretical (Baseline & Marginal)
(Conditional)
0.2
0.4
0.1
0.2
0.0
0.0
0.0
0.5
1.0
1.5
2.0
0.0
0.2
0.4
0.6
0.8
1.0
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 varies the logistic regression Lipschitz constant, and varies the correlation between the two synthetic features I1, T2.
A.2 Synthetic data with model perturbation
Here, construct logistic regression classifier as well 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 more robust against model perturbation, as suggested by Theorem
Setup. The synthetic data generation follows that of Appendix A.1 A logistic regression classifier constructed with coefficients [5, 3,3,77 1], and its perturbed counterpart f' has coefficients [o, 5, 1] Because the synthetic features 1 are identical with the default value p = 1, we have the functional distance If oolloo.xx = O, so our theoretical bound suggests that [lo(f,z) of', x) = 0 for conditional feature removal. When for conditional rather than marginal removal due to Lemma 4 because significant probability mass lies in a subdomain where f and f' are functionally close (near the hyperplane where Φ = x22). Feature attributions are computed for 100 explicands with the same feature removal settings as in Appendix A.I]
consider smaller values for p. we still expect better robustness
Results. As suggested by our theoretical bound model perturbation when p 1, the Shapley attribution difference between f' is near zero when features are removed using their conditional distribution, whereas the attribution difference can be large with baseline and marginal feature removal 5). When we decrease p in the synthetic data generation, the models f, become
16

---

p 1.0
0.8
0.4
0.0
10
alx, I)l
Shapley
0.6
||
0.4
0.01
Morginal
Marginal
Marginal
Marginal
Figure 5: Shapley value attribution differences between the logistic regression classifiers f. as constructed in Appendix A.2 with baseline, marginal, and conditional feature removal, and varying correlation p.
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 [[] 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 unified analysis for removal- attribution methods instead of focusing on the Shapley value, and we consider cases where models disagree even on the data manifold (Lemma
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 D suggesting that training a model with better Lipschitz regularity can help obtain more robust removal-based feature attributions. Weight decay one way achieve this, as 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.
consists 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, use MNIST 44 CIFAR-10 [2], and ten classes of ImageNet [19] (i.e., Imagenette)) From the official training set of MNIST and CIFAR-I0, 10, 10,000 images are randomly chosen as the validation
Datasets. The UCI white wine quality dataset [4 24
Setup. Fully connected networks (FCNs) and convolutional networks (CNNs) are trained with range of weight decay values for the wine quality dataset and MNIST, respectively, where use the values {0, 0.001, 0.0025, 0.005, 0.01). ResNet- networks [33] are trained for CIFAR-10 10, with weight decay values (0.0075, 0.01, 0.025, 0.05, 0.075). ResNet-50 networks pre- 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 2I). Details on model architectures and other hyperparameters are in Appendix | 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 x 4 as features in MNIST and CIFAR- 10, and of size 16 x 16 as features in Imagenette (see Appendix for discussion of how our theory extends the use of feature groups). For the wine quality dataset, removal- bassd attributions are computed exactly with all 2 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, 10,000 random feature sets, and Shapley and Banzhaf use least squares regression similar KernelSHAP 45
17

---

For CIFAR-10 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.
||x-x|||> 0.1
||X- XX|>> 0.5
x 2.0
0.020
0.00
0.150
0.25
Marginal
Basslim
Shap ey
0.323
0.20
Conzit onal
0.015
0.05
5,200
0.04
0.873
0.000
0.02
0.000
0.10
000 0.022 0SE
F006 0 968 0.000
1000 902 0S4
0:006 0 368 0.030
1.000 9 302
1.000 102 dos 006 a 505 0.030 0:006 o 508 0.030
Weight Decay
Weight Decay
Weight Decay
Weight Decay
[b]
|| x
| x 0.5
| Ix-x =.0
x |x-xxla =
110 461 xllls
Baseline Baseline
Shapley
0.08-
0.08
0,500 8,002 0,004 0.006 0.009 9,910 0,800 9,602 0,004 Weight Decay
6.006 0.009 9,836 0'000 0,502 6,004 0,006 0.008 0,530 0.000 0,502 5,504 6,006 0,008 0.510 Weight Decay
Wel Weight ght Decay
Weight Decay
Figure 6: Identical to Figure 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, 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 pairs.
95% confidence intervals across
[|x Xll2 -
|K x112
X -
[
X -
9.0375
8,0018
0.000
0.018
0.030
Dcdusioo
0.016.
9.6925
5.0300
3,0034
0.007
0.814
3.0250
0.005
0.012
3.0225
0.0900
MC x01ls
3.0075
3.5275
0.0275
3.0275
0.0230
0.0225
0.0225
0.0223
0.0255
19
0,0000
3.0200
0.0200
5.0000
3.115
5.0175
8.0150
3.0150
0.0150
5.0150
0.125
0.125
0.100
0.120
0.120
0.120
Shapley
0.115
0.115
0.315
0.155
0.115
0.105
0.105
0135
0.100
0.100
0.100
01:0
0.55
0.00
0.095
0.095
0.44
6,04
Weight Decay
Weight Decay
Weight Decay
Weight Decay
Figure 7: Attribution difference for ResNet-18 networks trained on CIFAR-10 with increasing weight decay, under input perturbations with perturbation norms 10.1, 0.5, 1, 2). Error bars show the means and 95% confidence intervals across applicand-ertuuraation pairs.
18

---

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-I0. 10, and with zeros (identical to the training set means due to normalization) for Imagenette. The wine quality training set used for removing features with their marginal distributions, while a conditional variational autoen- coder (CVAE) is trained following the procedure in Frye et al. ( to approximate feature removal with the conditional distribution (details in Appendix (]. 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 with 10 random perturbations applied per explicand and with perturbation norms {0.1, 0.5, 1. 2.2] for the wine quality dataset. MNIST, and CIFAR-10. For Imagenette, 5 random perturbations with norms (50, 100, 200, } are applied per explicand to obtain similar perturbation strength per input.
generated as in Appendix A.
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 measure explanation similarity using Pearson Spearman rank correlation rather than &g distance. For ResNet- 18 networks trained on CIFAR- 10 and ResNet-50 networks trained Imagenette, similar trends observed for Occlusion, RISE, and Shapley attributions (Figure and Figure 8).
||x xx']>= 50.0
|| xx|2 100.0
|| x122 - 200.0
| x1b
400.0
8.6275
0.000
0.300
9,0250
0.055
0.255
Occluson
0.050
0.12
0.200
8,0000
0.10
0.25
6,0175
0.00
0.025
0.08
0.030
0.150
0.025
0.00
0.000
5,0025
0.055
0.40
1 Att.
6.0o0n
0.050
0.35
RISE
3,6275
0.055
0.12
0.30
0,0250
0.000
3,0000
0.00
0.20
0.14
0.16
Shapley
elr,
0.2
0.12
0.03 Weight Decay
0.03 0.5
001
0.02
0.03
0.35
0,88
0.0
0.02
0.03
0.35
Weight Decay
Weight Decay
Weight Decay
Figure& 8: Attribution difference for 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
pairs.
A.4 Implication II: sanity checks for removal- -base 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 to sequentially randomize its parameters from the output to input layer, procedure known as cascading randomization introduced by Adebayo
19

---

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.
Baseline
Marginal
Conditional
1.0
1.0
1.0
Peanson Corslation
0.8
0.8
0.8
6.6
0.6-
0.4.
0.4
0.4
0.2
0.0
0.0
0.0
1.0
Rank Condittion
0.8
0.8
266
0.6-
0.6
2.4
0.4
6.2
0.2
0.2
6.0
0.0
Lnbir
L
Linar
original
Libir
Linar
Liaarz
Liwer O
Driginal
Lhear
Liner
Linarr O
Occlusion
Banzhar RUSE
LIME
Vanilia
IntGrad
Grad
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.
1.0
1.0
Pearson Correaaion
0.8
0.8
Rank Correlation
0.6-
0.6
0.4
0.4
0.2
0.0
0.0
Original
Unner
Linear
Conv
Convo
original
Linea 1
Uneer D
Occlusion (Meen Baseline) Banzhaf (Mean Baseline)
Shepley (Mean Baseline)
LIEE (Mean Baseline)
IntGrad
RISE (Mean Baseline)
Vanilla Grad
Grad Input
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 Appendix
Setup. The best performing networks from Appendix A.3 are used here. The settings for removal- based attributions are the same in Appendix A.3 except 50,000 random feature sets are used for more stable Shapley estimations in Imagenette. Parameters for each network are randomized in 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
20

---

3 times for Imagenette. Gradient bradient-based attribution methods such as Vanilla Gradients [56], Integrated Gradients [60], and Grad 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 2}. This shows that removal-based attributions reliably pass the model randomization sanity check. In contrast, we observe that for Integrated Gradients and Grad X 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 Input (Figure [T0). Similar trends are observed for the ResNet- trained on CIFAR-10 and ResNet-50 trained on Imagenette (Figure 11
Figure 12. and Figure 29)
Pearson Correlation
0.8
Rank Carrelation
0.8
0.6
0.6
0.4
0.4
0.2
0.0
Origiaal
Unear
Layer4.1 1. Cono22
Layers,1 1. Convi
Conv2.
Layer4 O Corv1
Layer3.1 1. Corv1
eyyer3.0.Convz
w..r3...coovv Convl
Layer2. 1.oov2
1. Convi
Layer2. Convv..
Layerz 0. Conv1
Layer1 1. Conv2
Layer1 1.ver.1..oom Convi
Layeri 0. Conva.
Conv1:
Original.
Linar-
Layer4, 1. Conv2
Layer4, Convi
Layerd O Corv2
Layer4.0 Convi
Layer3.1 1.1 Conv2.
Layer3.0 0.0.3.0..O0 Conv2
Layer3.0. C.novl
Layer2.1 1. Conv2
Syver.1.Coovv
Layer2.0 2...r....oo2222 Conv2.
Layer1 1. Cono2
Layer1 bayer1.1.oom 1orrv
avver....com Cono2
Layeri .. Convi
Layer3,
Layera
Layerz.
Layors.
Layer2.0
Occlusion (Mean Baseline)
RISE (Mean Baseline)
IntGrad
Grad Input
Shapley
Baseline)
Vanilla Grad
Figure 11: Sanity checks for attributions using cascading randomization for the ResNet-18 trained on CIFAR-10. Features removed by replacing them with training set means. Attribution similarity measured by Pearson correlation and Spearman rank correlation. We show the means and 95% confidence intervals across 10 random seeds. Missing points correspond undefined correlations due to constant attributions under model randomization.
Pearson Correlation
1.0
Rank Comrelation
0.8
0.6-
0.2
0.2
0.0
Original
LOR
ESRV
ESRV
CORV
Conv
Comve
Conv
Conv
Cond
Original
Comv2
Layerr 22
Occlusion (Zero Baseline)
RISE iZero Baseline)
IntGrad
Grad Input
Shapley (Zero Baseline)
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 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.
21

---

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 (as) for any feature set g L--ipschitz continuous:
lf(xs) f(r's) - L 1|xs dlslla V 2ssr's E R
Proof. The key insight for this result is that removing features using either the baseline or marginal approach is equivalent to using distribution gzz)) that does not depend on the observed features. Because of that, we have:
-(ss)) - s(''s) / (as, (()s)))) /
((d's
< / [M (J)) - s(r's k5) ((zg)des
$ L. 1 z's s1l12 / q(xg)d23
= L. ||xs - z's 12
Lemma 2. When removing features using the conditional approach, the prediction function f(s) for a feature set æi satisfies
((ss) - [(ds)) L.
L. 1|55 d'sll2 + 2B drv (plars (rs),p(ag z'ss):
where the total variation distance defined via the L functional distance as
drv (p(xs p(as) a's)) p(æs xg) p(æs a''s)l;
1
Proof. In this case, we must be careful to account for the different distributions used when marginal- izing out the removed features.
/ f(as. 2z)p(a3 zs)drş /
If(=s) (iss)-f((sa)) f(x's)]
=g)s,, d's)drs
(((xs)z))((z æs) f( xs)p(rs s's) dzs
$ (((5522)p(as æs) f(d's x p xs) (dzs
+ ( as)) / ( 2 ))((s (ds)
$ / [((s) (s)ss))))) p(as as))do
+ / [ 1p(ss æz) - p(as x's) drs
< L. ||æ - x'slle / plas us)dx + B / |plss æs) p(zs Is's)drss
= L- saa - L'sll + 2B . drv (p(rs as). a's)).
22

---

Example 1. For Gaussian random variable a ~ (H,,E) with mean E R and covariance E E gaxd, Assumption holds with M =1 Ammaax((2 (2-1)~ Amin (2-1). If is assumed to be standardized, this captures the case where independent features yield M means that if one dimension is roughly linear combination of the others, or if Xmax is large, the conditional distribution can change quickly as a function of the conditioning variables.
Intuitively, it also
Proof. Consider a Gaussian random variable æ with mean μ E IRd and positive definite covariance E Raxd We can partition the variable into two parts, denoted æ1 and æ2, and then partition the variable and parameters as follows:
11
E1 21
Ei E2
x
Σ
222
122
We can then consider the distribution for æn conditioned on specific 22 value. The conditional distribution is æι 12 ( Σ) with the parameters defined as
p= p= p1 E12 (32 12)
Σ En E1 E2
Our goal is to understand the total variation distance when conditioning different values 22, 2/2 Pinsker's inequality gives us the following upper bound based on the KI. divergence:
1 DKL p(z1 (2) |l p(æ1 x²))
drv (plz) 22). p(ø1 x22) -
The KL divergence between Gaussian distributions has closed-form solution, and in the case with equal covariance matrices reduces to the following:
DKL. (p(zr) (x)) || p(z) x52)
1 --((2 (x2
E2 (En1 Di1
E12 (x2 x12)
To provide upper bound this term based ||اا x/2 21.2. it suffices 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:
mー1=
(->
EEn - En E21)
E12 E21 Ej2Sigh (Enr
-E E2 (En1-
E12 Ezg E21)
-1
En
Egg
E21
(En1
Σ Eng
I-IIa1 -1
Ez1)
En
We therefore have the following, because the matrix we're interested in is principal minor:
Amax
Mnnl Mnn 221 (En Ei22igh Ei122ng s Amxx
Weyl's inequality also gives us:
Amin (Egz') + Amax EmnEnn Ev1 EEn En Ee1)
EizEzd
Amax
Eng
E21
EEn
Σ12 Dn E21)
Combining this with the previous inequality, arrive the following:
23

---

A Amax
Eng E21
EE1
S2 E21)
E12
S x mmxx (
Amiu (Emz2)
For the subtractive term, we can write:
Amin (E22)
Ami
(m (--1).
Amax (2122
Amax
We therefore have the following upper bound on the KL divergence:
DKL p(21 x2) ll p(æ1
x24)
1 -
| 2'212 (Amax() - Amin Aain(--)) 1)
2
Combining this with Pinsker's inequality, arrive the final result:
drv (p(21 (22).p(a)) x221) S || 2/212 v Amaax
Amiin ( (E-1).
This means that we have M A Amiin (E-1).
Notice that this does imply M 0 when the features æ are independent: it only captures this for the case of isotropic variance. However, if we assume that æ is standardized, so that each dimension has variance equal to M 0 is implied for independent features because have Σ I.
Lemma Under Assumption 3 the prediction function fxs) defined using the conditional ap-
proach is Lipschitz continuous with constant + 2BM for any feature set
Proof. The result follows directly from Lemma and Assumption
Is(zs) - f(''s)) < L ||ES x 2B drv (p(2s) xs). a's))
< L. ||x5 - 23312+2BM 2BM | z'sll2
= (BM) ||05- r'sll2.
24

---

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' IR + IR and a subdomain X' C TRd, the prediction functions
f(+s), I (xs) for any feature set x5 satisfy
(f(z)) f' (zs) Il x' Qzs (2') + 2B (1-Qag(X"))
where Qzs (X') is the probability of imputed samples lying in X' based on the distribution (<g)):
Qay (X') = Byy(s)) [T ((xs, a. E xx
Proof. The main idea here is that when integrating over the difference in function outputs, we can separate the integral over (Is, Øg) values falling in the domains X' and IR \ X'. Because we have single value for the observed variables sg, the distribution g( is identical for both terms regardless of the removal technique.
/ fass ((5s)))) / f (ass (((ss)))
If(zs) - f(s)) f
Σ / [[(rs - (xs, g(zg)drs
/ [((ss - (xss, ((ss. E X²) g g(xg)drs
+ [(ss, 25) 8' [zs, m3) ( as) & x') ((g))dr
$ f1 f'loo,x / 1(TS) 25) E X3} a(rg)das
+ +BB / /1((FS,0)) $ x') g(zs)drg
|| f'loo, Qag (x) + 2B (1 Qus (2'))
Lemma 5. When removing features using the conditional approach, the prediction functions fl (zs) for two models f and any feature set 25 satisfy
(f(xs) ('(xs) S 1 floo.
Proof. The result follows directly from Lemma 4 when we set X' X. Note that when we remove features using the conditional distribution, we have Qus (X) = 1 for all xs lying on the data manifold, or where there exists gg such that p(as, xs)) > 0.
Lemma 6. When removing features using the baseline or marginal approach, the functionns f( (xs) for two models f, f' and any feature set satisfy
If(xs) - (xs) S llf f lla
Proof. The result follows directly from Lemma Qzy (I) for all 2g IR|S
when we set X" Rd. because we have
25

---

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 linear operator
A Rx vector IR? representing the predictions with each each feature set,
off. xff.z)= = Av,
where the linear operator A for each method is listed in Table 3 and v is defined as f(xs)) for each _ [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 generates attribution scores. Assuming fixed input where f(xs) denotes the prediction with feature set, each method's 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 (f(f.x) ) is defined as
=( = - ((
We refer to this summary technique as "leave-one-out" in Table
Shapley value 45 54 : this method calculates the impact of removing a single feature, and then
averages this over possible preceding subsets. The attribution score @i (f, 2) is defined as
1 Σ d ()(rrss))))))) -((s))
(((,,))=
SEIVOO
Banzhaf value 9 13:
this method similar to the Shapley value, but it applies a uniform weighting when averaging over all preceding subsets. The attribution score (f. is defined as
E (f(Esu(O) -((ss)))-((ss). SCOMA
of(f.z)
2d-1
RISE [51]: this method is similar to the Banzhaf value, but it does not subtract the value for subsets where feature is not included. The attribution score (f x) defined as
Φ1 x)=
Σ f(
2d-1
SCIVNE)
We refer to this summary technique as "mean when included" in Table
LIME this method more flexible than the others and generally involves solving a regression problem that treats the features inputs and the predictions labels. LIME can be viewed as a removal-based explanation when the features are represented a binary fashion, indicating whether they provided the model not [17]. There is also choice of weighting kernel (SS) 2 0, of regularization term, and of whether to fit an intercept term.
We assume that the regularization is omitted and that we fit intercept term. Denoting the model's coefficients as (Ba....,,B, Ba),
the problem is the following:
Σ (S)
Σ 3, 3 f iES
min So.....9
SC(4]
We can rewrite this problem in a simpler fashion to derive its solution. First, we let W e [0, 1 be diagonal matrix containing the weights w(S) for each subset. Next, we let
26

---

the predictions be represented by vector e E IR? where vy = f(zs), and we let the matrix Z {0, 1 enumerate all subsets binary fashion. Finally, to accommodate the intercept term, let = (So...., ,aa E Rd+! represent the parameters, and we prepend columns of ones to Z, resulting in Z' E {0, 1) xd++1, The problem then becomes:
......,,).
min (Z - v) W(Z33 0') = Z' Ilw.
Assuming that Z W'' is invertible, which guaranteed when w(S) for all S [d]. this problem has closed-form solution:
8 =- -Z WZ)
LIME's attribution scores given by (z(f.z) fr for i 1, ,..... d, where discard the intercept term. We can therefore conclude that each attribution a linear function of u with coefficients given by all but the first row of (Z' W W 77)-11 Z W. If denote all but the first of (Z' W W
(Z" W Z') [1:j then have
(Z" W Z') [1 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
= (Z W ' W. We can similarly allow for ridge regularization: given a penalty parameter
0, we then have A (Z T WZ Xlaa) W.
These coefficients depend on the specific choice of weighting kernel w wSS). and previous work has shown that leave-one-out, Shapley and Banzhaf values all correspond specific weighting kernels [17. Our analysis of those methods therefore pertains to LIME, at least in certain special cases. However, LIME often uses heuristic exponential kernel in practice, which 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
We therefore that in each method, (f( is linear function of the predictions for each S [d]. Across all methods, we can let the predictions be represented by a vector v E IR." and the coefficients by matrix A E _dx22".
Lemma 7. The difference in allributions given the same summary technique A and different model outputs v, v' can be bounded as
||Av Ao'll ||Al22 0'112
||A - Au'll2 S [|A]]1,00 111 of lloos
where ||A||2 is the spectral norm, and the operator norm
||
1,00 is the square root of the sum of
squared row I-norms, with values for each given Table
Proof. We first derive the two matrix inequalities. The first inequality follows from the definition of the spectral norm, which is defined the maximum singular value:
|uull |ul|2
||||| = sup
T
V Amax (AAT).
a
We prove the second bound using Holder's incquality as follows, where A, denotes the ith row:
d
Auull
Σ (AT
Σ AAR |ul 1|///oo
Σ | = |||| I|A
In the special case where all of A's row 1-norms are equal, or where we have ||::|| ||Ay|| i for all the simplified expression ||] ||A1]|1 Vd.
i,j j [d].
27

---

Next, we derive the specific norm values for each linear operator in Table 3
Beginning with the operator norm || remark that several methods have equal 1-norms across all rows: following their interpretation as semivalues [46], leave-one-out, Shapley and Banzhaf for all [d]. This yields the
all have | 2 for all E [d]. Similarly, RISE has |A:Il
norms ||||,, 2 for leave-one-our Banzhaf and Shapley, and ||AA||i,co v for RISE.
For the spectral norm, we consider each method separately. We find that for all methods, is simplest calculate AA and then Amax AAax entries of the matrix are determined by the inner product (AAT)s ) A. Aj. For leave-one-out, we derive AA as follows:
find the maximum singular value maxx
Fiax(A), (A). Note that the
[2 if (AA
otherwise.
The matrix AAT is therefore the of rank-one component and identity term, and have AA La + Ia. Recognizing that 1 is the leading eigenvector, we get the following result:
(Lale + Ta)
Lad + 1a
rmax (A) / Xmmax (AAT)
+d+1. 1.
|all2
|all
For the Banzhaf value, we derive AAT as follows:
11/22d-2 if
(AAT): )g
10 otherwise.
The max eigenvalue is therefore Amax (AAT) 1/20-2, which yields |||l2 1/2 For the mean when included technique used by RISE, we derive AAT as follows:
[1/2d-1 fi= 11/2d otherwise.
(AAT) )ig
The matrix AAT is therefore the of rank-one component and identity term, have AAT = 1/la/2d + Ta/2d. Recognizing that 1 the leading eigenvector, we get the following result:
Tmax (A) = V Amax (AAT)
/ (1ola/2d 22 + La/2d)1a
1ad/2d + 1a/2d
d d+1
V
1
|al|2
2d
For the Shapley value, we derive AAT by first considering the diagonal entries:
28

---

(IS) (15) 1)1(d (S|)!
|S|!(d= [S| 1)! 2
)a==== AAA Σ
+
d!
d!
S:iES
SitES
d k ( ( (
1)((d (k-1)1(d- 1:)1
d (d-1) (
[bldd - 1s - 1)1 27
d!
k
d!
1)1(d= - k:)!
d-1
k!(dd
k 1)!
d!
d
d!
k=
(k:
1)((d - k- - 1)!
2
d!
d
k=
d-1 2
-1
Σ
+
d(d - 1) k=
k- 1
dc
We next consider the off-diagonal entries:
(IS) 1)!(d (S])! (IS) 1)((d |SD
(AAT): Daj ATA,
Σ
d!
(IS) 1)l(d IID! s( - -1S) [S] 1)!
Σ SHES
d!
dp
|S| |S| 1)! (15) 1)1(d IS))
Σ SaES,jES
d!
|S|!(d |S| 1)! sSldd [8| 1)!
d!
dp
S:ESS
d 2
(k: 1)((d - k)! (k: - 1)((d - k))
k- 2)
d!
d!
dd -2-
2)
(k: 1)l(d= - k)) k!(d - k: - 1)
k 1
d!
dp
d
(d - 2)
k!(d - 1)1 kl(d k 1)
Σ K=0 k
d!
d!
(k 1) (k: 1)!(d - k)!
Σ
d(d 1)
d!
d-1 1
klld - K)!
d(d- - 1) d!
k=1
d-2
(d - k 11) k!(d - k: - 1)!
k= d(d - 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 & values simplify the above:
29

---

kk: 1)1(d- - k 1)!
2(d 1)
2
A
d(d 1)
(d - 1)!
km 0,kmd d
1
d-2
(k 1)l(d - k - 1)!
2
2
did 1)
1d=
d
d(d 1)2
k=2
d-2
d - 2
-1
2
Σ
+
d(d - 1)3
h: -
d
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
1- (d-2) kkml k-1
2
i otherwise.
(AAT); Dog
d(d d-2
+
dd- 1)2
From this, we can see that AAA is the sum of a rank-one component and an identity term, or that we have AAT 1a1ab + Ia(e b) with values for b (AAT), )ij and = (AA) )a given above. This matrix has two eigenvalues, C - and db + (c b), and we must consider which of the two is larger. Under our assumption of > we can observe that
2
d-2
2
-1
b b= (AA
)oj =
d
ddd - 1)2
1)
d(d - 1)2
k=
2
d-2
d-2
-1
2
d²
d(d - 1)2 km
1
d(d 1)2
2
d-3
2
d d(d - 1) 22) d(d 1)2
2(d - (d 2) - d(d 3) 2d - 2)
d2 - () 2)
It can be verified that the numerator of the above fraction is positive for d > 2. This implies that
db- (c b) the larger eigenvalue. We can therefore calculate Amax (AA as follows:
Amax (AAT) = db + (c - b) = c (d 1)b
2
2(d -
2
d² + d(d - 1)
d2
d(d - 1)
k=1,k=d-1
2(d 1) + 2 1)2
d3 1)
d
On the other hand, when we have d 2 have the following entries for AAT.
(AAT0a = {1
ifi=j otherwise.
30

---

This yields Amax (AAT) Finally, this let us conclude that the spectral norm for the Shapley value the following:
) = 1, which coincides with the formula Amax (AAT ) = 2/d derived above.
|| V Amax (AAT) 12.
Lemma 8. When the linear operator A satisfies the boundedness property, we have ||A||1,0 2v
Proof. Following Theorem in Monderer and Samet [46], solution concepts satisfying linearity and the boundedness property (also known the Milnor axiom) are probabilistic values. This means that each attribution is the average feature's marginal contributions, or that for each feature E [d] there exists distribution pa (S) with support the power set of [d] {i}. such that the attributions are given by
= Ep((s) (((SS)))))))))) f(x)]]
E ((()((()))))))) SEIA(()
f(ss))
(D.0.1)
The corresponding matrix Rdx2" therefore has rows A, with entries A,s 2 0 when E S and Ais & 0 when & S, where the positive and negative entries sum to and -1 respectively. We can therefore conclude that each row satisfies ||A,l|a 2, and that || 2vd.
Lemma 9. When the linear operator A satisfies the boundedness and symmetry properties, the spectral norm is bounded as follows,
24/2-1 11AA|22 S Vd+1,
with 1/20/2-1 achieved by the Banzhaf value and d achieved by leave-one-out.
Proof. Following Monderer and Samet [6). solution concepts satisfying linearity, boundedness and symmetry are semivalues. This means that they are probabilistic values for which the probability distribution PS(S) depends only on the cardinality [S] (see the definition of probabilistic values above). As result, we can parameterize the summary technique by how much weight places on each cardinality.
1,...,d, d. We write the Let Rdd denote probability distribution the cardinalities k 1,..., attributions generated by any semivalue as follows:
off.z) Σ an Σ
d k (--1)
(f(rsu)))) f((s)).
SCI 1
Note that this resembles eq. (DO.1). 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 E Rdx22 with each semivalue, similar to Proposition Furthermore, we can let be linear combination basis elements corresponding to each cardinality. We define the basis matrix (A) Rdx22" for each cardinality 1,..., as follows:
d
if E S and |S| = k
[S]
1
A
d-
if & S and |S| k: otherwise.
0
31

---

With this, we can write the matrix associated with any semivalue as a function of the cardinality weights
d A(a) = -e ag (4)) (4)
1=1
An important insight is that by composing linear function with the convex spectral norm func- tion [T], we can see that [|A(a)) ll2 is convex function of a. This lets find both the upper and lower bound of [[a)ll subject to being a valid probability distribution.
For the lower bound, refer readers Theorem C.3 from Wang and Jia 664 who prove that
the Banzhaf value achieves the minimum with a² such that oj (6=1)/2d-1 and spectral norm ||A(a Ill 1/2d/2-1, Notice that this is a valid probability distribution because we have EES (2=1) 2 Also note that the proof in (64 parameterizes slightly differently.
For the upper bound, we can leverage convexity to make the following argument. The fact that spectral norm is convex with a representing probabilities lets us write the following for any a value:
0 max (A(a)) = max
au Σ k=1
A(k)
Σ as max (A S max { (mmax (A(h))) }.
&
Finding the upper bound therefore reduces to finding the basis element with the largest spectral norm. For each basis element A(A), we can derive the spectral norm by first finding the maximum eigenvalue of (k) ( T E Bdddd In deriving this matrix, we first consider the diagonal entries:
A(())
E
|S| 111-1)
d
E
(d 1
-2
|S|
S=ES,1S
|=
( (d-1)
d (
+ d ( (1-11111
-2
k -
Next, we consider the off-diagonal entries:
(A())A(A) ),, = Σ
d (111
Σ
d 1
[S] 1-|--1
|8|
S.ES.JES,(SI =
Sad S = k-1
dd 2
d d-11
-2
( (d-2) 2)
(d (d-1)
-2
2
1-1
1: - 1/
4 4-1
d
d- k d 1
+
d k-1
d - k 1-11 1
To summarize, the entries of this matrix are the following:
(2(=)) ifi=j
((0)A(0)))))) ",, =
otherwise.
32

---

We therefore see that A(k) is equal a rank-one matrix plus an identity matrix, or ( A(R) bl,1j (c ba with values for b = (A(k) A(k) )ig and (A(k) A )x+ given above.
Recognizing that b > and therefore that 1a is the leading eigenvector, have the following expression for the maximum eigenvalue:
Amax (A)))
cc (d 1}b
(A(4)) ) (d-1), (d(1))((4()))
(d-1) -1
= (d + 1)
k: -
Finally, we must consider how to maximize this as a function of k. Setting k: = or k: = d achieves the same value of d which yields the following spectral norm:
|| 1|2
V Amax A(d)) = Vd+I.
The case with = corresponds to leave-one-out (see Lemma and = 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.
22
A|ll2 2 c0s
d+1
with 2 achieved by the Shapley value.
Proof. Following Monderer 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 notation,
to denote ordering over [d]. and T(i) to denote the position of index within the ordering. We also let Pre(m, i) denote the set of all elements preceding in the ordering, or Pre(m, i) (j : (s) r()). Similarly, we let Pre( i, j, be indicator of whether directly precedes in the ordering, or Pre( I, i, [[((()) = (i) - 1}. Finally, we use to denote the set of all orderings, where II] d!.
Now, given a distribution over orderings p(T), we can write the attributions for any random-order value as follows:
dx(fx)) = Σ (())) [f(an(((ommt)))ttt))) - (spre(t.a))
*E
Σ
Σ
p(T)) (n))) ( s(as))
SCIdM
T:Pre(T.i)=S
Based on this formulation, we can associate matrix Rdx2 with each random- value based on the coefficients applied to each prediction. The entries are defined as follows:
S0 P(T) if S
A
-
Emde((s)ss otherwise.
Furthermore, view A as linear combination of basis elements ( E 2 corresponding to each individual ordering. We define the basis matrix A (*) associated with each Il as follows:
33

---

Pre( i) - S {i} if Pro(T,, i) S
A
otherwise.
With this, write the linear operator corresponding distribution p(π)) the following linear combination:
A(p) = Σ p(n) A WEll
By composing a linear operation with the convex spectral norm function [11. we can see that [|A(p)] ll2 a convex function of the probabilities p(T)) assigned to each ordering. This will allow us to derive form for both the upper and lower bound on [pA(p)]
Beginning with the upper bound, can use convexity to write the following:
max (A(p)) max
ΣP(T)
p(T) A(
ΣP(T) o max (A( I max maxx (A(r)) }.
xell
TE
TENl
Finding the upper bound therefore reduces to finding the basis element with the largest spectral norm. For an arbitrary ordering æ,, we can derive the spectral norm via the matrix A(x) A(r) whose entries given by:
2
if jj
(A())
if j and Pre(x, 4, if and Pre(z, otherwise.
is
Within these matrices, the diagonal entries are always equal to 2, most rows contain two -1's, and two rows (corresponding the first and last element in the ordering α) contain a single -1. We can that the matrices A (m)) are identical up permutation of the rows and columns, that they are similar matrices. They therefore share the same eigenvalues.
Without loss of generality, we choose to focus on the ordering x' 'hat places the indices in their original order, or N (i) = for E [d]. The corresponding matrix A has the following form:
1
0
0
0
-1
-1
0
0
-1 2
-1
0
-1
-1 2
This is tridiagonal Toeplitz matrix, a type of matrix whose eigenvalues have a known closed-form The maximum eigenvalue is the following:
solution 49
Amax (1()) ((()))(=T))
= 2 (-1) (-1)
COs
d
= -2+2 2 cos
d d+1
34

---

We therefore have the following upper bound on the spectral norm for a random-order value:
||A||2 S // 2 2- cos
d
This value tends asymptotically towards 2, and is achieved by any random-orded value that places all its probability mass single ordering. Such single-order value approach has not been discussed the literature, but is a special case of proposed variant of the Shapley value [26].
Next, we consider the lower bound. We now exploit the convexity of [|A(p))lz to argue that specific distribution achieves lower spectral norm than any other distribution.
Consider an arbitrary distribution p(T)) over orderings, with the associated linear operator A(p) and spectral norm AAp 212- Our proof technique is to find other distributions that result in different linear operators A that share the same spectral norm. For any ordering x' E 11, we can imagine permuting the distribution p(TT)) in the following sense: we generate new distribution ps (TT) such that par (T) = p( (T<((')), where π (Tx') is a modified ordering with (( (() π (T)(i)). That is, the probability mass reassigned based on a permutation of the indices.
Crucially, the linear operator A(px) is different from A(p) but shares the same spectral norm. This is because ((p=)A(p=>)) is similar to A(p) A(p) T, or they are equal up to a permutation of the rows and columns. We therefore have [|A(pr)) 2 [|A(p) 2l,, and due to convexity we have
1 |A(pp22 - 5 1114(p>))22 -1A())lla ||A())]l
This does not imply that p/2 + Px. 2 is minimizer, but we can now generalize this technique to identify a lower bound. Rather than generating single permuted distribution, all possible permuted distributions. That is, we consider the set of all d! permutations n' E IT and generate the corresponding distributions {p-' x IT]. Notice that when take the mean of these distributions, the result identical regardless of the original distributions (T),, and that it assigns uniform probability mass to all orderings, because have:
propose generating
Pr
1 d! r'ell
Σ
p()))
d! E'Ell
d!
We can then exploit convexity to write the following inequality:
| Σ px pdt |1 S [|A())]l.
T'Ell
We can make this argument for any original distribution so this provides a global lower bound. The mandom-order value with uniform weighting is the Shapley value. whose spectral norm derived in Lemma 7 We therefore conclude with the following lower bound for random-order values, which achieved by the Shapley value:
2 || 2 V d'
Lemma 11.
When the linear operator A satisfies the boundedness. symmetry and efficiency proper- ties, the spectral norm is ||Al|/ 22/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 [66].
35

---

E Proofs: main results
We now re-state and prove our results from Section
Theorem 1. The robustness of removal-basee explanations to input perturbations given by the following meta-formula,
ll(f. off S g(removal) h(summary) || x' 2211
where the factors for each method are defined follows:
g(removal) = [L
if removal baseline or marginal
L + 2BM if removal = conditional
(2vd
summary Shapley, Banzhaf, or leave-one-out Vd if summary mean when included.
h(summary)
Proof. Given two inputs æ and x', the results Lemma | and Lemma Soainaaiing imply the following bound for any feature set æs.
((a's) S L' || z's112.
where have I' L when using the baseline marginal approach and L' L + 2BM when using the conditional approach. This bound depends on the specific subset, we also have a bound across all subsets because || ='sll2 s |ær x' for all S [d]. In terms of the corresponding
vectors v E IR="". this implies the following bound:
|| - o' 10x
max
Its
d'sl
< L' ||2 2'112.
sc
Next, Lemma 7 that we have the following bound on the difference in attributions:
((f=)))((f ||AAu ||Ac S ||A] 1| olo
where the norm 1 depends on the chosen summary technique. Combining these bounds, arrive at the following result:
[ - +(Sx'')la S A| L', || 221
Substituting in the specific values for L' and || completes the proof.
Theorem 2. The robustness of removal-based explanations to model perturbations is given by the following meta-formula,
[Kffx)) - 6(f' x))l S h(summary) Ilf
- f'll.
where the functional distance and factor associated with the summary technique are specified as follows:
[ [1. - f -
if removal baseline or marginal if removal conditional,
If f'll
(2vd summary Shapley, Banzhaf. or leave-one-out "summary mean when included.
h(summary)
Vd
36

---

Proof. Given two models f and f', the results in Lemmas and iimply the following bound for any feature set as.
(f(as) - f' S 11 f'Ils
where the specific norm is ff when using the baseline marginal approach and I f' when using the conditional approach. In terms of the corresponding vectors v, E R²", this implies the following bound:
|/ 10
max Us Scld
u'sl - fl
-f'll
Next, Lemma showed that we have the following bound on the difference in attributions:
lØ(f. 2)) o(f', = ||Av - AA - ||A||1,0 11v0~
11u-ed e'lloo
where the norm |, depends on the chosen summary technique. Combining these bounds, we arrive at the following result:
[ ((5",))1: S || 1'I1-
Substituting in the specific values for |l - fll and || 1,o0 completes the proof.
Next, we present 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 is easy to analyze given our previous analysis.
Corollary 1. The robustness of removal-based explanations to simultaneous input and model pertur- bations is given by the following meta-formula,
[((f,x) - 6(f', x') 222 & g(removal) h(summary) || - 1'112 + h(summary) I f'll.
where the functional distance and factors associated with the removal and summary technique are given in Theorems andd2
Proof. This result follows from triangle inequality:
[Ø(5,x) "('',""''' S Il((,x)) *(S, x' 122 + [0((5x')) &(', ((',"'') x2)1
The bounds for the first and second terms are given by Theorem and Theorem 2 respectively. D
37

---

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 g(wz) for removed features (see eq. 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(rs) with a subset of features and its estimator f(xs) with independent samples from g(æg).
[( x5 (1) be the empirical estimator of ((rs) in eq. where Lemma 12. Let f(rs) = :. I
(6) Lid. g(zg). For any å E (0,1) and > 0. if the sample size satisfies 2 2B3 log(2/0) then Is
(1((x)) f(xs) S 2 8
Proof. With Assumption 2. have B
B. Also, B[f(=))] f(zz) because I'g
n
m
sampled from g(zz). Applying Hoeffding's 's inequality
(35] obtain the probabilistic bound
2ea
P(I((es) f(rs)) s 2
(2B)
Lim1
m
2BB² log(2/6)
e2
ms²
= 2exp
S 2exp
= 5,
2B
2B²
and hence for the complement event we have
P(I((zs) f(xs)) 5c)
Lemma 12 can be applied to bound the difference between the original and perturbed attributions
when the integral in eq. (T) is estimated via sampling. For example, under imput perturbations, the difference [f(zs) f(z's) decomposes into !f(xs)
+ > - f(x'g))+ ((a))) 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. (e.g..
Lemma or Lemma 3. Intuitively, by combining probabilistic bounds over all the feature subsets, is possible to derive final probabilistic upper bound with the following meta-formula
[l(ff - Φ(f, f x x2) I gremovall) hssummary) || x' + Ismpling (removal), (F.1.)
where Isumpling (removal)
function of 8, E, m when removal marginal or conditional, and Ssmpling (baseline) = 0 because no sampling is needed. Similarly, for our bound under model perturbations, the term Ssampling (removal) can be included account for sampling in feature removal.
F.2 Empirical observations
In eq. F.I.II. the sampling-dependenl term Isampling (removal) does not depend on the input pertur- bation strength ||= - xx') 1|2, so we empirically estimate Isampling (removal) with the intercept of the empirical upper bound (as illustrated in Figure 3 when || - z [| O. As shown in Figure 13 (left), Isampling (baseline) is indeed zero, whereas Isampling (marginal) and Isampting (conditional) empirically decrease with increasing sample size for feature removal as expected. We observe that the sampling enror tends to be lower for the conditional approach compared to the marginal approach, potentially due to lower variances in the conditional versus marginal distributions. Finally, Figure [3 (right) confirms that the cntturaaiiondeppndeen term g(removal) h fssummaryy || - 221l does not depend on the removal sample size.
38

---

(Dependert on Iriuu Perturbation)
Thecratical
0.005
Ma rginal Conditional
(Beseline Marginal)
Intercept of Upper Bound
(Dependent on Sampling)
Stopa of Upper Bound
Marginal
0.25
Conditioral
Occluston
0.004
0.20
6.003
0.15
0.00
0.001
0.5
2500
5200
7503
10000
325.00
19503
17500
20330
0.0
2500
5263
7500
15030
12500
15003
17500
20003
(Dependent on inpuut Perturbation)
0.005
0.30-
Intercept of Upper Bound
(Dependent on Sampling!
Stope of Upper Bound
0.25
Banzhal
0.004
0.0-
0.15
0.10
0.001
0.5-
2500
5200
7503
10000
325.00
15003
17500 20020
0.00
2500
5203
7500
10030
32300
15003
37500
20000
0.007
0.55
(Deperident on Iriuit PerturBation:
0.006
0.30-
Intercept of Upper Bound
(Dependentt on Sempling)
Stope of Upper Bound
025-
Shappey
5.001
0.00
0.15
0.10
7503
10030
12500 15303
10030
0.007
0.35-
(Deperdeent on Input Perturbattonn
0.30
Intercept of Upper Bound (Dependent on Sampling)
Slope of Upper Bound
0.5
RISE
5.004
0.20
0.002
7503
13303
13003
5.007
0.35
Slope Upper Bound (Dependent on Inpit Perturbatton)
0.30
Intercept of Upper Bound (Dependent on Sempling))
0.003
025-
LIME
0.0
0.002
0.000
Samples
Number Samples
Figure 13: Intercepts and slopes of empirical and theoretical upper bounds with respect to the input perturbation norm, for Occlusion, Banzhaf, Shapley, RISE, and LIME attribution differences in the synthetic experiment in Appendix A.I and with varying sample size for estimating feature removal.
39

---

G
Feature grouping
Rather than generate attribution scores for each feature z;, another option partition the feature set
and generate group-level scores. This is common in practice for images [ 16 53 where removing individual pixels unlikely to produce meaningful prediction differences, and reducing the number of features makes attributions less computationally intensive. For example, LIME pre-computes superpixels using a basic segmentation algorithm [53). and ViT Shapley uses grid-shaped patches defined by vision transformer [16].
When computing group-level attributions, {G....., Gg J. where each group satisfies G2 C [d]. We assume that the groups are non-overlapping and cover all the indices, so that we have G2 nG, approach, the attributions are vector of the form off, æ) E Rs, and we compute them by querying
require
of feature index groups that we denote as
O for i + j and G₁ UU..UG = [d]. Under this
the model only with subsets that are unions of the groups. For this, we let T C [g] denote subset of
group indices and denote the corresponding union of groups as Gr = Uier G3 S [d]. Using this
notation and taking the Shapley value as an example, the attribution for the ith group G, is defined as
d((,,)==
1 Σ 9 ()
((mi)
9
(J(ECrop)) ((xr)),
where fIG,) is defined based on the chosen feature removal approach.
To generalize our earlier results this setting, we first remark that querying the model with each subset induces vector E ? containing the predictions for each union GT with T S [g]. The attributions are then calculated according to ø(f, x) Av. where A E R2x2 is the linear operator associated
with each summary technique (see Table 3 Our approach for bounding the attribution difference
under either input or model perturbation is similar here, because we can use the same inequalities
from Lemma 7 or ||A(r 12)l < || l22 |/ 0'112 and ||A(v - -2' < 1, ||u - o'llac.
depending which norm is available for the vectors of predictions |/e w'll-
Following the same proof techniques used previously, we can guarantee the following bounds in the feature grouping setting.
Input perturbation. Lemma 1 and Lemma guarantee that f(x)) is Lipschitz continuous for any subset S C [d]. These results automatically apply to the specific subsets induced by feature grouping. If fxx)) is L'-Lipschitz continuous for all S C [d]. we therefore have the following inequality,
(f(#or) f('c,) )1 s L' ||=Gr 2lc,
for all group subsets T S [g]. As a result, we can guarantee that the induced vectors v, e IR2 for two inputs x, x' have the following bound:
|| - o'llo S L'' 1| 2'112-
Finally, we see that Theorem I only needs to be changed by replacing with in the h(summary) term. This bound may even be tighter when we use feature groups, because the Lipschitz constant may be smaller when we restrict our attention to a smaller number of subsets.
Model perturbation. Lemma 5 and Lemma guarantee that the predictions f(xs) and f'(xs) for two models f' can differ by no more than the functional difference ||ʃ -''l, where the specific norm depends on the feature removal approach. These results apply to all subsets S C [d]. so they automatically apply to those induced by feature grouping. We therefore have the following inequality,
If(rr,) f((,,) S 1lf-f'1l f'll.
for all x and all group subsets T C [g]. For the induced vectors v, v' E RR" we also have ||v d' /loc S Ilf - f'll-
Finally, we that Theorem applles as well, where the only change required is replacing d with g the (suummary) term.
40

---

H
LIME weighted least squares linear operator
In this section, we discuss the linear operator A E Rux22 associated with LIME 53 theoretical characterization that focuses on LIME's default implementation choices, and we then corroborate these points with empirical results.
We begin with
H.1 Theory
The observation that LIME is equivalent to linear operator A E IRdx2" is discussed in Appendix D We pointed out there that LIME's specific operator depends several implementation choices, including the choice of weighting kernel, regularization and intercept term. Here, we consider that no regularization used, that an intercept term is fit, and we focus on the default weighting kernel in LIME's popular implementation The default is the following exponential kernel,
11 (1-VIsI/d
w(S) exp
(H.1.1)
σ²
where o² is hyperparameter. The
set this parameter arbitrarily, but LIME also considers default values for each data type: o² = 0.25 for images, o² = 25 for language, and x² = 0.75.d d for tabular data. We find that these choices approach two limiting cases, o² - oc and a² 0, where LIME reduces other summary techniques; namely, we find that these limiting cases respectively recover the Banzhaf value [9] and leave-one-out [67] summary techniques. We argue this first from theoretical perspective, and Appendix H.2 then provides further empirical evidence.
Limiting case o² + oc. In this case, is easy to see from eq. (H.1.1) that we have w(S) - for all S S [d]. The weighted least squares objective for LIME then reduces to unweighted least squares problem, which according to results from Hammer and Holzman [32] recovers the Banzhaf value. We might expect that the default settings of x² approach this outcome in practice (e.g., o² 25), and we verify this in our empirical results.
Limiting case 2 - 0. As we let o² approach O, increasing weight is applied to w(S) when the cardinality [S] is large. Similar to a softmax activation, the differences between cardinalities become exaggerated with low temperatures; for very low temperatures, the weight applied to |S = d is much larger than the weight applied to |S| = d 1, which much larger than the weight for S d etc. In this regime, applying weight only to the subset S [d] leaves the problem undetermined, so the small residual weight on ||] = d plays an important role. The problem begins to resemble one where only subsets with |S| d are taken into account, which yields unique solution: following Covert et al. [17], this reduces LIME to leave-one-out (Occlusion). More formally, if we denote ALIM(( the linear operator resulting from specific o² value, argue that lim ALIME (a²) the operator associated with leave-one-out. We do not prove this result due to the difficulty of deriving AL) for specific o² value, but Appendix H.2 provides empirical evidence for this claim.
H.2 Empirical analysis
Here, we denote the LIME linear operator with the exponential kernel as ALME instead of AL2) for brevity. Recall from Appendix D that the LIME linear operator has the form ALIME = (2 W 2' T W, where Z' = {0, 1)2d x(d+1) x(d+1) is a binary matrix with all ones in the first column and the other columns indicating whether each feature is in each subset, and W E 122x2" diagonal matrix with Ws.s w(S) for all S [d] and zeros elsewhere.
To empirically compare the linear operator of LIME to those of the Banzhaf value and leave-one-out (Occlusion), first derive closed-fon expressions of z' W Z' and Z'T W as in the following proposition.
https://ggthab..co/aaarottttt//
41

---

Proposition 2. Define
d
(-(1- -k/vv//) 2 100).
a0
exp
k)
d
- 1
- (-(1-- vk/aa 102) and
k -1)
exp
k=1
d - 2
(1-vVk/a) 10²).
02
exp
k 2)
k=2
Then (Z" WZ') is (d + 1) x (d + 1) matrix with the following form:
(Z wz')
6)
B
where B = (A1...., a1) E Rixd, and C is d d matrix with a1 on the diagonal entries and a2 the off-diagonal entries. Furthermore, for all [d]
(S)).
fi==1
() is ii E S)w(S)
otherwise.
Proof. We first find an expression for each element in T W Z', Consider the first diagonal entry:
wz') z wzj 12d W124
w(S)
SC
exp (-(1- vk/d)² 102) ag.
k=
For the rest of the diagonal entries, we have:
(z (Z'T 2 =
Σ w(S)
S:-1E8
Σ
exp (-(1- 11--vk/d) 102 = 01
1=1
for 2...,, d For the off-diagonal entries, we have:
(z WZ') (Z' Wz') 1 Σ w(S) a1
S4-1Es
for = 2.... d +1 1. Finally, consider the rest of the off-diagonal entries:
(z WZ'),g Z,TZZZ Wz}-
Σ
w(S)
S:i 1681
Σ
d
d-
2 k カ-2 exp ((1-VE/d)
/02) = 022
for i,j 2..... d + such that + j. Taken all these together, Z'T W Z' has the closed form expression shown in the proposition. We now proceed to find an expression for each element 2T W, which has the first row entries:
(Z' W) = 1.dd Ws w(S),
for all S S [d]. whereas the rest of the rows have entries:
W2'T = z Ws = 1(i 1 S)w(S),
for 2, d + 1. Hence the closed form expression of z"T W in the proposition follows.
42

---

Number of features = 10
Number of features 15
Number of features 20
Number of features 25
33
144
||A Asverall
32-
Tabular
Tabubr
0.5
Taburare
2.5
1.5
0.5
0A
0.1
Figure 14: Frobenius norm of the difference between the LIME linear operator ALIME and the Banzhaf linear operator ABazharf with varying LIME exponential kernel hyperparameter o² and number of features. We include reference lines for the default LIME hyperparameters: o² 0.25 for images, d² 0.75vd for tabular data, and d² 25 for text data.
With closed form expressions for z WZ' and Z' we numerically compute the LIME linear operator ALIME = (
W (1 T' W and compare it to the linear operators of the Banzhaf value (ABanhat) and Occlusion (Aocclusion) with both large and small o² values.
As shown in Figure 14 as a² increases, the Frobenius norm difference between the LIME and Banzhaf linear operators approaches zero, verifying theoretical argument for the limiting case of o² CO. Furthermore, we observe that ALIME and ABanzhaf have near zero difference for all default o² values when there are {10, 15, 20, 25} features. As o² approaches zero, the Frobenius norm difference between the LIME and Occlusion linear operators empirically approaches zero (Figure [5 confirming our intuition for the limiting
where x² -> 0. Based on computational results, LIME with the default の² values indeed has similar spectral norm I|-1|2 and the operator norm ||:11,,0c compared to the Banzhaf value (Figure 16
Number features 10
Number features
Number of features
Number of features
|Aue Aocamoll
20:7
20
13
107
10
10
10
10
107
10.1
10
10
10
107
30
10
15"
10:1
Figure 15: Frobenius norm of the difference between the LIME linear operator ALME and the Occlusion linear operator AOcclusion: with varying LIME parameter o² and number of features.
10
Occlusion
Banzhaf
Shapley
4-
RISE
LIME (Tabular)
[|A]].
3-
(Image)
LIME (Text)
2-
1-
Number of Features
Number of Features
Figure 16: Computed norms of summary technique linear operators. Note that all lines except for RISE overlap on the left-hand-side showing the operator norm.
43

---

I Additional results for synthetic experiments
Baseline
Marginal
Conditional
0.16
Thawedia
Thouettel UisertBoutt
0.14
Drafice Upper
Engriial Uppar Zound
0.12
Occlusion
0.10
0.08
0.06
0.04
0.02
0.00
Banzhaf
RISE
0.10
LIME
0.08
0.08
0.6
000
0.
3.5
2.0
15
25
100
1.
2.0
3.0
0.5
25
25
20
||x x/122
||x x'Ila
x
Figure 17: Occlusion, Banzhaf, RISE, and LIME 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 empirical upper bounds. Theoretical upper bounds are included for baseline and marginal feature removal, and omitted for conditional feature removal because is infinite for p
44

---

0.4
Marginal
Marginal
Slope of Upper Bound
0.8
Theoretical
Theoretical
0.3
(Baseline Marginal)
(Conditional)
Occlusion
0.6
Theuretical
(Baseline Marginal)
0.4
0.0
0.5
1.0
1.5
2.0
0.0
0.2
0.4
0.6
1.0
1.0
0.4
Slope of Upper Bound
0.3
Banzhaf
0.2
0.4
0.1
0.2
0.0
0.0
0.0
0.5
1.5
2.0
0.0
0.2
0.4
0.6
0.8
1.0
1.0
0.20
Slope of Upper Bound
0.8
0.15
0.6
RISE
0.10
0.05
0.0
0.0
0.5
1.0
1.5
2.0
0.0
0.2
0.4
0.6
0.8
1.0
1.0
0.4
Slope of Upper Bound
0.3
0.6
0.2
0.4
0.1
0.2
0.0
0.0
0.5
1.5
2.0
0.2
0.6
1.0
a
P
Figure 18: Slopes of empirical and theoretical upper bounds with respect the input perturbation norm, for Ocelusion, Banzhaf, RISE, and LIME attribution differences in the synthetic experiment. The parameter varies the logistic regression Lipschitz constant, and p varies the correlation between the two synthetic features æ1, Z2.
45

---

1.0
0.8
0.4
0.0
Filp
1.2
Occlusion
1.0
||
0.4
Baseine Morginal Candibional
Baseline Marginal Conditional
Baseline Marginal Conditional
Baseline Marginal Conditiona
1.2
Banzhaf
0.6
||
0.4
Baseline
Marginal Candibional
Baséline
Marginaa Conditional
Baseline Marginal Conditinal
Baseline Manginal Conditional
Fill:
RISE
1| 0
0.6
Baseinee
Marginal Condibiona
Baseline Marginal Conditional
Baseline Marginal Conditonal
Baseline Marginal Conditonnl
1.44
LIME
0.8
0-
0.6-
0.0
Baséline
Marginal Conditional
Baseline
Marginal Conditional
Baseline Marginal Conditinnl
Baseline Manginal Conditinna
Figure 19: Occlusion, Banzhaf, RISE, and LIME attribution differences between the logistic regres- sion classifiers f.f f' as constructed in Appendix A.2 with baseline, marginal, and conditional feature removal, and varying correlation
Implementation details
Fully connected networks. For the wine quality dataset, we trained FCNs with varying levels of weight decay. The architecture consists of hidden layers of width 128 with ReLU activations and was trained with Adam 50 epochs with learning rate 0.001. We also trained a conditional VAE (CVAE) following the procedure in Frye et al. 27 to approximate conditional feature removal. The CVAE encoder and decoder both consist of 2 hidden layers of width 64 with ReLU activations and latent dimension 8. We trained the CVAE using prior regularization strength of 0.5 for at most 500 epochs with Adam and learning rate 0.0001, with early stopping if the validation mean squared error did not improve for 100 consecutive epochs. All models were trained with NVIDIA GeForce RTX 2080 Ti GPUs with 11GB memory.
Convolutional networks. We trained CNNs for MNIST with varying levels of weight decay. Following the LeNet-like architecture Adebayo al. 22. our network has the following structure: Input + Conv(5 x 5. 32) - MaxPool(2 x 2) Conv(5 x 5, 64) - MaxPool(2 2) - Flatten + Linear(1022) 1024) - Linear(1024, 10). We trained the model with Adam for 200 epochs with learning rate 0.001. The same GPUs for FCN training were used to train the CNNs.
46

---

ResNets. We trained ResNet-18 networks for CIFAR-10 from scratch with Adam for 500 epochs with learning rate 0.00005. The same GPUs for FCN training were used to train the ResNet- networks. We used ResNet-50 network pre-trained on ImageNet and fine-tuned it for Imagenette with Adam for epochs with learning rate 0.0005. The ResNet-50 networks were trained with NVIDIA Quadro RTXX000 GPUs with 24GB memory.
K Additional results for practical implications
Wine Quality Test Accuracy
0.95
MNIST Test Accuracy
0.98
0.90
0.96
0.85
0.94
0.80
0.92
0.75
0.90
0.000
0.002
0.004
0.006
0.008
0.010
0.000
0.002
0.004
0.006
0.008
0.010
Weight Decay
Weight Decay
Figure 20: Test accuracy of FCN trained on the wine quality dataset and CNN trained on MNIST with varying degree of weight decay.
0.84-
Imagenette Test Accuracy
0.75
CIFAR-10 Test Accuracy
0.70
0.65
0.80
0.60
0.78
0.50
0.76
0.45
0.40
0.01
0.02
0.07
0.03
0.05
Weight Decay
Weight Decay
Figure 21: Test accuracy of ResNet-18 trained on CIFAR-10 and ResNet-50 fine-tuned on Imagenette with varying degree of weight decay.
47

---

[|x - = 0.1
x x X 0.5
||x x 1.0
| | 2.0
1.2
Baseinc
0.30
Margnal
1.0
Credk
Occlusion
0.8
0.06
0.6
0.2
0.2
0.07
0.255
0.034
ll0fr eGG,**lls
0,200
8.012
0.75
9.68
0.008
8.006
0.56
0.100
0.54
0,075
0.050
0.0000
(MC 660 Mt
0.0175
0.38
0.100
0.25
0.125
0.0125
0.26
0.0000
0.24
0.975
0.0000
0.10
0.92
0.07
0.25
9.012
0.56
9.32
0,200
LIME
3.000
0.25
0.30
9.68
0.125
8.606
3.56
0.100
3.008
0.075
2.0000 5.5075 0,0050 0,0055 6.0003
5.0000 3:025 3,5050 0.8555 8.100
0,0030 0.025 0.055 0.0755 0.6600
93000
0.803
0.0355 8.075 0.050
Weight Decay
Weight Decay
Weight Decay
Weight Decay
Figure 22: Attribution difference for networks trained with increasing weight decay, under input perturbations with perturbation norms (0.1, 0.5, 1,2). The results include the wine quality dataset with FCNs and baseline (replacing with training set minimums), marginal, and conditional feature removal. Error bars show the mean and 95% confidence intervals across explicand-perturaatio perturbation pairs.
48

---

||-x|| = 0.1
||x ||<xx|| 0.5
||<-x||> =11x=====1000 =1.0
| - x = 2.0
1.00
Corrilation
0.45
0.5
0.5
Occlusion
0,7
0.0
0.25
0.97
0.6
0.4
Bise ins Marg Conditonal
0.2
1,000
Correaaion
0.22
0.9
0.598
0.8
8.997
0.90
0.7
0.940-3
0.95
0.890
0.995
0.94
0.825
1.0000
1.00
0.95
0.2298
0.44
0.993
Banzhal
0.3290
Pcarson
0.9985
0.97
0.900
0.2200
0.96
0.859
0.575
0.9945
0.22
0.9990
0.96
0.925
Pearson
0.6985
0.75
0.998
0.96
0.873
0.79
0.65
0.69
1.0000
1.00
0.93
Pearson Corrciation
0.9995
0.93
0.950
0.9900
0.44
0.85
0.83
0.9555
0.97
0.903
0.75
0.850
0.95
0.0500
5.0100
0.0000
Weight Decay
Weight Decay
Decay
Weight Decay
Figure 23: Pearson correlation of attributions for networks trained with increasing weight decay, un- der input perturbations with perturbation norms {0.1, 0.5, 1, The results include the wine quality dataset with FCNs and baseline (replacing with training minimums), marginal, and conditional fea- ture removal. Error bars show the means and 95% confidence intervals across explicand- explicand-perturbation perturbation pairs.
49

---

||x 1 0.1
[x - x 0.5
[ =10
||x x 2.0
0.5
0.35
0.0
0.8
Occlusion
0.90
0.8
0.6
05
0.94
HIHH
Basciin Mangral
0.6
0.3
Conditional
0.92
Rank Corclatioo
0.995
0.93
0.95
0.5
Shapley
0.8
0.990
5.85
0.985
0.92
2.50
8.75
0.998
0.998
Banzhal
0.95
5.00
0.8
0.80
0.
2.55
0.0
0.395
0.994
Rank
0.992
0.94
8.55
0.7
0.990
0.80
0.988
0.992
Correlation
0.
0.996
0.994:
0.95
3.50
0.8
0.92
0.000 0.202 0.044 0.006 9,008 0.010
0.000 0.002 0.004 0:506 0.088 6.000
0.000 0.022 9.04 0:006 0.16 sille
0,000 0.022 0.004 0,005 9.000 0.016
Weight Decay
Weight Decay
Weight Decay
Weight Decay
Figure 24: Spearman rank correlation of attributions for networks trained with increasing weight decay, under input perturbations with perturbation norms {0.1, 0.5, 1, 2). The results include the wine quality dataset with FCNs and baseline (replacing with training set minimums), marginal, and conditional feature removal. Error bars show the means and 95% confidence intervals across explicand-perturaation perturbation pairs.
50

---

- x = 0.1
||x - x = 0.5
1 - x = 1.0
| xx 2.0 base tie (Mian) Base the (2er))
0.23
3.50
Occlusion
0.003
0.013
0.03
0.01
0.02
3.24
9.02
0.06
0.05
0.07
9.20
0.09
Banzhaf
0.05
0.05
10.0
3.06
0.04
9.05
0.01
0.04
0.02
0,066
0.065
0016
0.800
0.20
0,017
0.055
0.040
0.08
0.0255
0.250
0.00
0.006
0.0000
0.045
9.67
0.00
0.0355
0.02
0.033
0.043
0.25
8,6250
0.090
0.040
6,6225
0.023
0.535
9.06
5,0000
0.030
0,0175
0,025
3.24
8,0125
0.015
0,020
9,0100
0,0025 00100
0,0020 0.0025
8:0055 Weight Decay
0.050
Decay
Weight Decay
Weight: Decay
Figure 25: Attribution difference for networks trained with increasing weight decay, under input perturbations with perturbation norms {0.1, 0.5, 1, The results include MNIST with CNNs and baseline feature removal with either training set means or zeros. Error bars show the means and 95% confidence intervals across
pairs.
51

---

X 0.1
[ |x-x|| = 0.5
[ -x1b = 1.0
[ |
10203
1.009
3.00
1.0
Pearson Correlation
0.59
3.6
Orclusion
0.903
238
9.54
257
0.0
0.9385
538
3.8
0.203
0.9389
5.5
0.66
(Mean) (eer)
0.943
0.259
Pearson Carrela
0.985
0.955
0.955
0.30
Shapley
0.983
0.955
0.915
0.15
0.910
0.39
0.50
0.955
0.855
0.055
0.998
8.995
0.908
0.950
8.991
0.224
0.992
Pearson
8.992
0.522
0.990
0.88
0.950
8.990-
0.88
0.05
0.56
0.91
9.54
0.92
3.92
2.02
RUSE
0.00
5.50
0.86
0.88
2.84
Pearson Correlatin
0.993
0.329
0.092
0.950
0.955
0.944
8.945
0.944
0.900
0.975
0.993
0.0000 Decay
0.0030 o.oors Weight Decay
Weight Decay
Weight Decay
Figure 26: Pearson correlation of attributions for networks trained with increasing weight decay, under input perturbations with perturbation norms {0.1, 1,2]. The results include MNIST with CNNs and baseline feature removal with either training set means zeros. Error bars show the means and 95% confidence intervals across explicand-perturbation pairs.
52

---

1x 0.1
|K - = 0.5
||x 1 1.0
[ - =200
10
0.0
0.0
Corclatiion
0.
0.7
Occlusion
0.6
93
0.4
0.4
0.3
Boseline
Baceline (Zem)
0.2
0.80
Corretation
0.75
0.75
0.75
0.75
Shapley
0.70
0.00
0.02
0.55
0.55
0.65
0.85
0.05
Corrlaatio
0.80
0.00
0.82
Banzhal
0.75
0.70
0.75
0.05
0.65
0.00
0.60
0.725
0.723
0.725
0.725
3.700
0.700
Rank Corrdation
0.500
0.90
Correlttio
0.37
0.97
0.96
0.42
LIME
0.36
0.02
0.34
0.34
0.04
0.38
0,53
0.90
0.010
0.019
0,000 Weight Decay
Weight Decay
Weight Decay
Weight Decay
Figure 27: Spearman rank correlation of attributions for networks trained with increasing weight decay, under input perturbations with perturbation norms (0.1, 0.5, 1, 2). The results include MNIST with CNNs and baseline feature removal with either training set means or zeros. Error bars show the means and 95% confidence intervals across peglicand-peruurbationn pairs.
53

---

1.0
1.0
Pearson Correlation
0.8
Rank Correlation
0.6
0.6
0.4
0.4
Original
Linear 1
Linear
Conv
Conv 0
Original
Linear
Linear
Conv
Conv 0
Occlusion (Zero Baseline) Banzhaf (Zero Baseline)
Shapley (Zero Baseline)
LIME (Zero Baseline)
IntGrad
RISE (Zero Baseline)
Vanilla Grad
Grad x Input
Sanity checks for attributions using cascading randomization for the CNN trained on MNIST. Features are removed by replacing them with zeros. Autribution similarity is measured by Pearson correlation and Spearman rank correlation. We show the means and 95% confidence intervals
Figure 28:
across 10 random seeds.
1.0
1.0
Pearson Correlation
Rank Correlation
0.8
0.6
0.4
0.2
0.0
0.0
Original-
Linear
Layer4,1 1. Conv2.
Layer4 1......Co. Convi
Layer4. Conv2.
Layer4. 0. c.novl
Layer3,1 Convv
Layer3.1.Com/ Convi
Layer3, Layer3.0.convz 0.Co22.
Layer3 C.Cavvl
Layer2.11 Conv2-
Layer2.1 1...1.... Convi
Layer2.0 0.0.Co2
Layer2.0 C.Covll
Layer1. 1. Conv2.
Layer1 w.w.222
Layer1 1.c.1..ccov
Convl
Original
Linear
Layer4,1 Conv2.
Layer Conovvl
Layer4.0.Com/ 0.n.C22
Layer4.0m ..o.4.
Layer3.1.Come Cono2
Layer3.1 1. Convi
Layer3. Cano22
Layer3 C.Covvl
Layer2.1.Com Conov
Layer2.1.Coms 1.Covvl
Layer2. 0 Conv2.
Layer2 0 Convl-
Layer1.1 Conv2.
Layer1. Conovl
Layer 0.c.vo2
Layer1.0.Com convv..
Convl.
Layer1.
Occlusion (Zero Baseline)
RISE (Zero Baseline)
IntGrad
Grad Input
Shapley (Zero Baseline)
Vanilla Grad
Figure 29: Sanity checks for attributions using cascading randomization for the ResNet-18 trained on CIFAR-10. Features are removed by replacing them with zeros. 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 undefined correlations due to constant attributions under model randomization.
54


Reviews:
[{"review_id": "DopUx9tRXZ", "reviewer": "Reviewer_aeSn", "summary": "Previous research has primarily focused on the robustness of gradient-based feature attributions, but the robustness properties of removal-based attribution methods are not well understood. To fill this gap, the authors of the paragraph aim to theoretically analyze and characterize the robustness of removal-based feature attributions. They provide a unified analysis of such methods and prove upper bounds for the difference between intact and perturbed attributions under various settings of input and model perturbation. The authors validate their theoretical findings through empirical experiments on synthetic and real-world data, showcasing their practical implications.", "questions": "Please refer to the Weakness.", "limitations": "The authors have adequately addressed the limitations of this work in the last section.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. This paper addresses a crucial XAI problem: Interpretation robustness concerning both input and model perturbations. The study may benefit the derivation of the theoretical impact of some important tricks in XAI, such as baseline value selection and marginal distribution approximation.\n\n2. The derivation of robustness considers input perturbation and model perturbation comprehensively.\n\n3. Section 4 provides valuable insights, revealing that interpretation robustness relies on both the Lipschitz continuity constant and feature number for input perturbation, while it depends on the ∞-norm and feature number in the case of model perturbation.\n\n4. The experiment observations in Section 5 are also enlightening, highlighting the advantages of removal-based attribution over gradient-based attribution.", "weaknesses": "The paper should be re-arranged by including experiment results in the main content and moving some theoretical results to the appendix.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> “The paper should be re-arranged by including experiment results in the main content and moving some theoretical results to the appendix.”\n\nWe thank the reviewer for this suggestion to improve our paper. To address this suggestion, we will move Corollary 1 to the Appendix to make room for a longer discussion of the experiments. Also, with the extra page allowed in the camera-ready version, we can move Figures 4 & 5 to the main text to highlight the empirical implications of our theory. While we agree that it would be ideal to include our full experiments section in the main text, we believe that a focus on theoretical insights along a concise description of experiments will allow the average reader to best understand the findings from our work. "}, {"review_id": "3mcnBOvo0H", "reviewer": "Reviewer_isMK", "summary": "This paper studies the robustness properties of an explanation to small perturbations in the input space (i.e. like an adversarial example) and also to the model parameters. The authors use a number of Lipschitz-style bounds to then derive overall limits on how much explanations can change. \n\nUpdate: As the reviewers have addressed the main concern on the Lipschitz constant, I have no major reason to reject this paper.", "questions": "Can you measure or estimate the Lipschitz constants and show that these are meaningfully non-vacuous assumptions? \n\nIt would be better to highlight some kind of experimental result in the main paper. \n\nIs there a reason or limiting factor preventing the work from applying to higher dimensional work such as text or image settings?", "limitations": "The limitations discuss some conservativeness of the bounds, which could be addressed with work on certified robustness.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The bounds are principled and rigorous for explanations, as opposed to heuristic. \n\nThe paper considers multiple removal techniques and summary statistics to show how the bounds change based on the methods. \n\nThese bounds can sometimes lead to significant asymptotic differences.", "weaknesses": "The elephant in the room is that the results all need some kind of Lipschitz or Lipschitz-like bound, and that the final bound depends on the constant. This constant could be significant and not necessarily ignored, but there does not seem to be any evidence that the actual Lipschitz constant is at all close to being small enough to be useful. \n\nThere are barely any experimental results in the main paper. Almost the entirety is deferred to the appendix. \n\nThe results seem to focus mainly on synthetic settings and linear settings where Lipschitz constants and other constants can be directly computed. They also focus on fairly simple datasets such as UCI wine and MNIST.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> “The elephant in the room is that the results all need some kind of Lipschitz or Lipschitz-like bound […] there does not seem to be any evidence that the actual Lipschitz constant is at all close to being small enough to be useful.”\n\nFirst, it’s worth emphasizing that our work aims to understand these attribution methods, not necessarily advocate for them, so it’s not a weakness that our analysis highlights the dependence on a quantity that is difficult to measure. In settings where $L$ can be large, our theory explains their lack of robustness; however, more positively, it also highlights the opportunity to ensure robustness via methods for controlling $L$, such as Lipschitz-constrained networks. Ultimately, this doesn't seem like a reason to reject our work, as we have correctly characterized this family of attribution methods. Note also that earlier works on gradient-based methods like Dombrowski et al. [R1] are based on properties of the network’s Hessian, which are more difficult to measure than the Lipschitz constant.\n\nNonetheless, we can also provide a brief description of the current literature on estimating $L$. Exact computation of the Lipschitz constant, even for a two-layer neural network, is NP-hard [R2]. Even for simple networks on MNIST, current methods for estimating their Lipschitz constants give large estimates [R2-R4], which can be conservative and obscure meaningful conclusions with our theoretical results. That said, our theoretical analysis suggests empirical techniques for improving the robustness (e.g., increasing weight decay to shrink an upper bound on the Lipschitz constant), and these show the intended effect in our experiments. Therefore, our theoretical results regarding Lipschitz continuity are indeed useful in practice.\n\n> “It would be better to highlight some kind of experimental result in the main paper.”\n\nWe thank the reviewer for this suggestion to improve our paper. To address this suggestion, we propose to move Corollary 1 to the Appendix to make room for a longer discussion of the experiments. Also, with the extra page allowed in the camera-ready version, we can move Figures 4 & 5 to the main text to highlight the empirical implications of our theory. We believe that a focus on theoretical insights and a concise description of experiments in the main text will allow the average reader to best understand the contributions and findings from our paper. \n\n> “Is there a reason or limiting factor preventing the work from applying to higher dimensional work such as text or image settings?”\n\nWe thank the reviewer for this suggestion to improve the impact of our paper. To address this point, we ran additional experiments with (i) ResNet-18 trained on CIFAR-10, and (ii) ResNet-50 trained on 10 classes from ImageNet (i.e., Imagenette). The figures are included in the attached file in the “global” response. We again observe empirical results consistent with our theoretical insights. In particular, ResNet-18 and ResNet-50 trained with stronger weight decay have more robust attributions under input perturbations. Also, under cascading model randomization (a form of increasing model perturbation) [R5], the unperturbed and perturbed ResNets indeed have increasingly dissimilar attributions (and as in our previous experiments, more so for removal-based methods than gradient-based ones).\n\n"}, {"review_id": "KcnaQq0PqA", "reviewer": "Reviewer_2wCy", "summary": "This paper theoretically analyzes the robustness of removal-based feature attributions against input perturbation, and model perturbation with different summary techniques. Empirical experiments on synthetic datasets support their theoretical analyses such as conditional sampling is more robust to model perturbations compared to baseline or marginal samplings. Experiments on real-world datasets (UCI wine quality and MNIST) give more insights regarding the robustness under different model training settings and comparison to gradient-based explanation methods.", "questions": "The overall robustness proposed at the end of Section 4, which combines theorems 1 and 2, seems to be synthesized and is hard to adapt into a realistic scenario. For instance, the input perturbation robustness measures the consistency of explanations generated by one model, while model perturbation robustness analyzes the consistency given different (similar) target models. If use the overall robustness, it does not provide a clear message about different removal-based attributions.", "limitations": "Broader impact of the proposed framework should be discussed. More examples to use proposed robustness would be necessary to make the contribution stronger.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "(1) The proposed theoretical analysis of removal-based feature attributions is technically sound. The input perturbation and model perturbation together with the summary technique cover the different aspects of removal-based attributions. \n\n(2) As there are more explanation techniques proposed, this work gives a good starting point to analyze the robustness limitations from a theoretical perspective, which can enable solid explanation evaluation on robustness among different explanation methods.\n\n(3) Besides providing theory proofs, empirical experiments also validate and support the findings from theoretical analyses. \n\n(4) Messages from the proposed theories are clear. The paper is well written and presented.", "weaknesses": "(1) The analyses are only limited to several removal-based feature attribution methods. Extending the analyses to other explanation techniques such as gradient-based explanations would make the contribution stronger. In fact, gradient-based explanations are more popularly used. \n\n(2) The project is inspiring and gives insights into different removal-based attributions. However, the impact of the analyses can be broader if the authors can propose a technically sound robustness evaluation benchmark based on the theory. \n\n(3) To analyze the robustness is computationally costly if the input data is high-dimensional or the model is huge. The author should consider extending the proposed analyses on complex datasets and some experimental results would make the paper stronger.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> “Extending the analyses to other explanation techniques such as gradient-based explanations would make the contribution stronger. In fact, gradient-based explanations are more popularly used.”\n\nOther works have focused on the robustness of gradient-based methods [R1, R6], but the feature attribution literature lacks comprehensive theoretical analyses on the robustness of removal-based methods, which are also highly popular (e.g., the SHAP GitHub repository has ~19.9k stars, and the LIME GitHub repository has ~10.8k stars). Our paper makes a distinct contribution by addressing this gap. In the Discussion, we draw connections between the robustness of removal- and gradient-based methods. Primarily, the robustness of removal-based methods under input perturbations is related to a model’s Lipschitz continuity, whereas robustness of gradient-based methods is often related to the Lipschitzness of model gradient (i.e., Lipschitz smoothness). It would be interesting to pursue a unified, theoretical analysis of the robustness of both classes of methods, but we leave this topic to future work.\n\n> “The impact of the analyses can be broader if the authors can propose a technically sound robustness evaluation benchmark based on the theory.”\n\nThe goal of a robustness benchmark is to provide an objective comparison between attribution methods, perhaps in the context of a specific model and dataset, but ideally with more generalizable conclusions. Several existing works address this topic empirically, and our work intends to provide a complementary theoretical perspective. We represent the essential properties of the model and dataset through our assumptions (namely $L$, $B$ and $M$), and this allows us to compare removal-based explanations in a more generalized sense by focusing on their implementation choices.\n\nThat said, there is also a clear connection between our analysis and possible empirical benchmarks, and this can be seen in our experiments. For robustness to input perturbations, we verify that our theoretical findings are demonstrated in practice (Theorem 1) by testing the stability of attributions when adding noise to the input (e.g., Figure 4). These experiments resemble the empirical analysis conducted by Ghorbani et al. for gradient-based methods [R7], and they are also similar to metrics used in an existing benchmark library [R8]. In other words, our theory provides a way to characterize properties that other works study empirically.\n\nRegarding model perturbation, two types of benchmarks could include (i) making small perturbations to verify that the attributions don’t change much, and (ii) making large perturbations to verify that the attributions change significantly. In terms of existing works, Anders et al. explore the first direction [R6], and Adebayo et al. explore the latter [R5]. We incorporated the sanity checks from Adebayo et al. into our experiments (Figures 5-6), and Figure 2 includes an experiment with logistic regression similar to (i). There is certainly room to create new benchmarks along these lines, but a difficult design choice is how to define the small perturbations. Giving this choice proper consideration and generating results across a wide class of methods is best left as a topic for future work. Again, it is perhaps best to view our work as studying theoretically what other works have studied empirically, including not only [R5] and [R6] but also Slack et al. [R9].\n\n> “To analyze the robustness is computationally costly if the input data is high-dimensional or the model is huge. The author should consider extending the proposed analyses on complex datasets and some experimental results would make the paper stronger.”\n\nWe thank the reviewer for this suggestion to improve the impact of our paper. We ran additional experiments with (i) ResNet-18 trained on CIFAR-10, and (ii) ResNet-50 trained on 10 classes from ImageNet (i.e., Imagenette). The figures are included in the attached file in the “global” response. We again observe empirical results consistent with our theoretical insights. In particular, ResNet-18 and ResNet-50 trained with stronger weight decay have more robust attributions under input perturbations. Also, under cascading model randomization (a form of increasing model perturbation) [R5], the unperturbed and perturbed ResNets indeed have increasingly dissimilar attributions (and as in our previous experiments, more so for removal-based methods than gradient-based ones).\n\n> “The overall robustness proposed at the end of Section 4, which combines theorems 1 and 2, seems to be synthesized and is hard to adapt into a realistic scenario.”\n\nThe overall robustness in Corollary 1 is included for the sake of completeness and to show that our theory can account for this case. The result can be useful in situations where we want to interpret a model $f$ on a sample $x$, but the system is subject to simultaneous perturbations, such that we only have access to a perturbed model $f’$ and perturbed sample $x’$ (e.g., hiding a black-box model’s undesirable dependency on certain features with simultaneous, small model and input perturbations). However, we agree that such simultaneous perturbations may be uncommon in practice, so Corollary 1 can be deferred to the Appendix (thus making more room for our experiments).\n"}, {"review_id": "HwFDpj5HOT", "reviewer": "Reviewer_dzDD", "summary": "This paper derives robustness results regarding a general class of feature attribution methods referred to as “removal-based feature attributions”, which includes occlusion methods, but also Shapley values and LIME explanations. The authors study Lipschitz-continuity properties of these methods with respect to changes in the modeled function and the inputs. They find that, under Lipschitz assumptions on the classifiers output, the removal-based attributions are also Lipschitz-continuous w.r.t to deviation in the inputs and the model.", "questions": "I would appreciate if the authors could elaborate on the following points:\n\n* How can the results be operationalized in practice? Can something constructive be derived from them?\n* How would authors quantify the degree of robustness that is desirable (e.g., not invariant, not noisy)\n* Do the results bear any connections to the benchmarking of feature attributions with is frequently done through removals (e.g., Tomsett et al., 2020; Rong et al, 2022)? \n\nA discussion of these points in the paper could help to strengthen this work and I would be willing to reconsider my overall assessment based on the authors' replies.\n", "limitations": "I think the technical limitations are sufficiently clear, and mostly concern the Lipschitz-continuity assumptions and the bounded domains.", "rating": 6, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 2, "strengths": "**Generality.** The class of feature attributions studied in this work is kept quite general and it thus covers many relevant implementations such as Occlusion, LIME, and Shapley values.\n\n**Technical soundness.** I checked the proofs of some main results and was not able to identify any flaws This makes this work is a rigorous technical contribution in my view. \n\n**Good exposition and presentation.** The formalization and notation were introduced well, and the results were good to follow. The proofs that I checked were sufficiently detailed to convince me of the correctness of the theoretical claims. The tables help to get an overview over the various results.\n \n**Interesting results on robust summarization techniques.** I particularly like the results on the robustness of the summarization techniques and the statements of the most robust aggregation techniques under each set of axioms (Lemma 8 - Lemma 11). The do not rely on any assumptions and can be operationalized directly, for instance for robust data valuation as done in the recent work mentioned (Wang & Jia, 2023). \n\n", "weaknesses": "**How robust should it be?** While the results are interesting, they are hard to interpret, because there is not optimal robustness. As the authors correctly state, too small robustness can be considered as invariance as in the sense of Adebayo et al. On the contrary, low robustness allows for adversarial attributions. Therefore, the computation of the robustness score allows for no clear interpretation.\n\n**Given the formalization done in prior work, the main results are not unexpected.** First, I would like to underline that central parts of the formalization in this work were transferred from Covert et al., 2021. For instance, the distinction between the removal strategy and the summarization technique and the representation of the methods in this scheme was introduced in this work. I think this should be pronounced more clearly. Further, we assume the model to be Lipschitz-continuous (Assumption 1), the probability divergence when changing the set of present features to be Lipschitz as well (Assumption 3), and we consider only attributions that can be represented as linear operator of the model outputs, it is not very surprising that we arrive at Lipschitz-continuity of the attributions w.r.t. to inputs and the model. While I see that intuitive results also need to be rigorously proven, I am uncertain about how many new insights are added. \n\n**Implications and Experiments.** The practical implications of the results still remain unclear to some extent. Taking a practitioner’s perspective, I wonder about what novel insights can be gained from this work. The finding that regularized networks have more robust (or less noisy) attributions, has already been established (e.g., by Shah et al, 2021, Dombrowski et al.). On a sidenote, the Lipschitz-continuity does not imply that the method needs to pass the parameter randomization test, as a constant function would still be Lipschitz-continuous. I would be more interested in some constructive way to operationalize the findings, for instance to prevent “adversarial” explanations as hinted in the text? Can they help us construct better techniques? I think these questions are interesting and to this end, to me, it is seems unfortunate that the experiments already conducted and some further experimentation towards this direction did not make it into the main paper.\n\n**Minor Points**\n\nA line of related work uses removal strategies to benchmark feature attributions (e.g., Tomsett et al., 2020; Rong et al., 2022), which could be discussed as potential a related work. I wonder if the results in this work bear some connection with attribution evaluation problem, which could be a possible application.\n\nTable 2 seems not to be referenced in the text and thus appears a bit context-less.\n\n--------------------------------------\n\n**Summary.**\nA rigorously executed paper with general theoretic results on the robustness of attribution methods. However, the practical implications remain a bit unclear and there are no experiments in the main paper. Overall, the paper seems just above the bar of acceptance to me. I would be willing to further increase my score if the authors can convince me of the practical relevance and impact of their results.\n\n\n**References**\n\nIan C Covert, Scott Lundberg, and Su-In Lee. Explaining by removing: A uniﬁed framework for model explanation. The Journal of Machine Learning Research, 22(1):9477–9566, 2021.\n\nShah, Harshay, Prateek Jain, and Praneeth Netrapalli. \"Do input gradients highlight discriminative features?\" Advances in Neural Information Processing Systems 34 (2021): 2046-2059.\n\nTomsett, Richard, et al. \"Sanity checks for saliency metrics.\" Proceedings of the AAAI conference on artificial intelligence. Vol. 34. No. 04. 2020.\n\nRong, Yao, et al. A Consistent and Efficient Evaluation Strategy for Attribution Methods. In International Conference on Machine Learning (pp. 18770-18795). PMLR, 2022\n\nWang, Jiachen T., and Ruoxi Jia. \"Data banzhaf: A robust data valuation framework for machine learning.\" International Conference on Artificial Intelligence and Statistics. PMLR, 2023.\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> “Given the formalization done in prior work, the main results are not unexpected.”\n\nRegarding the “removal-based explanation” formalization from Covert et al., 2021, we do not intend to frame this as a contribution, hence its exposition in the Background (Section 2.2). Regarding the assumptions required for our analysis and whether the results are surprising, our goal is to understand these methods, so it is not a weakness in our analysis to have identified these necessary conditions. Viewed differently, these assumptions perhaps reflect how difficult it is to achieve robust attributions. To highlight several insights that are worthy and perhaps “surprising” technical contributions, these include: identifying the role of Lipschitz-like continuity in the conditional distribution and its role in making the conditional removal approach less robust; bridging two notions of implementation invariance (weak and strong) through the lens of robustness to model perturbation; and deriving the most robust summary techniques within classes of game-theoretic options.\n\n> “The finding that regularized networks have more robust (or less noisy) attributions, has already been established (e.g., by Shah et al., 2021, Dombrowski et al.).”\n\nThe novelty lies in the different reasons why network regularization works, which are implied by our theoretical analysis, and which should be a topic of sufficient interest for readers. The previous works demonstrating this effect specifically study gradient-based methods [R1, R10], where weight decay plays a role in bounding a network’s Hessian. Here, we use weight decay because of its role in bounding the network’s Lipschitz constant. Furthermore, the role of Lipschitz continuity in our analysis implies that we can use other techniques to improve the robustness: beyond weight decay regularization, we can also use recent techniques for Lipschitz-constrained networks [R11-R13]. To help readers appreciate these points, we will be sure to clarify them in our revised Discussion.\n\n> “The Lipschitz continuity does not imply that the method needs to pass the parameter randomization test.”\n\nThe parameter randomization test is a form of model perturbation, and Lipschitz continuity has no role in our analysis of robustness to model perturbations.\n\n> “How can the results be operationalized in practice? Can something constructive be derived from them?”\n\nFrom a practical perspective, our theory suggests (i) how to choose an attribution with strong robustness, and (ii) how to train a model that achieves more robust explanations. Regarding the latter point and specifically robustness to input perturbations, by largely reducing this problem to the model's Lipschitz constant, our work shows how to leverage a range of solutions developed in the robust learning literature, such as regularizing or constraining the network’s Lipschitz constant.\n\nThere has been an ongoing debate in the literature about whether features should be removed using baseline values/marginal distribution or conditional distribution, as discussed and summarized in Chen et al. [R14]. Our analysis provides constructive guidelines for this choice from a robustness perspective: that is, if input perturbation is a concern, then baseline or marginal removal provides superior robustness. However, if model perturbation is a concern (e.g., “fairwashing” a model’s dependency on features that encode social biases), then conditional removal may be a better choice. \n\n> “How would authors quantify the degree of robustness that is desirable (e.g., not invariant, not noisy)?”\n\nThis is an interesting question that we think is somewhat open-ended in terms of how to address it. The perspective we explored is how to choose an attribution that’s as robust as possible while maintaining a specific sense of being meaningful, and our results in Section 3.3 provide several options for how \"meaningful\" can be defined (see Definition 2). Briefly, we found that different combinations of constraints on the attribution lead to different optimal choices, for example the Banzhaf or Shapley value, which are both known in the literature. However, we acknowledge that there could be other definitions for “meaningful” explanations, and exploring alternative perspectives on this question is an interesting topic for future work.\n\n\n> “Do the results bear any connections to the benchmarking of feature attributions which is frequently done through removals (e.g. Tomsett et al. 2020; Rong et al., 2022)?”\n\nThank you for raising this question, we are happy to cite several papers from that related line of work. After carefully considering the subject, we do not think there’s an interesting connection between these research directions; instead, we think they are best viewed as distinct, yet important notions of good behavior from an attribution method. And in fact, they may represent competing notions of good behavior, because as we allude to in Section 3.3, perfect robustness can only be achieved by designing an uninformative attribution method (which would necessarily perform badly on any removal/ablation-based metric). Thus, future work on new attribution methods should be sure to account for their performance from the perspective of both robustness and ablation-based metrics.\n\n> “Table 2 seems not to be referenced in the text and thus appears a bit context-less.”\n\nThank you for reading our paper carefully. We will include a short comment after Lemmas 5 & 6 to reference Table 2."}]

Author rebuttal:
We thank all the reviewers for their insightful and constructive comments, which have helped us further improve our paper. The specific issues raised by each reviewer are addressed in the individual responses below. We hope you will consider raising the scores if we have adequately addressed your comments. References relevant to all responses are included below.

References

[R1] Towards robust explanations for deep neural networks - Dombrowski et al. 2022

[R2] Lipschitz regularity of deep neural networks: analysis and efficient estimation - Scaman and Virmaux 2018

[R3] Efficient and Accurate Estimation of Lipschitz Constants for Deep Neural Networks - Fazlyab et al. 2019 

[R4] Efficiently Computing Local Lipschitz Constants of Neural Networks via Bound Propagation - Shi et al. 2022

[R5] Sanity Checks for Saliency Maps - Adebayo et al. 2018

[R6] Fairwashing Explanations with Off-Manifold Detergent - Anders et al. 2020

[R7] Interpretation of Neural Networks is Fragile - Ghorbani and Abid et al. 2019

[R8] OpenXAI: Towards a Transparent Evaluation of Post hoc Model Explanations - Agarwal et al. 2022

[R9] Fooling LIME and SHAP: Adversarial Attacks on Post hoc Explanation Methods - Slack and Hilgard et al. 2020

[R10] Do Input Gradients Highlight Discriminative Features? - Shah et al. 2021

[R11] Regularisation of Neural Networks by Enforcing Lipschitz Continuity - Gouk et al. 2018

[R12] The Singular Values of Convolutional Layers - Sedghi et al. 2018

[R13] Lipschitz-Margin Training: Scalable Certification of Perturbation Invariance for Deep Neural Networks - Tsuzuku et al. 2018

[R14] True to the Model or True to the Data? - Chen and Janizek et al. 2020
 

Meta-review: This paper presents a theoretical analysis of the robustness of removal-based feature attributions, deriving bounds for how much attributions may change under input or model perturbations. The bounds are based on Lipschitz continuity of the model, with factors that depend on the type of removal distribution and summary method. This is primarily a theory paper, with the main contribution being the bounds given in Section 4, as well as some discussion of the implications of these bounds and how to improve them.

Reviews were positive, agreeing that robustness of attributions are an important problem in XAI and that the paper makes a solid theoretical contribution. Weaknesses include the limited experimental results which are almost entirely in the appendix, although reviewers appreciated these and the authors have said they will move some to the main paper for camera-ready. One reviewer also expressed concern that it may be difficult to obtain a Lipschitz constant low enough to provide a useful bound; however, this concern was allayed by the discussion and the score was raised significantly.
Decision: Accept (poster)