{'fraction_non_alphanumeric': 0.04030799991491715, 'fraction_numerical': 0.012294471741858635, 'mean_word_length': 4.460394889663182, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 4, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Making Sense of Dependence: Efficient Black-box Explanations Using Dependence Measure

This paper presents a new efficient black-box attribution method based on Hilbert-Schmidt Independence Criterion (HSIC), a dependence measure based on Reproducing Kernel Hilbert Spaces (RKHS). HSIC measures the dependence between regions of an input image and the output of a model based on kernel embeddings of distributions. It thus provides explanations enriched by RKHS representation capabilities. HSIC can be estimated very efficiently, significantly reducing the computational cost compared to other black-box attribution methods. Our experiments show that HSIC is up to 8 times faster than the previous best black-box attribution methods while being as faithful. Indeed, we improve or match the state-of-the-art of both black-box and white-box attribution methods for several fidelity metrics on Imagenet with various recent model architectures. Importantly, we show that these advances can be transposed to efficiently and faithfully explain object detection models such as YOLOv4. Finally, we extend the traditional attribution methods by proposing a new kernel enabling an ANOVA-like orthogonal decomposition of importance scores based on HSIC, allowing us to evaluate not only the importance of each image patch but also the importance of their pairwise interactions. Our implementation is available at https://github.com/paulnovello/HSIC-Attribution-Method.


Introduction
Artificial Intelligence has established itself as the reference technique for tackling many real-world automation tasks. Consequently, the diversity of its applications is growing and reaching fields where its outputs can contribute to critical decision-making. In such cases, it is essential to be able to provide explanations for each link of the decision chain, including AI algorithms. Over the past decade, many techniques have emerged to explain the predictions of these algorithms [48,40,45,17,58,36,35,30,29,42], marking the birth of a new field called Explainable Artificial Intelligence (XAI). The tools developed in this research field, mostly designed to explain neural networks, have already proven helpful. For instance, it has been used in model debugging, identification of new development strategies for practitioners, and failure understanding.
Initial approaches are based on analyzing the internal state of neural networks during inference, often relying on input gradients or activation values of hidden layers [48,45,17,28]. However, the gradient only reflects the model's operation in an infinitesimal neighborhood around an input and can therefore be misleading [20]. Furthermore, their applicability is limited to the case where the final user has access to the implementation of the model. Therefore, such methods cannot be applied in the most common use cases, e.g. when models are made available by third parties through 36th Conference on Neural Information Processing Systems (NeurIPS 2022). Figure 1: HSIC Explainability method. We sample random binary masks M that we use to perturb the input image X. We obtain a perturbed output and measure the dependence between the distribution of each patch P M of the binary mask and that of the output P y . We use a dependence measure, Hilbert Schmidt Independence Criterion (HSIC), based on the kernel embedding of this distribution in a Reproducing Kernel Hilbert Space (RKHS). Each patch is then assigned the value of this measure: the more independent a patch is from y, the less important it is to explain it.
API calls or specialized hardware. In order to address this issue, some black-box approaches have been recently proposed, relying on the perturbation of the input and the observation of its effect on the output [61,40,36,13]. One challenge of such perturbation methods is assessing this effect together with taking into account complex interactions inherent to deep neural networks. To account for these characteristics, black-box methods resort to complex Monte Carlo methods that require a high number of model forward passes, which can be expensive for recent neural networks that are growing larger.
In [13,32], the authors propose methods to reduce the required number of forward passes, but the obtained performance improvements still do not make them close to white-box methods. In this work, similarly to [13], we cast perturbation studies as Global Sensitivity Analysis (GSA) [38]. However, we rely on a whole different approach based on dependence measure rather than analysis of variance. We measure the dependence between patch-wise perturbations of an input image and the model's output by comparing the distribution of perturbed inputs and outputs embedded in a Reproducing Kernel Hilbert Space (RKHS). More specifically, we use Hilbert-Schmidt Independence Criterion (HSIC), a dependence measure based on the Hilbert-Schmidt norm of the empirical cross-covariance operator evaluated between the represented distribution. HSIC leverages the rich theory of RKHS, thereby capturing more diverse information than variance-based indices such as Sobol. In addition, it can be estimated more efficiently, even bridging the performance gap between black-box and white-box methods.
Our contributions are as follows: (1) we introduce a new efficient black-box attribution method relying on HSIC; (2) we derive a new kernel that confers an Analysis of Variance (ANOVA)-like orthogonal decomposition property, allowing us to go beyond usual attribution methods and evaluate interactions between patches of the image; (3) we conduct experiments to assess the fidelity of our method on ImageNet and show that it improves or matches the state-of-the art for different metrics while bridging the computational gap between black-box and white-box attribution methods; (4) we demonstrate its versatility and its potential by successfully applying it to a less common test case: explanations of object detection; and a new test case: evaluation of pairwise interactions between patches of the input image.

Related work
Our work builds on prior efforts aiming to develop attribution methods in order to explain the prediction of a deep neural network by pointing to input variables that support the predictiontypically pixels or groups of pixels, i.e. patches in the image -which lead to importance maps.
Attribution methods for white-box models A large number of attribution methods have been developed relying on the gradient of the decision studied. The first method was introduced in [4] and improved in [48,61,56] and consists of explaining the decisions of a convolution model by back-propagating the gradient from the output to the input, indicating which pixels affect the decision score the most. However, this family of methods is limited because they focus on the influence of individual pixels in an infinitesimal neighborhood around the input image in the image space. For instance, it has been shown that gradients often vanish when the prediction score to be explained is near the maximum value [58]. Integrated Gradient [58] and SmoothGrad [50] partially address this issue by accumulating gradients. Another family of attribution methods relies on the neural network's activations. Popular examples include CAM [62], which computes an attribution score based on a weighted sum of feature channel activities -right before the classification layer. GradCAM [45] extends CAM via the use of gradients, reweighting each feature channel to take into account their importance for the predicted class. Nevertheless, the choice of the layer has a huge impact on the quality of the explanation. In comparison, our proposed approach is model-agnostic and hence does not require access to internal computations.
Attribution methods for black-box models In this paper, we extend the problem by restricting it to a black-box model: the analytical form and potential internal states of the model are unknown. Several methods compute influence scores for each individual pixel or group of pixels.
The first method, Occlusion [61], masks individual image regions -one at a time -with an occluding mask set to a baseline value and assigns the corresponding prediction scores to all pixels within the occluded region. Then the explanation is given by these prediction scores and can be easily interpreted. However, occlusion fails to account for the joint (higher-order) interactions between multiple image regions. For instance, occluding two image regions -one at a time -may only decrease the model's prediction minimally (say a single eye or mouth component on a face), while occluding these two regions together may yield a substantial change in the model's prediction if these two regions interact non-linearly as is expected for a deep neural network. Our work, together with related methods such as LIME [40], RISE [36] and more recently Sobol [13] addresses this problem by randomly perturbating the input image in multiple regions at a time.
Surprisingly, RISE [36] and Sobol Attribution [13] have recently shown that black-box attribution methods can rival and even surpass the white-box methods commonly used without recourse to internal states. However, despite the efforts in [13] to limit their computational overhead, black-box methods remain far from white-box methods in terms of execution time. In this work, we show that it is possible to match or even surpass the performances of current black-box methods while reaching computation times lower than some white-box methods by using dependence measure-based Global Sensitivity Analysis (GSA).
Global sensitivity analysis using dependence measures Our attribution method builds on the GSA framework. The approach was introduced in the 70s [9] and was popularized with variancebased sensitivity analysis and Sobol indices [51]. It consists of evaluating the sensitivity of a model's output of interest to some input design variables. GSA is currently used in many fields, especially for the study of physical phenomena [27,38]. Recently, dependence measure-based sensitivity analysis was introduced in [10] and was shown to circumvent some practical issues of variancebased sensitivity analysis. In particular, by relying on the representation capabilities of RKHS, the dependence measure that we use in this work, HSIC [22], captures more diverse information than traditional variance-based indices for far fewer model evaluations.

Explanations using Hilbert-Schmidt Independence Criterion
In this section, we describe sensitivity analysis-based attribution methods, define Hilbert-Schmidt Independence Criterion (HSIC) [22] and explain how we can use it and adapt it to design a new black-box attribution method whose efficiency competes with white-box methods. We also explain the theoretical advantages of HSIC that we can build on to go beyond traditional attribution methods.

Sensitivity analysis of perturbed black-box models
Let f : X = X 1 , ..., X n → Y be the model under study, x i ∈ X i the input variables and y = f (x 1 , ..., x n ) ∈ Y the output value of the model f . GSA studies the sensitivity of y to each input x i by considering them as iid (independent and identically distributed) random variables and assessing the link between their distribution and that of the output after an initial input sampling. Given an input vector X = (x 1 , ..., x n ), a prediction y = f (X) can thus be explained using sensitivity analysis by applying random perturbations x = (x 1 , ..., x n ), x i ∼ P Xi of the original X and analyzing the importance of each x i for explaining the variations of y -which is considered a random variable, y ∼ P y .
For image data, the inputs x i are pixels. However, pixel perturbations would only emphasize low level explanations. To obtain high level and meaningful explanations, we rather consider a random perturbation mask M = (M 1 , ..., M d ) ∈ [0, 1] d . We upsample this mask using a Nearest Neighbor interpolation method to obtain u(M ) ∈ [0, 1] n , a patch-perturbated vector that we apply on the input image X using a mask operator τ : X × [0, 1] d → X . More specifically, we use the inpainting operator defined by τ (X, M ) = X u(M ) + (1 − u(M ))µ, with the Hadamard product and µ a baseline value (here, µ is a black image with all pixels' value = 0 [40,61]). Hence, the mask M aggregates the patch-wise random perturbations M i that are sampled independently for each patch (M i are iid). In practice, the perturbations contained in the mask are binary perturbations, to simulate whether the information contained in the patch is kept in the image or not.. We thereby assess the effect of each image patch, represented by M i , on the output. .., f (τ (X, M (p) ))} and (4) statistically assessing the effect of each mask M i on y by estimating a sensitivity measure between each P Mi and P y from the previous sampling. In this paper, we consider that the more independent M i is from y, the less important the corresponding image patch is to explain it [10]. In the following, we describe HSIC, a dependence measure, and how to use it in practice.

Hilbert-Schmidt Independence Criterion
Let x and y be two random variables of probability distribution P x and P y defined on X and Y. HSIC measures the dependence between P x and P y based on their embedding in Reproducing Kernel Hilbert Space (RKHS). Let ϕ : X → F and ψ : Y → G two continuous feature mapping between X , Y, and two RKHS F, G, such that the inner product between the feature embeddings of x, x ∈ X in F is given by the kernel k(x, x ) = ϕ(x), ϕ(x ) (and l(y, y ) = ψ(y), ψ(y ) for y, y ∈ Y). The cross-covariance operator C xy : G → F between the random variables x and y is defined in [19] and can be written.
where µ x = E x [ϕ(x)] and µ y = E y [ψ(y)] are the mean embedding of x and y in F and G. When F and G are universal RKHS on the compact domains X and Y, then C xy HS = 0 if and only if x and y are independent, where · HS denotes the Hilbert Schmidt norm (see [23]). In [22], the authors define HSIC as C xy 2 HS , which can be written: where x, x and y, y are pairwise iid. HSIC can also be expressed in terms of Maximum Mean Discrepancy (MMD), which is a distance between mean embeddings defined in an RKHS [57]. More specifically, let the product RKHS P = F × G with kernel v((x, y), (x , y )) = k(x, x )l(y, y ). Then, HSIC(x, y) = γ 2 v (P x P y , P x,y ), where γ v is the MMD operator on P. HSIC thus measures the distance between P xy and P y P x embedded in P [10]. Since x ⊥ y ⇒ P xy = P y P x , the closer these distributions are, in the sense of the MMD, the more independent they are.
Thus, given a set of inputs {x 1 , ..., x p } and the associated outputs {y 1 , ..., y p }, [22] shows that HSIC can be estimated by where H, L, K ∈ R p×p and K ij = k(x i , x j ), L i,j = l(y i , y j ) and H ij = δ(i = j) − p −1 . Using this formula, H p x,y can be computed with a O(p 2 ) complexity. In this work, the input variables M i are the patch perturbations. Therefore, we compute H p Mi,y i.e. the estimation of the HSIC between a patch M i and the output y, for each patch (i ∈ {1, ..., d}), see Algorithm 1. We denote H p i := H p Mi,y for clarity. In the next section, we discuss the kernels k and l and show that we can obtain a valuable ANOVA-like orthogonal decomposition property for HSIC-based indices, allowing to assess interactions between input variables.

Orthogonalisation of HSIC to enable evaluation of interactions
One question of interest in explainability is the measurement of the importance of a specific group of variables. Indeed, it is notorious that neural networks are highly non-linear, and as it has been demonstrated in several works [16,13], the effects of the groups of variables are far from being additive. Concretely, some areas of the image may only be important in interaction with other areas, affecting the output only when both areas are perturbed at the same time -as we shall see in Section 4.4 (for instance, for mustaches of a puma).
Let {x 1 , ..., x n } ∈ X n be a set of n univariate random variables. For any subset A = {l 1 , ..., l |A| } ⊆ {1, ..., n}, we denote x A = (x l1 , ..., x l |A| ) the vector of input variables with indices in A. When using Sobol indices based GSA, it is possible to measure the interactions between variables using ANOVA decomposition. For HSIC, the corresponding decomposition property (which is not ANOVA since HSIC does not measure the variance) can be stated as follows: Property 1 (Orthogonal decomposition property). The orthogonal decomposition property is fulfilled if: HSIC(x, y) = A⊆{1,...,n} where each term HSIC A is given by and HSIC(x B , y) is defined as in equation (1) with kernels l and k A .
In Appendix A, we introduce an example to illustrate why this property is necessary in order to properly assess the importance of interactions. This property was lacking for dependence measurebased sensitivity analysis until the work of [11], which shows that when using HSIC, a specific choice of kernel k can enable this decomposition. For any choice of l and characteristic univariate kernel k, it is possible to obtain the orthogonal decomposition property by defining the input kernel k A such that These conditions can be stringent, especially the right one, which implies to compute integrals (analytically or empirically). It can be non trivial for continuous input distributions p x and classical kernel choices such as Radial Basis Function (RBF) of bandwidth σ, k(x, x ) = exp( x − x /2σ 2 ). This condition is alleviated when using discrete input variables of known densities, for which the integrals can be computed analytically. In particular, in this work, we rely on Proposition 1: Proposition 1. Let x i be a Bernoulli variable of parameter p = 1 2 , and δ(x = x ) the Dirac kernel. Then the following kernel k A satisfies the decomposition property: The proof is in Appendix A. As a practical consequence, if we sample binary masks from a Bernoulli variable of parameter p = 1/2, i.e. M i ∼ B(p) for i ∈ {1, ..., d}, and use the kernel defined in equation (5), we can assess not only the importance of each patch in the image but also the importance of the interactions between patches. It allows to go beyond classical attribution methods and reveal areas of the image that are only important in interaction with other areas, i.e. that affect the output only when both areas are perturbated at the same time. Concretely, for two image patches indexed by i and j, the interaction HSIC, H i×j , can be obtained with [11]: We insist on the fact that if the decomposition property is not valid, subtracting H p Mi,y and H p Mj ,y to H p (Mi,Mj ),y does not ensure that we assess the importance of the interactions only. Traces of the independent importance of M i and M j may remain in the obtained metric. Some qualitative benefits of such a property are illustrated in Section 4.4. This decomposition holds for any choice of kernel l : Y → G. Therefore, in the following, we use the RBF kernel, with the common practice of choosing the bandwidth as the median of the output [10,22,54].

Sample efficiency of HSIC estimator
Several types of metrics are classically used for sensitivity analysis. The most famous one, Sobol indices [52] and its variants [18,6] classically require p 2 model evaluations [25] to reach an estimation error of O( 1 √ p ). Recent design of experiments managed to reduce this requirement to p × (d + 2) [43] (Theorem 1), but with the increase in complexity of state-of-the-art architectures and the high dimensionality of inputs (although mitigated by the use of d patches instead of all pixels), it can still be cumbersome. Despite this drawback, [13] uses Sobol indices to obtain explanations and still achieves execution time improvement compared to other state-of-the-art attribution methods, which are even less efficient in terms of samples requirements.
HSIC is much less expensive to estimate than Sobol indices: for a same estimation error of O( 1 √ p ), p forward passes are needed instead of p × (d + 2) [22] (Theorem 1). This allows using far fewer samples to obtain relevant explanations, thereby dramatically increasing the efficiency of the method compared to previous black-box approaches. This huge advantage is empirically illustrated in Section 4.2, where we demonstrate that our method defines a new standard in terms of efficiency for black-box attribution methods. It even bridges the efficiency gap between black-box and white-box approaches.

Implementation of the method
A summary of the whole attribution method is provided in Algorithm 1. The computation of H p i is O(p 2 ) but it is possible to vectorize it using any library optimized for tensor operations (e.g. tensorflow). As a result, the computation time of H p i is negligible compared to that of the p forward passes. Furthermore, we implemented a sampling based on Latin Hypercube Sampling [33], a Quasi-Monte Carlo (QMC) method designed to efficiently fill the input space in Monte Carlo integration. Once the grid of H p i is obtained, we use a bilinear upsampling to be able to apply it to the image. Compute H p i using equation (2) and assign this value to the i-th patch of the input image.

Experiments
This section showcases the benefits of our approach compared to other attribution methods. These benefits are threefold. First, the computational cost of HSIC attribution method is significantly lower than previous state-of-the-art methods, even bridging the performance gap between black-box and white-box methods. Second, our method improves state-of-the-art for several fidelity metrics, for black-box as well as white-box methods. Finally, the orthogonal decomposition property of HSIC allows to go beyond usual attribution methods and assess interactions between image patches.
In Section 4.1, we compute explanations of the predictions in the ILSVRC-2012 [12] classification task (ImageNet), for four common architectures, namely MobileNet [44], ResNet50 [24], EfficientNet [59] and VGG16 [49]. Then, we compare those explanations with these of other state-of-the-art black-box and white-box attribution methods in terms of fidelity and efficiency. In Section 4.2, we investigate the convergence of our method by measuring its correlation with a high sample estimator  [13]. In the remaining sections, we conduct additional experiments that show the versatility of our method. In Section 4.3, we evaluate HSIC attribution method to explain object detection on COCO dataset [31] with YOlOv4 [39]. We conclude the experiments with Section 4.4, where we showcase the use of the HSIC orthogonal decomposition property to assess interactions between image patches.

Fidelity of classification explanations
In this section, we evaluate the fidelity of the explanations with three fidelity metrics. The first, Deletion [36], assumes that the more faithful an explanation is, the quicker the prediction score should drop when pixels that are considered important are shut down. The second one, Insertion [36], instead adds pixels on a baseline image, starting with pixels that are associated with the highest importance scores of the explanation. Finally, µFidelity [5] creates random pixels subsets which are assigned a baseline value and measure the correlation between the drop in the score and the importance of the explanation. Those metrics are further described in Appendix F.
In Table 1, we report the results of several different attribution methods for explaining the classification of MobileNet [44], ResNet50 [24], EfficientNet [59] and VGG16 [49] on 1000 images sampled from the ImageNet validation dataset. The models used for the experiments have been accessed from tensorflow [1] with the keras API [8]. We introduce two variants of our method, H p i eff. and H p i acc. The words "eff" and "acc" stand for efficient and accurate because we use p = 764 and p = 1536 samples, respectively. We use our method with a grid size of 7 × 7 (d = 49). To evaluate the different methods, we use the Xplique [15], a library dedicated to explainability. For black-box and white-box methods, we bold the best result and underline the second. When the differences between some methods are not statistically significant, we highlight both. Note that for µFidelity, the estimation variance is high (typically about 20%), so we only use the bold notation and leave the Table in Appendix B. The exact error bars are also left in the Appendix to make the presentation lighter.
For the Deletion metric, our method obtains the best results among the tested black-box methods for all the architectures in both its efficiency and accurate variants. Except for EfficientNet, we even beat all tested white-box methods. RISE is still the best of black-box methods for the Insertion metric, but the accurate variant is systematically second. Besides, our methods are among the best of both black-box and white-box methods.
While HSIC is systematically better in Deletion, we can note that RISE overshadows it in Insertion. This could be explained by how these metrics are constructed. Deletion and Insertion metrics consist in measuring Area Under the Curve (AUC) of scores that respectively decrease and increase when deleting and adding patches, starting from a baseline image (see Appendix F for a detailed definition). Since Deletion measures a drop in the score starting from the original image, the faster the score drops, the better the metric. Hence, Deletion will favor methods that sharply identify important regions. On the contrary, since Insertion starts from an arbitrary baseline image, if the explanation map is more spread out, more relevant secondary information will be added, so the score will be better. As we can see in the maps of Appendix C, RISE saliency maps are way more spread out than HSIC's, which are sharper. It may explain why RISE is better in the Insertion benchmark and why HSIC attribution method dominates the Deletion benchmark. We provide additional quantitative examples to illustrate this link in Appendix F. Note that even if RISE dominates Insertion, it is far behind in Deletion. This is not the case for HSIC, which is still competitive in Insertion while dominating Deletion.
These results, as such, are already satisfactory. But it goes even further: we obtained these results with far fewer forward passes than other state-of-the-art black-box attribution methods. With the efficient variants, competitive results are obtained more than 8 times faster than RISE, the current standard of black-box attribution methods. It improves on Sobol, a recent and promising attribution method that was already branded as more efficient, by a factor 5. The time improvement factors for the accurate variant of our method are still very appealing (5 and 3). We even beat the execution time of Integrated Gradients white-box method [58], a popular and successful white-box method. The efficiency of HSIC attribution method is investigated more thoroughly in the next section. The advantage of black-box attribution methods lies in providing explanations without access to the model's internal state or the gradients. However, this advantage comes at a cost since many forward passes are needed to obtain meaningful explanations. This cost is all the more constraining since recent architectures are increasingly heavy in terms of computational time.

Estimator efficiency
Therefore, it is critical for such attribution methods to use as few forward passes as possible. Results reported in Section 4.1 attest that our approach based on H p i shines in that regard, and we refer to this section for more comments. It motivates us to study the efficiency of our method further. To that end, we compute an "asymptotical" 1 explanation with 13, 000 forward passes, for HSIC, Sobol and RISE attribution methods, on 100 images of ImageNet validation set.
We apply the three methods on EfficientNet with d = 7 × 7 masks and image patches, like in [36] and [13]. We then compute explanations for an increasing number of forward passes and compare the obtained explanation with the baseline "asymptotical" explanation. We use Spearman rank correlation [55] like theoretically and empirically argued in [21,2,60,14]. This experiment allows comparing the convergence speed of our method with RISE and Sobol. The curves are plotted in Figure 2 and show that our method, HSIC converges much faster than RISE and Sobol.

Explanation of object detection
Explaining model's predictions is more challenging for object detection than for image classification. Indeed, recent object detection models usually predict three pieces of information: localization (bounding box corners), objectness score (probability that a bounding box contains an object of any class), classification information (probability of each different possible class). Recently, it has been demonstrated that it was possible to use attribution methods to explain object detection by constructing a score aggregating for the previous information. In [37] the authors combine intersection over union for localization, cosine similarity for the class probability, and focus on high objectness areas. As a result, they can use RISE to explain the object detection using this score as the output of the model.
In this section, we test our method for explaining the object detections of YOLOv4 [39] on COCO dataset [31] compared to the approach presented in [37], D-RISE. We also compare the explanations with these of Kernel Shap [32], another black-box attribution method. The explanations for 1, 000 validation images are evaluated with the Deletion, Insertion, and µFidelity metrics. This experiment is time-consuming, so we use H p i eff. and 5000 samples for D-RISE and Kernel Shap.  Even if our method is not the best for Deletion, it is for Insertion and µFidelity, and more importantly, it is 5 and 6 times faster than D-RISE and Kernel Shap. Figure 3 displays visualizations of object detection explanations. While the first images show a standard detection explanation, the rightmost one is more interesting since it emphasizes an error of the object detector. The model identifies a zebra as a cat, and our method manages to explain this error by emphasizing the cat at the bottom right corner of the image. Note that we did not obtain such an explanation with D-RISE, even with a high sample number and different grid sizes -visualizations can be found in Appendix C).

Finding spacial interactions in the model
Usual attribution methods provide explanations in the form of heat maps that assign each pixel (or patch) an importance score. However, the scope of such explanations is limited since the reason for a prediction may not be explained only by the single importance of independent patches. In [13], the authors use Sobol total indices that account for the importance of a patch in interaction with all other patches, but they cannot localize the interactions. Thanks to its orthogonal decomposition property, our HSIC-based attribution method is able not only to assess the importance of each patch, but also the importance of interactions between specific patches by computing the HSIC of the joint patches and subtracting the contribution of each patch taken independently 2 . It is then possible to identify regions of the image that affect the output only when both areas are perturbated jointly. We illustrate this property in Figure 4. For each image, we computed all the possible interactions between patches and reported the most important ones. Note that not all images exhibit significant interactions, so we performed a qualitative selection for pedagogical purposes. On the upper row, we plot usual heatmaps that traduce the importance of each patch taken independently. On the bottom row, important interactions are represented in red. We can see that the middle part of the snake body is not the only important element in the picture, that one part of the mountain interacts with one of the skiers, so do the mustaches of the puma and two corners of the roof. Note that for each image, the maximum values of H p i are 40.6, 11.4, 8.2 and 18.8 and the maximum values of H p i×j are 19.6, 6.6, 9.2 and 6.3 respectively. Thanks to the orthogonal decomposition property, these metrics can be compared and we can deduce that some interactions are as significant as some important independent patches. For the image with a puma, the interactions between the mustaches are even more important than the eye / forehead for identifying the animal (H p i×j = 9.2 when i and j are the right and left mustaches and H p if f = 8.2 when i is the eye / forehead).

Conclusion
We have introduced a new attribution method based on a dependence measure, Hilbert-Schmidt Independence Criterion, which leverages representation capabilities of Reproducing Kernel Hilbert Spaces, thus being able to capture complex information. This attribution method is black-box, so it is applicable even when the implementation of the neural network to explain is not available. Nonetheless, it alleviates the computational burden of traditional black-box methods, improving on the state-of-the-art of both black-box and white-box attribution methods while being closer to the latter than the former in terms of efficiency. In addition, we showed how the rich framework of RKHS could be used to assess and localize interactions between pairs of patches of the input image that are relevant for explaining the output. We hope that the introduced framework will open up research avenues for attribution methods beyond traditional pixel-wise or patch-wise explanations.

A On the Orthogonal Decomposition Property 1
In this part, we state the orthogonal decomposition Property, motivate its importance with a pedagogical example, and finally prove Proposition 1, which enables the decomposition property in the context of HSIC attribution method.

A.1 Orthogonal Decomposition Property
Let x = {x 1 , ..., x n } ∈ X n be a set of n univariate random input variables. For any subset A = {l 1 , ..., l |A| } ⊆ {1, ..., n}, we denote x A = (x l1 , ..., x l |A| ) the vector of input variables with indices in A. Let y the random output variable defined by y = f (x), F the RKHS defined by the kernel k A : X |A| → R and G the RKHS defined by the kernel l : Y → R.
In [11], the author shows that for any choice of kernel l, if we respect some constraints on the kernel k A , we can construct indices HSIC(x A , y) that satisfy the following decomposition property. Property 2 (Decomposition property). For any kernel l, the kernel k A satisfies the decomposition property if: where each term HSIC A is given by and HSIC(x B , y) is defined as in equation (1) with kernels l and k A .
The constraints on the kernel k A are detailed in the main document and in the last section of this appendix. Before describing these constraints and how to fulfill them with Proposition 1, let us illustrate the importance of the property with a motivating, pedagogical example.

A.2 Motivating example
In this section, we introduce a pedagogical example to motivate the interest in assessing the interactions and the importance of the Orthogonal Decomposition Property in that regard. Let f : [0, 2] 3 → {0, 1} such that The relation between x 2 , x 3 and the output of function f is illustrated on Figure 7. Here x i is analogous to M i . In that case, it is clear that x 1 is important to explain the output. However, assessing the effect of x 2 and x 3 is more difficult. Given the definition of f ; they are important to explain the output y, but it can be shown theoretically that HSIC(x 2 , y) = 0 and HSIC(x 3 , y) = 0 [34]. This motivates to assess the interactions between input variables. One way to retrieve the information that x 2 and x 3 are important is to assess HSIC(x 2,3 , y), where x 2,3 = (x 2 , x 3 ) One could assess HSIC(x 2,3 , y) without any constraints on the kernel k, and the obtained value for HSIC(x 2,3 , y) would indeed be > 0 . However, by doing so, we would also obtain that HSIC(x 1,2 , y) > 0 and HSIC(x 1,3 , y) > 0, whereas x 1 does not interact with x 2 and x 3 , only because of the individual effect of x 1 . We empirically illustrate this by assessing these metrics using the estimator of Eq. 2 with p = 10000, and kernels k, l chosen as the Radial Basis Function (RBF). The results are found in Table 1 below and show that: • HSIC(x 2 , y) ≈ HSIC(x 3 , y) ≈ 0 • HSIC(x 1,2 , y) ≈ HSIC(x 1,3 , y) > HSIC(x 2,3 , y) HSIC(x 1 , y) HSIC(x 2 , y) HSIC(x 3 , y) HSIC(x 1,3 , y) HSIC(x 1,2 , y) HSIC(x 2,3 , y) 1.79 × 10 −2 2.28 × 10 −6 9.63 × 10 −6 1.36 × 10 −2 1.36 × 10 −2 2.92 × 10 −3 Table 1: HSIC metrics with k taken as RBF In order to correctly assess the pairwise interactions of input variables x 1 and x 2 , one has to remove the individual effect of each variable from the HSIC(x 1,2 , y). The orthogonal decomposition property [11] allows to do so by simply computing HSIC inter (x 1,2 , y) as If the decomposition property does not hold, we are not guaranteed to fully remove the individual effect of x 1 and x 2 using the previous formula. We estimate HSIC inter (x 1,2 , y) when the kernel k satisfies the decomposition property (in that case we choose a Sobolev kernel as in [1]), and when it does not, and show that the correct information of HSIC inter (x 1,2 , y) is only retrieved when the decomposition property is satisfied. As previously, this is illustrated in the experiment, whose results are found in Table 2.

A.3 Proof of Proposition 1
To benefit from Property 2, the kernel k A must satisfy the following assumption [11]: where k 0 (x, x ) = k(x, x ) − k(x, t)dP (t) k(x , t)dP (t) k(s, t)dP (s)dP (t) .
We now recall and prove the introduced Proposition 1 defined in Section 3.3. Proposition 1. Let x a Bernoulli variable of parameter p = 1/2, and δ(x = x ) the dirac kernel such that δ(x = x ) = 1 if x = x and 0 otherwise. Let k 0 be defined as in equation (4). Then, the kernel k A satisfies the decomposition property (Property 1) if it is defined according to Assumption 1, with Proof. Let s and t be two iid random Bernoulli variables of parameter p with probability density functions p s and p t . We have that dP (s) = p s (s)ds = pδ(s = 1) + (1 − p)δ(s = 0) ds Now, let's consider two Bernoulli variables x and x , two samples x ∼ x and x ∼ x , and a kernel k such that k(x, x ) = δ(x = x ).

E Attribution methods
In the following section, we give the formulation of the different attribution methods used in this work. The library used to generate the attribution maps is Xplique [15]. By simplification of notation, we define f (x) the logit score (before softmax) for the class of interest (we omit c). We recall that an attribution method provides an importance score for each input variable x i . We will denote the explanation functionnal mapping an input of interest x = (x 1 , ..., x d ) as g(x).
Saliency [48] is a visualization technique based on the gradient of a class score relative to the input, indicating in an infinitesimal neighborhood which pixels must be modified to most affect the score of the class of interest.
Gradient Input [46] is based on the gradient of a class score relative to the input, element-wise with the input, it was introduced to improve the sharpness of the attribution maps. A theoretical analysis conducted by [3] showed that Gradient Input is equivalent to -LRP and DeepLIFT [46] methods under certain conditions -using a baseline of zero and with all biases to zero.
Integrated Gradients [58] consists of summing the gradient values along the path from a baseline state to the current value. The baseline x 0 used is zero. This integral can be approximated with a set of m points at regular intervals between the baseline and the point of interest. In order to approximate from a finite number of steps, we use a Trapezoidal rule and not a left-Riemann summation, which allows for more accurate results and improved performance (see [53] for a comparison).
SmoothGrad [50] is also a gradient-based explanation method, which, as the name suggests, averages the gradient at several points corresponding to small perturbations (drawn i.i.d from an isotropic normal distribution of standard deviation σ) around the point of interest. The smoothing effect induced by the average help reducing the visual noise, and hence improve the explanations. The attribution is obtained by averaging after sampling m points. For all the experiments, we took m = 80 and σ = 0.2 × (x max − x min ) where (x min , x max ) being the input range of the dataset.
VarGrad [26] is similar to SmoothGrad as it employs the same methodology to construct the attribution maps: using a set of m noisy inputs, it aggregate the gradients using the variance rather than the mean. For the experiment, m and σ are the same as Smoothgrad. Formally: Grad-CAM [45] can only be used on Convolutional Neural Network (CNN). Thus we couldn't use it for the MNIST dataset. The method uses the gradient and the feature maps A k of the last convolution layer. More precisely, to obtain the localization map for a class, we need to compute the weights α k c associated to each of the feature map activation A k , with k the number of filters and Z the number of features in each feature map, with α c k = 1 As the size of the explanation depends on the size (width, height) of the last feature map, a bilinear interpolation is performed in order to find the same dimensions as the input. For all the experiments, we used the last convolutional layer of each model to compute the explanation.
Grad-CAM++ (G+) [7] is an extension of Grad-CAM combining the positive partial derivatives of feature maps of a convolutional layer with a weighted special class score. The weights α (k) c associated with each feature map are computed as follows: Occlusion [61] is a sensitivity method that sweeps a patch that occludes pixels over the images using a baseline state and uses the variations of the model prediction to deduce critical areas. For all the experiments, we took a patch size and a patch stride of 1 7 of the image size. Moreover, the baseline state x 0 was zero. [36] is a black-box method that consists of probing the model with N randomly masked versions of the input image to deduce the importance of each pixel using the corresponding outputs. The masks m ∼ M are generated randomly in a subspace of the input space. For all the experiments, we use a subspace of size 7 × 7 and E(M) = 0.5.

F Evaluation
For the purpose of the experiments, three fidelity metrics have been chosen. For the whole set of metrics, f (x) score is the score after softmax of the models. We first describe these metrics and then discuss the trade-off between Deletion and Insertion scores.

F.1 Definitions
Deletion. [36] The first metric is Deletion, it consists in measuring the drop in the score when the important variables are set to a baseline state. Intuitively, a sharper drop indicates that the explanation method has well identified the important variables for the decision. The operation is repeated on the whole image until all the pixels are at a baseline state. Formally, at step k, with u the most important variables according to an attribution method, the Deletion (k) score is given by: We then measure the AUC of the Deletion scores. For all the experiments, the baseline state is fixed at x 0 = 0.
Insertion. [36] Insertion consists in performing the inverse of Deletion, starting with an image in a baseline state and then progressively adding the most important variables. Formally, at step k, with u the most important variables according to an attribution method, the Insertion (k) score is given by: We then measure the AUC of the Deletion scores. The baseline is the same as for Deletion.
µFidelity [5] consists in measuring the correlation between the fall of the score when variables are put at a baseline state and the importance of these variables. Formally: For all experiments, k is equal to 20% of the total number of variables and the baseline is the same as the one used by Deletion.

F.2 Trade-off between Insertion and Deletion
Deletion and Insertion metrics consist in measuring AUC of scores that respectively decrease and increase when deleting and adding patches, starting from a baseline image. Since the patches deleted/added are those that are the most important (in the sense of the tested attribution method), most of the score will come from the first patch deletions/additions. Using those different methods has two important consequences, detailed below.
Deletion is preferable There is a key difference between those two evaluations that makes Deletion more suited to explanation evaluation than Insertion. In Deletion, since we start from the original image and sequentially delete patches, the score is tested in a region of the input image space that is close to the input image. On the contrary, Insertion starts from an arbitrary baseline (here, pure black image), which is far from the input image. It is likely that the value of the baseline has an undesired impact on the score for Insertion. That is why we tend to favor Deletion over Insertion.
Some methods are more suited to Deletion or Insertion Since Deletion measures a drop in the score, the faster the score drops, the better the metric. Hence, Deletion will favor methods that sharply identify important regions. On the contrary, since Insertion starts from an arbitrary baseline image, if the explanation map is more spread out, more relevant secondary information will be added, so the score will be better. To illustrate this observation, in table 5 we show the value of Insertion and Deletion metrics for HSIC method and for different grid sizes, obtained after a grid search for MobileNetV2 on 1000 ImageNet validation images. The metrics are averaged over 27 runs (with a different number of samples and different samplers). Table 5 gives an idea of the trend of the evolution of Insertion and Deletion with respect to the grid size. As we can see, Deletion improves when the grid size increases, i.e. when the explanation map becomes sharper, and Insertion improves when the grid size decreases, i.e. when the map becomes more spread out.  This trend also explains why RISE shines in the Insertion benchmark and why our HSIC attribution method dominates the Deletion benchmark. Indeed, as we can see in the maps of Appendix C, RISE saliency maps are way more spread out than HSIC's, which are sharper.

G Additional experiments on stability
In this section, we report the evolution of the Deletion score for HSIC, RISE, and Sobol with respect to the number of forward passes, with a Resnet50 on 100 Imagenet validation images. The scores for Sobol and RISE are less stable than for HSIC, which corroborates that HSIC attribution method can be used with fewer forward passes.