Interpretability (2).ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Interpretability.ipynb",
      "version": "0.3.2",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "J2yvF3O3XGyQ",
        "colab_type": "text"
      },
      "source": [
        "<h1><center>Interpretability of Convolutional Neural Networks</center></h1>\n",
        "<h3><center>by Sílvia Casacuberta Puig</center></h3>\n",
        "<h4><center>supervised by Dr. Seth Flaxman</center></h4>\n",
        "<h4><center>Imperial College London</center></h4>\n",
        "<h4><center>August 2019</center></h4>\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kFPM-IYbo62k",
        "colab_type": "text"
      },
      "source": [
        "## 1. Introduction\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "G-Ata476y7Ui",
        "colab_type": "text"
      },
      "source": [
        "In recent years, use of **convolutional neural networks** (CNN) has become widespread due to their success at performing tasks such as image classification or speech recognition. CNNs have achieved outstanding results at the ImageNet challenge and popularized the use of such architectures (Yosinski et al., 2015). Remarkably, these same networks that outperform at the ImageNet challenge are also successful when used for other tasks such as object detection on the PASCAL VOC dataset (Zhou et al., 2015). Large training datasets, powerful GUPs and the implementation of regularization techniques such as dropout or regularization have also helped boost the performance of CNN architectures (Zeiler & Fergus, 2014). "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DaR8SxSNPfv7",
        "colab_type": "text"
      },
      "source": [
        "While CNNs continue to excel at countless tasks and competitions, their internal procedures remain a mystery and there is very little insight into how these architectures achieve such outstanding results. We effectively treat these neural networks as **black boxes** and it is extremely challenging for humans to try to understand how they operate, due to their complexity and large number of interacting parts. However complex these black boxes might be, it is vital to acquire a deeper understanding of how they work. There are several reasons for encouraging research in this area. From a scientific perspective, understanding how deep neural networks work is interesting on its own, but mastering their inner mechanisms should also allow us to improve their results and accuracy. Without any insight on how these black boxes work, their development into better models can only be achieved by trial-and-error. From a social perspective, we should not be allowing a system that applies untransparent and unexplainable algorithms to make decisions that govern health care, banking, or politics. Once the decisions taken by deep neural networks shift from innocuously classifying hand-written digits to deciding whether someone has a particular disease or is eligible for a bank loan, it becomes imperative to advocate for a **right to explanation** (Goodman & Flaxman, 2016). "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "X1K5g79nPh1K",
        "colab_type": "text"
      },
      "source": [
        "This tutorial arises from a research project that tackles the following question: **how can we rank the hidden units of a convolutional layer** in order of importance towards the final class classification? We discuss this question thoroughly in sections 4 and 5, and in the latter we propose a novel statistical method called **PCACE** that ranks the hidden units of a convolutional layer according to their relevance to the final class score. We also analyze our results obtained with PCACE on the MNSIT dataset to point out the flaws of the activation maximization method. Moreover, in previous sections we provide different methods to **evaluate, visualize, and interpret** what is happening in the hidden layers and units of the convolutional layers, and examine how they relate to the output of the network. We provide **Python code** for all of these methods to encourage the use of interpretability tools. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HISO_gHtPi-x",
        "colab_type": "text"
      },
      "source": [
        "To make explanations and codes clearer, we have chosen a simple CNN architecture for the **MNIST** digit classification. Before we delve into our interpretation tools, in section 2 we carry out a bibliographical survey of recent papers on the area of interpretability of neural networks. In section 3 we discuss the CNN architecture that we are going to use throughout the tutorial and apply some basic interpretability tools available in the Keras documentation. In section 4 we apply some of the ideas found in the literature to numerically quantify the importance of each hidden unit. However, a statistical method appears to be more reliable, and so we turn to our PCACE algorithm to rank the hidden units of a convolutional layer. This algorithm is based on two statistical methods: Principal Component Analysis (PCA) and Alternating Conditional Expectations (ACE). Finally, we analyze our results to detect weaknesses of the activation maximization method."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wcLmxaX4MRJR",
        "colab_type": "text"
      },
      "source": [
        "Note: in this tutorial we will use the words *channel*, *unit*, and *neuron* indistinctly. We did not find a consensus in the literature to decide which term should be prioritized. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "z7lyW7W7uFm-",
        "colab_type": "text"
      },
      "source": [
        "## 2. Bibliographical Research and Recent Papers on Interpretability\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "U-ne0_YGlICb",
        "colab_type": "text"
      },
      "source": [
        "After the rapid success of deep neural networks, researchers have recently started to acknowledge the fact that we have very limited understanding of how these architectures work and how they are able to achieve such impressive results. New methods for trying to visualize and interpret neural networks have been developed in recent years, and in this section we highlight some of these methods. Some of the first visualization tools that were proposed were **saliency maps** (Simonyan et al., 2014), which indicate the areas of the input image that are discriminative and most important with respect to the given class by using the intensity of the pixels. That is, given an input image, the saliency method ranks its pixels based on their influence on the final class score by using derivatives and backpropagation. A similar idea is the **Grad-CAM** method (Selvaraju et al., 2017), which employs the gradients of any target concept to produce a localization map that highlights the important regions in the image that predict the concept. Grad-CAM can also be combined with existing pixel-space visualizations (Guided Grad-CAM) to achieve visualizations that are both high-resolution and class-discriminative. Another way to study the importance of the pixels in the input image is to **perturb** it by occluding patches of the image and see how the classification score drops (Zeiler & Fergus, 2014). Other authors follow a similar idea by directly **deleting some parts** of the image in order to find the part of the image that makes the final class score drop the most (Fong & Vedaldi, 2017). A bit differently, Koh and Liang use influence functions to gain understanding on which training points are more relevant to the final classification without having to retrain the network. Their method essentially shows how the model parameters change as we upweight a training point by an infinitesimal amount (Koh & Liang, 2017)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZBj673DAldbh",
        "colab_type": "text"
      },
      "source": [
        "A different paper (Erhan et al., 2009) focuses on creating an artificial input image that maximizes the activation of a particular layer or unit instead of highlighting the pixels in the input image. This is called **activation maximization**, and the artificial image is created by optimization with gradient ascent. In (Zeiler & Fergus, 2014), the authors try to understand CNNs backwards by creating a **Deconvolutional Network**. A DeConvNet can be thought as a ConvNet model with the same components but in reverse, so that it maps features to pixels. A similar idea is followed by (Dosovitskiy & Brox, 2016), who study image representations by **inverting** them with an up-convolutional neural network."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JB8mK_Qtlum_",
        "colab_type": "text"
      },
      "source": [
        "Some others focus on developing human-friendly concepts to help understand the machine. (Kim et al., 2018) introduces Concept Activation Vectors to provide an interpretation of the internal state of a network in **human-friendly concepts**, and Deep Dream and Lucid are Google projects that intend to *humanize* what the hidden layers and neurons see in the input images (Olah et al., 2017). Their paper *Feature Visualization* develops outstanding methods for allowing humans to visualize what each neuron is focusing on. In the conclusion, the authors point out that one of the issues that still stands out in network interpretability is finding which units are most meaningful for understanding neural net activations, which is precisely our research goal. Other papers also engage in trying to bridge human concepts and neural networks; for example, (Zhou et al., 2015) investigate how transferable are features in deep neural networks by differentiating between general and specific features learnt by the architecture. Moreover, (Alsallakh et al., 2017) asks if CNNs learn class hierarchy, and (Bau et al., 2017) study the semantic concepts learnt by the units, such as colors, scenes, and textures. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mE1aoftM2L65",
        "colab_type": "text"
      },
      "source": [
        "Other papers bring up **criticism** to some of the methods we just described. In (Kindermans et al., 2017a) the authors argue that saliency methods lack reliability when the explanation is sensitive to factors that do not contribute to the model prediction, and (Kindermans et al., 2017b) argues that DeConvNets and Guided Backpropagation do not produce the theoretically correct explanations for a linear model, and so even less for a multi-layer network with millions of parameters. Finally, (Fong & Veldadi, 2018) and (Nguyen et al., 2016) both argue that neurons do not encode single concepts and that they are in fact **multifaceted**, and also that some concepts are encoded by a group of neurons rather than by a sole neuron by itself."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YnQenC-YmohE",
        "colab_type": "text"
      },
      "source": [
        "In this literature background search we have not encountered any paper that tries to *quantify* the importance of each hidden neuron rather than studying the specific features and concepts that it is learning. We have also noted that most papers focus on analyzing the input images instead of the neurons themselves. In this tutorial we develop the PCACE algorithm, which ranks neurons in the network in order of importance, and to our knowledge there is no other algorithm or method in the literature that attempts to do so. Neither have we found in the literature the ideas that we described in the first paragraph of this section but applied to the hidden neurons instead of the input images. The aforementioned papers developed methods for ranking the importance of the pixels in the input image in relation to the final class score of the image; why not apply these methods to the elements of the weight or activation matrix of a hidden neuron to also rank their importance in relation to the final class score? In this tutorial we will apply similar ideas to the saliency, perturbing, and occlusion methods but to the hidden neurons instead of the pixels of the input image in order to rank them according to how relevant they are towards the final class score. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "w17_91-ryiHk",
        "colab_type": "text"
      },
      "source": [
        "## 3. MNIST Architecture and Basic Interpretability Tools"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LEWJl9TJmDU5",
        "colab_type": "text"
      },
      "source": [
        "### 3.1. Building a CNN Architecture"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UiRX7BnHyopy",
        "colab_type": "text"
      },
      "source": [
        "To present our methods we first need to define the architecture that we are going to be working with. As suggested in the introduction, for the sake of simplicity and efficiency we are going to present our ideas with the MNIST architecture that can be found in the Keras documentation. However, we will add some explanations to it that can be helpful to new users. Finally, we will also provide code for some basic tools for interpretability as introduced in the previous section. \n",
        "\n",
        "Note: we assume that the user has set up the necessary environment for coding deep neural networks. \n",
        "\n",
        "First of all, we need to download the necessary libraries for our code. For this tutorial we need Keras, Numpy, Matplotlib, Sklearn, ACE, Pandas, Seaborn and Scipy. From each library we are going to use the following **packages** throughout the tutorial:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "PAvICOfOzgeO",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import numpy as np\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "import matplotlib.image as img\n",
        "\n",
        "from keras import backend as K\n",
        "from keras import activations\n",
        "from keras.models import Sequential, Model\n",
        "from keras.utils import np_utils \n",
        "from keras.datasets import mnist  \n",
        "from keras.preprocessing.image import load_img\n",
        "from keras.layers import Dense, Dropout, Activation, Flatten, \n",
        "Convolution2D, MaxPooling2D, Input, Dense, load_model\n",
        "\n",
        "from sklearn import preprocessing\n",
        "from sklearn.preprocessing import StandardScaler\n",
        "from sklearn.decomposition import PCA\n",
        "from sklearn.neighbors.kde import KernelDensity\n",
        "from sklearn.metrics import confusion_matrix\n",
        "\n",
        "from vis.visualization import visualize_activation\n",
        "from vis.utils import utils\n",
        "\n",
        "from ace import ace\n",
        "from ace import model\n",
        "\n",
        "from matplotlib.pyplot import figure\n",
        "from string import ascii_uppercase\n",
        "from pandas import DataFrame\n",
        "import seaborn as sn\n",
        "from scipy import stats\n",
        "\n",
        "K.set_image_dim_ordering('th')"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "C58zLLJKzh07",
        "colab_type": "text"
      },
      "source": [
        "The MNIST dataset is readily available from Keras datasets, so we just need to import it. Recall that the MNIST dataset consists of 60,000 28x28 images of hand-written digits with the corresponding correct labels. You can check that these numbers are right by running\n",
        "```\n",
        "print(np.shape(X_train))\n",
        "```\n",
        "which outputs (60000, 1, 28, 28). Before we start building the network we need to **preprocess** the data by normalizing it and splitting it between training and testing sets:\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9b_pM1RQzk_G",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "(X_train, y_train), (X_test, y_test) = mnist.load_data() \n",
        "\n",
        "X_train = X_train.reshape(X_train.shape[0], 1, 28, 28) \n",
        "X_test = X_test.reshape(X_test.shape[0], 1, 28, 28)\n",
        "\n",
        "X_train = X_train.astype('float32') \n",
        "X_test = X_test.astype('float32')\n",
        "X_train /= 255  \n",
        "X_test /= 255\n",
        "\n",
        "Y_train = np_utils.to_categorical(y_train, 10)\n",
        "Y_test = np_utils.to_categorical(y_test, 10) "
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Jd_g8CuOAOOV",
        "colab_type": "text"
      },
      "source": [
        "We are now ready to build the network's architecture. We are going to use the same architecture as the one found in the Keras documentation: \n",
        "\n",
        "* 2D convolutional layer with 3x3 kernel size\n",
        "* ReLu activation function\n",
        "* 2D convolutional layer with 3x3 kernel size\n",
        "* ReLu activation function\n",
        "* 2D Max pooling with 2x2 pooling size\n",
        "* Dropout with 0.25 probability\n",
        "* Flatten\n",
        "* Dense layer\n",
        "* ReLu activation function\n",
        "* Dense layer\n",
        "* Softmax function\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "okU7sxLpBdz8",
        "colab_type": "text"
      },
      "source": [
        "In our further analysis of this network we will mainly be concerned with the first convolutional layer with its activation function and the final dense layer with the softmax function. However, we need the max pooling and dropout layers for regularization purposes. Coding this architecture is straightforward using the Keras sequential model:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "v3w5aWhtCNtA",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "model = Sequential()  \n",
        "\n",
        "model.add(Convolution2D(32, 3, 3, input_shape=(1,28,28)))\n",
        "model.add(Activation('relu'))\n",
        "\n",
        "model.add(Convolution2D(32, 3, 3))\n",
        "model.add(Activation('relu'))\n",
        "\n",
        "model.add(MaxPooling2D(pool_size=(2,2)))\n",
        "model.add(Dropout(0.25))  \n",
        "model.add(Flatten()) \n",
        "\n",
        "model.add(Dense(128)) \n",
        "model.add(Activation('relu'))\n",
        "\n",
        "model.add(Dropout(0.5)) \n",
        "model.add(Dense(10))  \n",
        "model.add(Activation('softmax'))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OCGAwGNdCQIi",
        "colab_type": "text"
      },
      "source": [
        "We can obtain a summary of the different layers by running\n",
        "```\n",
        "model.summary()\n",
        "```\n",
        "Finally, we choose a loss function and an optimizer to **compile** the model, and then we **fit** the model with our training data. Once this is done, we save the model in an h5 file. In this case we only run one epoch because it already achieves high accuracy. However, one can run more epochs to improve the performance of the model (but should be aware of overfitting problems). \n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_qd04az4CciC",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "model.compile(loss='categorical_crossentropy',  \n",
        "              optimizer='adam', \n",
        "              metrics=['accuracy']) \n",
        "\n",
        "model.fit(X_train, Y_train, \n",
        "          batch_size=32, nb_epoch=1, verbose=1)\n",
        "\n",
        "model.save('mnist_meu.h5')"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jhCEsPZIDD3p",
        "colab_type": "text"
      },
      "source": [
        "Once the model has been trained, we can **evaluate** its loss and accuracy with the testing data (that the model has never seen). We can also output the number of correct and incorrect classifications. Together they should add up to 1,000 because this is the size of the testing data."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "KRIua2gtDNDe",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "mnist_model = load_model('mnist.h5')\n",
        "predicted_classes = mnist_model.predict_classes(X_test)\n",
        "\n",
        "correct_indices = np.nonzero(predicted_classes == y_test)[0]  \n",
        "incorrect_indices = np.nonzero(predicted_classes != y_test)[0]\n",
        "\n",
        "print(len(correct_indices),\" classified correctly\")\n",
        "print(len(incorrect_indices),\" classified incorrectly\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qoa-HibIDQlo",
        "colab_type": "text"
      },
      "source": [
        "We can play around a bit more with the correctly and incorrectly classified images. You can run the following code to plot a correctly classified image and an incorrect one, including the false label it obtained. It is important to remember to reshape the images to size (28, 28) to be able to plot them with Matplotlib:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YvZdVN1lDsK_",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "num1 = np.random.randint(0, len(correct))\n",
        "num2 = np.random.randint(0, len(incorrect))\n",
        "\n",
        "plt.imshow(X_test[correct[num1]].reshape(28,28))\n",
        "plt.title(\n",
        "    \"Predicted: {}, Truth: {}\".format(predicted_classes[correct[num1]],\n",
        "                                    y_test[correct[num1]]))\n",
        "plt.figure()\n",
        "\n",
        "\n",
        "plt.imshow(X_test[incorrect[num2]].reshape(28,28))\n",
        "plt.title(\n",
        "  \"Predicted: {}, Truth: {}\".format(predicted_classes[incorrect[num2]],\n",
        "                                    y_test[incorrect[num2]]))\n",
        "plt.figure()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pd9PAh98KJbl",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1yYonww08MJkNz7-wwo5IOzrgosnhUqb4' width=400px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_cZX9NaeEe8F",
        "colab_type": "text"
      },
      "source": [
        "A simple but very visual and informative tool that we can use to understand why our network is classifying incorrectly some of the input images is a **confusion matrix**. Each row of the matrix represents the instances in a predicted class while each column represents the instances in an actual class. \n",
        "\n",
        "With the following code we can print a regular confusion matrix and then another one that contains the same information but looks nicer:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "TCqOKxSxE6bY",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "model = load_model('mnist.h5')\n",
        "\n",
        "test_predictions = model.predict_classes(X_test)\n",
        "cm = confusion_matrix(y_true = y_test, y_pred = test_predictions) \n",
        "print(cm)\n",
        "\n",
        "columns = ['class %s' %(i) for i in list(ascii_uppercase)\n",
        "           [0:len(np.unique(y_test))]]\n",
        "confm = confusion_matrix(y_test, test_predictions)\n",
        "df_cm = DataFrame(confm, index=columns, columns=columns)\n",
        "\n",
        "plt.figure(figsize = (15,5))\n",
        "sn.heatmap(df_cm, annot=True, cmap='Blues', linewidths=.5, fmt='g')"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YwzIlOR2HpA6",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1OYgRY1hrCC1OnoCfxYBPAIpGBEBzAGoY' width=700px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "sWiY3-PlmS6d",
        "colab_type": "text"
      },
      "source": [
        "### 3.2. Weights and Activations"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pXYdsiRoFiji",
        "colab_type": "text"
      },
      "source": [
        "Once we have built, trained, and tested our network, we can start using some basic interpretability tools that are part of the Keras documentation. First of all, there are two main aspects of the network that we can study and visualize: the learnt **weights** of each layer during the training process and the **activations** that these produce when we feed an input image into the network. \n",
        "\n",
        "To visualize the **weights,** we can use Keras' function ```get_weights```. To do so we need to indicate the layer that we are interested in. The layers are simply indexed in the order in which we built them. So for example, in our network, layer 0 corresponds to the first convolutional layer, and layer 2 to the second convolutional layer. So if we wanted to obtain the weights of the first convolutional layer we would run:\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YHGK90BS370h",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "weights = model.layers[0].get_weights()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LwTj9VWl4ApE",
        "colab_type": "text"
      },
      "source": [
        "If we then printed  ```np.shape(weights[0])```, we would obtain (3, 3, 1, 32) because there are 32 channels in the first conv layer and we applied a 3x3 kernel size (the 1 appears because the MNIST dataset is black and white; this would be a 3 in the case of RGB images). Later in the tutorial we will study these weights a bit more, but now we turn to visualizing the **activations**. We can print the 32 3x3 matrices (one for each channel) in matrix form or by using the function ```imshow``` to print the weights in a more visual way. This is the code that will output the 32 numeric matrices and the mean weight of each channel (the mean of the 9 elements in the matrix):\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "IS9fctujocJ_",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "mean = []\n",
        "for i in range(0, 32):\n",
        "    print('Channel ', i)\n",
        "    print(weights[0][:, :, 0, i])\n",
        "    print(np.mean(weights[0][:, :, 0, i]))\n",
        "    mean.append(np.mean(weights[0][:, :, 0, i]))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-O4SHtOOJWH5",
        "colab_type": "text"
      },
      "source": [
        "![alt text](https://drive.google.com/uc?id=1eiWlvCgwnVJbhxH0o1cFszv8dL7GVaSs)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3eT5RYfRooHQ",
        "colab_type": "text"
      },
      "source": [
        "Then we can plot the 32 mean values to see how the weights vary across the first layer (the x-axis corresponds to the channel index and the y-axis to the mean activation of the channels):"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "NOFX5_2boxhN",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "plt.plot(mean, 'bo')\n",
        "plt.title('Mean of the channel weights in the first layer')\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NubzBOM-KbDA",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=12rEuP02EKXOldR9p4RBTBySUz2gnZriM' width=350px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MLA25zOrzIXJ",
        "colab_type": "text"
      },
      "source": [
        "If we analyze the mean weight for each channel, we observe that there are two types of channels: the ones with positive mean weight and the ones with negative mean weight. Let us visually separate these two types of channels by manually adding a red line on the previous plot:"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0KENrkr8zE9A",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1R-9b9gZnuI73DGLFDYrCASt-cm8cSt68' width=350px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VNjlpuPy0kj6",
        "colab_type": "text"
      },
      "source": [
        "As we can see, the channels that have **negative mean weight** are: 5, 10, 12, 13, 15, 22, 24, 25, and 26. Later in the tutorial, when we examine the activations of all the channels we will see that these exact same channels that have negative mean weight are the ones that will focus on specific parts of the input digit instead of the global shape. We will refer to these channels as **reverse channels** and they will play a very important role in our research findings. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "X1J8_Scno6zm",
        "colab_type": "text"
      },
      "source": [
        "On a different note, we can also print the 32 3x3 matrices in a more visual way by using Matplotlib. We should unify the scale of the color bar so that all 32 channels use the same colors for the same values by using the ```vmin``` and ```vmax``` options in ```plt.imshow```. In this particular example we use -0.4 and 0.3 because we computed the maximum and minimum values in all the 32 3x3 matrices. For other particular examples, these values can be computed by simply using ```np.max``` and ```np.min``` functions from the Numpy library."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "t0vhNU43pJNu",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "for i in range(0, 32):\n",
        "    plt.imshow(weights[0][:, :, 0, i], vmin=-0.4, vmax=0.3)\n",
        "    plt.title('Channel {}'.format(i))\n",
        "    plt.colorbar()\n",
        "    plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GGN2XxujLJyk",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1rd8OBR5b_3GccsCZGRS3X8n5mVZGMWCK' width=1000px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6-OHZKZVoLMu",
        "colab_type": "text"
      },
      "source": [
        "To visualize the activations we need to use Keras' functional API Model. We first run:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6Dwd7mQ347Px",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "model = load_model('mnist.h5')\n",
        "layer_outputs = [layer.output for layer in model.layers] \n",
        "activation_model = Model(inputs=model.input, outputs=layer_outputs)  "
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DjG02xvU5W_B",
        "colab_type": "text"
      },
      "source": [
        "Recall that, contrary to the weights, which are intrinsic to the network after training, **activations** are obtained when an input image is fed into the network. In this example, we will feed into the network the first image in the testing dataset, which turns out to be a 7. We must not forget to reshape the image into the dimensions that our Model is expecting:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "il60hz8X59QJ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "activations = activation_model.predict(X_test[0].reshape(1, 1, 28, 28)) "
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wk_fvw1x5_Nd",
        "colab_type": "text"
      },
      "source": [
        "Now we can visualize the activations as images. Similarly to the visualization of weights, we need to indicate which layers' activations we want to show by using their index. For example, ```activations[1]``` corresponds to the activations of the first layer with the input image after the ReLu activation function. If we run ```print(np.shape(activations[1])```, we are going to obtain (1, 32, 26, 26) because the first layer has 32 channels and each outputs a 26x26 image of the activation of that particular channel in the first layer. To make this clearer, let us visualize the activations of the tenth channel of the first layer when we feed the first image of the testing dataset into the network:\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DwVcrOpG6_Hy",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "plt.matshow(activations[1][0, 10, :, :], cmap='viridis') \n",
        "plt.title('Layer 1, Channel 10')\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NYRJqEe4MZBq",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1nPsxOdQJKtm4u37b7bxRIBaObYEjgMte' width=230px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BCh3rES8mZ4i",
        "colab_type": "text"
      },
      "source": [
        "### 3.3. Saliency Maps and Grad-CAM"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ojdlFPXn7P-r",
        "colab_type": "text"
      },
      "source": [
        "Once we know how to visualize the weights and the activations of each layer, let us turn to more sophisticated methods that we introduced in the previous bibliographical and literature research part of this tutorial. These methods try to offer an insight on *how* the network is classifying. Similar to the distinction between weights and activations (where the weights are static and instrinsic to the trained network whereas each activation corresponds to a particular input image), we can either feed an image into the network and extract information on which features of this particular image are given more importance by the network when classifying it, or we can extract information on the mechanisms of the network without having to feed any input image. For the first case, we can use **saliency maps** and **grad-CAM** methods; for the latter, we can use **activation maximization**. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3DVsvoshQTne",
        "colab_type": "text"
      },
      "source": [
        "Both saliency maps and grad-CAM methods arise from the need of understanding *what* is it in the input that makes the network classify the image into a certain class. For example, given the image of a cat, is the network focusing on its ears, tail, fur, or is it focusing on something irrelevant to the class (like the background sofa of the image)? As mentioned in section 2, saliency maps were introduced in the paper *Deep Inside Convolutional Networks: Visualising Image Classification Models and Saliency Maps* (Simonyan et al., 2014). The idea behind saliency maps is to compute the gradient of the output category with respect to the input image. This indicates how the output values change when the input image varies a bit. Then these gradients are used to highlight the pixels in the input image that contribute the most towards the final class classification. Keras makes it very handy to visualize saliency maps with its visualization toolkit."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WgoN5dcx_COo",
        "colab_type": "text"
      },
      "source": [
        "For the following example on how to produce a **saliency map**, we are going to input into the network the second image of the MNIST testing dataset, which happens to be a 2. We also need to specify the layer whose activations we want to maximize in the input image. In this case, we want to generate a heatmap that shows us the parts of the input image that were more important for the final class classification, so this is why we need to use the layer ```dense_10``` (the last one). As a small trick, instead of counting the index of the specific layer ourselves, we can just have Python look it up for us by using its name: "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Q1ToToq8A8b7",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "model = load_model('mnist.h5')\n",
        "layer_idx = utils.find_layer_idx(model, 'dense_10') \n",
        "model.layers[layer_idx].activation = activations.linear\n",
        "model = utils.apply_modifications(model) "
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4wI-wz8PBDm7",
        "colab_type": "text"
      },
      "source": [
        "We can now produce the saliency map. We are going to obtain different results depending on which backprop modifier we use: vanilla gradient, guided backprop or rectified saliency (as explained in the Keras documentation). Here is the code for the three of them:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "b4_ou680BgiX",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# original\n",
        "plt.imshow(X_test[1].reshape(28,28))\n",
        "plt.title('original')\n",
        "plt.show()\n",
        "\n",
        "# vanilla gradient\n",
        "grads = visualize_saliency(model, layer_idx, filter_indices=0, \n",
        "                      seed_input=X_test[1])\n",
        "plt.imshow(grads, cmap='jet')\n",
        "plt.title('vanilla')\n",
        "plt.show()\n",
        "\n",
        "# guided backprop\n",
        "grads = visualize_saliency(model, layer_idx, filter_indices=0, \n",
        "                       seed_input=X_test[1], backprop_modifier='guided')\n",
        "plt.imshow(grads, cmap='jet')\n",
        "plt.title('guided')\n",
        "plt.show()\n",
        "\n",
        "# rectified saliency or deconv saliency\n",
        "grads = visualize_saliency(model, layer_idx, filter_indices=0, \n",
        "                      seed_input=X_test[1], backprop_modifier='relu')\n",
        "plt.imshow(grads, cmap='jet')\n",
        "plt.title('relu')\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kRU86M6SN-QM",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1RleQqrTb2oOtT0FnsvKnpZb90VlLrdqx' width=750px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vVOKQgC0Bs8b",
        "colab_type": "text"
      },
      "source": [
        "Using a similar code, it is also straightforward to produce similar heatmaps but using the **grad-CAM** method instead of the saliency one. As found in the Keras documentation, grad-CAM is another way of visualizing attention over the input image. Contrary to the saliency method, which uses gradients with respect to the output, grad-CAM uses the penultimate convolutional layer (the one before the last dense layer) to take into account the spatial information that is lost in the dense layer.\n",
        "\n",
        "In this case we can also employ three different backprop modifiers:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "j6uqBCEeB8Ic",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# original\n",
        "plt.imshow(X_test[1].reshape(28,28))\n",
        "plt.title('original')\n",
        "plt.show()\n",
        "\n",
        "# vanilla\n",
        "grads = visualize_cam(model, layer_idx, filter_indices=0, \n",
        "                      seed_input=X_test[imgnum], \n",
        "                      backprop_modifier=None)  \n",
        "plt.imshow(grads, cmap='jet')\n",
        "plt.title('vanilla')\n",
        "plt.show()\n",
        "\n",
        "# guided\n",
        "grads = visualize_cam(model, layer_idx, filter_indices=0, \n",
        "                     seed_input=X_test[imgnum], \n",
        "                      backprop_modifier='guided')  \n",
        "plt.title('guided')\n",
        "plt.imshow(grads, cmap='jet')\n",
        "plt.show()\n",
        "\n",
        "# relu\n",
        "grads = visualize_cam(model, layer_idx, filter_indices=0, \n",
        "                      seed_input=X_test[imgnum], \n",
        "                      backprop_modifier='relu')  \n",
        "plt.title('relu')\n",
        "plt.imshow(grads, cmap='jet')\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JVgmw9RuOlQz",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1anWQNDjT0F6_qgGkJ3IZLVgc3gXAuP8W' width=750px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "B3ZNfctgTjew",
        "colab_type": "text"
      },
      "source": [
        "As it is clear from the previous explanation of saliency maps and grad-CAM, both methods rank the pixels in the input in order of importance towards the final classification of the image. This is exactly our research goal, but we \n",
        "intend to rank the hidden units instead of the pixels in the input image. However, we will apply the idea behind saliency maps and grad-CAM to the hidden units in order to rank them: how does a small change in the hidden unit affect the final class score? Similar to the saliency and grad-CAM methods, the assumption is that the more a pixel or hidden unit changes the final class score when altered a tiny bit, the more important it is."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "q_vu4DY9mfbd",
        "colab_type": "text"
      },
      "source": [
        "### 3.4. Activation Maximization"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ys9dXa6FCHE6",
        "colab_type": "text"
      },
      "source": [
        "Finally, we turn to the **activation maximization** method. This is one that we are going to use more in our later analysis of interpretability of neural networks. The idea of the activation maximization method is to produce an artificial input image that maximizes the activation of a particular layer. This way, if we produce an artificial image that maximizes the activations of the last dense layer, we are in fact creating an artificial image that maximizes its final class score. As found in the Keras documentation, the activation maximization method generates an input image that maximizes the filter output activations. This allows us to visualize which inputs activate the most a particular layer. Therefore, for each of the 10 class scores (digit 0, ..., digit 9), we can create the *ideal* image for this number because it is the one that obtains the highest class score. Later on we will keep referring to them as *ideal* images.\n",
        "\n",
        "The following code will create and let us visualize the **ideal image** of each of the ten classes (which should resemble the digits from 0 to 9):"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_SBXuDDLEo9B",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "layer_idx = utils.find_layer_idx(model, 'activation_16')  \n",
        "\n",
        "model.layers[layer_idx].activation = activations.linear\n",
        "model = utils.apply_modifications(model) \n",
        "\n",
        "for i in range(0,10):\n",
        "    img = visualize_activation(model, layer_idx, filter_indices=i, \n",
        "                               input_range=(0., 1.))\n",
        "    plt.imshow(img[..., 0])\n",
        "    plt.title('number {}'.format(i))\n",
        "    plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TJgCQU8RQLPF",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1j3zT4X-YcqqMVFMkJjWZ_3nAbIQAkeZQ' width=750px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RSU84Lw5E5yY",
        "colab_type": "text"
      },
      "source": [
        "Now comes an idea that we develop to understand a network better: first we create the artificial  input image that maximizes a particular class score and then we feed this image into the network to study the activations of each layer. In theory, we would expect this image to fire the neurons *ideally*. However, the results that we obtain reveal some flaws of the activation maximization method. Refer to section 4.1. to visualize the activations of the hidden units of the convolutional layers when the input of the network is the ideal image of each digit."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KWL6VX7Kt8e3",
        "colab_type": "text"
      },
      "source": [
        "### 3.5. Importance of Each Layer for Image Classification"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "sFJI2QBit-Sf",
        "colab_type": "text"
      },
      "source": [
        "Once we have presented the basic tools for interpretability that can be found in the Keras documentation, we will start presenting new ideas and methods that can be used building up on the ones we have just explained. This tutorial has been created following a research question on interpretability: which hidden units are more important for the final class classification? As always, by important we mean how relevant they are to the final output of the network (the assigned correct class score to the input digit)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9RhoEgdKlqjZ",
        "colab_type": "text"
      },
      "source": [
        "In the following part of this tutorial we are going to develop this idea and present statistical and variation tools to rank the hidden channels in order of importance, but in this section we will first focus on the layers as a whole without breaking them up into channels in order to introduce these ideas. Therefore, we first ask ourselves: which of the layers of the network matters the most? How can we quantify this importance? We propose the following method: for each layer of the network, we get the weights *w* of that layer via ```get_weights``` and then add a small quantity to them to obtain the new value *w'*. Statistically, we found 0.1 times the standard deviation of the original weights to be a good choice for the small quantity that we add, which in this particular case turns out to be 0.005. Then we build a new network in which the weights of the layer are *w'* instead of *w* by using Keras' function ```set_weights```, and finally we evaluate the new test accuracy of the network, which will be lower than the original one. We repeat this process several times (keep adding 0.005 to the weights of the layer) and see how the test accuracy drops. As you can see, this method attempts to apply the idea of the **gradient:** how does a small variation in the input affect the final output? The bigger the final change is, the more relevant that layer appears to be."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "l1nuOeWquTgs",
        "colab_type": "text"
      },
      "source": [
        "Now let us apply this method to our MNIST case. First of all, there are only 3 layers to compare: first conv layer, second conv layer, and final dense layer (the others serve for regularization purposes only). We are going to add 0.005 to the original weights 50 times, and for each time we will evaluate the network with the new weights on the testing data. Here is the code for the first convolutional layer:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4VA_ZyTqvEkR",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "accuracy = []\n",
        "\n",
        "for i in range(0, 50):\n",
        "    weights = model.layers[0].get_weights()\n",
        "    change = weights\n",
        "    change = [x+0.005 for x in weights]\n",
        "    \n",
        "    model.layers[0].set_weights(change)  \n",
        "    weights_nous = model.layers[0].get_weights()\n",
        "    \n",
        "    evaluating = model.evaluate(X_test, Y_test, verbose=2)\n",
        "    print(\"Test Loss\", evaluating[0])\n",
        "    print(\"Test Accuracy\", evaluating[1])\n",
        "    \n",
        "    accuracy.append(evaluating[1])"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6YoLgduWvR79",
        "colab_type": "text"
      },
      "source": [
        "Once the program is finalized, we can plot the evolution of the test accuracy across the 50 iterations:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "QU6POfrpvWyf",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "plt.plot(accuracy)\n",
        "plt.title('Test accuracy change when perturbing layer 0')\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "F6De2TwpQ2ro",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1h1tMmM7x09WRnGrMKM3VN6V1Sw_dQQmc' width=350px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "r2h8ZCZdvmth",
        "colab_type": "text"
      },
      "source": [
        "Recall that the test accuracy is always a number between 0 and 1 (if one wants the percentual test accuracy, simply multiply by 100). By repeating the exact same procedure with the second convolutional layer and the final dense layer, we can put these 3 plots together to see how we can rank the 3 layers in order of importance:"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ehV4fMJtQ9zW",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1CAYlhxM9Vwx2O6wXYdUZsC5IYXMRtt7V' width=350px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rKIAcxcqzl1B",
        "colab_type": "text"
      },
      "source": [
        "It is also very interesting to compute the confusion matrix that results from perturbing each layer to see how each layer relates not only to the test accuracy but to each digit class in particular. We are not going to analyze these confusion matrices because we present them as only a useful tool for the user to know, but for example we can see that perturbing layer 5 makes the network mistake almost all digits 1 for digit 9 whereas it does not change too much the classification of digit 2. Another observation that we make is that the first layers of the network (the two convolutional ones) do not mistake a particular class with another one (the mistakes are fairly split across different classes), whereas the last dense layer mixes up two very particular classes, and digits from most classes are mistaken to be a 9. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "J0A9gSO9SN2H",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=11JK1uSNH63cckgd6X_Vi8k2J_NqPRUWd' width=700px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QNlIdMTuvvxX",
        "colab_type": "text"
      },
      "source": [
        "Going back to the previous plot of the accuracy loss, we can observe how the layers are ranked in inverse order: the final layer seems to be the most important for the final classification and the first layer the least important. From these results, one might hypothesize that in any network the last layer will be the most important because it is the closest one to the final softmax function; then the second to last, etc., with the first layer being the least important by the same logic. To test this hypothesis, which we found to be false, we applied the exact same procedure to the Cifar-10 architecture. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zDnV5lzAwaHX",
        "colab_type": "text"
      },
      "source": [
        "**Cifar-10** is a dataset that contains 60,000 color images in 10 different classes. After the MNIST dataset, this is the simplest one to work with. We can build a CNN architecture to train this dataset. To do so we are going to follow the architecture found in the Keras documentation: \n",
        "\n",
        "* Two convolutional layers with their respective activation functions\n",
        "* One max pooling and one dropout layers\n",
        "* Again, two convolutional layers with their respective activation functions\n",
        "* Again, one max pooling and one dropout layers\n",
        "* After a flattening, a dense layer with an activation function\n",
        "* One dropout layer\n",
        "* Final dense layer and softmax function"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gt3YtMQ-wbiM",
        "colab_type": "text"
      },
      "source": [
        "As we just outlined, the Cifar-10 architecture is more complex than the MNIST one, and there are 5 layers to consider: 4 convolutional layers and 1 final dense layer. The code to build and train the model is essentially the same one we presented for the MNIST architecture, but loading the Cifar-10 data from Keras instead (```(x_train, y_train), (x_test, y_test) = cifar10.load_data()```), and running more epochs because it is harder to achieve a high test accuracy in this case. To apply the gradient idea to the Cifar-10 layers to rank them in order of importance, as we just did with the MNIST architecture, we also use essentially the same code but with different layer indices. \n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rLC1NZs4SwAk",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1MiGDI-7RpcVjXRljWyP9Z-E5pG-U2Zrx' width=400px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jL6tBcSKERpG",
        "colab_type": "text"
      },
      "source": [
        "Since the plot is not as clear as before and the final test accuracies after the 50 iterations are closer to each other, these are the precise final test accuracies that we obtain with the code for each layer of the Cifar-10 architecture:\n",
        "\n",
        "* First convolutional layer: 0.0923\n",
        "* Second convolutional layer: 0.1089\n",
        "* Third convolutional layer: 0.1173\n",
        "* Fourth convolutional layer: 0.1163\n",
        "* Final dense layer: 0.1001"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wSl6LLU6F406",
        "colab_type": "text"
      },
      "source": [
        "In the Cifar-10 case, as it can be seen, although the final test accuracies are much closer to one another than in the MNIST case, the importance of the layers is not ranked in inverse order. The first convolutional layer appears to be the one that causes the largest variation on the final classifications when changed a little bit, therefore indicating that it is the layer that matters the most (if we assume this method properly quantifies the importance of a layer), whereas in the MNIST case the final layer was the one that was ranked as the most important. This is an example of how each architecture behaves differently and generalizations should be treated carefully.\n",
        "\n",
        "As a final remark, this method that we used to rank the layers can easily be applied to any known and already trained architecture. Keras offers the possibility to quickly obtain models that have already been trained (and therefore we are able to directly obtain their weights). Keras has the following available models: Xception, VGG16, VGG19, ResNet, ResNetV2, ResNeXt, InceptionV3, InceptionResNetV2, MobileNet, MobileNetV2, DenseNet and NASNet. \n",
        "\n",
        "For example, to extract the weights learnt by the VGG16 architecture to experiment, it suffices to import the model in the following way:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "BT5wzONQyVAL",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "from keras.applications.vgg16 import VGG16\n",
        "\n",
        "model = VGG16()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OLZSLV7oz9WW",
        "colab_type": "text"
      },
      "source": [
        "## 4. Basic Interpretability Tools Applied to the Hidden Units"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XqgKS2QW0PRO",
        "colab_type": "text"
      },
      "source": [
        "After introducing the basic CNN architectures along with Keras' interpretability tools, we are going to present our own tools and ideas for interpreting a neural network. As we have mentioned, our research question is clear: how can we rank the hidden units in order of importance?"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Jh0BhL1VaKEj",
        "colab_type": "text"
      },
      "source": [
        "### 4.1. Visualizing the Hidden Units"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GgY78ka20CRz",
        "colab_type": "text"
      },
      "source": [
        "First of all, it is important to visualize these hidden units to understand what we are talking about. The best way to visualize them is by plotting their activations when an input image is fed into the network. Now the question is: which image should be fed into the network?\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_Z5T9PV6Soym",
        "colab_type": "text"
      },
      "source": [
        "#### 4.1.1. One Training Image"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Onk0qicCaVhj",
        "colab_type": "text"
      },
      "source": [
        "One possibility is to input a random **image from the training set** for each digit. So let us plot the activations of the 32 channels of the first convolutional layer when we input a random image from the training set for each of the 10 digits. It is important to remark that from now on we are going to show the code and the plots for the channels of the first convolutional layer only for the sake of simplicity, but this can easily be reproduced with any other layer of the network by just changing the layer index that we are using. For the purposes of analysis and study of these activations, while plotting these images we are also going to compute the mean activation of each of the 32 channels. Recall that by the MNIST architecture that we built before, the first convolutional layer contains 32 channels whose activation matrix has dimensions 26x26 so with the following code we should expect to obtain 320 26x26 images showing the activations of each channel. \n",
        "\n",
        "First we are going to create an array that will contain all the images of the training set that correspond to one particular digit. So for example, for the digit 3, we check all the labels in ```y_train``` and we append the corresponding ```X_train``` image to our array ```digit``` if the label is a 3. Here is the code for digit 3:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZVY28TOA0BiU",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "n = 3\n",
        "digit = []\n",
        "\n",
        "for i in range (0, len(y_train)):\n",
        "    if y_train[i] == n: \n",
        "        digit.append(X_train[i])"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6m5_CCmm24cb",
        "colab_type": "text"
      },
      "source": [
        "but of course ```n``` can be any integer value between 0 and 9. Now we are going to plot the activations of each of the 32 channels when the first digit of the array ```digit``` is fed into the network and also compute the **mean activation:**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "FUAYKvP_3uxe",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "for n in range(0, 9):\n",
        "  digit = []\n",
        "\n",
        "  for i in range (0, len(y_train)):\n",
        "      if y_train[i] == n: \n",
        "          digit.append(X_train[i])\n",
        "\n",
        "  plt.imshow(digit[0].reshape(28,28))\n",
        "  plt.title('input: random image')\n",
        "  plt.show()\n",
        "\n",
        "  mean = []\n",
        "\n",
        "  for z in range(0, 32): \n",
        "      img = digit[0]\n",
        "      activations = activation_model.predict(img.reshape(1, 1, 28, 28))\n",
        "\n",
        "      plt.imshow(activations[1][0, z, :, :])\n",
        "      plt.title(\"Digit {}, Channel {}\".format(n, z))\n",
        "      plt.colorbar()\n",
        "      plt.show()\n",
        "\n",
        "      mean.append(np.mean(activations[1][0, z, :, :]))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jzy1V4UYSftI",
        "colab_type": "text"
      },
      "source": [
        "#### Show "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zDW4PtFiWs4u",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1QJD6OQwum7EtvyXR7F2eGwZETGMR5zXB' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6ynI0dZ0XWcL",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1qA04qDKeyh1fcwYcQZbBw3bIvfw7gFB7' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cWfnWHw5X0NO",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1CFbWsB7feNlhg4f9bwS7-XVHSMhisyf3' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DJ7SGfHjX4Sl",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1AQc7qUKqulgTqTukodRL2GFqW8tAEwuV' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "72FtmYZlYW5c",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1lgkrcVR5RpZN5UnKBYViK0JOqQRKjp-Q' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "78YInMy9b3JL",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1SKOot_cf-17zHltz5ps0hoyfph9DFYuG' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fYhiYiRzcPlb",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1Bo6K-5VctiDv6vmLSgWt22YmrxXmhc-K' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vh23PZCwcxR0",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1gvrAmCm2bEmmBG2Z0-s70qpH5xi_6xjN' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BkCQrsCPdTg_",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=12sC4LNljrUXXWIm0CYaxn_RTGIwDYntK' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "V-sDdxGEdW1C",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=10G03QQvU3TPItXVzBtPnmDjsF4DtpRD7' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "b4a-Oka-4IF4",
        "colab_type": "text"
      },
      "source": [
        "Since we have stored the mean activations for each channel and for each digit, if we wanted to plot them we would simply run"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jI38oXiH4J6A",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "plt.plot(mean, 'bo')\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JWe0E9hIM03e",
        "colab_type": "text"
      },
      "source": [
        "Just by plotting the activations of each channel with a random image from the training set, we already see that there are some channels that appear to be different. The majority of the channels see the digit in its pretty much original shape in their activations. However, for all digits, channels 5, 10, 12, 13, 15, 22, 24, 25 and 26 seem to be focusing on something different: they all focus on the **surroundings of the digit** (edges, corners, background image). Some activate only at some corners and edges (like channel 15) and some also at the entire background image that is *not* the digit (like channel 5). Therefore, these channels focus on the specific details of the digit **rather than the entire bigger picture** of the digit. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1eijmOnp4QFd",
        "colab_type": "text"
      },
      "source": [
        "Even more fascinating is the fact that these are *exactly* the channels that have a **mean negative weight** and that we called **reverse channels** in section 3.2. Why so? We believe that the negative elements of the weight matrix of the reverse channels invert what the channel is focusing on and therefore activates with the pixels of the image that are not the digit itself. Moreover, in general, we see that the more negative the mean weight is, the less the channel focuses on the main shape of the digit (the *more reverse* it is). For example, we classified channel 10 as a reverse channel because its mean is negative; however, its weight mean is very close to 0. At the same time, we can see that for all digits, while channel 10 does focus on the edges and corners of the digit as the other channels do, it also activates with the main shape of the digit much more than the other reverse channels. However, channels 5, 15, or 22, which have a more negative mean weight, do not activate at all at the parts of the image that correspond to the main body of the image. We will discuss the impact of the negative weights of the reverse channels and their impact on the activation function a bit more in the conclusion. These *reverse* channels are going to be the main focus of our research study from now on."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "sUDSvb_LSzVQ",
        "colab_type": "text"
      },
      "source": [
        "#### 4.1.2. Mean Real Image"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "n_pAtgBe5tgC",
        "colab_type": "text"
      },
      "source": [
        "But this is not the only type of input image we can feed into the network. Since a random image of the training set may not be very representative of the digit class, we can also compute the *mean* image of a digit by computing the mean of all the images of that particular digit in the training dataset (elementwise). Nevertheless, we must be aware that this method works very well for the MNIST dataset because all the hand-written digits in the training set have approximately the same size and are centered. In other types of data sets, computing the mean image may not make as much sense. We shall call the final image after computing the mean the *mean real image* from now on. \n",
        "\n",
        "Since the training set might be very large, we can also decide how many training images of a particular digit in the dataset we want to use to build the **mean real image**. In the following code, we output the main real image for digit 5 by using 500 hand-written images of a 5 from the training dataset:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "uOriWWfiYp0i",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "iterations = 500  \n",
        "resultats = []\n",
        "\n",
        "n = 5  \n",
        "digit = []\n",
        "\n",
        "for i in range (0, len(y_train)):\n",
        "    if y_train[i] == n: \n",
        "        digit.append(X_train[i])\n",
        "\n",
        "digit2 = []\n",
        "for i in range(0, iterations):\n",
        "    digit2.append(digit[i])\n",
        "\n",
        "w, h = 28, 28;  \n",
        "mean_image = [[0 for x in range(w)] for y in range(h)] \n",
        "\n",
        "for k in range(0, len(digit2)):\n",
        "    for i in range(0, 28):\n",
        "        for j in range(0, 28):\n",
        "            mean_image[i][j] += digit2[k][0][i][j]\n",
        "\n",
        "for i in range(0, 28):\n",
        "    for j in range(0, 28):\n",
        "        mean_image[i][j] /= len(digit2)     \n",
        "\n",
        "plt.imshow(mean_image)\n",
        "plt.title('Mean Real Image')\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KowMntzD10NJ",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1EJ9nyFTl-QIEtIhZ2_Zi-54IKUGXzvbq' width=200px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Ou2UY4gOZR9-",
        "colab_type": "text"
      },
      "source": [
        "Once we have generated the mean real image for each digit, we can now feed it into the network by using the same code as before and print the activations for each of the 32 channels of the first layer. We can interpret these results as the mean activations that we get by feeding all (or a subset of) the training dataset images."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3HF5CALgfB6R",
        "colab_type": "text"
      },
      "source": [
        "#### Show"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "b4p8B5fJfCwl",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1wVuUO1Dyu-v9GRpr7KaClHUlspALg06l' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3D3XRFt_fZe7",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=17klnbmzGQl2Tc5zoS9-nandA7QKBUNv0' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yOfgieyPguwC",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1CFnUGaVoZBm5gzGSmfVEzz1ldj8pMoyw' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "r05TNzNLhOXR",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1S7TxKxgOVDOpFHfODH3bY3Mow5xhhVG2' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MXgl26lCiHOm",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1kUyL2GmkAbLpYS7PnzbZ2pUWLT2GNAsU' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nV5y6qb2ieKP",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1QevfdDUreApaDQgUC54xc7jceF2LLa6B' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "IgrFfE4tjAn1",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1NBZuXV4Ihs4BMhpjCEpojZ6i62Di2KG8' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3uqj4DIwjr1D",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1ck8-2uPtohHv7ki2yhwZ4YdGm2fQk7Os' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Udhb4VzNkBxd",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1lSgIfY-tYo4NCEy8yVYxqWB8VCrHwcHV' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0jmfga1slGSD",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1cQoRn0_P5AVn0VwVib8OaPhCzqp8MMzh' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RKFZBvflaDup",
        "colab_type": "text"
      },
      "source": [
        "As we can observe, we continue to have the same *reverse* channels we obtained when feeding a random image from the training set, but they look sharper and more different. These visualizations make these differences even more evident than they were before. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bVMijWwkS8o8",
        "colab_type": "text"
      },
      "source": [
        "#### 4.1.3. Ideal Digit"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "snkg31psbLDg",
        "colab_type": "text"
      },
      "source": [
        "Finally, there is another input image that we can use to visualize the activations of each channel. Instead of using real images from the training dataset, we can use the activation maximization method that we introduced before to create an artificial image that maximizes the final class score for each digit, that is, what we called the *ideal* image for each digit. Then we feed the **ideal digit** into the network and plot the activations of the hidden channels. We would expect to find the *ideal* activations for each channel as well."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7ilgItkYljR5",
        "colab_type": "text"
      },
      "source": [
        "#### Show"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "u2YAIFI7lkAp",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1ys9NzsoljubhopAPrqLpicbs2Cx9R-Lx' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "543s3biRl9nX",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1utfMiMT2sguDlJc4K9gpoMoYX_O1KCTU' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tbH2wm_BmYKx",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1c9GvByLQC3JBp06gDeS9xWrnJoOb2gqr' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gjRnP9KLmtb2",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1Diyv6PHBf8neezzEx3hSKItrlXOdi2TY' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GzjpZ966nPqS",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1fRe7ildvCDc69-LVrTCzSDgcwLzyH67J' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TJRYuCzlnmO1",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1Wg8QlxPmJpnFl4bLiXrGKbtIAv0COwb4' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TRcKHzknn8Cq",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1ekd6SftUazrMIlXM56soptCo3bwjPwc9' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "o5IwLiDQoRjk",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1QlnvW3LChcV4EVwZp2Z8D4jxTUYD5Mq2' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Uhrdz3gMoyi4",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1BAg9AE8y6V7CHCOMMA2L9ZgGtOzBdoC1' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2Y5u2j4VpUrB",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1l5JAqrL9-GR-eLRIgmKyW_ngjLYXcC6B' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "t3a_xEJGb1hC",
        "colab_type": "text"
      },
      "source": [
        "However, these results turn out to be quite surprising. For every digit, the activations of the reverse channels are completely 0 (dead neurons, as they are called sometimes in the literature) or nearly 0. Recall that the reverse channels as we defined them before (the ones that focus on the surroundings of the digit and that have negative mean weight) are 5, 10, 12, 13, 15, 22, 24, 25 and 26. Now, if we focus closely on the activations that each ideal digit causes on the channels of the first layer, we can observe that for all the ideal digits, channels 5, 15, and 22 are completely **dead or mostly dead**, and channels 12, 24, 25, 26 are also nearly dead in all cases. Therefore, we seem to observe that the activation maximization method produces ideal images that disregard the channels that focus on the details of the input. The digits are produced to maximize the part that corresponds to the global big picture but the the surroundings are not, and therefore the **reverse channels are not activated**."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ihOiR8rNM6j7",
        "colab_type": "text"
      },
      "source": [
        "In order to show more clearly how the reverse channels have very low activations for each of the ideal digits in comparison to the other channels, we have plotted the mean activation for the 32 channels for each of the 10 digits. We have manually added a red line that separates the channels with the lowest activations (below the line; closer to 0) and thus it is very easy to see that for all the 10 digits, the channels below the red line (the nearly dead ones) are always reverse channels. More precisely, channels 5, 12, 15, 22, 24, 25, 26 are the lowest ones and the ones below the red line, but the other reverse channels are also low in comparison to the other channels. The x-axis of the plots corresponds to the index of the channel and the y-axis to the mean activation of that channel when the ideal digit is fed into the network:"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "u_Dg1imfQUlb",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1ohDCdztD9orFleELNTgjljmYxgaz5DPg' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PH9VmCw7T3YA",
        "colab_type": "text"
      },
      "source": [
        "From the results that we are obtaining, one fear that might arise is the following: is any of the channels focusing on a particular location of the input images due to overfitting or some other issue? This is a fear that we might have with the reverse channels, since they focus on edges and corners and the surroundings of the digits. This would not make much sense in a convolutional layer since CNNs are meant to be **location invariant** and therefore focus on aspects of the image independently of their location. To overcome this doubt, we can visualize the activations of the channels in yet another way: by creating a **noisy** 28x28 random image and feeding it into the network. We are going to create a random input image with Numpy's function ```np.random.uniform```. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JowsaidiUu4K",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "fake = X_train[0]  # just to get the right size\n",
        "\n",
        "for i in range(0, 28):\n",
        "    fake[0, i, :] = np.random.uniform(low=0, high=1, size=(28,))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LIDcFP1FU8vP",
        "colab_type": "text"
      },
      "source": [
        "Then, using the same code as before, we can plot the 32 activations that we obtain for each digit. This is what we obtain:"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vSJfI1nGpfEi",
        "colab_type": "text"
      },
      "source": [
        "#### Show"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YxeDnXW-pgEs",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1qIuKDdbxyG9vN2REJfola-gEcclIlnYL' width=950px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gIiKSlx1OF_T",
        "colab_type": "text"
      },
      "source": [
        "We can observe that the reverse channels do not activate with any fixed particular location of the image, and therefore we have double checked that our architecture is not overfitting and that CNN layers are indeed location invariant."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "N0OP229eVROQ",
        "colab_type": "text"
      },
      "source": [
        "Lastly, another way to visualize activations with some randomness would be to **add noise** to an image from the training dataset to perturb it a little bit, feed it into the network, and then visualize the activations of the channels. We call ```img``` the image that we obtained from the training dataset. By using again Numpy's ```np.random.normal``` we add noise to the initial image by creating random numbers which are part of a Normal Gaussian distribution, centered at the mean of the channel activation, and with standard deviation that of the activation matrix as well. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "znX003D4VzNY",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "amean = np.mean(img[0, :, :])\n",
        "astd = np.std(img[0, :, :])\n",
        "for i in range(0, 28):\n",
        "    noise = np.random.normal(amean, astd, 28)\n",
        "    img[0, i, :] += noise"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5fKQHji63q5Y",
        "colab_type": "text"
      },
      "source": [
        "We obtain the following activations which again evidentiate the reverse channels."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TZ6i10SvsnYs",
        "colab_type": "text"
      },
      "source": [
        "#### Show"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LO9C_akSsoTq",
        "colab_type": "text"
      },
      "source": [
        "![alt text](https://drive.google.com/uc?id=1lUM15uPSrF3wC_DSoXfi5BnGDDCC0yDY)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UxMLhJtI3ASY",
        "colab_type": "text"
      },
      "source": [
        "### 4.2. Variation Methods"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "InUhG7mzcUpM",
        "colab_type": "text"
      },
      "source": [
        "Once we have visualized the activations of each neurons by using three different types of input images, we turn back to our main research question: how can we **rank these channels in order of importance**? In this section we will attempt to apply some versions of the ideas that we saw in the literature that ranked the pixels of the input images but to the hidden units instead. We will explain and provide the code needed for three of these methods: **gradient, killing, and noise methods**, as we call them. Each method can be applied to either the weights or the activations of the layer."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "W3HtUKqsM_a5",
        "colab_type": "text"
      },
      "source": [
        "However, while these ideas carry useful insight into understanding the role of each hidden neuron, we believe that they are not the best method to purposely rank them in order of importance. The variations that we observe do not seem to be significant enough and might neglect some features of the network when trying to compare the variations across the hidden neurons. Therefore, the gradient, killing and noise methods can be effectively used to study the impact of a hidden unit towards the final class score, but they should be more carefully considered when comparing the results across the different units. To rank the neurons in order of importance, we find a much more reliable statistical method that allows a nuanced and standardized comparison between the importance of the different units: PCACE. This method is presented and discussed in section 5."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Ltdr_VRYdkY6",
        "colab_type": "text"
      },
      "source": [
        "#### 4.2.1. Gradient"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CLWtj0VZdsWo",
        "colab_type": "text"
      },
      "source": [
        "This method follows in spirit the idea behind the **saliency and grad-CAM methods** that we explained in section 3.3. (Simonyan et al., 2014; Selvaraju et al., 2017) and that we applied to the entire layers in section 3.5.: how is the final class score affected when we apply a small change to the neuron? Instead of using **derivatives**, we will use the intuitive idea behind them: look up the original class score of the particular digit, add a small quantity to the neuron, and look up the new class score to compute the variation with respect to the original class score. As a remark, we take the class score of the digit before applying the softmax function (that is, from the penultimate layer instead of the last one), and we recommend doing so because the softmax function changes the relationship between the 10 class scores and modifies them so that they are all between 0 and 1, which alters the work previously done by the network (Simonyan et al., 2014). This is why we call this method *gradient*, because despite not using calculus explicitly, we apply the idea of the gradient: we **add a small amount** to the neuron and see how much this impacts the final class score. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "S29pyrk8NEs4",
        "colab_type": "text"
      },
      "source": [
        "As mentioned in the introduction, this method can be applied to either **the weights or the activations** of each neuron, depending on what we wish to study. The difference is that if we apply the variation to the weights, then we only need to do so once and then we can feed input images without changing anything else. When applied to the activations, it takes a bit more work: for *every* input image to the network, we let it run through until it activates the neuron that we are focusing on, stop, modify its activations, build a new model, and let it finish run through to obtain the final class score. Another problem that we need to consider is the following: what quantity should we add to the weight or activation matrix? If we add a fixed quantity, then we are not taking into account the values of the matrices themselves. If we add the same fixed quantity, a neuron with high magnitude weights or activations will induce less change to the final class score than a neuron with lower weights or activations. Therefore, we should take into account somehow the distribution of the weights or activations. From a statistical standpoint, we have decided to add 0.1 times the standard deviation of the elements in the weight or activation matrix of the channel."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "iYsD5LvHNITv",
        "colab_type": "text"
      },
      "source": [
        "We will first provide the code to apply this method to the **weights**. First of all, we need to look up the original class score that the input image ```img``` obtains when fed into the network. We argued we should take the class score of digit ```n``` before the softmax function is applied, so this corresponds to layer 10 as per how we built the MNIST architecture."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "h7JfNeZD3QlS",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "layer_outputs = [layer.output for layer in model.layers] \n",
        "activation_model = Model(inputs=model.input, outputs=layer_outputs)\n",
        "activations = activation_model.predict(img.reshape(1, 1, 28, 28\n",
        "                                                   \n",
        "original_score = activations[10][0][n]\n",
        "print('Original score: ', original_score)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_IShxCqd3eST",
        "colab_type": "text"
      },
      "source": [
        "Now we are going to modify the weights of channel ```z``` of the first convolutional layer by adding 0.1 times its standard deviation. We use ```i``` and ```j``` up until 3 because the kernel size that we used for this convolutional layer is 3x3. Take into account that it is necessary to first examine the shape of the weights that one wants to modify in order to avoid mistakes. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "skKTjUjm3wdQ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "weights = model.layers[0].get_weights()\n",
        "\n",
        "desvi = np.std(activations[1][0, z, :, :])\n",
        "for i in range(0, 3):\n",
        "    for j in range(0, 3):\n",
        "        weights[0][i, j, 0, z] += 0.1*desvi"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OFRLyFDf3-nw",
        "colab_type": "text"
      },
      "source": [
        "We set up these new weights and compute the new class score for image ```img```:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rLHxYQp94kko",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "model.layers[0].set_weights(weights) \n",
        "\n",
        "layer_outputs = [layer.output for layer in model.layers] \n",
        "activation_model = Model(inputs=model.input, outputs=layer_outputs)\n",
        "activations = activation_model.predict(img.reshape(1, 1, 28, 28))\n",
        "\n",
        "new_score = activations[10][0][n]\n",
        "print('New score: ', new_score)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FKIk8Z5w436e",
        "colab_type": "text"
      },
      "source": [
        "Finally, the change that this modification had on the final class score is simply the magnitude of the change of the two scores:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lxKEomNm4-zl",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "change = abs(new_score - original_score)\n",
        "print('Change: ', change)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MqHWnT3p5E9q",
        "colab_type": "text"
      },
      "source": [
        "If you run this code in a loop, remember to load the original model before every time. Otherwise, since we are using the function ```set_weights```, we would be keeping the weight changes that we applied to the previous neurons. Finally, one can also compute the test accuracy of the network instead of the change in the final class scores (as we did with the layers in section 3.5.), but most likely a small variation in one particular hidden unit will not be enough to change the classification between classes, and so the test accuracy will not change."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "F9ZKjioz5_Sx",
        "colab_type": "text"
      },
      "source": [
        "If we want to apply the gradient method to the **activations** instead of the weights we will need a different code. However, it will take much more time to obtain results with many input images using the activations than using the weights. Why? For the weights, we need to change them once, and then the network is forever fixed. If we wanted to record the change between the original and the modified class score for a large group of input images, we still would only need to set the weights only once: feed all the input images into the original network, then change the weights, and feed them all again. However, recall that activations are not intrinsic to the network – each image produces its corresponding activations. For this reason, for *every* input image we need to stop the process inside the network, modify the activations, and finalize. Moreover, the only way to do so is to build a new Keras functional Model every time, as we will now code."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yKXxuSgRiz5K",
        "colab_type": "text"
      },
      "source": [
        "First, we feed an input image ```img``` into the network and look at the corresponding class score of digit ```n``` before the softmax function (layer 10 in this case):"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6L1fBWuL3Rcn",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "activations = activation_model.predict(img.reshape(1, 1, 28, 28))\n",
        "score = activations[10][0][n] \n",
        "print('Original score: ', activations[10][0][n])"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cUgKRQCZjZSK",
        "colab_type": "text"
      },
      "source": [
        "Following this, we apply a small change to the activations that the image has ignited in channel ```z```:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "wyABJ1qLkZ0m",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "activations = activation_model.predict(img.reshape(1, 1, 28, 28))\n",
        "desvi = np.std(activations[1][0, z, :, :])\n",
        "for i in range(0, 26):\n",
        "    for j in range(0, 26):\n",
        "        activations[1][0, z, i, j] += desvi*0.1"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1nNjW4AhkYi0",
        "colab_type": "text"
      },
      "source": [
        "We then have to build a new **Keras Model** to fix these new activations. To do so, we will create a model in which the input is not an image, but the first convolutional layer of our original model (which is the index ```idx``` refers to). It is vital to keep the correct input shapes that we used to build the network, and so we must specify that the new model will have inputs with shape (32, 26, 26), which is the shape of the first convolutional layer of the original model. The (0) that appears in the second line of the code is needed for structure purposes (the so-called nodes of the layer). As explained in the documentation of Keras' functional API, when a layer on some input is called, a new tensor is created and a *node* is added to the layer to link the input tensor to the output tensor, which is why we need to add the 0."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "elxg85lMNP3q",
        "colab_type": "text"
      },
      "source": [
        "This new model will have the same layers as before except without the final softmax layer. We stop at the second to last layer because we want this new model to output the final class score of the digit before the softmax function. This is why in the code we add all the layers in the new network up until and not including the 11th layer. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "U1lppFJNknHm",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "idx = 1 \n",
        "input_shape = model.layers[idx].get_input_shape_at(0) \n",
        "\n",
        "layer_input = Input(shape=(32, 26, 26)) \n",
        "\n",
        "x = layer_input\n",
        "for layer in model.layers[1:11]: \n",
        "    x = layer(x) \n",
        "\n",
        "new_model = Model(layer_input, x)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "G9WzHT42mSQi",
        "colab_type": "text"
      },
      "source": [
        "We then feed the new activations of the first layer into the network (where the activations of 31 channels remain the same ones except the one channel that we are studying and that therefore has been added 0.1 times its standard deviation) and we look up the new class score of the particular digit:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "iwbphyHPnEOB",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "new_activations = new_model.predict(activations[1]) \n",
        "print('New score: ', new_activations[0][n])"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Y8ohpUiF3R2I",
        "colab_type": "text"
      },
      "source": [
        "Finally, as in the case of the weights, we compute the change that this variation on the neuron has induced by computing the absolute value between the original and the final class scores. Moreover, this gradient method applied either to the weights or the activations leaves room for a lot of creativity. For example, one could analyze how a variation on a particular neuron affects each of the 10 final class scores to examine if a neuron is more relevant to certain class scores than others. However, this method has some flaws that we will examine at the end of the section."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_yTSR7l-n7O6",
        "colab_type": "text"
      },
      "source": [
        "#### 4.2.2. Killing"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ckfWj3Zg3TcN",
        "colab_type": "text"
      },
      "source": [
        "The second variation method that we present is the **killing method**. Instead of adding a small quantity to a particular channel to see how much it affects the final class score of the digit, we directly kill this particular channel. By killing we mean that we set the activation matrix to 0 (so all the elements in the 26x26 matrix are set to 0), therefore effectively turning this channel into a **dead neuron**, and then look up the new class score of the digit to see how much it varied. In spirit, this method follows the idea of **occlusion** that some papers applied to the input image to see how important each pixel was (Fong & Vedaldi, 2017): compute the class score of the original input image, occlude a small subset of the pixels, and compute again the class score to see how much it varied. \n",
        "\n",
        "This killing method can also be applied to either the weights or the activations – however, in this case, the results should be the same. In both cases, the code is the same as the one for the gradient method but instead of adding 0.1 times the standard deviation to the matrix, we do:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "A9bYcZcn3Ufc",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "for i in range(0, 26):\n",
        "    for j in range(0, 26):\n",
        "        activations[1][0, z, i, j] = 0"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Zc9pUqPtprS_",
        "colab_type": "text"
      },
      "source": [
        "#### 4.2.3. Noise"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9O1CbKN93aCi",
        "colab_type": "text"
      },
      "source": [
        "Finally, the last variation method is the **noise method**. This one is similar to the gradient method in spirit, but we use the idea of **perturbing** that some papers use on the input images (Zeiler & Fergus, 2014), although we apply these perturbations to the hidden channels instead. So instead of adding a small quantity to each particular neuron, we perturb them. Perturbing means that we add random noise to the channel. There is a Numpy function that we can use to generate this random noise, ```np.random.normal```, and then we simply add it to the activation matrix of the particular channel that we are studying. This noise can we computed in different ways, but statistically we recommend using random values that are part of a Normal (Gaussian) distribution. This is why we use ```np.random.normal``` instead of ```np.random```. The mean and standard deviation for this Normal distribution are going to be the mean and standard deviation of the particular activation matrix. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "l0bf15CF3a4W",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "amean = np.mean(activations[1][0, z, :, :])\n",
        "astd = np.std(activations[1][0, z, :, :])\n",
        "for i in range(0, 26):\n",
        "    noise = np.random.normal(amean, astd, 26)\n",
        "    activations[1][0, z, i, :] += noise"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KlA0erJZrK30",
        "colab_type": "text"
      },
      "source": [
        "After this step, similar to the gradient and killing methods, we use the same code to create a new model and get the new final class score to compute the difference between the original and the new one."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Q_pYcmNo_1XG",
        "colab_type": "text"
      },
      "source": [
        "#### 4.2.4. Some Considerations"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3K2-hitFNkka",
        "colab_type": "text"
      },
      "source": [
        "As we suggested in the introduction, we believe that the gradient, killing, and noise methods are valuable when studying the relationship between a particular hidden unit and the final class score of a digit. They offer insight into the behavior of each hidden unit of the network and are able to quantify their functioning. However, while they might be used to study a particular unit, we do not recommend using these methods to explicitly rank the hidden units and draw comparisons between them. We point out some of the reasons."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3rgX_546MgpS",
        "colab_type": "text"
      },
      "source": [
        "In the gradient method, different rankings may emerge depending on the small quantity that we add. A percentual quantity that depends on the magnitude of the original weights or activations appears to be more logical than a fixed quantity, but the choice of such quantity might have a great impact on the final results. For future work, we believe it would be better to actually use derivatives and back propagation for the gradient method (as it is done to produce saliency maps). A similar problem occurs with the noise method, for which we have to choose the range of values from where we draw the random numbers that we will add to the neuron. Moreover, the gradient and killing method give very little information about neurons that have low original activations and weights, and in general the variation of the final class score tends to be quite small and might not be significant enough to draw conclusions. Finally, the magnitude of the variation of the final class score is unbounded and can turn out to be any number. This **non-standardization** makes it hard to compare units of different layers, given that (as we saw before) layers also have varying importance between them and therefore the output changes cannot be compared across. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HWDY3qo-NmvK",
        "colab_type": "text"
      },
      "source": [
        "In the next section we propose PCACE, a reliable statistical method that can be used to both study particular neurons and, most importantly, properly rank the hidden units in order of importance, as we wanted to achieve with our original research question. PCACE assigns a number between -1 and 1 to each hidden unit, which corresponds to the maximal correlation coefficient between that hidden unit and the corresponding final class score. Therefore, it allows us to make comparisons between the PCACE values between different units and of different layers. Furthermore, the use of PCA allows the use of many input images without the algorithm becoming slow and inefficient, which was not possible in the gradient, killing, and noise methods when using the activations instead of the weights. The following section digs into the theory behind PCACE and its applications. However, as indicated before, the gradient, killing, and noise methods are still valuable tools for the study of the interpretability of neural networks and are backed-up by the literature since these ideas have been applied to the pixels of the input images."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Ik3G0pcz3eIY",
        "colab_type": "text"
      },
      "source": [
        "## 5. Statistical Methods and PCACE"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QTwGcCBcxiOO",
        "colab_type": "text"
      },
      "source": [
        "### 5.1. ACE"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fIJnY_RMNp1l",
        "colab_type": "text"
      },
      "source": [
        "This section introduces and develops a novel statistical method that can rank the hidden units of any neural network in order of importance towards the final class score: **PCACE**. The name PCACE results from the fusion of **PCA** (Principal Component Analysis), that we use for dimensionality reduction, and **ACE** (Alternating Condition Expectation), that we use to compute the maximal correlation coefficient between a dependent variable (the final class score) and multiple independent variables (each element of the activation matrix of a particular hidden unit). Therefore, we will **rank the units based on the strength of the non-linear relationship to the final class score**. We will first explain the ACE algorithm and then build PCA on top of it. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "39jD-t-xUDwX",
        "colab_type": "text"
      },
      "source": [
        "The **Alternating Conditional Expectation** (ACE) algorithm is a non-parametric approach for estimating the transformations that lead to the **maximal multiple correlation** of a response and a set of independent variables in regression and correlation analysis. This iterative algorithm was published in 1985 by Leo Breiman and Jerome H. Friedman (Breiman & Friedman, 1985), and we can apply it on our code by installing the ACE library that is readily available for Python. ACE estimates optimal transformations for multiple regression that find the maximal correlation between multiple independent variables *X* and a dependent variable *Y*. These transformations minimize the unexplained variance of a linear relationship between the transformed response variable and the sum of the transformed predictor variables (Wang & Murphy, 2005). Moreover, ACE does not require any assumptions of the response or predictor variables. \n",
        "\n",
        "Since ACE will assign a **number between -1 and 1** for each hidden unit, it allows us to rank all the hidden units of the network in a reliable, normalized, and clear way. The higher the ACE value is, the stronger is the relationship between that channel and the class score of a particular digit. This is exactly what we were looking for."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DL999UViXzyv",
        "colab_type": "text"
      },
      "source": [
        "Let us now get more concrete by applying ACE to our MNIST architecture. We will explain and code how to compute the maximal correlation coefficient between one concrete hidden neuron in the network and the output result (the class score of a particular digit). Therefore, we first pick a hidden unit ```u``` and one of the 10 possible digits ```d```. Recall that each input image that we feed into the network is going to produce an activation matrix for each of the hidden neurons and a final class score (the one that is assigned to digit ```d```). Each of the elements in the activation matrix is a **predictor variable** *X*, and the final class score of the digit is the **response variable** *Y*. In the MNIST architecture, the size of the activation matrix of the hidden units of the first layer is 26x26, which means that we have 676 predictor variables. Therefore, for each input image, we record the activation value that each of the elements in the activation matrix achieves, store it, and record the final class score. This means that if we feed ```m``` input images into the network, we will end up with 676 vectors, each of length ```m```, for the predictor variables, and one vector of length ```m``` for the response variable. However, per the implementation of the ACE library in Python, we need to input a matrix *X* in which each row is a predictor variable, which means that we will end up with a 676x```m``` matrix (we just put the 676 vectors together in a matrix as columns) and a vector of length ```m```. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "kPZKufOSaaZE"
      },
      "source": [
        "Let us try to implement this method. We are going to use 500 input images of each particular digit from the training dataset. Per the ACE Python documentation, we need to put the 676 vectors for each predictor variable into a single matrix *X* of dimensions 676x500 (that is, the number of rows in the matrix corresponds to the number of predictor variables, and the number of columns in the matrix corresponds to the number of input images), and *Y* is a vector of length 500. We will now provide the code to implement ACE to rank the channels of the first layer of the MNIST architecture by order of importance. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "C8e34g23NsmZ",
        "colab_type": "text"
      },
      "source": [
        "By previous code, recall that ```digit2``` contains 500 training images of a particular digit ```n```. For each of the 500 input images that we will feed into the network, the matrix  ```channelact``` is going to store the activations of the 676 activations of a particular channel ```i``` and ```classscore``` is a vector that is going to store the final class score of each input image. Recall that the matrix *X* must have the number of rows equal to the number of predictor variables, meaning that for each input image to the network we are going to fill another column of the matrix. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "H2O5IQVd4nJo",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "w, h = iterations, 676;  \n",
        "channelact = [[0 for x in range(w)] for y in range(h)] \n",
        "classscore = []\n",
        "\n",
        "for j in range(0, len(digit2)):\n",
        "    img = digit2[j]\n",
        "    activations = activation_model.predict(img.reshape(1, 1, 28, 28))  \n",
        "    classscore.append(activations[10][0][n])\n",
        "\n",
        "    s = 0\n",
        "    for k in range(0, 26):\n",
        "        for z in range(0, 26):\n",
        "            channelact[s][j] = activations[1][0, i, k, z]\n",
        "            s += 1"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DYmGbqsV5dSV",
        "colab_type": "text"
      },
      "source": [
        "Next we use the ACE Python library to compute the maximal correlation coefficient between channel ```i``` and digit ```n```. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8rqBRQB46BFV",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "x = channelact\n",
        "y = classscore\n",
        "\n",
        "myace = ace.ACESolver()\n",
        "myace.specify_data_set(x, y)\n",
        "myace.solve()\n",
        "\n",
        "val = np.corrcoef(myace.x_transforms[0], myace.y_transform)[0,1]\n",
        "print('Maximal correlation coefficient: ', val)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gqXF05nn6Uwz",
        "colab_type": "text"
      },
      "source": [
        "However, this method has two flaws: first, the **dimensions** of both the matrix *X* and the vector *Y* might be too large for ACE to be able to compute the maximal correlation coefficient in a reasonable amount of time. In our case the channel's activations have size 26x26 but with architectures with large input images this might be quite larger. The second problem has to do with  **```nan``` values** encountered while running the ACE algorithm. ACE is going to face problems if the values in matrix *X* are too small (because at one point of the iteration it attempts to divide by a small number and thus obtains a  ```nan```). One way to fix this is to standardize each row of the matrix *X* (so standardize each of the 676 predictor variable vectors): center to the mean and componentwise scale to unit variance. We can use the function ```preprocess.scale``` form the Sklearn library to do exactly this. However, we still face problems if the standard deviation of any of the 676 predictor variables is 0, because ACE will encounter again some  ```nan``` problems. Therefore, we are going to delete all rows whose standard deviation is 0. Here is the code to modify the matrix *X* to avoid the  ```nan``` problems (```channelact``` is the 676x200 matrix that we built in the previous code snippet):"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "NTYqgrKz80H8",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "final_act = []\n",
        "for r in range(0, 676):\n",
        "    act_normalized[r] = preprocessing.scale(channelact[r])\n",
        "    if (np.std(act_normalized[r]) != 0):\n",
        "        final_act.append(actnormalized[r])"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4pKJOMB_8_Tb",
        "colab_type": "text"
      },
      "source": [
        "Finally, we call ACE with ```x = final_act``` instead of ```channelact```. Nevertheless, while we applied these pre-processing methods to try to avoid ```nan``` problems when running ACE, it might be that this solution does not suffice and one still encounters ```nan``` issues. This will be quite evident to the user because the ACE library for Python will show a warning message. We recommend trying to use other ways of pre-processing the data to avoid such problems. Or, even better, use the PCACE method that we present in the next section, which solves both the problem of large dimensionalities and the problem of ```nan``` values."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BQ8zK8UICixs",
        "colab_type": "text"
      },
      "source": [
        "### 5.2. PCACE"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gVvKfPDGNxYb",
        "colab_type": "text"
      },
      "source": [
        "This section presents the **statistical method** that we encourage to use in order to rank channels in order of importance: **PCACE**. The PCACE method consists of the using the **Principal Component Analysis** statistical procedure for **dimensionality reduction** on top of the **ACE algorithm** to reduce the dimensionality of the *X* matrix while maintaining most of the relevant information. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4mT_huY5AXom",
        "colab_type": "text"
      },
      "source": [
        "Recall that our main goal with using PCA before ACE is to reduce the dimensionality of the predictor variable matrix *X*. Without PCA, the dimensions of *X* are huge: the number of rows corresponds to the number of elements in the activation matrix of the neuron and the number of columns to the number of input images that we fed into the network. In the case of the MNIST architecture, the activation matrices have size 26x26 (therefore having 676 elements), but this is quite small and other networks have way larger input images. When the number of classes is also larger than 10, we need to increase the number of input images in order to obtain significant results. It is unfeasible to apply the ACE algorithm to such large matrices in a way that is **computationally efficient**. This is why we turn to the Principal Component Analysis method: to reduce the dimensionality of the *X* predictor variables matrix while keeping as much information as possible. PCA performs a linear mapping of the data in the original matrix to a lower-dimensional space such that the variance of the data in the low-dimensional representation is maximized by using eigenvalue decomposition (Jolliffe, 2014). Therefore, after applying PCA, our new matrix *X* will have smaller dimensions but will have inherited the maximum possible variance of the original data. Once the dimension of the original matrix *X* has been reduced, we can now apply ACE in a computationally efficient way. \n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OpwN4WwWA_vr",
        "colab_type": "text"
      },
      "source": [
        "We believe that PCACE is a method that can be effectively used to rank hidden channels for several reasons:\n",
        "\n",
        "* It brings together two beautiful, effective, reliable, and powerful statistical methods: the Alternating Conditional Expectation (ACE) and the Principal Component Analysis (PCA). \n",
        "* By using ACE on top of PCA, we can significantly reduce the dimensions of matrix *X*, therefore making the ranking computation much more time efficient.\n",
        "* Since the number of predictor variables is significantly reduced, we can increment the number of input images and PCACE will still be computationally efficient.\n",
        "* The ```nan``` problems that we discussed in the previous section are also solved by using PCA before applying ACE because PCA already standardizes the original data. \n",
        "* The PCACE values are all between -1 and 1 and therefore it is a better method to make comparisons across different neurons and from different layers.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CY6Q57zrD3JQ",
        "colab_type": "text"
      },
      "source": [
        "Therefore, in order to rank the hidden channels of a layer of a neural network, we propose the following **PCACE algorithm:**\n",
        "```\n",
        "> For each channel c in the layer, do:\n",
        "      > For each input image of the training set fed into the network, do:\n",
        "           > Store the activations of channel c into matrix X\n",
        "           > Store the final class score into vector Y\n",
        "           \n",
        "      > Apply PCA to reduce the dimensionality of matrix X\n",
        "      > Apply ACE to compute the maximal correlation coefficient between the newly reduced matrix X' and vector Y\n",
        "\n",
        "> Sort all the ACE values\n",
        "```\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "U2nXahp0Fbtx",
        "colab_type": "text"
      },
      "source": [
        "Of course, we do not have to focus on the hidden neurons in only one particular layer and we can compute the ACE value of all the neurons in the network to then rank them. In our case, as we discussed before, our original *X* matrix has size 676x500 (we use 500 input images), and we are going to use PCA to reduce this matrix down to dimension 10x500. It is important to remark that when using PCA for the purposes of PCACE we want to reduce the number of rows (the number of predictor variables) and not the columns (the number of input images), even if it is not necessarily 10. We could reduce the dimensions a bit more if we find that ACE if still not efficient after this reduction."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Kt1FcT_BGZh8",
        "colab_type": "text"
      },
      "source": [
        "Finally, we are going to provide the **code to apply PCACE** to our neural network. Even if some parts of the following code have been introduced, we wanted to provide the full PCACE code to obtain the ranking of the 32 channels of the first layer of the MNIST architecture for every one of the 10 digits. Therefore, the following code is going to produce 10 rankings of the channels: one for each of the 10 digits. One final but important remark of the PCACE code has to do with one small technical incongruence between the PCA algorithm provided by the Sklearn library and the ACE algorithm provided by its corresponding Python library: when applying PCA we need to transpose the *X* matrix (because otherwise it is going to reduce the number of columns instead of the number of rows), and after applying PCA we need to transpose it again before applying ACE or else the dimensions between *X* and *Y* will not match. Also, note that the original data is standardized before applying PCA."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "oYQFKobzHhys",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "for n in range(0, 10):\n",
        "    iterations = 500\n",
        "    results = []\n",
        "    digit = []\n",
        "\n",
        "    for i in range (0, len(y_train)):\n",
        "        if y_train[i] == n: \n",
        "            digit.append(X_train[i])\n",
        "\n",
        "    digit2 = []\n",
        "    for i in range(0, iterations):\n",
        "        digit2.append(digit[i])\n",
        "\n",
        "    for i in range(0, 32):\n",
        "        w, h = iterations, 676;  \n",
        "        channelact = [[0 for x in range(w)] for y in range(h)] \n",
        "        classscore = []\n",
        "\n",
        "        for j in range(0, len(digit2)):\n",
        "            img = digit2[j]\n",
        "            activations = activation_model.predict(img.reshape\n",
        "                                                   (1, 1, 28, 28))  \n",
        "            classscore.append(activations[10][0][n])\n",
        "\n",
        "            s = 0\n",
        "            for k in range(0, 26):\n",
        "                for z in range(0, 26):\n",
        "                    channelact[s][j] = activations[1][0, i, k, z]\n",
        "                    s += 1\n",
        "\n",
        "        x = channelact\n",
        "        y = classscore\n",
        "\n",
        "        s = np.transpose(x)\n",
        "        c = s.tolist()\n",
        "        x_std = StandardScaler().fit_transform(c)\n",
        "\n",
        "        pca = PCA(n_components = 10) \n",
        "        nou_x = pca.fit_transform(x_std) \n",
        "        s = np.transpose(nou_x)\n",
        "        c = s.tolist()\n",
        "\n",
        "        x = c\n",
        "        y = classscore\n",
        "\n",
        "        myace = ace.ACESolver()\n",
        "        myace.specify_data_set(x, y)\n",
        "        myace.solve()\n",
        "\n",
        "        val = np.corrcoef(myace.x_transforms[0], myace.y_transform)[0,1]\n",
        "        results.append(val)\n",
        "\n",
        "    print('Digit ', n)    \n",
        "    print(np.argsort(results))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CyarfXNEtV5Z",
        "colab_type": "text"
      },
      "source": [
        "![alt text](https://drive.google.com/uc?id=1BnOi3ABc0TeAMf08PpufSjsdvKQzfdmp)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oKjCF7Z15Ghr",
        "colab_type": "text"
      },
      "source": [
        "### 5.3. Analysis and Comparisons"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MpJdw4UAN4Cn",
        "colab_type": "text"
      },
      "source": [
        "To end this statistical section, we are going to present two effective ways to statistically **compare the rankings** that we have obtained with PCACE for each digit. An interesting analysis that we can do between the 10 rankings that we have obtained with PCACE (one for each digit) is to compare how similar each pair of two rankings is, in order to find out if there are certain digits that share a similar order of importance. We can compare two arrays in two different ways: we either compare the sorted arrays by displaying each of the ACE values for each channel or we compare the sorting with the channel index instead of its ACE value. In the first case we would use **Pearson's coefficient** to compute how similar the two arrays are. In the second case we would use **Spearman's rank correlation coefficient** to compute their similarity."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rV1v49S2N68E",
        "colab_type": "text"
      },
      "source": [
        "The widely used Pearson's coefficient measures the linear correlation between two variables *X* and *Y* (Benesty et al., 2009), which in our case are 32 ordered ACE values for one digit for *X* and 32 ordered ACE values for another digit for *Y*. Spearman's rank correlation coefficient is a nonparametric measure of rank correlation; that is, the statistical dependence between the rankings of two variables (Zar, 1972). In our case, these two variables are the indices of 32 neurons ordered by their ACE value for one digit and 32 neurons ordered by their ACE value for another digit. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9r_pYGhzIaCt",
        "colab_type": "text"
      },
      "source": [
        "Since one might want to compare rankings of the channels that have been obtained through different methods, we believe that in this case it is better to use Spearman's rank correlation and have each array display the index of the channel (in our case, the integers from 0 to 31), so that the values sorted by each method are not taken into account. However, it is true that we might be losing some information by not using the original values that we used to sort. For example, it surely is the case that the *distance* between different pairs of channels that are side to side in the ranking is different, and so they are not equally spread in terms of importance. However, it must be noted that the comparison between the values of Pearson's and Spearman's coefficients on the same sets of data is currently being studied and it is not clear which method is more suitable (Hauke & Kossowski, 2011). "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MXDF2DywN-Dc",
        "colab_type": "text"
      },
      "source": [
        "In any case, we can easily compute the Spearman's rank correlation coefficient between arrays ```a``` and ```b``` with the Scipy library:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lAI1dxq8MFEy",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "coef, p = stats.spearmanr(a, b)\n",
        "print('Spearman coefficient: ', coef)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zKUZzY8fMMEM",
        "colab_type": "text"
      },
      "source": [
        "We can use Spearman's rank correlation coefficient to display in a beautiful way how similar the 10 rankings obtained by the PCACE algorithm are by computing the correlation coefficient between all pairs of rankings and then displaying them using a heatmap. We have computed Spearman's rank correlation coefficient between every pair of the 10 rankings obtained with PCACE and stored them in the correct order in the matrix called ```uniform_data```. By correct order we mean that the value in element ```(i, j)``` and ```(j, i)``` of the matrix should be the same and equal to the correlation coefficient between the ranking of digit ```i``` and digit ```j```.  This is how we print the final **heatmap:**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "LHjx_cLYMks1",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "w, h = 10, 10;  \n",
        "uniform_data = [[0 for x in range(w)] for y in range(h)] \n",
        "    \n",
        "for i in range(0, 10):\n",
        "    for j in range(0, 10):\n",
        "        coef, p = stats.spearmanr(taula[i], taula[j])\n",
        "        uniform_data[i][j] = coef\n",
        "        \n",
        "ax = sns.heatmap(uniform_data, linewidth=0.2, cmap=\"YlGnBu\", annot=True)\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "l57qBe1_PPIO",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1Kxd7rqj_povRwrDtJE08dspGnUfb72uC' width=360px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qVP5e_LpPTU5",
        "colab_type": "text"
      },
      "source": [
        "Finally, another nice way to visualize the rankings that we have obtained with PCACE for the channels of the first layer is to pick the **5 top channels** for each ranking and visualize the activation of the corresponding digit in three different ways as we did in section 4.2.1. That is, for each digit we pick the 5 most important channels according to the PCACE classification and we visualize these 5 channels with **3 different input images:** the ideal image, the mean real image, and a random image of the training set for each digit. This is what we obtain: "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "El-CG1M2yR62",
        "colab_type": "text"
      },
      "source": [
        "![alt text](https://drive.google.com/uc?id=1usUV-lZMBStOFsm9PXJ0D2ku2ysv-xes)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "an9nrZqDyIkO",
        "colab_type": "text"
      },
      "source": [
        "![alt text](https://drive.google.com/uc?id=1OrxdI-lK5wAb0SLidU-xiIxyQRZqOObX)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KSuGaFYChw69",
        "colab_type": "text"
      },
      "source": [
        "If we analyze the PCACE rankings carefully, we can make the following observations:\n",
        "\n",
        "* **Many of the channels that PCACE ranks as the most important for each digit are precisely the ones that we characterized as *reverse* channels** in section 3.5.; that is, the channels that focus on the surroundings and details of the digit instead of the global shape. Recall that the reverse channels activate at the edges, corners, and surroundings of the digit when the input digit is a real image from the training set (or a mean of the images in the training set), but is nearly dead when the input image is an ideal digit created with the activation maximization method. As a side note, the reverse channels are numerically characterized by having a negative weight mean, and thus, as we saw in section 3.2., the reverse channels of the first layer in this MNIST architecture are 5, 10, 12, 13, 15, 22, 24, 25 and 26. However, as we already observed in section 3.5., when the input image corresponds to the **ideal image** of the digit, the reverse channels are **not activated** at all or barely activated. This indicates a flaw of the activation maximization method: the artificial image produced by the activation maximization does not activate most of the channels that are found to be the most important for the correct classification of the digit by the PCACE algorithm. This is because these channels focus on what is *not* the digit and are therefore not fired because the pixels that they focus on are not maximized by the backpropagation used in the activation maximization method.\n",
        "\n",
        "* If we examine the channels that are ranked as top 5 important by the PCACE algorithm, we see that all digits except 0 have one or several reverse channels in their top 5 classification. To make this observation more clear, we have put together the rank of the 32 channels for each digit as produced by the PCACE algorithm and we have marked in blue all the reverse channels and circled in purple the top 5 channels for each digit (recall that the array is written in ascending order of importance). This is the result:"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NOPothxsU_Ua",
        "colab_type": "text"
      },
      "source": [
        "<img src='https://drive.google.com/uc?id=1SuXRwJkzWXDTJL5lNSXypKNeY2BFAR1V' width=900px>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3JDNyQoRVNo7",
        "colab_type": "text"
      },
      "source": [
        "* Indeed, all digits except 0 have at least one reverse channel in the top 5 classification. Digit 6 has only one reverse channel; digits 4 and 7 have two; digits 8 and 9 have three; digits 1, 2, and 5 have four, and digit 3 has five reverse channels on the top 5 classification. Therefore, we can conclude that in general reverse channels obtain high importance scores (in particular, 5, 12 and 15 are the ones that are most often in the top 5 classification). \n",
        "\n",
        "* It can also be observed that the importance of each channel varies across the digits. For example, channel 22 is the least important for digit 9 according to the PCACE ranking, whereas it is the most important for digit 8. However, other channels have more consistent positions, and as we just mentioned channels 5, 12, and 15 tend to be in high positions (although this is not true for all digits). This indicates that the PCACE algorithm should be applied to each class separately, and that we should be wary if we try to obtain a **global ranking** of the hidden units instead of creating a separate ranking for each digit. We believe that the difference between the top 5 channels for each score will depend on each layer and architecture. A way to compute the difference between rankings is to use Pearson's or Spearman's coefficients, as we explained before in this section. To compute the global top 5 most important channels, we would then compute the mean PCACE value of each channel across all the different classes and then obtain the final global top 5 classification.\n",
        "\n",
        "* Another very important conclusion that we extract from applying the PCACE algorithm to the first layer of the MNIST architecture is that it highlights a **weakness of the activation maximization method**. As we explained before, the ideal digit for each class that is artificially created by using the activation maximization method does **not** fire the reverse channels. That is, when the ideal digit is fed into the network, the ideal channels are dead or nearly dead. From this fact we would infer that reverse channels are the least important channels towards the final class classification of the input digit, because when the digit that maximally activates its corresponding class score is fed into the network these channels are not activated and therefore do not contribute to the final score. However, and surprisingly, the PCACE classification shows that these reverse channels are in fact very important towards the final score of the digit and in most cases are the ones that demonstrate the strongest relationship between their activations and the final class score of the digit. Therefore, we conclude that this method fails in artificially maximizing parts of the ideal input image that will fire the reverse channels, and this is necessary since they are generally the ones that matter the most towards the final classification. We hypothesize that the current activation maximization method does not fire the reverse channels because it creates an artificial digit focusing on the big picture and shape of the digit. However, because the reverse channels focus on the **surroundings of the digit**, these are not taken care of when creating the ideal image and thus the reverse channels are not fired. Most probably, the parts of the digits that the reverse channels focus on are not activated by the activation maximization method because of the negative weights of the reverse channels. The activation maximization method tries to increase maximizations and thus does not try to maximize parts that would obtain a low maximization due to the negative sign of the weights of the reverse channels.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "03uBHQ7W6ota",
        "colab_type": "text"
      },
      "source": [
        "## 6. Conclusions"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oIf3cFIfOEPD",
        "colab_type": "text"
      },
      "source": [
        "In this tutorial and research work we have delved into the inner functioning of convolutional neural networks to try to open these black boxes a little bit. First, we made a **bibliographical search** of existing methods and tools for the interpretability of CNNs. We mentioned the current existence of saliency maps, grad-CAM, occluding and perturbing methods, deconvolutional networks, activation maximization method, and the development of human-understandable visualizations, among others, that researchers have created during the past years in order to help us understand and interpret what is going on inside neural networks. Then we provided the **code** for the MNIST architecture that we used throughout the tutorial along with some **basic visualization tools** that can be found in the Keras documentation: how to get the weights and activations of a layer or unit, how to produce a saliency map, the grad-CAM method and the activation maximization method. We also used the weights of the layers to investigate how much each layer influenced the final class score. Finally, we showed that this code was easily applicable to other networks such as the Cifar-10 architecture or any that can be found in the Keras documentation. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7IRaxyhtOLBS",
        "colab_type": "text"
      },
      "source": [
        "Next we focused on our research topic: the **importance of the hidden units** of the neural network. First, we **visualized the activation matrices** for all the 32 channels in the first convolutional layer by using three different input images for each digit: a random image of the training dataset, the mean image of the digit built from the training dataset, and the ideal image generated with the activation maximization method. We also fed some randomly generated images and training images with noise into the network to ensure that the channels of the CNN were indeed location invariant."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_2U4C9P3OyEL",
        "colab_type": "text"
      },
      "source": [
        "By visualizing the activations of the 32 channels, we observed that there were some **reverse channels** that **focused on what was *not* the digit**; either the edges of the digit or everything in the image *but* the digit. When visualizing the activations of the reverse channels with images of the training dataset, it was clear that they were activating with parts of the image that were not the digit itself. However, these channels appeared to be **dead** or nearly dead (not activated) when the input image was the one artificially **created with the activation maximization** method, pointing out that the activation maximization method was disregarding the necessary activation of the reverse channels. Moreover, we saw that the reverse channels were precisely the ones whose **mean weight** was negative. We found this to be the explanation to why the reverse channels were focusing on everything that was not the digit – the negative weights made them focus on the *negative* parts of the image, meaning the pixels in the image that were not the digit itself. Or, another way to put it, the negative weights **inverted** what the channels regularly focused on (the digit itself)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lfkYbZDDObJP",
        "colab_type": "text"
      },
      "source": [
        "We then used these tools and applied some of the ideas found in the literature to rank the pixels of the input image according to how much they matter for the final class score classification, but we applied these ideas to the hidden units instead. We explained and provided the code for the **gradient, killing, and noise methods**, all applied to both the weights and the activations of the channels. While these methods examine the relevance of a particular hidden unit for the final class score and are valuable interpretability tools, we believed it was better to use a statistical method to purposely rank the hidden units in order of importance."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_oFNyrmcOdTe",
        "colab_type": "text"
      },
      "source": [
        "This is how we got to the main and final part of the tutorial, by presenting the **PCACE algorithm**. PCACE is a statistical method that puts together the Principal Component Analysis for dimensionality reduction algorithm and the Alternating Conditional Expectations algorithm to find the maximal correlation coefficient between a hidden neuron and the final class score. PCACE returns a number between -1 and 1 that indicates the strength of the non-linear relationship between the multiple predictor variables (each element in the activation matrix of a particular channel when feeding multiple input images into the network) and the response variable (each correct final class score before the softmax function when feeding multiple input images into the network). We then use the different PCACE values to rank the hidden units in order of importance."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FeVJzPfbOgie",
        "colab_type": "text"
      },
      "source": [
        "We used the results we obtained from applying PCACE to the MNIST architecture to point out the **flaws of the activation maximization method:** PCACE indicated that the reverse channels were usually among the most important towards the correct class classification, yet they are not activated at all or barely activated by the ideal image digit that is artificially generated with the activation maximization method. The activation maximization method generates an image focusing on activating the pixels that correspond to the digit, but they neglect the parts of the image that are *not* the digit, and the activation visualizations that we made indicate that the reverse channels learn and focus on precisely the parts of the image that do not belong to the digit. We believe that the mean negative weights of the reverse channels explain why these channels were focusing on elements of the image that were not on the digit itself."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yXAG--usOjVY",
        "colab_type": "text"
      },
      "source": [
        "In this conclusion we must also acknowledge the **flaws of our methods** and suggest **future work** on this area. First, for all our codes, methods, and conclusions we have always used the MNIST architecture, which is the **simplest** one and therefore is surely neglecting several nuances of the network mechanism because we are working with a dataset and an architecture that are fairly simple. For example, all the digits in the MNIST training dataset are centered and have more or less the same shape. Therefore, for example, it might be much more complicated in other networks to define the reverse channels and to visualize what they are focusing on since there is much less uniformity and hegemony in the dataset. Thus, we believe that the focusing on the negative weights is a much more reliable way of identifying the reverse channels instead of doing it visually. It would be interesting to investigate if the reverse channels always focus on the surroundings of the important part of the image towards the classification in all architectures and datasets. However, it should also be noted that using the mean of the weights might not be too informative if the dimensions of the matrix are large – in that case, other parameters should be taken into account, such as the standard deviation. In our case, the weight matrix of the channels of the first layer had size 3x3, which meant it was a mean of only nine elements."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LE6TDK8uOleI",
        "colab_type": "text"
      },
      "source": [
        "Moreover, we have developed and tested our original methods on only the **first layer** of the MNIST architecture. In future work, we would replicate these methods on other layers to see if there are any substantial changes. Another variation that could be applied on the architecture to replicate our methods would be to change the activation function that we use for each layer. In our MNIST architecture we have always applied the **ReLu function**, which kills the inputs to the function that are negative (because it outputs *max(0, x)*, where *x* is the input). Since we have carefully examined the relationship between reverse channels, negative weights, activations, and dead neurons, it would be helpful to analyze the role of the ReLu function in all these dynamics. The choice of the activation function plays a fundamental role in the inner values of the neural network and thus in the future we should be able to better explain how the role of the reverse channels is affected by the choice of the activation function."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "51I8s3QaOnqp",
        "colab_type": "text"
      },
      "source": [
        "We could have also explored and used current existing visualization methods much more for our research purposes. For example, we showed how to obtain saliency and grad-CAM maps but we did not **incorporate these methods** into our research. We could have produced a saliency map of the ideal digit once it is fed into the network to produce a heatmap of the most important pixels and compare it to the one obtained with a real image from the training set. We could have also used the activation maximization method to create the ideal input image that would maximize a reverse channel and another to maximize a non-reverse channel to study how they differ. Many other ideas can be applied to further interpret and visualize what the hidden units are focusing on. Furthermore, we could also have studied in more depth the ideal digits produced with the activation maximization method and have manually explored the activation maximization method a bit more. Finally, we already analyzed the flaws of the gradient, killing, and noise methods in section 4.2.4."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Oj3afO1OOqOw",
        "colab_type": "text"
      },
      "source": [
        "It could be questioned ***how* we quantified the importance of a neuron** in the first place. Regarding PCACE, it is questionable whether the maximal correlation coefficient between a neuron and a class score is what should be the measure of importantness of the unit. One might argue that there are other ways of interpreting what the *importance* of a neuron means. Maybe a neuron is important for a certain class score but is at the same time harmful for another class score. Or maybe the fact that a neuron is not highly *correlated* with the final class score does not rule out the possibility that this neuron might be important in other ways. For example, just visually, we would probably have guessed that a channel that gets activated with the digit itself is clearly more important than a channel that focuses on the edges of the digit, and PCACE indicated that this is not true in most cases. Another possible criticism would be that we cannot break down the functioning of the network to each particular neuron because they have been found to also work ***in groups*** (Fong & Veldadi, 2018; Nguyen et al., 2016). Others might argue that while it is valid to compute the importantness of a neuron, it might not make the most sense to try to compare them and even less to compare neurons of different layers. Maybe the work done by neurons from distinct layers is just different in essence and cannot be compared across. Or, if we take this criticism even further, maybe it does not even make sense to try to interpret neural networks from a **human perspective** in the first place. Maybe neural networks are not intrinsically black boxes, but rather black boxes *to us* in the same way that humans do not have the cognitive capacity needed to understand other processes in the world."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VCeNCUVOOsmy",
        "colab_type": "text"
      },
      "source": [
        "While all these criticisms are perfectly valid, we believe that the work behind this research project and tutorial is still well-founded and justified. Most of the times humans try to quantify a phenomenon, concept, or process in the world around us there is no clear way of how this should be done, and there is debate on which is the most reliable way to do so. However, it is very important to still try to quantify the world around us to make it more understandable and interpretable. And in the case of neural networks it is even more crucial to be able to understand and **be able to explain** how these networks classify and reach certain conclusions because they are currently being used in companies, banks, and hospitals. Moreover, neural networks are not even an object that we found in nature that we are trying to decipher – *we* built these architectures. It is quite alarming to think that we currently need to be researching how these architectures work as if they were a living brain, when it is humans who came up with the idea of building neural networks in the first place, not nature. We built an architecture that works **without understanding how**."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4Pm1cIfSZZGz",
        "colab_type": "text"
      },
      "source": [
        "This is why it is important to keep doing research on this topic and to keep trying to disentangle the black boxes that are making decisions that **affect our daily life**. Some directions for future work include improving the PCACE methods; quantifying the importance of the hidden units in other ways; studying the relationship between each hidden neuron and all the final class scores instead of only the one that corresponds to the input digit; studying neurons in groups instead of individual units; making the interpretability methods more computationally efficient; further investigating why the activation maximization fails to fire the reverse channels, and applying Lucid or Deep Dream to the hidden neurons, among many others, for this is quite a new topic of study. We encourage users who have followed this tutorial to try themselves the methods and tools we described, for **there is still much to be done to open up these black boxes**.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JzA6qNlG6xj6",
        "colab_type": "text"
      },
      "source": [
        "## 7. Acknowledgements"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OfWy63ZYl6DD",
        "colab_type": "text"
      },
      "source": [
        "I would like to express my gratitude towards **Dr. Seth Flaxman** from Imperial College London for supervising this project. His guidance and vision have made these eight weeks a major enriching and fulfilling experience. I would also like to thank **Imperial College** and the **Mathematics Department** for organizing the Undergraduate Research Opportunities Program (UROP), which allows undergraduate students to carry out research during the summer. I am also grateful to **Dr. David Ham** for letting me participate in the seminars and communication trainings planned by Imperial College to encourage women pursue a PhD in Mathematics, and to **Dr. Lorin Crawford** from Brown University and **Dr. Wenjia Bai** from Imperial College for their time and comments. I would like to thank **Harvard College Office for Career Services** for sponsoring this summer research through the David Rockefeller International Experience Fund Grant. Finally, I am indebted to Isabel Diersen and Oliver Riskin-Kutz for great discussions and company."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "M2heU1Ivye7F",
        "colab_type": "text"
      },
      "source": [
        "## 8. References"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "sEyj1T_avxqd",
        "colab_type": "text"
      },
      "source": [
        "B. Alsallakh, A. Jourabloo, M. Ye, X. and Liu, L. Ren. Do Convolutional Neural Networks Learn Class Hierarchy? IEEE, 2018. [PDF](https://arxiv.org/pdf/1710.06501.pdf)\n",
        "\n",
        "D. Bau, B. Zhou, A. Khosla, A. Oliva, and A. Torralba. Network Dissection: Quantifying Interpretability of Deep Visual Representations. IEEE, 2017. [PDF](https://arxiv.org/pdf/1704.05796.pdf)\n",
        "\n",
        "J. Benesty, J. Chen, Y. Huang, I. Cohen. Pearson Correlation Coefficient. In: Noise Reduction in Speech Processing. Springer Topics in Signal Processing, vol 2. Springer, Berlin, Heidelberg, 2009. [LINK](https://link.springer.com/chapter/10.1007/978-3-642-00296-0_5)\n",
        "\n",
        "L. Breiman and J. H. Friedman. Estimating Optimal Transformations for Multiple Regression and Correlation. Journal of the American Statistical Association, Vol. 80, No. 391, pp. 580-598, 1985. [PDF](https://pdfs.semanticscholar.org/e1ee/012f12793b4021352bb953f2fe9a40c33cf2.pdf)\n",
        "\n",
        "A. Dosovitskiy and T. Brox. Inverting Visual Representations with Convolutional Networks. arXiv:1506.02753, 2015. [PDF](https://arxiv.org/pdf/1506.02753.pdf)\n",
        "\n",
        "D. Erhan, Y. Bengio, A. Courville, and V. Pascal. Visualizing Higher-Layer Features of a Deep Network. Technical report, University of Montreal, 2009. [PDF](https://pdfs.semanticscholar.org/65d9/94fb778a8d9e0f632659fb33a082949a50d3.pdf)\n",
        "\n",
        "R. Fong and A. Veldadi. Interpretable Explanations of Black Boxes by Meaningful Perturbation. IEEE, 2017. [PDF](https://arxiv.org/pdf/1704.03296.pdf)\n",
        "\n",
        "R. Fong and A. Vedaldi. Net2Vec: Quantifying and Explaining how Concepts are Encoded by FIlters in Deep Neural Networks. arXiv:1801.03454, 2018. [PDF](https://arxiv.org/pdf/1801.03454.pdf)\n",
        "\n",
        "B. Goodman and S. Flaxman. European Union regulations on algorithmic decision-making and a “right to explanation”. AI Magazine, Vol 38, No 3, 2017. [PDF](https://arxiv.org/pdf/1606.08813.pdf)\n",
        "\n",
        "J. Hauke, T. Kossowski. Comparison of Values of Pearson's and Spearman's Correlation Coefficients on the Same Sets of Data. Quaestiones geographicae 30, 2 (2011), 87. [LINK](https://content.sciendo.com/view/journals/quageo/30/2/article-p87.xml)\n",
        "\n",
        "I. Jolliffe. Principal Component Analysis. Springer, 2014. [LINK](https://www.springer.com/gp/book/9780387954424)\n",
        "\n",
        "B. Kim, M. Wattenberg, J. Gilmer, C. Cai, J. Wexler, F. Viegas, and R. Sayres. Interpretability Beyond Feature Attribution: Quantitative Testing with Concept Activation Vectors (TCAV). ICML, 2018. [LINK](https://arxiv.org/pdf/1711.11279.pdf)\n",
        "\n",
        "P. Kindermans, S. Hooker, J. Adebayo, M. Alber, K. T. Schütt, S. Dähne, D. Erhan, and B. Kim. The (Un)reliability of Saliency Methods. arXiv:1711.00867, 2017a. [PDF](https://arxiv.org/pdf/1711.00867.pdf)\n",
        "\n",
        "P. Kindermans, K. T. Schütt, M. Alber, K. Müller, D. Erhan, B. Kim, S. Dähne. Learning How to Explain Neural Networks: PatternNet and PatternAttribution. arXiv:1705.05598, 2017b. [PDF](https://arxiv.org/pdf/1705.05598.pdf)\n",
        "\n",
        "P. W. Koh and P. Liang. Understanding Black-box Predictions via Influence Functions. Proceedings of the 34 th International Conference on Machine Learning, Sydney, Australia, PMLR 70, 2017. [PDF](https://arxiv.org/pdf/1703.04730.pdf)\n",
        "\n",
        "A. Nguyen, J. Tosinski, and J. Clune. Multifaceted Feature Visualization: Uncovering the Different Types of Features Learned By Each Neuron in Deep Neural Networks. arXiv:1602.03616, 2016. [PDF](https://arxiv.org/pdf/1602.03616.pdf)\n",
        "\n",
        "C. Olah, A. Satyanarayan, I. Johnson, S. Carter, L. Schubert, K. Ye, A. Mordvintsev. Feature Visualization. Distill, 2018. [LINK](https://distill.pub/2017/feature-visualization/)\n",
        "\n",
        "R. R. Selvaraju, M. Cogswell, A. Das, R. Vedantam, D. Parikh, and D. Batra. Grad-CAM: Visual Explanations from Deep Networks via Gradient-based Localization. arXiv:1610.02391, 2016. [PDF](https://arxiv.org/pdf/1610.02391.pdf)\n",
        "\n",
        "K. Simonyan, A. Vedaldi, A. Zisserman. Deep Inside Convolutional Networks: Visualising Image Classification Models and Saliency Maps. arXiv:1409.1556, 2014. [PDF](https://arxiv.org/pdf/1312.6034.pdf)\n",
        "\n",
        "D. Wang and M. Murphy. Identifying Nonlinear Relationships in Regression using the ACE Algorithm. Journal of Applied Statistics, 32:3, 243-258, 2005. [LINK](https://www.tandfonline.com/doi/abs/10.1080/02664760500054517)\n",
        "\n",
        "J. Yosinski, J. Clune, A. Nguyen, T.  Fuchs, and H. Lipson. Understanding Neural Networks Through Deep Visualization. ICML, 2015. [PDF](https://arxiv.org/pdf/1506.06579.pdf)\n",
        "\n",
        "J. H. Zar. Significance Testing of the Spearman Rank Correlation Coefficient. Journal of the American Statistical Association, Vol. 67, No. 339, 1972. [LINK](https://www.jstor.org/stable/2284441?seq=1#metadata_info_tab_contents)\n",
        "\n",
        "M. D. Zeiler and R. Fergus. Visualizing and Understanding Convolutional Networks. ECCV, 2014. [PDF](https://arxiv.org/pdf/1311.2901.pdf)\n",
        "\n",
        "B. Zhou, A. Khosla, A. Lapedriza, A. Oliva, and A. Torralba. Object Detectors Emerge in Deep Scene CNNs. ICLR, 2015. [PDF](https://arxiv.org/pdf/1412.6856.pdf)"
      ]
    }
  ]
}