When Explanations Lie: Why Many Modified BP Attributions Fail

Leon Sixt, Maximilian Granz, Tim Landgraf

Introduction

Explainable AI (XAI) aims to improve the interpretability of machine learning models. For deep convolutional networks, attribution methods visualize the areas relevant for the prediction with so-called saliency maps. Various attribution methods have been proposed, but do they reflect the model behavior correctly?

(Adebayo et al., 2018) proposed a sanity check: if the parameters of the model are randomized and therefore the network output changes, do the saliency maps change too? Surprisingly, the saliency maps of GuidedBP (Springenberg et al., 2014) stay identical, when the last layer (fc3) is randomized (see Figure 1(a)). A method ignoring the last layer can not explain the network’s prediction faithfully.

In addition to (Adebayo et al., 2018), which only reported GuidedBP to fail, we found several modified backpropagation (BP) methods fail too: Layer-wise Relevance Propagation (LRP), Deep Taylor Decomposition (DTD), PatternAttribution, Excitation BP, Deconv, GuidedBP, and RectGrad (Bach et al., 2015; Montavon et al., 2017; Kindermans et al., 2018; Zhang et al., 2018; Zeiler & Fergus, 2014; Springenberg et al., 2014; Kim et al., 2019). The only tested modified BP method passing is DeepLIFT (Shrikumar et al., 2017).

Modified BP methods estimate relevant areas by backpropagating a custom relevance score instead of the gradient. For example, DTD only backpropagates positive relevance scores. Modified BP methods are popular with practitioners (Yang et al., 2018; Sturm et al., 2016; Eitel et al., 2019). For example, (Schiller et al., 2019) uses saliency maps to improve the classification of whale sounds or (Böhle et al., 2019) use LRPα1β0 to localize evidence for Alzheimer’s disease in brain MRIs.

Deep neural networks are composed of linear layers (dense, conv.) and non-linear activations. For each linear layer, the weight vector reflects the importance of each input directly. (Bach et al., 2015; Kindermans et al., 2018; Montavon et al., 2017) argue that aggregating explanations of each linear model can explain a a deep neural network. Why do these methods then fail the sanity check?

Empirically, we quantify the convergence to a rank-1 matrix using our novel cosine similarity convergence (CSC) metric. CSC allows to retrace, layer by layer, how modified BP methods lose information about previous layers. Using CSC, we observe that all analyzed modified BP methods, except for DeepLIFT, converge towards a rank-1 matrix on VGG-16 and ResNet-50. For sufficiently large values of α\alpha and β\beta, LRPαβ does not converge but also produces rather noisy saliency maps.

The paper focuses on modified BP methods, as other attribution methods do not suffer from the converges problem. They either rely on the gradient directly (Smilkov et al., 2017; Sundararajan et al., 2017), which does not converge or consider the model as a black-box (Ribeiro et al., 2016; Lundberg & Lee, 2017).

Our findings show that many modified BP methods are prone to class-insensitive explanations and provide saliency maps that rather highlight low-level features. Negative relevance scores are crucial to avoid the convergence to a rank-1 matrix — a possible future research direction.

Theoretical Analysis

For our theoretical analysis, we consider feed-forward neural networks with a ReLU activation function [x]+=max⁡(0,x)[{\bm{x}}]^{+}=\max(0,{\bm{x}}). The neural network f(x)f({\bm{x}}) contains nn layers, each with weight matrices WlW_{l}. The output of the ll-th layer is denoted by hl{\bm{h}}_{l}. We use [ij]{[ij]} to index the i,ji,j element in WlW_{l} as in Wl[ij]W_{l_{[ij]}}. To simplify notation, we absorb the bias terms into the weight matrix, and we omit the final softmax layer. We refer to the input with h0=x{\bm{h}}_{0}={\bm{x}} and to the output with hn=f(x){\bm{h}}_{n}=f({\bm{x}}). The output of the ll-th layer is given by:

All the results apply to convolutional neural networks as convolution can be expressed as matrix multiplication.

Gradient

The gradient of the kk-th output of the neural network w.r.t. the input x{\bm{x}} is given by:

where Ml=diag⁡(1hl>0)M_{l}=\operatorname{diag}(1_{h_{l}>0}) denotes the gradient mask of the ReLU operation. The last equality follows from recursive expansion. The vector vk{\bm{v}}_{k} is a one-hot vector to select the kk-th output.

The gradient of residual blocks is also a product of matrices. The gradient of hl+1=hl+g(hl){\bm{h}}_{l+1}={\bm{h}}_{l}+g({\bm{h}}_{l}) is:

where G∂g(hl)/∂hlG_{\partial g({\bm{h}}_{l})/\partial{\bm{h}}_{l}} denotes the derivation matrix of the residual block, and II is the identity matrix. For the gradient, the final saliency map is usually obtained by summing the absolute channel values of the relevance vector r0∇(x)r^{\nabla}_{0}({\bm{x}}) of the input layer.

The following methods modify the gradient definition and to distinguish the rules, we introduce the notation: rl∇(x)=∂f(x)∂hlr_{l}^{\nabla}({\bm{x}})=\frac{\partial f({\bm{x}})}{\partial{\bm{h}}_{l}} which denotes the relevance at layer ll for an input x{\bm{x}}.

Interpretability of Linear Models

The relevance of the input of a linear model can be calculated directly. Let y=wTxy={\bm{w}}^{T}{\bm{x}} be a linear model with a single output scalar. The relevance of the input x{\bm{x}} to the ii-th output y[i]{\bm{y}}_{[i]} is :

z^{+}-Rule The z+z^{+}-rule is used by DTD (Montavon et al., 2017), Excitation BP (Zhang et al., 2018) and also corresponds to the LRPα1β0 rule (Bach et al., 2015). The z+z^{+}-rule backpropagates positive relevance values, which are supposed to correspond to the positive evidence for the prediction. Let wijw_{ij} be an entry in the weight matrix WlW_{l}:

The relevance of multiple layers is computed by applying the z+z^{+}-rule to each of them. Similar to the gradient, we obtain a product of non-negative matrices: Ck=∏lkZl+C_{k}=\prod^{k}_{l}Z^{+}_{l}.

Let A1,A2,A3…A_{1},A_{2},A_{3}\dots be a sequence of non-negative matrices for which lim⁡n→∞An\lim_{n\to\infty}A_{n} exists. We exclude the cases where one column of lim⁡n→∞An\lim_{n\to\infty}A_{n} is the zero vector or two columns are orthogonal to each other. Then the product of all terms of the sequence converges to a rank-1 matrix Cˉ\bar{C}:

(Hajnal, 1976; Friedland, 2006) proved a similar result for squared matrices. In appendix A, we provide a rigorous proof of the theorem using the cosine similarity.

The geometric intuition of the proof is depicted in Figure 2. The column vectors of the first matrix are all non-negative and therefore in the positive quadrant. For the matrix multiplication AiAjA_{i}A_{j}, observe that AiakA_{i}{\bm{a}}_{k} is a non-negative linear combination of the column vectors of AiA_{i}, where ak{\bm{a}}_{k} is the kk-th column vector Aj[:k]A_{j_{[:k]}}. The result will remain in the convex cone of the column vectors of AiA_{i}. The conditions stated in the theorem ensure that the cone shrinks with every iteration and it converges towards a single vector. In the appendix B, we simulate different matrix properties and find non-negative matrices to converge exponentially fast.

The Z+Z^{+} matrices of dense layers fulfill the conditions of theorem 1. Convolutions can be written as matrix multiplications. For 1x1 convolutions, the kernels do not overlap and the row vectors corresponding to each location are orthogonal. In this case, the convergence happens only locally per feature map location. For convolutions with overlapping kernels, the global convergence is slower than for dense layers. In a ResNet-50 where the last convolutional stack has a size of (7x7), the overlapping of multiple (3x3) convolutions still induces a considerable global convergence (see LRPCMP on ResNet-50 in section 5).

If an attribution method converges, the contributions of the layers shrink by depth. In the worst-case scenario, when converged up to floating-point imprecision, the last layer can only change the scaling of the saliency map. However, the last layer is responsible for the network’s final prediction.

2 Modified BP algorithms

The LRPz rule of Layer-wise Relevance Propagation modifies the backpropagation rule as follows:

If only max-pooling, linear layers, and ReLU activations are used, it was shown that LRPz corresponds to gradient⊙\odotinput, i.e. r0z−LRP(x)=x⊙∂f(x)∂xr_{0}^{z-\text{LRP}}({\bm{x}})={\bm{x}}\odot\frac{\partial f({\bm{x}})}{\partial{\bm{x}}} (Shrikumar et al., 2016; Kindermans et al., 2016; Ancona et al., 2017). LRPz can be considered a gradient-based and not a modified BP method. The gradient is not converging to a rank-1 matrix and therefore gradient⊙\odotinput is also not converging.

LRPαβ

separates the positive and negative influences:

where Zl+Z_{l}^{+} and Zl−Z_{l}^{-} correspond to the positive and negative entries of the matrix ZZ. (Bach et al., 2015) propose to weight positives more: α≥1\alpha\geq 1 and α−β=1\alpha-\beta=1. For LRPα1β0, this rule corresponds to the z+z^{+}-rule, which converges. For α>1\alpha>1 and β>0\beta>0, the matrix Zl=αZl+−βZl−Z_{l}=\alpha{Z}_{l}^{+}-\beta Z_{l}^{-} can contain negative entries. Our empirical results show that LRPαβ still converges for the most commonly used parameters α=2,β=1\alpha=2,\beta=1 and even for a higher α=5\alpha=5 it converges considerable on the ResNet-50.

Deep Taylor Decomposition

PatternNet & PatternAttribution

takes into account that the input hl{\bm{h}}_{l} contains noise. If dl{\bm{d}}_{l} corresponds to the noise and sl{\bm{s}}_{l} to the signal, than hl=sl+dl{\bm{h}}_{l}={\bm{s}}_{l}+{\bm{d}}_{l}. To assign the relevance towards the signal direction, it is estimated using the following equation:

where aia_{i} is the estimated signal direction for the i−thi-th neuron with input h{\bm{h}} and weight vector wi=W[i:]{\bm{w}}_{i}=W_{[i:]}. PatternNet is designed to recover the relevant signal in the data. Let Al[i:]=aiA_{l{[i:]}}={\bm{a}}_{i} be the corresponding signal matrix to the weight matrix WlW_{l}, the rule for PatternNet is:

PatternNet is also prone to converge to a rank-1 matrix. To recover the relevant signal, it might be even desired to converge to the a single direction – the signal direction.

The convergence of PatternNet follows from the computation of the pattern vectors ai{\bm{a}}_{i} in equation 9. It is similar to a single step of the power iteration method vk+1=Cvk/∥Cvk∥{\bm{v}}_{k+1}=C{\bm{v}}_{k}/\left\lVert C{\bm{v}}_{k}\right\rVert. In appendix C, we provide details on the relationship to power iteration and also derive equation 9 from the equation given in (Kindermans et al., 2018). The power iteration method converges to the eigenvector with the largest eigenvalue exponentially fast.

All column vectors in A[i:]=aiA_{[i:]}={\bm{a}}_{i} underwent a single step of the power iteration and therefore tend to point towards the first eigenvector of cov⁡[h]\operatorname{cov}[{\bm{h}}]. This can also be verified empirically: the ratio of the first and second singular value σ1(A)/σ2(A)>6\sigma_{1}(A)/\sigma_{2}(A)>6 for almost all the VGG-16 patterns (see Figure 3(a)), indicating a strong convergence of the matrix chain towards a single direction.

The findings from PatternNet are hard to transfer to PatternAttribution. The rule for PatternAttribution uses the Hadamard product of AlA_{l} and WlW_{l}:

The Hadamard product complicates any analytic argument using the properties of AlA_{l} or WlW_{l}. The theoretical results available (Ando et al., 1987; Zhan, 1997) did not allow us to show that PatternAttribution converges to a rank-1 matrix necessarily.

We provide a mix of theoretical and empirical insights on why it converges. The conditions of convergence can be studied well on the singular value decomposition: (Wl⊙Al)T=UlΣlVl(W_{l}\odot A_{l})^{T}=U_{l}\Sigma_{l}V_{l}. Loosely speaking, the matrix chain will converge to a rank-1 matrix if the first σ1\sigma_{1} and second σ2\sigma_{2} singular values in Σl\Sigma_{l} differ and if VlV_{l} and Ul+1U_{l+1} are aligned such that higher singular values of Σl\Sigma_{l} and Σl+1\Sigma_{l+1} are multiplied together such that the ratio σ1/σ2\sigma_{1}/\sigma_{2} grows.

To see how well the per layer matrices align, we look at the inter-layer chain members: Tl=ΣlVlUl+1Σl+1T_{l}=\sqrt{\Sigma_{l}}V_{l}U_{l+1}\sqrt{\Sigma_{l+1}}. In Figure 3, we display the ratio between the first and second singular values σ1(Tl)/σ2(Tl)\sigma_{1}(T_{l})/\sigma_{2}(T_{l}). For W⊙AW\odot A, the first singular value is considerably larger than for the plain weights WW. Interestingly, the singular value ratio of inter-layer matrices shrinks for the plain WW matrix. Whereas for PatternAttribution, the ratio increases for some layers indicating that the Hadamard product leads to more alignment of the matrices.

DeepLIFT

is the only tested modified BP method which does not converge to a rank-1 matrix. It is an extension of the backpropagation algorithm to finite differences:

For the gradient, one would take the limit x0→x{\bm{x}}^{0}\to{\bm{x}}. DeepLIFT uses a so-called reference point for x0{\bm{x}}^{0} instead, such as zeros or for images a blurred version of x{\bm{x}}. The finite differences are backpropagated, similar to infinitesimal differences. The final relevance is the difference in the kk-th logit: rlDL(x)=fk(x)−fk(x0)r^{DL}_{l}({\bm{x}})=f_{k}(x)-f_{k}(x^{0}).

Additionally to the vanilla gradient, DeepLIFT separates positive and negative contributions. For ReLU activations, DeepLIFT uses either the RevealCancel or the Rescale rule. Please refer to (Shrikumar et al., 2017) for a description. The rule for linear layers is most interesting because it is the reason why DeepLIFT does not converge:

where the mask M>0M_{>0} selects the weight rows corresponding to positive deltas (0<Δhl=hl−hl00<\Delta{\bm{h}}_{l}={\bm{h}}_{l}-{\bm{h}}_{l}^{0}). For negative relevance rlDL−r^{DL-}_{l}, the rule is defined analogously. An interesting property of the rule (13) is that negative and positive relevance can influence each other.

If the intermixing is removed by only considering W+W^{+} for the positive rule and W−W^{-} for the negative rule, the two matrix chains become decoupled and converge. For the positive chain, this is clear. For the negative chain, observe that the multiplication of two non-positive matrices gives a non-negative matrix. Non-positive vectors b,c{\bm{b}},{\bm{c}} have an angle ≤90∘\leq 90^{\circ} and cTb=∥c∥∥b∥cos⁡(c,b)≥0{\bm{c}}^{T}{\bm{b}}=\left\lVert{\bm{c}}\right\rVert\left\lVert{\bm{b}}\right\rVert\cos({\bm{c}},{\bm{b}})\geq 0. In the evaluation, we included this variant as DeepLIFT Ablation, and as predicted by the theory, it converges.

Guided BP & Deconv & RectGrad

apply an additional ReLU to the gradient and it was shown to be invariant to the randomization of later layers previously in (Adebayo et al., 2018) and analyzed theoretically in (Nie et al., 2018):

Ml=diag⁡(1h1>0)M_{l}=\operatorname{diag}(1_{h_{1}>0}) denotes the gradient mask of the ReLU operation. For Deconv, the mask of the forward ReLU is omitted, and the gradients are rectified directly. RectGrad (Kim et al., 2019) is related to GuidedBP and set the lowest qq percentile of the gradient to zero. As recommended in the paper, we used q=98q=98.

As a ReLU operation is applied to the gradient, the backpropagation is no longer a linear function. The ReLU also results in a different failure than before. (Nie et al., 2018) provides a theoretical analysis for GuidedBP. Our results align with them.

Evaluation

We report results on a small network trained on CIFAR-10 (4x conv., 2x dense, see appendix D), a VGG-16 (Simonyan & Zisserman, 2014), and ResNet-50 (He et al., 2016). The last two are trained on the ImageNet dataset (Russakovsky et al., 2015), the standard dataset to evaluate attribution methods. The different networks cover different concepts: shallow vs. deep, forward vs. residual connections, multiple dense layers vs. a single one, using batch normalization. All results were computed on 200 images from the validation set. To justify the sample size, we show bootstrap confidence intervals in Figure 4(b) (Efron, 1979). We used the implementation from the innvestigate and deeplift package (Alber et al., 2019; Shrikumar et al., 2017) and added support for residual connections. The experiments were run on a single machine with two graphic cards and take about a day to complete.

Random Logit

We display the difference of saliency maps explaining the ground-truth and a random logit in Figure 4(a). As the logit value is responsible for the predicted class, the saliency maps should change. We use the SSIM metric (Wang et al., 2004) as in (Adebayo et al., 2018).

Sanity Check

We followed (Adebayo et al., 2018) and randomized the parameters starting from the last layer to the first layer. For DTD and LRPα1β0, randomizing the last layer flips the sign of the saliency map sometimes. We, therefore, compute the SSIM also between the inverted saliency map and report the maximum. In Figure 4(b), we report the SSIM between the saliency maps (see also Figure 1(a) and appendix G).For GuidedBP, we report different saliency maps than shown in Figure 2 of (Adebayo et al., 2018). We were able to confirm a bug in their implementation, resulting in saliency maps of GuidedBP and Guided-GradCAM to remain identical for early layers.

Cosine Similarity Convergence Metric (CSC)

An alternative way to measure convergence would have been to construct the derivation matrix Ck=∏l=1kZlC_{k}=\prod^{k}_{l=1}Z_{l} and measure the ratio σ1(Ck)/σ2(Ck)\sigma_{1}(C_{k})/\sigma_{2}(C_{k}) of the first to the second-largest singular value of CkC_{k}. Although this approach is well motivated theoretically, it has some performance downsides. CkC_{k} would be large and computing the singular values costly.

We use five different random vectors per sample – in total 1000 convergence paths. As the vectors are sampled randomly, it is unlikely to miss a region of non-convergence (Bergstra & Bengio, 2012).

For convolution layers, we compute the cosine similarity per feature map location. For a shape of (h,w,c)(h,w,c), we obtain h⋅wh\cdot w values. The jump in cosine similarity for the input is a result of the input’s low dimension of 3 channels. In Figure 5, we plot the median cosine similarity for different networks and attribution methods (see appendix F for additional Figures). We also report the histogram of the CSC at the first convolutional layer in Figures 5(e)-5(g).

Results

Our random logit analysis reveals that converging methods produce almost identical saliency maps, independently of the output logit (SSIM very close to 1). The rest of the field (SSIM between 0.4 and 0.8) produces saliency maps different from the ground-truth logit’s map (see Figure 4(a)).

We observe the same distribution in the sanity check results (see Figure 4(b)). One group of methods produces similar saliency maps even when convolutional layers are randomized (SSIM close to 1). Again, the rest of the field is sensitive to parameter randomization. The same clustering can be observed for ResNet-50 (appendix E, Figure 8).

Our CSC analysis confirms that random relevance vectors align throughout the backpropagation steps (see Figure 5). Except for LRPz and DeepLIFT, all methods show convergence up to at least 0.99 cosine similarity. LRPα5β4 converges less strongly for VGG-16. Among the converging methods, the rate of convergence varies. LRPα1β0, PatternNet, the ablation of DeepLIFT converges fastest. PatternAttribution has a slower convergence rate – still exponential. For DeepLIFT Ablation, numerical instabilities result in a cosine similarity of 0 for the first layers of the ResNet-50. Even on the small 6-layer network, the median CSC is greater than 1-1e-6 for LRPα1β0 (see Figure 5(d)).

Discussion

When many modified BP methods do not explain the network faithfully, why was this not widely noticed before? First, it is easy to blame the network for unreasonable explanations – no ground truth exists. Second, MNIST, CIFAR, and ImageNet contain only a single object class per image – not revealing the class insensitivity. Finally, it might not be too problematic for some applications if the saliency maps are independent of the later network’s layers. For example, to explain Alzheimer’s disease (Böhle et al., 2019), local low-level features are sufficient as they are predictive for the disease and the data lacks conflicting evidences (i.e. the whole brain is affected).

When noticed, different ways to address the issue were proposed and an improved class sensitivity was reported (Kohlbrenner et al., 2019; Gu et al., 2018; Zhang et al., 2018). We find that the underlying convergence problem remains unchanged and discuss the methods below.

(Kohlbrenner et al., 2019; Lapuschkin et al., 2017) use LRPz for the final dense layers and LRPαβ for the convolutional layer. We report results for α=1,2\alpha=1,2 as in (Kohlbrenner et al., 2019) in Figure 6(a).

For VGG-16, the saliency maps change when the network parameters are randomized. However, structurally, the underlying image structure seems to be scaled only locally (see Figure 6(a)). Inspecting the CSC path of the two LRPCMP variants in Figure 6(c), we can see why. For dense layers, both methods do not converge as LRPz is used, but the convergence start when LRPαβ is applied. The relevance vectors of the dense layer can change the coarse local scaling. However, they cannot alter the direction of the relevance vectors of earlier layers to highlight different details.

In the backward-pass of the ResNet-50, the global-averaging layer assigns the identical gradient vector to each location of the last convolutional layer. Furthermore, the later convolutional layers operate on (7x7), where even a few 3x3 convolutions have a dense field-of-view. LRPCMP does not resolve the global convergence for the ResNet-50.

Contrastive LRP

(Gu et al., 2018) noted the lack of class sensitivity and proposed to increase it by subtracting two saliency maps. The first saliency map explains only the logit yk=y⊙mk{\bm{y}}_{k}={\bm{y}}\odot{\bm{m}}_{k}, where mk{\bm{m}}_{k} is a one-hot vector and the second explains the opposite y¬k=y⊙(1−mk){\bm{y}}_{\neg k}={\bm{y}}\odot(1-{\bm{m}}_{k}):

n(.)n(.) normalizes each saliency map by its sum. The results of Contrastive LRP are similar to Figre 1(e), no max⁡\max is applied. The underlying convergence problem is not resolved.

Contrastive Excitation BP

The lack of class sensitivity of the z+z^{+}-rule was noted in (Zhang et al., 2018) and to increase it, they proposed to change the backpropagation rule of the final fully-connected layer to:

where mk{\bm{m}}_{k} is a one-hot vector selecting the explained class. The added Nfinal fc+N^{+}_{\text{final fc}} is computed as the Zfinal fc+Z^{+}_{\text{final fc}} but on the negative weights −Wfinal fc-W_{\text{final fc}}. Note that the combination of the two matrices introduces negative entries. Class sensitivity is increased. It does also not resolve the underlying convergence problem. If, for example, more fully-connected layers would be used, the saliency maps would become globally class insensitive again.

Texture vs. Contours

(Geirhos et al., 2019) found that deep convolutional networks are more sensitive towards texture and not the shape of the object. For example, the shape of a cat filled with an elephant texture will be wrongly classified as an elephant. However, modified BP methods highlight the contours of objects rather.

Recurrent Neural Networks

Modified BP methods are focused on convolutional neural networks and are mostly applied on vision tasks. The innvestigate package does not yet support recurrent models. To our knowledge, (Arras et al., 2017) is the only work that applied modified BP rules to RNNs (LRPz for LSTMs). Training, and applying modified backpropagation rules to RNNs, involves unrolling the network, essentially transforming it to a feed-forward architecture. Due to our theoretical results, modified BP rules that yield positive relevance matrices (e.g. z+z^{+}-rule) will converge. However, further work would be needed to measure how RNN architectures (LSTM, GRU) differ in their specific convergence behavior.

Not Converging Attribution Methods

Besides modified BP attribution methods, there also exist gradient averaging and black-box methods. SmoothGrad (Smilkov et al., 2017) and Integrated Gradients (Sundararajan et al., 2017) average the gradient. CAM and Grad-CAM (Zhou et al., 2016; Selvaraju et al., 2017) determine important areas by the activation of the last convolutional layer. Black-box attribution methods only modify the model’s input but do not rely on the gradient or other model internals. The most prominent black-box methods are Occlusion, LIME, SHAP (Zeiler & Fergus, 2014; Ribeiro et al., 2016; Lundberg & Lee, 2017). IBA (Schulz et al., 2020) applies an information bottleneck to remove unimportant information. TCAV (Kim et al., 2018) explains models using higher-level concepts.

All here mentioned attribution methods do not converge, as they either rely on the gradient or treat the model as black-box. Only when the BP algorithm is modified, the convergence problem can occur. The here mentioned algorithms might still suffer from other limitations.

Limitations

Also, we tried to include most modified BP attribution methods, we left some out for our evaluation (Nam et al., 2019; Wang et al., 2019; Huber et al., 2019). In our theoretical analysis of PatternAttribution, we based our argument on why it converges on empirical observations performed on a single set of pattern matrices.

Related Work

The limitations of explanation methods were studied before. (Viering et al., 2019) alter the explanations of Grad-CAM arbitrarily by modifying the model architecture only slightly. Similarly, (Slack et al., 2020) construct a biased classifier that can hide its biases from LIME and SHAP. The theoretic analysis (Nie et al., 2018) indicates that GuidedBP tends to reconstruct the input instead of explaining the network’s decision. (Adebayo et al., 2018) showed GuidedBP to be independent of later layers’ parameters. (Atrey et al., 2020) tested saliency methods in a reinforcement learning setting.

(Kindermans et al., 2018) show that LRP, GuidedBP, and Deconv produce incorrect explanations for linear models if the input contains noise. (Rieger, 2017; Zhang et al., 2018; Gu et al., 2018; Kohlbrenner et al., 2019; Montavon et al., 2019; Tsunakawa et al., 2019) noted the class-insensitivity of different modified BP methods, but they rather proposed ways to improve the class sensitivity than to provide correct reasons why modified BP methods are class insensitive. Other than argued in (Gu et al., 2018), the class insensitivity is not caused by missing ReLU masks and Pooling switches. To the best of our knowledge, we are the first to identify the reason why many modified BP methods do not explain the decision of deep neural networks faithfully.

Evaluation metrics for attribution

As no ground-truth data exists for feature importance, different proxy tasks were proposed to measure the performance of attribution algorithms. One approach is to test how much relevance falls into ground-truth bounding boxes (Schulz et al., 2020; Zhang et al., 2018).

The MoRF and LeRF evaluation removes the most and least relevant input features and measures the change in model performance (Samek et al., 2016). The relevant image parts are masked usually to zero. On these modified samples, the model might not be reliable. The ROAR score improves it by retraining the model from scratch (Hooker et al., 2018). While computationally expensive, it ensures the model performance does not drop due to out-of-distribution samples. The ROAR performance of Int.Grad. and GuidedBP is equally bad, worse than a random baseline (see Figre 4 in (Hooker et al., 2018)). Thus, ROAR does not separate converging from non-converging methods.

Our CSC measure has some similarities with the work (Balduzzi et al., 2017), which analyzes the effect of skip connections on the gradient. They measure the convergence between the gradient vector from different samples using the effective rank (Vershynin, 2012). The CSC metric applies to modified BP methods and is an efficient tool to trace the degree of convergence.

A different approach to verify attribution methods is to measure how helpful they are for humans (Alqaraawi et al., 2020; Doshi-Velez & Kim, 2017; Lage et al., 2018).

Conclusion

In our paper, we analyzed modified BP methods, which aim to explain the predictions of deep neural networks. Our analysis revealed that most of these attribution methods have theoretical properties contrary to their goal. PatternAttribution and LRP cite Deep Taylor Decomposition as the theoretical motivation. In the light of our results, revisiting the theoretical derivation of Deep Taylor Decomposition may prove insightful. Our theoretical analysis stresses the importance of negative relevance values. A possible way to increase class-sensitivity and resolve the convergence problem could be to backpropagate negative relevance similar to DeepLIFT, the only method passing our test.

Acknowledgements

We are grateful to the comments by our reviewers, which help to improve the manuscript further. We thank Benjamin Wild and David Dormagen for stimulating discussions. We also thank Avanti Shrikumar for answering our questions and helping us with the DeepLIFT implementation. The comments by Agathe Balayn, Karl Schulz, and Julian Stastny improved the manuscript. A special thanks goes to the anonymous reviewer 1 of our paper (Schulz et al., 2020), who encouraged us to report results on the sanity checks — the starting point of this paper. The Elsa-Neumann-Scholarship by the state of Berlin supported LS. We are also grateful to Nvidia for a Titan Xp and to ZEDAT for access their HPC system.

References

Appendix A Proof of Theorem 1

In (Friedland, 2006) the theorem is proven for square matrices. In fact Theorem 1 can be deduced from this case by the following argument:

The matrices Aˉ1,Aˉ2,Aˉ3,...\bar{A}_{1},\bar{A}_{2},\bar{A}_{3},... define a sequence of non-negative square matrices that fulfill the conditions in (Friedland, 2006) and therefore converge to a rank-1 matrix.

Our proof requires no knowledge of algebraic geometry and uses the cosine similarity to show convergence. First, we outline the conditions on the matrix sequence AnA_{n}. Then, we state the theorem again and sketch our proof to give the reader a better overview. Finally, we prove the theorem in 5 steps.

Theorem 1. Let A1,A2,A3…A_{1},A_{2},A_{3}\dots be a sequence of non-negative matrices as described above such that lim⁡n→∞An\lim_{n\to\infty}A_{n} exists. We exclude the cases where one column of lim⁡n→∞An\lim_{n\to\infty}A_{n} is the zero vector or two columns are orthogonal to each other. Then the product of all terms of the sequence converges to a rank-1 matrix Cˉ\bar{C}:

Matrices of this form are excluded for lim⁡n→∞An\lim_{n\to\infty}A_{n}:

Proof sketch

To show that ∏i∞Ai\prod_{i}^{\infty}A_{i} converges to a rank-1 matrix, we do the following steps:

We define a sequence sns_{n} as the cosine of the maximum angle between the column vectors of Mn:=∏i=1nAiM_{n}:=\prod^{n}_{i=1}A_{i}.

We show that the sequence sns_{n} is monotonic and bounded and therefore converging.

We assume lim⁡n→∞sn≠1\lim_{n\to\infty}s_{n}\neq 1 and analyze two cases where we do not get a contradiction. Each case yields an equation on lim⁡n→∞An\lim_{n\to\infty}A_{n}.

In both cases, we find lower bounds on sns_{n}: αnsn−1\alpha_{n}s_{n-1} and αn′sn−1\alpha^{\prime}_{n}s_{n-1} that are becoming infinitely large, unless we have lim⁡n→∞αn=1\lim_{n\to\infty}\alpha_{n}=1 (case 1) or lim⁡n→∞αn′=1\lim_{n\to\infty}\alpha^{\prime}_{n}=1 (case 2).

The lower bounds lead to equations on lim⁡n→∞An\lim_{n\to\infty}A_{n} for non-convergence. The only solutions, we obtain for lim⁡n→∞An\lim_{n\to\infty}A_{n}, are those explicitly excluded in the theorem. We still get a contradiction and lim⁡n→∞sn=1⇒∏i∞Ai=cˉγT\lim_{n\to\infty}s_{n}=1\Rightarrow\prod_{i}^{\infty}A_{i}=\bar{c}\gamma^{T}.

Proof (1) Let Mn:=∏i=1nAiM_{n}:=\prod_{i=1}^{n}A_{i} be the product of the matrices A1⋅…⋅AnA_{1}\cdot\ldots\cdot A_{n}. We define a sequence on the angles of column vectors of MnM_{n} using the cosine similarity. Let v1(n),...,vk(n)(n){\bm{v}}_{1}(n),...,{\bm{v}}_{k(n)}(n) be the column vectors of MnM_{n}. Note, the angles are well defined between the columns of MnM_{n}. The columns of MnM_{n} cannot be a zero vector as we required AnA_{n} to have no zero columns. Let sns_{n} be the cosine of the maximal angle between the columns of MnM_{n}:

where ⟨⋅,⋅⟩\langle\cdot,\cdot\rangle denotes the dot product. We show that the maximal angle converges to 0 as lim⁡n→∞sn=1\lim_{n\to\infty}s_{n}=1, which is equivalent to MnM_{n} converging to a rank-1 matrix. In the following, we take a look at two consecutive elements of the sequence sns_{n} and check by how much the sequence increases.

(2) We show that the sequence sns_{n} is monotonic and bounded and therefore converging. Assume an+1{\bm{a}}_{n+1} and bn+1{\bm{b}}_{n+1} are the two columns of An+1A_{n+1} which produce the columns vm(n+1){\bm{v}}_{m}(n+1) and vm′(n+1){\bm{v}}_{m^{\prime}}(n+1) of Mn+1M_{n+1} with the maximum angle:

We also assume that ∥vi(n)∥=1\|{\bm{v}}_{i}(n)\|=1 for all ii, since the angle is independent of length. To declutter notation, we write vi(n)=:vi{\bm{v}}_{i}(n)=:{\bm{v}}_{i}, an=a=(a1,...,ak)T{\bm{a}}_{n}={\bm{a}}=(a_{1},...,a_{k})^{T} and bn=b=(b1,...,bk)T{\bm{b}}_{n}={\bm{b}}=(b_{1},...,b_{k})^{T}. We now show that sns_{n} is monotonic and use the definition of the cosine similarity:

Using the triangle inequality ∥∑iaivi∥≤∑iai∥vi∥\|\sum_{i}a_{i}{\bm{v}}_{i}\|\leq\sum_{i}a_{i}\|{\bm{v}}_{i}\| we get:

As we assumed that the ∥vi∥=1\left\lVert{\bm{v}}_{i}\right\rVert=1, we know that ⟨vi,vj⟩=scos⁡(vi,vj)\langle{\bm{v}}_{i},{\bm{v}}_{j}\rangle=s_{\cos}({\bm{v}}_{i},{\bm{v}}_{j}) which must be greater than the smallest cosine similarity sns_{n}:

Therefore sns_{n} is monotonically increasing and upper-bounded by 11 as the cosine. Due to the monotone convergence theorem, it will converge. The rest of the proof investigates if the sequence sns_{n} converges to 1 and if so, under which conditions.

(3) We look at two consecutive sequence elements and measure the factor α\alpha by which they increase: sn+1≥αsns_{n+1}\geq\alpha s_{n}. We are using proof by contradiction and assume that sns_{n} does not converge to 1. We get two cases, each with a lower bound on the factor of increase. For both cases, we find a lower bound for α>1\alpha>1 for all nn which would mean that sns_{n} is diverging to ∞\infty – a contradiction. Under certain conditions on lim⁡n→∞An\lim_{n\to\infty}A_{n}, we do not find a lower bound α>1\alpha>1. We find that these conditions correspond to the conditions explicitly excluded in the theorem and therefore MnM_{n} converges to a rank-1 matrix.

We multiply the first lower bound of equation 23 by 1=snsn1=\frac{s_{n}}{s_{n}} and get:

We will now pull terms corresponding to the pair (l,m)(l,m) out of the sum and for all terms in the sum, we lower bound ⟨vi,vj⟩sn≥1\frac{\langle{\bm{v}}_{i},{\bm{v}}_{j}\rangle}{s_{n}}\geq 1 by one. Let the set I:={(i,j)∣(i,j)≠(l,m),(m,l)}I:=\{(i,j)\mid(i,j)\neq(l,m),(m,l)\} index all other terms:

We know that ⟨vl(nk),vm(nk)⟩≥(1+ε)snk\langle{\bm{v}}_{l}(n_{k}),{\bm{v}}_{m}(n_{k})\rangle\geq(1+\varepsilon)s_{n_{k}}:

We absorb the mm, ll factors back into the sum:

where qˉ\bar{q} is an upper bound on ∑ijaibj\sum_{ij}a_{i}b_{j} which exists since lim⁡n→∞An\lim_{n\to\infty}A_{n} exists, which is also why lim⁡k→∞rnk\lim_{k\to\infty}r_{n_{k}} exists.

(4) Case 1: Define rn=al(n)bm(n)+am(n)bl(n)r_{n}=a_{l}(n)b_{m}(n)+a_{m}(n)b_{l}(n) analogous to rnkr_{n_{k}}. So if lim⁡n→∞rn=lim⁡k→∞rnk≠0\lim_{n\to\infty}r_{n}=\lim_{k\to\infty}r_{n_{k}}\neq 0, the factor by which sns_{n} increases would be greater that one by a constant – a contradiction:

where nk−n′n_{k}-n^{\prime} is the number of cases where rnk=0r_{n_{k}}=0 and c>0c>0 is a lower bound on the set {εrnkqˉ≠0}\{\frac{\varepsilon r_{n_{k}}}{\bar{q}}\neq 0\}. As we assumed lim⁡n→∞rn≠0\lim_{n\to\infty}r_{n}\neq 0, c>0c>0 for an infinite number of cases and therefore n′→∞n^{\prime}\to\infty when k→∞k\to\infty.

To end case 1, we have to ensure that the first sequence element is greater than zero: sn1≥s1>0s_{n_{1}}\geq s_{1}>0. This is not the case if the first NN matrices have two orthogonal columns, s1 ⁣= ⁣... ⁣= ⁣sN ⁣= ⁣0s_{1}\!=\!...\!=\!s_{N}\!=\!0. We can then skip the first NN matrices and define s1s_{1} on AN+1A_{N+1} (set Ai=AN+iA_{i}=A_{N+i}). We know NN has to be finite, as lim⁡n→∞An\lim_{n\to\infty}A_{n} has no two columns that are orthogonal.

where we used ⟨vi,vj⟩≥sn\langle{\bm{v}}_{i},{\bm{v}}_{j}\rangle\geq s_{n}. We now find a lower bound for the square of this factor. The steps are similar to case 1:

where q2=(∑ijaibj)2q^{2}=(\sum_{ij}a_{i}b_{j})^{2} and qˉ2\bar{q}^{2} is an upper bound on q2q^{2} for all nn. Note that q2−ε′2rn′>0q^{2}-\varepsilon^{\prime 2}r_{n}^{\prime}>0, since q2q^{2} has all terms that rn′r_{n}^{\prime} has but more.

(4) Case 2: So rn′r_{n}^{\prime} is a sequence that converges to zero. Otherwise, the factor by which sns_{n} increases would be greater than one by at least a constant for infinitely many nn. As in the previous case 1, this would lead to a contradiction.

(5) Case 1 and case 2 are complements from which we obtain two equations for lim⁡n→∞An\lim_{n\to\infty}A_{n}. Let ak=(a1,...,ak)T{\bm{a}}_{k}=(a_{1},...,a_{k})^{T} and bk=(b1,...,bk)T{\bm{b}}_{k}=(b_{1},...,b_{k})^{T} be columns of lim⁡n→∞An\lim_{n\to\infty}A_{n}. We get one equation per case. For all (i,j)(i,j) with i<ji<j we have:

where the first equation comes from lim⁡n→∞rn=0\lim_{n\to\infty}r_{n}=0 and the second from lim⁡n→∞rn′=0\lim_{n\to\infty}r_{n}^{\prime}=0.

For equation 36 to be true, the following set of equations have to hold:

This is equivalent to the matrix being rank one already or one of the columns is the zero vector or two are orthogonal to each other. To show why this statement holds, we are using induction on kk. For k=2k=2, we have the following set of solutions:

Case one provides rank 1 matrices only and case two gives orthogonal columns. So, the statement holds for k=2k=2.

Next assume we solved the problem for columns with kk entries and want to deduce the case where we have k+1k+1 entries (i.e. they satisfy the equations in S(k+1)S(k+1)). The pair ak+1,bk+1a_{k+1},b_{k+1} satisfies either one of the three equations in line equation 37: ak+1=bk+1=0a_{k+1}=b_{k+1}=0, ak+1=0≠bk+1a_{k+1}=0\neq b_{k+1} or ak+1≠0=bk+1a_{k+1}\neq 0=b_{k+1}. If ak+1=bk+1=0a_{k+1}=b_{k+1}=0, the rest of the non-trivial equations will be the same set of equations that one will get in the case of kk entries. If ak+1=0≠bk+1a_{k+1}=0\neq b_{k+1}, we would be left with equation aibk+1=0a_{i}b_{k+1}=0 (Case 1) or bibk+1=0b_{i}b_{k+1}=0 (Case 2) which would mean that for all i≤ki\leq k either ai=0a_{i}=0 or bi=0b_{i}=0, which will satisfy the equations in S(k+1)S(k+1). We get an analogous argument in the case of ak+1≠0=bk+1a_{k+1}\neq 0=b_{k+1}.

The other possibility is that ak+1≠0≠bk+1a_{k+1}\neq 0\neq b_{k+1}. But in this case both equations from case one and two ak+1bi=aibk+1=0a_{k+1}b_{i}=a_{i}b_{k+1}=0 and ak+1ai=bk+1bi=0a_{k+1}a_{i}=b_{k+1}b_{i}=0 lead to ai=bi=0a_{i}=b_{i}=0 for all i≤ki\leq k and this satisfies S(k+1)S(k+1) in line equation 38, concluding the induction.

This completes the proof. Since lim⁡n→∞sn≠1\lim_{n\to\infty}s_{n}\neq 1 only if lim⁡n→∞An\lim_{n\to\infty}A_{n} has a column that is the zero vector, a multiple of a standard basis vector or it has two columns that are orthogonal to each. Exactly, the conditions excluded in the theorem. For all other cases, we get a contradiction: therefore lim⁡n→∞sn=1\lim_{n\to\infty}s_{n}=1 and MnM_{n} converges to a rank-1 matrix. ∎

Appendix B Convergence Speed & Simulation of Matrix Convergences

We proved that Mn=∏inAiM_{n}=\prod_{i}^{n}A_{i} converges to a rank-1 matrix for n→∞n\to\infty, but which practical implications has this for a 16 weight-layered network? How quickly is the convergence for matrices considered in neural networks?

We know that sns_{n} increases by a factor (1+c)(1+c) greater than 1 (c>0c>0):

Each iteration yields such a factor and we get a chain of factors:

Although the multiplication chain of cnc_{n} has some similarities to an exponential form γn\gamma^{n}, sns_{n} does not have to converge exponentially as the individual cnc_{n} have to decrease (sns_{n} bounded by 11). We investigated the convergence speed using a simulation of random matrices and find that non-negative matrices decay exponentially fast towards 1.

We report the converging behavior for matrix chains which resembles a VGG-16. As in the backward pass, we start from the last layer. The convolutional kernels are considered to be 1x1, e.g. for a kernel of size (3, 3, 256, 128), we use a matrix of size (256, 128).

We test out the effect of different matrix properties. For vanilla, we sample the matrix entries from a normal distribution. Next, we apply a ReLU operation after each multiplication. For ReLU learned, we used the corresponding learned VGG parameters. We generate non-negative matrices containing 50% zeros by clipping random matrices to [0,∞][0,\infty]. And positive matrices by taking the absolute value. We report the median cosine similarity between the column vectors of the matrix.

The y-axis of Figure 7(a) has a logarithmic scale. We observe that the positive, stochastic, and non-negative matrices yield a linear path, indicating an exponential decay of the form: 1−exp⁡(−λn)1-\exp(-\lambda n). The 50% zeros in the non-negative matrices only result in a bit lower convergence slope. After 7 iterations, they converged to a single vector up to floating point imprecision.

We also investigated how a slightly negative matrix influences the convergence. In Figure 7(b), we show the converges of matrices: αW++βW−\alpha W^{+}+\beta W^{-} where W+=max⁡(0,W),W−=min⁡(0,W)W^{+}=\max(0,W),W^{-}=\min(0,W) and W∼N(0,I)W\sim\mathcal{N}(0,I). We find that for small enough β<4\beta<4 values the matrix chains still converge. This simulation motivated us to include LRPα5β4 in our evaluation which show less convergence on VGG-16, but its saliency maps also contain more noise.

Appendix C Pattern Attribution

We derive equation 9 from the original equation given in (Kindermans et al., 2018). We will use the notation from the original paper and denote a weight vector with w=Wl[i,:]{\bm{w}}=W_{l_{[i,:]}} and the corresponding pattern with a=Al[i,:]{\bm{a}}=A_{l_{[i,:]}}. The output is y=wTxy={\bm{w}}^{T}{\bm{x}}.

For the positive patterns of the two-component estimator Sa+−S_{{\bm{a}}+-}, the expectation is taken only over {x∣wTx>0}\{{\bm{x}}|{\bm{w}}^{T}{\bm{x}}>0\}. We only show it for the positive patterns a+{\bm{a}}_{+}. As our derivation is independent of the subset of x{\bm{x}} considered, it would work analogously for negative patterns or the linear estimator SaS_{\bm{a}}.

The formula to compute the pattern a+{\bm{a}}_{+} is given by:

Using the notation cov⁡[h]=cov⁡[x,x]\operatorname{cov}[{\bm{h}}]=\operatorname{cov}[{\bm{x}},{\bm{x}}] gives equation 9.

Connection to power iteration

A step of the power iteration is given by:

The denominator in equation 9 is wTcov⁡[h]w{\bm{w}}^{T}\operatorname{cov}[{\bm{h}}]{\bm{w}}. Using the symmetry of cov⁡[h]\operatorname{cov}[{\bm{h}}], we have:

This should be similar to the norm ∥cov⁡[h]w∥\left\lVert\operatorname{cov}[h]{\bm{w}}\right\rVert. As only a single step of the power iteration is performed, the scaling should not matter that much. The purpose of the scaling in the power-iteration algorithm is to keep the vector vkv_{k} from exploding or converging to zero.

Appendix D CIFAR-10 Network Architecture

Appendix E Results on ResNet-50

Appendix F Additional Cosine Similarity Figures

Appendix G Saliency maps for Sanity Checks

For visualization, we normalized the saliency maps to be in ifthemethodproduceonlypositiverelevance.Ifthemethodalsoestimatesnegativerelevance,thanitisnormalizedtoif the method produce only positive relevance. If the method also estimates negative relevance, than it is normalized to. The negative and positive values are scaled equally by the absolute maximum. For the sanity checks, we scale all saliency maps to be in $$.