Open rodrigobdz opened 2 years ago
Hey Rodrigo, thank you for raising this issue! Just a short message for now, I will investigate this more thoroughly later.
I have seen this behavior before with the Gamma rule used in combination with the ResNet-Canonization (replacing the residual connections and weighting them by contribution) in different implementations unrelated to Zennit. I think this might be a problem with this setup specifically. Therefore it was rather low on my priority list, as I perceived this more as an issue of the compatibility with EpsilonGammaBox and ResNet.
I know for a fact that the same instability also happens if you choose the ZBox bounds too low, much lower than the actual bounds of the data.
Do you have a specific setup for VGG16? At which epsilon did this happen? For VGG16 I have not seen it before, except for the aforementioned bounds.
Could you also list the specific Composites for which you have seen this behaviour? I have not seen this happening for anything that did not use the Gamma rule (see the example heatmaps for different composites for ResNet50 in the README).
Thank you for the prompt reply, Christopher!
The snippet above reproduces the instability with the following setup:
Could you also list the specific Composites for which you have seen this behaviour?
The following snippet uses the following setup:
lrp-tutorial
—i.e., the one from Grégoire's repo: gmontavon/lrp-tutorial, which we discussed in https://github.com/chr5tphr/zennit/issues/76#issue-1076622341```python import cv2 import numpy import torch from torch.nn import AvgPool2d, Conv2d, Linear from torchvision.models import resnet50 from zennit.composites import EpsilonGammaBox, NameMapComposite from zennit.core import BasicHook, collect_leaves, stabilize from zennit.rules import Epsilon, Gamma, ZBox from zennit.torchvision import ResNetCanonizer from matplotlib import pyplot as plt from zennit.image import imgify # the LRP-Epsilon from the tutorial class GMontavonEpsilon(BasicHook): def __init__(self, stabilize_epsilon=1e-6, epsilon=0.25): super().__init__( input_modifiers=[lambda input: input], param_modifiers=[lambda param, _: param], output_modifiers=[lambda output: output], gradient_mapper=(lambda out_grad, outputs: out_grad / stabilize( outputs[0] + epsilon * (outputs[0] ** 2).mean() ** .5, stabilize_epsilon)), reducer=(lambda inputs, gradients: inputs[0] * gradients[0]) ) # use the gpu if requested and available, else use the cpu device = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu') class BatchNormalize: def __init__(self, mean, std, device=None): self.mean = torch.tensor(mean, device=device)[None, :, None, None] self.std = torch.tensor(std, device=device)[None, :, None, None] def __call__(self, tensor): return (tensor - self.mean) / self.std # mean and std of ILSVRC2012 as computed for the torchvision models norm_fn = BatchNormalize((0.485, 0.456, 0.406), (0.229, 0.224, 0.225), device=device) batch_size = 1 # the maximal input shape, needed for the ZBox rule shape = (batch_size, 3, 224, 224) # the highest and lowest pixel values for the ZBox rule low = norm_fn(torch.zeros(*shape, device=device)) high = norm_fn(torch.ones(*shape, device=device)) model = resnet50(pretrained=True) model.eval() # only these get rules, linear layers will be attributed by the gradient alone # target_types = (Conv2d, AvgPool2d) target_types = (Conv2d, AvgPool2d, Linear) # lookup module -> name child_name = {module: name for name, module in model.named_modules()} # the layers in sequential order without any containers etc. layers = list(enumerate(collect_leaves(model))) # list of tuples [([names..], rule)] as used by NameMapComposite name_map = [ ([child_name[module] for n, module in layers if n == 0 and isinstance(module, target_types)], ZBox(low=low, high=high)), ([child_name[module] for n, module in layers if 1 <= n <= 16 and isinstance(module, target_types)], Gamma(0.25)), ([child_name[module] for n, module in layers if 17 <= n <= 30 and isinstance(module, target_types)], GMontavonEpsilon(stabilize_epsilon=0, epsilon=0.25)), ([child_name[module] for n, module in layers if 31 <= n and isinstance(module, target_types)], Epsilon(0)), ] # create the composite from the name map composite = NameMapComposite(name_map, canonizers=[ResNetCanonizer()]) R = None with composite.context(model) as modified_model: # compute attribution # Returns a numpy array in BGR color space, not RGB img = cv2.imread('castle.jpg') # Convert from BGR to RGB color space img = img[..., ::-1] # img.shape is (224, 224, 3), where 3 corresponds to RGB channels # Divide by 255 (max. RGB value) to normalize pixel values to [0,1] img = img/255.0 data = norm_fn( torch.FloatTensor( img[numpy.newaxis].transpose([0, 3, 1, 2])*1 ) ) data.requires_grad = True output = modified_model(data) output[0].max().backward() # print absolute sum of attribution print(data.grad.abs().sum().item()) R = data.grad heatmap = imgify( R.detach().cpu().sum(1), symmetric=True, grid=True, cmap='seismic', ) plt.imshow(heatmap) ``` **Input(s):** - Input image: ![castle.jpg](https://user-images.githubusercontent.com/14152377/173832410-88bcc3a0-0f27-46c3-a958-20d6e161b1ac.jpg) **Outputs:** - Text: ```python 2.66506890811329e+25 ``` - Heatmap: ![output](https://user-images.githubusercontent.com/14152377/173869994-68fe185c-1aa3-402c-8f5f-8a357127d6a9.png)
From my observations, I've narrowed down the issue to the denominators in the equations below.
Generic LRP rule:
> Montavon, Grégoire, Alexander Binder, Sebastian Lapuschkin, Wojciech Samek, and Klaus-Robert Müller. "Layer-wise relevance propagation: an overview." Explainable AI: interpreting, explaining and visualizing deep learning (2019): 193-209.
> Samek, Wojciech, Grégoire Montavon, Sebastian Lapuschkin, Christopher J. Anders, and Klaus-Robert Müller. "Explaining deep neural networks and beyond: A review of methods and applications." Proceedings of the IEEE 109, no. 3 (2021): 247-278.
The issue arises depending on the implementation of the stabilize
method and it is input-dependent.
I have tested the following heuristics:
```python epsilon: float = 0.1 dividend: torch.Tensor = torch.Tensor([-epsilon, 5, -5, -10]) # tensor([ -0.1000, 5.0000, -5.0000, -10.0000]) ``` 1. **Heuristic from zennit:** Add epsilon to the absolute value of the dividend conserving the sign: ```python dividend + ((dividend == 0.).to(dividend) + dividend.sign()) * epsilon ``` Example: ```python dividend + ((dividend == 0.).to(dividend) + dividend.sign()) * epsilon # tensor([ -0.2000, 5.1000, -5.1000, -10.1000]) ``` 2. **Heuristic from `lrp-tutorial`:** Scale epsilon according to dividend's magnitude using quadratic mean ```python dividend + epsilon * (dividend**2).mean()**.5 + 1e-9 ``` Example: ```python dividend + epsilon * (dividend**2).mean()**.5 + 1e-9 # tensor([ 0.5124, 5.6124, -4.3876, -9.3876]) ``` 3. **Vanilla:** Add epsilon to dividend without heuristics ```python dividend + epsilon ``` Example: ```python dividend + epsilon # tensor([ 0.0000, 5.1000, -4.9000, -9.9000]) ```
Disclaimer: The following images have been generated with my own implementation of LRP. Nevertheless, the error can also be reproduced using zennit by modifying the heuristics in the stabilize
function.
It is worth noting that heatmaps that exhibit only few visual relevance concentrations, as shown in the ResNet50 heatmaps in the first comment, do have non-zero relevance scores elsewhere, their visibility strongly depend on the plotting settings—see example below.
Heatmap of relevance scores with numerical instability:
Same relevance scores with numerical instability, adjusted plotting settings:
Heatmap with input in the same plot:
![image](https://user-images.githubusercontent.com/14152377/173880561-322bbd88-f83f-4b4b-9c63-07de88d24481.png)
stabilize
functionsAt which epsilon did this happen?
I used epsilon=1e-6
in the stabilize
function. See the following hyperparameter grid search I conducted with the composite from lrp-tutorial
:
![image](https://user-images.githubusercontent.com/14152377/173878030-a663e7bb-8833-486a-8d94-79707ccf8acc.png)
Hey Rodrigo,
thank you for your insights! After investigating the instability, I have verified my previous expectation that this heavily dependent on the gamma rule, which can lead to vanishing contributions.
This is also why it is so dependent on the heuristic used in the stabilize function.
In my analysis I have also found that the ResNet residual connections can also lead to vanishing contributions when they are attributed with Norm()
, which also requires the stabilizer to be changed accordingly.
In my conclusion, there is no global solution that could be implemented in Zennit.
However, I think exposing more control over the stabilizer in the built-in Composites, especially the epsilon value for Norm()
, will greatly help the user. Gamma is otherwise already exposed in the built-in composites. I have implemented multiple heuristics for stabilize
, however, the general gist seems to be "lower gamma" -> "higher epsilon" (e.g. through heuristic). I will think a little whether it may be better to pass a function, or the parameters for a feature-rich stabilize function, maybe making both possible would give the most freedom.
As a side note, I have also noticed that for VGG the necessary epsilon grows when the bias is not included, as I am currently preparing the 'no-bias' feature. The obvious hypothesis here is that when the bias is positive, it acts like a stabilizer itself, similar to a higher epsilon.
Here is also the snippet I used for my analysis. I am thinking of also adding something similiar as a tutorial, to also debug conservativity.
Do you agree with my observation?
I would then proceed to provide stronger stabilizer customization for Epsilon
and Norm
, maybe through a Stabilizer
class, which can be instantiated with arbitrary parameters, something along the lines of
class Stabilizer:
def __init__(self, epsilon=1e-6, clip=False, mean_scale=False, dim=None):
self.epsilon = epsilon
self.clip = clip
self.mean_scale = mean_scale
self.dim = dim
def __call__(self, input):
sign = ((input == 0.).to(input) + input.sign())
epsilon = self.epsilon
if self.mean_scale:
dim = self.dim
if self.dim is None:
dim = tuple(range(1, input.ndim))
epsilon = epsilon * ((input ** 2).mean(dim=dim) ** .5)
if self.clip:
return sign * input.abs().clip(min=epsilon)
return input + sign * epsilon
And then I will maybe just check inside the Epsilon
and Norm
rules whether epsilon
is callable or a float, and act accordingly.
Hey,
Thank you Rodrigo for putting so much effort into solving the instability and thanks to Christopher for taking the time to answer all the issues in depth. Although I can not offer a solution, I'd like to add that it would also be interesting to observe the intermediate attribution by accessing the .grad attribute of hidden layer outputs. I also observed that the instability in heatmaps already begins at hidden layers and continues to grow. Unfortunately I can not experiment now with the code as I have no access to my PC.
Best, Reduan
As community service, here is the output for resnet18 with the standard values, but bias deactivated. The first column after the layer name is the distance from the maximum to the minimum value, and the second value is the sum (should stay 1 if conservative)
Hi Chris,
Thanks that's interesting. It confirms the correlation between heatmap quality and hidden concept attribution.
I might add that this could be a problem of the gamma-rule, which does not fit (as it is) to the skip-connections of the ResNet model. Usually, the gamma rule stabilizes the relevance propagation, as it favors positive contributions (by 1+gamma) and thereby enlarges the denominator (if the denominator is positive). Exploding relevances due to a small denominator are less likely.
Formula:
However, if the denominator is negative, it becomes more positive with larger gamma. Then, the denominator can become smaller in magnitude and relevances can explode. This usually does not happen, as negative denominator means negative pre-activation, which is set to zero with ReLU non-linearities. With skip-connections (where pre-activations are added together), a negative pre-activation can receive relevance when it is added to a larger positive pre-activation.
So the problem is, that the gamma-rule assumes, that only neurons with positive outputs have relevance/contributed in the forward pass. One way to fix it would be to make the gamma rule more generic, by (1) checking the sign of the output and then either favoring positive (if output pos.) or negative contributions (if output neg.). This way, the denominator always becomes larger.
A fix like that would probably also stabilize relevance propagation if other non-linearities are used, where negative outputs receive relevance/contribute in the forward pass. As I think of it, this could also be interesting for other rules, e.g. alpha-beta. Here, for alpha1_beta0 all positive contributions receive relevance if the output is negative. Would it not be more sensible if all negative contributions receive relevance?
I have run a small experiment and making the gamma-rule symmetric (for pos. and neg. output, indicated by *) seems to result in much more reasonable results. The total relevance at the input stays bounded and heatmaps seem sensible.
``` class GammaSymmetric(BasicHook): def __init__(self, gamma=0.25): self.gamma = gamma def gradient_map(out_grad, outputs): P = outputs[0] + outputs[1] # positive contributions N = outputs[2] + outputs[3] # negative contributions alpha = 1 + gamma * (P >= N.abs()) # gamma for positive contr. beta = 1 + gamma * (P < N.abs()) # gamma for negative contr. out = P * alpha + N * beta factors = [alpha, alpha, beta, beta] return [fac * out_grad / stabilize(out) for fac in factors] super().__init__( input_modifiers=[ lambda input: input.clamp(min=0), lambda input: input.clamp(max=0), lambda input: input.clamp(min=0), lambda input: input.clamp(max=0), ], param_modifiers=[ lambda param, _: param.clamp(min=0), lambda param, name: param.clamp(max=0) if name != 'bias' else torch.zeros_like(param), lambda param, _: param.clamp(max=0), lambda param, name: param.clamp(min=0) if name != 'bias' else torch.zeros_like(param), ], output_modifiers=[lambda output: output] * 4, gradient_mapper=(lambda out_grad, outputs: gradient_map(out_grad, outputs)), reducer=( lambda inputs, gradients: ( (inputs[0] * gradients[0] + inputs[1] * gradients[1]) + (inputs[2] * gradients[2] + inputs[3] * gradients[3]) ) ) ) ```
Thank you all for taking the time to look into this issue so thoroughly. I am yet to go in detail through the snippets and some comments but I am really inspired by the fruitful conversation we have going on right now.
A few thoughts came into my mind while reading the observations made so far:
Vanishing contributions might be a different issue, because whenever the numerical instability happens, the relevances elsewhere do not vanish as it might seem from looking at the initial heatmap. Rather, the relevances elsewhere are significantly low in magnitude relative to the scores blown out of proportion.
See Comparison section of Additional insights in https://github.com/chr5tphr/zennit/issues/148#issuecomment-1156701533
Are we sure that the numerical instability affects exclusively the gamma rule?
From an implementation perspective, the stabilize function is widely used across zennit's codebase.
From a theoretical standpoint, the implementation of the generic LRP rule from section 10.2.2 in [1] affects LRP-zero/epsilon/gamma rules.
An initial step towards a solution could be to warn the user that the heatmap might not be correct whenever the relevances blow out of proportion—i.e., when conservativity is violated. E.g., when sum is not (roughly) equal to one, as shown by @chr5tphr in https://github.com/chr5tphr/zennit/issues/148#issuecomment-1159082693.
The warning itself certainly does not solve the actual issue but it would be more user-friendly. A step further would be to provide recommendations together with the warning 💯.
The in-depth explanation of the possible root cause of the numerical instability sounds plausible. I am wary about the following suggestion, as I fear that it would go against the intrinsic purpose of the Gamma rule, which is to mainly favor positive contributions:
One way to fix it would be to make the gamma rule more generic, by (1) checking the sign of the output and then either favoring positive (if output pos.) or negative contributions (if output neg.). Source: https://github.com/chr5tphr/zennit/issues/148#issuecomment-1159448340
The parameter γ controls by how much positive contributions are favored. As γ increases, negative contributions start to disappear.
Source: [1]
Nevertheless, the plots by @maxdreyer in https://github.com/chr5tphr/zennit/issues/148#issuecomment-1159653216 seem to correctly favor positive over negative contributions. 🤔
Food for thought
These insights seem to be actionable to implement a solution:
for stabilize, however, the general gist seems to be "lower gamma" -> "higher epsilon"
for VGG the necessary epsilon grows when the bias is not included
Source: https://github.com/chr5tphr/zennit/issues/148#issuecomment-1159026371
@rachtibat Could you elaborate on this?
correlation between heatmap quality and hidden concept attribution. Source: https://github.com/chr5tphr/zennit/issues/148#issuecomment-1159178438
References:
[1] Montavon, Grégoire, Alexander Binder, Sebastian Lapuschkin, Wojciech Samek, and Klaus-Robert Müller. "Layer-wise relevance propagation: an overview." Explainable AI: interpreting, explaining and visualizing deep learning (2019): 193-209.
Specifically for ResNet50, there is also a pretty nasty artifact in modified backprop explanations caused by the 1x1conv stride 2 downsampling shortcuts. Short of excluding those shortcuts, there is not really a way to get around that. I suspect that this may further contribute to the concentrated attributions you are reporting (in addition to the issues with the gamma-rule). This artifact should also appear for other composites, but not for other models.
Below are some example heatmaps showing that artifact for PascalVOC.
Hey everyone,
coming back to @maxdreyer's suggestion, it is true that the Gamme rule is only defined on positive inputs. This previously was also the case for the ZPlus rule, which we decided to make general in #5. As it turns out, there is a generalized gamma version (Supplement C.2). Given that the original Gamma should not be used for negative inputs, and that the generalized version reduces to the original version when there are no negative inputs, I think I will modify the Gamma rule to the generalized version, as we did for the ZPlus rule.
For reference only:
I just made public my LRP implementation lrp-pf-auc
and I'd like to share these Jupyter notebooks (1 and 2) with additional insights on the numerical instability.
Calculating the relevance on ResNet50 seems to be prone to a numerical instability, producing heatmaps where all attribution is concentrated in a few spots because the values in those spots have become larger than the rest. See heatmap in bug reproduction section. I can confirm that this unexpected behavior also happens using different composites and different images.
I have also seen this issue on VGG16 in my own LRP implementation depending on the heuristic used in the stabilize function.
Bug reproduction
Code based on snippet provided in https://github.com/chr5tphr/zennit/issues/76#issuecomment-992581003. **Minimal reproducible example:** ```python import cv2 import numpy import torch from matplotlib import pyplot as plt from torchvision.models import resnet50 from zennit.composites import EpsilonGammaBox from zennit.image import imgify from zennit.torchvision import ResNetCanonizer # use the gpu if requested and available, else use the cpu device = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu') class BatchNormalize: def __init__(self, mean, std, device=None): self.mean = torch.tensor(mean, device=device)[None, :, None, None] self.std = torch.tensor(std, device=device)[None, :, None, None] def __call__(self, tensor): return (tensor - self.mean) / self.std # mean and std of ILSVRC2012 as computed for the torchvision models norm_fn = BatchNormalize((0.485, 0.456, 0.406), (0.229, 0.224, 0.225), device=device) batch_size = 1 # the maximal input shape, needed for the ZBox rule shape = (batch_size, 3, 224, 224) # the highest and lowest pixel values for the ZBox rule low = norm_fn(torch.zeros(*shape, device=device)) high = norm_fn(torch.ones(*shape, device=device)) model = resnet50(pretrained=True) model.eval() # create the composite from the name map composite = EpsilonGammaBox(low=-high, high=high, canonizers=[ResNetCanonizer()]) R = None with composite.context(model) as modified_model: # compute attribution # Returns a numpy array in BGR color space, not RGB img = cv2.imread('castle.jpg') # Convert from BGR to RGB color space img = img[..., ::-1] # img.shape is (224, 224, 3), where 3 corresponds to RGB channels # Divide by 255 (max. RGB value) to normalize pixel values to [0,1] img = img/255.0 data = norm_fn( torch.FloatTensor( img[numpy.newaxis].transpose([0, 3, 1, 2])*1 ) ) data.requires_grad = True output = modified_model(data) output[0].max().backward() # print absolute sum of attribution print(data.grad.abs().sum().item()) # relevance scores R = data.grad # show maximum and minimum attribution print(torch.aminmax(R)) heatmap = imgify( R.detach().cpu().sum(1), symmetric=True, grid=True, cmap='seismic', ) plt.imshow(heatmap) ``` **Input(s):** - Input image: ![castle.jpg](https://user-images.githubusercontent.com/14152377/173832410-88bcc3a0-0f27-46c3-a958-20d6e161b1ac.jpg) **Outputs:** - Text: ```python 755226.0625 torch.return_types.aminmax( min=tensor(-7540.4922), max=tensor(2985.0886)) ``` - Heatmap: ![resnet50-heatmap](https://user-images.githubusercontent.com/14152377/173833505-5fe6d378-ede9-4f91-9e71-7d00f9e98624.png)
Additional information
The bug is not limited to the `castle.jpg` image, it can also be reproduced using the following image. See the corresponding heatmap below. - Input image: ![castle2](https://user-images.githubusercontent.com/14152377/173835000-81ef1ce1-4678-4549-aae1-dfdecc2ed34d.jpg) - Heatmap: ![resnet50-castle2-heatmap](https://user-images.githubusercontent.com/14152377/173834567-5f318957-9de1-4934-be2f-2d0bdc280016.png)