{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2.3 Layer-wise Relevance Propagation Part 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "**This section corresponds to Sections 4.3 and 5.1 of the original paper.**\n",
    "\n",
    "Recall how in Section 2.1 we defined a relevance function $R$ that takes in an ith pixel $x_i$ of image $x$ as an argument and returns a scalar value $R(x_i)$ which indicates how much positive or negative contribution $x_i$ has had on the final decision $f(x)$. Then as we progressed throughout Sections 2.1 and 2.2, we added several constraints to $R(x_i)$ such that it better described the rationale behind the DNN's decision. I will restate the constraints below so we can refer back to it later.\n",
    "\n",
    "##### Constraint 1. $R(x_i) > 0$ should indicate positive contribution and $R(x_i) < 0$ negative contribution. If the range of $R$ is limited to nonnegative numbers, $R(x_i) = 0$ indicates no contribution and $R(x_i) > 0$ indicates positive contribution (Section 2.1).\n",
    "\n",
    "##### Constraint 2. Relevance Conservation Constraint (Section 2.2): \n",
    "\n",
    "\\begin{equation}\n",
    "f(x) = \\sum^V_{i=1} R(x_i)\n",
    "\\end{equation}\n",
    "\n",
    "By adding *Constraint (1)* in the Sensitivity Analysis framework, we quantified the amount of influence individual pixel $x_i$ has on $f(x)$. However, sensitivity analysis did not provide an straightforward explanation of the function value  $f(x)$, but rather suggestions. By introducing *Contraint (2)* in the Simple Taylor Decomposition framework, we created an explicit relationship between $R(x_i)$ and $f(x)$. We aim to improve $R$ by attaching a third constraint that exploits the feed-forward graph structure of the DNN.\n",
    "\n",
    "**Notation Warning: **\n",
    "so far, we defined the relevance score $R(x_i)$ as a function that takes in pixels only. However, we are now going to generalize the concept of relevance score such that it can be applied to all layers of a DNN. The $l$th layer of a DNN is modeled as a vector of pre-activations $z = (z^{(l)}_d)^{V(l)}_{d = 1}$ with dimensionality $V(l)$. That is, $z = (z^{(l)}_1, z^{(l)}_2, z^{(l)}_3, ..., z^{(l)}_{V(l)})$ where $V(l)$ is the number of nodes/neurons in the $l$th layer. Then we can define a relevance score $R^{(l)}_d$ for each pre-activation $z^{(l)}_d$ of the vector $z$ at layer $l$. See that we have now incorporated both the layer index and the neuron index into the notation of the relevance score. With this new notation, we can indicate the relevance of any neuron at any layer, even including the input layer! For example, the previous notation $R(x_i)$ is now equivalent to $R_i^{(1)}$ for $i = 1, 2, 3, ..., V(1)$. They both indicate the relevance of $i$th neuron at the first layer.\n",
    "\n",
    "## Layer-wise Relevance Propagation\n",
    "\n",
    "As its name implies, LRP makes explicit use of the feed-forward graph structure of a DNN. The first layer are the inputs, the pixels of the image, and the last layer is the real-valued prediction output of the classifier $f$. LRP assumes that we have a relevance score $R^{(l+1)}_d$ for each pre-activation $z^{(l+1)}_d$ of the vector $z$ at layer $l + 1$. The idea is to find a relevance score $R^{(l)}_d$ for each dimension $z^{(l+1)}_d$ of the vector $z$ at the next layer $l$ which is close to the input layer such that the following holds:\n",
    "\n",
    "\\begin{equation}\n",
    "f(x) = \\cdots = \\sum_{d = 1}^{V(l+1)} R^{(l+1)}_d = \\sum_{d = 1}^{V(l)} R^{(l)}_d = \\cdots = \\sum_{i = 1}^{V(1)} R^{(1)}_d\n",
    "\\end{equation}\n",
    "\n",
    "Iterating this equation from the last layer which is the classifier output $f(x)$ down to the input layer $x$ consisting of image pixels naturally leads to *Constraint (2)*. With this, we can refine *Relevance Conservation Constraint* into *Layer-wise Relevance Conservation Constraint*.\n",
    "\n",
    "##### Constraint 3. Layer-wise Relevance Conservation Constraint: total relevance must be preserved throughout layers.\n",
    "\n",
    "Given a DNN where $j$ and $k$ are indices for neurons at two successive layers $l$ and $(l+1)$. We can define $R_{j \\leftarrow k}^{(l,l+1)}$ as the portion of relevance that flows from neuron $k$ to neuron $j$. The portion is determined by the amount of contribution of neuron $j$ to $R_{k}^{(l+1)}$, subject to the *Layer-wise Relevance Conservation Constraint*:\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_j R_{j \\leftarrow k}^{(l,l+1)} = R_{k}^{(l+1)}\n",
    "\\end{equation}\n",
    "\n",
    "Same can be said for $R_{j}^{(l)}$:\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_k R_{j \\leftarrow k}^{(l,l+1)} = R_{j}^{(l)}\n",
    "\\end{equation}\n",
    "\n",
    "## Propagation Rules for DNNs\n",
    "\n",
    "Let the neurons of the DNN be described by the equation\n",
    "\n",
    "\\begin{equation}\n",
    "a_k = \\sigma \\left( \\sum_j a_j w_{jk} + b_k \\right)\n",
    "\\end{equation}\n",
    "\n",
    "with $a_k$ the neuron activation, $a_j$ the activations from the previous layer, $w_{jk}$ the weight and $b_k$ the bias parameters of the neuron. The activation function $\\sigma$ is a positive and monotonically increasing activation function (e.g. tanh, ReLU).\n",
    "\n",
    "One propagation rule that fulfills the three constraints mentioned above is the $\\alpha \\beta$-rule. Let $z_{jk}^{+} = a_j w_{jk}^{+}$ and $z_{k}^{+} = \\sum_j a_j w_{jk}^{+} + b_{k}^{+} = \\sum_j z_{jk}^{+} + b_{k}^{+}$. $()^{+}$ and $()^{-}$ denote the positive and negative parts respectively. Same applies to $z_{jk}^{-}$ and $z_{k}^{-}$. Then the $\\alpha \\beta$-rule is given by\n",
    "\n",
    "\\begin{equation}\n",
    "R_{j \\leftarrow k}^{(l,l+1)} = R_{k}^{(l+1)} \\cdot \\left(\\alpha \\cdot \\frac{z_{jk}^{+}}{z_{k}^{+}} - \\beta \\cdot \\frac{z_{jk}^{-}}{z_{k}^{-}} \\right)\n",
    "\\end{equation}\n",
    "\n",
    "where the parameters $\\alpha$ and $\\beta$ satisfy the constraint $\\alpha - \\beta = 1$ and $\\beta \\geq 0$. Then, the *Layer-wise Relevance Conservation Constraint* becomes:\n",
    "\n",
    "\\begin{align}\n",
    "\\sum_j R_{j \\leftarrow k}^{(l,l+1)} & = \\sum_j R_{k}^{(l+1)} \\cdot \\left( \\alpha \\cdot \\frac{z_{jk}^{+}}{z_{k}^{+}} - \\beta \\cdot \\frac{z_{jk}^{-}}{z_{k}^{-}} \\right) \\\\\n",
    "& = R_{k}^{(l+1)} \\cdot \\left( \\alpha \\cdot \\frac{\\sum_j z_{jk}^{+}}{z_{k}^{+}} - \\beta \\cdot \\frac{\\sum_j z_{jk}^{-}}{z_{k}^{-}} \\right) \\\\\n",
    "& = R_{k}^{(l+1)} \\cdot \\left( \\alpha \\cdot \\frac{z_{k}^{+} - b_{k}^{+}}{z_{k}^{+}} - \\beta \\cdot \\frac{z_{k}^{-} - b_{k}^{-}}{z_{k}^{-}} \\right) \\\\\n",
    "& = R_{k}^{(l+1)} \\cdot \\left( (\\alpha - \\beta) - \\alpha \\cdot \\frac{b_{k}^{+}}{z_{k}^{+}} + \\beta \\cdot \\frac{b_{k}^{-}}{z_{k}^{-}} \\right) \\\\\n",
    "& = R_{k}^{(l+1)} \\cdot \\left( 1 - \\alpha \\cdot \\frac{b_{k}^{+}}{z_{k}^{+}} + \\beta \\cdot \\frac{b_{k}^{-}}{z_{k}^{-}} \\right) \\\\\n",
    "\\end{align}\n",
    "\n",
    "For different combinations of $\\alpha$ and $\\beta$, we name the corresponding propagation rule by subscripting $\\alpha$ and $\\beta$ with corresponding values. For example, choosing the parameters $\\alpha = 2$ and $\\beta = 1$ will give us the LRP-$\\alpha_2 \\beta_1$ rule.\n",
    "\n",
    "Choosing LRP-$\\alpha_1 \\beta_0$ simplifies the propagation rule to and the *Layer-wise Relevance Conservation Constraint* to:\n",
    "\n",
    "\\begin{equation}\n",
    "R_{j \\leftarrow k}^{(l,l+1)} = \\frac{z_{jk}^{+}}{z_{k}^{+}} \\cdot R_{k}^{(l+1)}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_j R_{j \\leftarrow k}^{(l,l+1)} = R_{k}^{(l+1)} \\cdot \\left( 1 - \\frac{b_{k}^{+}}{z_{k}^{+}} \\right)\n",
    "\\end{equation}\n",
    "\n",
    "In addition, if all the biases are constrained or set to be $0$, the propagation rule and the *Layer-wise Relevance Conservation Constraint* can be further simplified into:\n",
    "\n",
    "\\begin{equation}\n",
    "R_{j \\leftarrow k}^{(l,l+1)} = R_{k}^{(l+1)} \\cdot \\frac{a_j w_{jk}^{+}}{\\sum_j a_j w_{jk}^{+}}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_j R_{j \\leftarrow k}^{(l,l+1)} = R_{k}^{(l+1)}\n",
    "\\end{equation}\n",
    "\n",
    "In this tutorial, we are going to apply the above rule to a ReLU network without bias (which is basically equivalent to a network with $0$ bias)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tensorflow Implementation Details\n",
    "\n",
    "Consider the LRP-$\\alpha_1 \\beta_0$ propagation rule of *Eq. (12)*:\n",
    "\n",
    "\\begin{equation}\n",
    "R_{j}^{(l)} = a_j \\sum_k \\frac{w_{jk}^{+}}{\\sum_j a_j w_{jk}^{+} + b_{k}^{+}} R_{k}^{(l+1)}\n",
    "\\end{equation}\n",
    "\n",
    "This rule can be written as four elementary computations, all of which can also be expressed in vector form:\n",
    "\n",
    "##### Element-wise\n",
    "\n",
    "\\begin{align*}\n",
    "z_k & \\leftarrow \\sum_j a_j w_{jk}^{+} \\\\\n",
    "s_k & \\leftarrow R_k / z_k \\\\\n",
    "c_j & \\leftarrow \\sum_k w_{jk}^{+} s_k \\\\\n",
    "R_j & \\leftarrow a_j c_j\n",
    "\\end{align*}\n",
    "\n",
    "##### Vector Form\n",
    "\n",
    "\\begin{align*}\n",
    "\\mathbf{z} & \\leftarrow W_{+}^{\\top} \\cdot \\mathbf{a} \\\\\n",
    "\\mathbf{s} & \\leftarrow \\mathbf{R} \\oslash \\mathbf{z} \\\\\n",
    "\\mathbf{c} & \\leftarrow W_{+} \\cdot \\mathbf{s} \\\\\n",
    "\\mathbf{R} & \\leftarrow \\mathbf{a} \\odot \\mathbf{c}\n",
    "\\end{align*}\n",
    "\n",
    "By applying the same procedure to negative parts, we obtain the LRP implementation for all cases of $\\alpha$ and $\\beta$. In fully-connected dense layers, LRP can be implemented by the following sequence of Tensorflow operations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def backprop_dense(activation, kernel, bias, relevance):\n",
    "    W_p = tf.maximum(0., kernel)\n",
    "    b_p = tf.maximum(0., bias)\n",
    "    z_p = tf.matmul(activation, W_p) + b_p\n",
    "    s_p = relevance / z_p\n",
    "    c_p = tf.matmul(s_p, tf.transpose(W_p))\n",
    "\n",
    "    W_n = tf.minimum(0., kernel)\n",
    "    b_n = tf.minimum(0., bias)\n",
    "    z_n = tf.matmul(activation, W_n) + b_n\n",
    "    s_n = relevance / z_n\n",
    "    c_n = tf.matmul(s_n, tf.transpose(W_n))\n",
    "\n",
    "    return activation * (self.alpha * c_p + (1 - self.alpha) * c_n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In convolution layers, the matrix-vector multiplications can be more efficiently implemented by `backprop` methods used for gradient propagation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def backprop_conv(self, activation, kernel, bias, relevance, strides, padding='SAME'):\n",
    "    W_p = tf.maximum(0., kernel)\n",
    "    b_p = tf.maximum(0., bias)\n",
    "    z_p = nn_ops.conv2d(activation, W_p, strides, padding) + b_p\n",
    "    s_p = relevance / z_p\n",
    "    c_p = nn_ops.conv2d_backprop_input(tf.shape(activation), W_p, s_p, strides, padding)\n",
    "\n",
    "    W_n = tf.minimum(0., kernel)\n",
    "    b_n = tf.minimum(0., bias)\n",
    "    z_n = nn_ops.conv2d(activation, W_n, strides, padding) + b_n\n",
    "    s_n = relevance / z_n\n",
    "    c_n = nn_ops.conv2d_backprop_input(tf.shape(activation), W_n, s_n, strides, padding)\n",
    "\n",
    "    return activation * (self.alpha * c_p + (1 - self.alpha) * c_n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "In max-pooling layers, the original paper by [Bach et al.](http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0130140) uses a winner-take-all redistribution policy where all relevance goes to the most activated neuron in the pool. On the other hand, the paper by [Montavon et al.](https://www.sciencedirect.com/science/article/pii/S0031320316303582) uses the proportional redistribution rule where the redistribution is proportional to neuron activations in the pool:\n",
    "\n",
    "\\begin{equation}\n",
    "R_{j}^{(l)} = \\frac{x_j}{\\sum_j x_j} R_{k}^{(l+1)}\n",
    "\\end{equation}\n",
    "\n",
    "Like the convolution layers, redistribution in pooling layers can also be efficiently implemend by `backprop` methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Bach et al.'s redistribution rule\n",
    "def backprop_pool(self, activation, relevance, ksize, strides, pooling_type, padding='SAME'):\n",
    "    z = nn_ops.max_pool(activation, ksize, strides, padding) + 1e-10\n",
    "    s = relevance / z\n",
    "    c = gen_nn_ops._max_pool_grad(activation, z, s, ksize, strides, padding)\n",
    "    return activation * c\n",
    "\n",
    "\n",
    "# Montavon et al.'s redistribution rule\n",
    "def backprop_pool(self, activation, relevance, ksize, strides, pooling_type, padding='SAME'):\n",
    "    z = nn_ops.avg_pool(activation, ksize, strides, padding) + 1e-10\n",
    "    s = relevance / z\n",
    "    c = gen_nn_ops._avg_pool_grad(tf.shape(activation), s, ksize, strides, padding)\n",
    "    return activation * c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have all the tools for LRP-$\\alpha \\beta$, we will see its application to a DNN trained on MNIST in the next part of the tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}