Lectures/Lecture29/3dasm_Lecture29.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<img src=../figures/Brown_logo.svg width=50%>\n",
    "\n",
    "## Data-Driven Design & Analyses of Structures & Materials (3dasm)\n",
    "\n",
    "## Lecture 29\n",
    "\n",
    "### Miguel A. Bessa | <a href = \"mailto: miguel_bessa@brown.edu\">miguel_bessa@brown.edu</a>  | Associate Professor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "**What:** A lecture of the \"3dasm\" course\n",
    "\n",
    "**Where:** This notebook comes from this [repository](https://github.com/bessagroup/3dasm_course)\n",
    "\n",
    "**Reference for entire course:** Murphy, Kevin P. *Probabilistic machine learning: an introduction*. MIT press, 2022. Available online [here](https://probml.github.io/pml-book/book1.html)\n",
    "\n",
    "**How:** We try to follow Murphy's book closely, but the sequence of Chapters and Sections is different. The intention is to use notebooks as an introduction to the topic and Murphy's book as a resource.\n",
    "* If working offline: Go through this notebook and read the book.\n",
    "* If attending class in person: listen to me (!) but also go through the notebook in your laptop at the same time. Read the book.\n",
    "* If attending lectures remotely: listen to me (!) via Zoom and (ideally) use two screens where you have the notebook open in 1 screen and you see the lectures on the other. Read the book."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "**Optional reference (the \"bible\" by the \"bishop\"... pun intended 😆) :** Bishop, Christopher M. *Pattern recognition and machine learning*. Springer Verlag, 2006.\n",
    "\n",
    "**References/resources to create this notebook:**\n",
    "* This simple tutorial is still based on a script I created for this article: https://imechanica.org/node/23957\n",
    "* It follows from some examples provided by the scikit-learn user guide, which seem to have originated from Mathieu Blondel, Jake Vanderplas, Vincent Dubourg, and Jan Hendrik Metzen.\n",
    "\n",
    "Apologies in advance if I missed some reference used in this notebook. Please contact me if that is the case, and I will gladly include it here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "## **OPTION 1**. Run this notebook **locally in your computer**:\n",
    "1. Confirm that you have the '3dasm' mamba (or conda) environment (see Lecture 1).\n",
    "2. Go to the 3dasm_course folder in your computer and pull the last updates of the [repository](https://github.com/bessagroup/3dasm_course):\n",
    "```\n",
    "git pull\n",
    "```\n",
    "    - Note: if you can't pull the repo due to conflicts (and you can't handle these conflicts), use this command (with **caution**!) and your repo becomes the same as the one online:\n",
    "```\n",
    "git reset --hard origin/main\n",
    "```\n",
    "3. Open command window and load jupyter notebook (it will open in your internet browser):\n",
    "```\n",
    "jupyter notebook\n",
    "```\n",
    "5. Open notebook of this Lecture and choose the '3dasm' kernel."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "## **OPTION 2**. Use **Google's Colab** (no installation required, but times out if idle):\n",
    "\n",
    "1. go to https://colab.research.google.com\n",
    "2. login\n",
    "3. File > Open notebook\n",
    "4. click on Github (no need to login or authorize anything)\n",
    "5. paste the git link: https://github.com/bessagroup/3dasm_course\n",
    "6. click search and then click on the notebook for this Lecture."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "scrolled": true,
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "outputs": [],
   "source": [
    "# Basic plotting tools needed in Python.\n",
    "\n",
    "import matplotlib.pyplot as plt # import plotting tools to create figures\n",
    "import numpy as np # import numpy to handle a lot of things!\n",
    "from IPython.display import display, Math # to print with Latex math\n",
    "\n",
    "%config InlineBackend.figure_format = \"retina\" # render higher resolution images in the notebook\n",
    "plt.rcParams[\"figure.figsize\"] = (8,4) # rescale figure size appropriately for slides"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Outline for today\n",
    "\n",
    "* Deriving the PPD of Artificial Neural Networks\n",
    "* What is Backpropagation? Why do we need it?\n",
    "\n",
    "**Reading material**: This notebook + (ANNs in Chapter 13)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Last lecture we gained (some!) intuition about what Artificial Neural Networks are doing by playing with https://playground.tensorflow.org\n",
    "\n",
    "Today, we are going to derive the PPD of ANNs. It's actually quite easy!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Defining ANNs from a Bayesian perspective\n",
    "\n",
    "Although ANNs are usually viewed in the literature from a deterministic perspective, this does not need to be the case.\n",
    "\n",
    "* For classification, they are a very trivial generalization from Logistic regression\n",
    "\n",
    "* For regression, they are a very trivial generalization from Linear regression."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "ANN for **classification** problems is defined similarly to what we did for Logistic regression, but using a nonlinear function to model the mean:\n",
    "\n",
    "1. Categorical observation distribution: $p(\\mathbf{y}|\\mathbf{x}, \\mathbf{z}) = \\text{Cat}\\left(\\mathbf{y}| \\boldsymbol{\\rho} = \\mathbf{f}(\\mathbf{x};\\, \\mathbf{z}) \\right)$\n",
    "\n",
    "    where $\\mathbf{z} \\equiv \\{\\mathbf{W}, \\mathbf{b}\\}$ are all the unknown model parameters (weights and biases of every layer) and where $\\mathbf{f}(\\mathbf{x};\\, \\mathbf{z})$ is nonlinear in the weights due to the choice of activation functions (the activation function in the **output layer** is the **softmax activation**, just like in logistic regression).\n",
    "\n",
    "2. Uniform prior distribution for each hidden rv $\\mathbf{z}$: $p(\\mathbf{z}) \\propto 1$\n",
    "    \n",
    "    * Note: we can also use other priors! For example a Gaussian prior ($\\ell_2$ regularization) or a Laplace prior ($\\ell_1$ regularization).\n",
    "\n",
    "3. MLE point estimate for posterior: $\\hat{\\mathbf{z}}_{\\text{mle}} = \\underset{z}{\\mathrm{argmin}}\\left[-\\sum_{i=1}^{N}\\log{ p(\\mathbf{y}=\\mathbf{y}_i|\\mathbf{x}=\\mathbf{x}_i, \\mathbf{z})}\\right]$\n",
    "    * Note: if we use a different prior, then it becomes the MAP estimate.\n",
    "\n",
    "Final prediction is given by the <font color='orange'>PPD</font>: $\\require{color}\n",
    "{\\color{orange}p(\\mathbf{y}^*|\\mathbf{x}^*, \\mathcal{D})} = \\int p(\\mathbf{y}^*|\\mathbf{x}^*,\\mathbf{z}) \\delta(\\mathbf{z}-\\hat{\\mathbf{z}}) d\\mathbf{z} = p(\\mathbf{y}^*|\\mathbf{x}^*, \\mathbf{z}=\\hat{\\mathbf{z}})$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "ANN for **regression** problems is defined similarly to what we did for Linear regression, but using a nonlinear function to model the mean:\n",
    "\n",
    "1. Gaussian observation distribution (**variance is a constant hyperparameter**): $p(\\mathbf{y}|\\mathbf{x}, \\mathbf{z}) = \\mathcal{N}(\\mathbf{y}| \\boldsymbol{\\mu}_{y|z} = \\mathbf{f}(\\mathbf{x};\\, \\mathbf{z}), \\sigma_{y|z}^2 = \\sigma^2)$\n",
    "\n",
    "    where $\\mathbf{z} \\equiv \\{\\mathbf{W}, \\mathbf{b}\\}$ are all the unknown model parameters (weights and biases of every layer) and where $\\mathbf{f}(\\mathbf{x};\\, \\mathbf{z})$ is nonlinear in the weights due to the choice of activation functions (the activation function in the **output layer** is the **linear activation**, just like in linear regression).\n",
    "\n",
    "2. Uniform prior distribution for each hidden rv $\\mathbf{z}$: $p(\\mathbf{z}) \\propto 1$\n",
    "    \n",
    "    * Note: we can also use other priors! For example a Gaussian prior ($\\ell_2$ regularization) or a Laplace prior ($\\ell_1$ regularization).\n",
    "\n",
    "3. MLE point estimate for posterior: $\\hat{\\mathbf{z}}_{\\text{mle}} = \\underset{z}{\\mathrm{argmin}}\\left[-\\sum_{i=1}^{N}\\log{ p(\\mathbf{y}=\\mathbf{y}_i|\\mathbf{x}=\\mathbf{x}_i, \\mathbf{z})}\\right]$\n",
    "    * Note: if we use a different prior, then it becomes the MAP estimate.\n",
    "\n",
    "Final prediction is given by the <font color='orange'>PPD</font>: $\\require{color}\n",
    "{\\color{orange}p(\\mathbf{y}^*|\\mathbf{x}^*, \\mathcal{D})} = \\int p(\\mathbf{y}^*|\\mathbf{x}^*,\\mathbf{z}) \\delta(\\mathbf{z}-\\hat{\\mathbf{z}}) d\\mathbf{z} = p(\\mathbf{y}^*|\\mathbf{x}^*, \\mathbf{z}=\\hat{\\mathbf{z}})$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "In both cases (classification and regression), we need to calculate the point estimate of the unknowns $\\mathbf{z}$.\n",
    "\n",
    "Without loss of generality, let's start by focusing on the **regression** case and assuming uniform prior, i.e. no \"regularization\" (as we did in Lecture 11 for Linear Least Squares):\n",
    "\n",
    "$$\\begin{align}\n",
    "\\hat{\\mathbf{z}}_{\\text{mle}} &= \\underset{z}{\\mathrm{argmin}}\\left[\\text{NLL}(\\mathbf{z})\\right]\n",
    "\\\\\n",
    "&= \\underset{z}{\\mathrm{argmin}}\\left[-\\sum_{n=1}^{N}\\log{ p(\\mathbf{y}=\\mathbf{y}_n|\\mathbf{x}=\\mathbf{x}_n, \\mathbf{z})}\\right]\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "\\hat{\\mathbf{z}}_{\\text{mle}} &= \\underset{z}{\\mathrm{argmin}}\\left[-\\sum_{n=1}^{N}\\log{ p(\\mathbf{y}=\\mathbf{y}_n|\\mathbf{x}=\\mathbf{x}_n, \\mathbf{z})}\\right] \\\\\n",
    "&= \\underset{z}{\\mathrm{argmin}}\\left[-\\sum_{n=1}^{N}\\log{\\left( \\frac{1}{\\sqrt{2\\pi \\sigma^2}} \\exp\\left\\{ -\\frac{1}{2\\sigma^2}\\left[\\mathbf{y}_n-\\mathbf{f}(\\mathbf{x}_n;\\, \\mathbf{z})\\right]^T\\left[\\mathbf{y}_n-\\mathbf{f}(\\mathbf{x}_n;\\, \\mathbf{z})\\right]\\right\\}\\right)}\\right]\\\\\n",
    "&= \\underset{z}{\\mathrm{argmin}}\\left[\\frac{N}{2}\\log{\\left(2\\pi \\sigma^2\\right)}+\\frac{1}{2 \\sigma^2}\\sum_{n=1}^{N}\\left[\\mathbf{y}_n-\\mathbf{f}(\\mathbf{x}_n;\\, \\mathbf{z})\\right]^T\\left[\\mathbf{y}_n-\\mathbf{f}(\\mathbf{x}_n;\\, \\mathbf{z})\\right] \\right]\\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where we recall that the unknowns are $\\mathbf{z}$ (the weights and biases of the network), and where we can ignore the variance term (because we assume it to be a constant hyperparameter, i.e. we will not solve for it)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Therefore, the point estimate of the ANN for regression becomes:\n",
    "\n",
    "$$\n",
    "\\hat{\\mathbf{z}}_{\\text{mle}} = \\underset{z}{\\mathrm{argmin}}\\left[\\frac{1}{2}\\sum_{n=1}^{N}\\left[\\mathbf{y}_n-\\mathbf{f}(\\mathbf{x}_n;\\, \\mathbf{z})\\right]^T\\left[\\mathbf{y}_n-\\mathbf{f}(\\mathbf{x}_n;\\, \\mathbf{z})\\right] \\right]\n",
    "$$\n",
    "\n",
    "where the argument is called the $\\ell_2$ loss function, i.e. the square of the $\\ell_2$ norm (equivalently, we could also use mean squared error loss function if we divide the $\\ell_2$ loss by the number of samples $N$).\n",
    "\n",
    "As we have seen previously, we can rewrite the $\\ell_2$ loss as follows:\n",
    "\n",
    "$$\n",
    "\\mathcal{L}\\left[(\\mathbf{x},\\mathbf{y}), \\mathbf{z} \\right] = \\frac{1}{2} ||\\mathbf{Y}-\\mathbf{f}(\\mathbf{X},\\mathbf{z}) ||_2^2\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* Note 1: Similarly, we can  find the loss function if we consider a nonuniform prior."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* Note 2: The loss function that you obtain for the **classification** case is different! It becomes the **cross-entropy** loss function, as we derived in Lecture 22, but where the prediction now depends on a nonlinear function: $\\mathbf{f}(\\mathbf{x}_n;\\, \\mathbf{z})$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Explanation of Note 1: ANN for regression with Gaussian prior\n",
    "\n",
    "As we showed in Lecture 12 (Ridge regression), if we assume a nonuniform prior and the MAP point estimate (Lecture 8):\n",
    "\n",
    "$\\require{color}\n",
    "\\hat{\\mathbf{z}}_{\\text{map}} = \\underset{z}{\\mathrm{argmin}}\\left[-\\log{p(\\mathbf{y}=\\mathbf{y}_n|\\mathbf{x}=\\mathbf{x}_n, \\mathbf{z})} - \\log{{\\color{red}p(z)}}\\right]$\n",
    "\n",
    "For example, if the prior $p(z)$ is a Gaussian distribution on the weights $\\mathbf{w}$ (we can also use a Gaussian distribution for the biases):\n",
    "\n",
    "$p(\\mathbf{w}) = \\mathcal{N}(\\mathbf{w}| \\mathbf{0}, \\overset{\\scriptscriptstyle <}{\\sigma}_w^2 \\mathbf{I})$\n",
    "\n",
    "where $\\overset{\\scriptscriptstyle <}{\\sigma}_w$ is specified by us (hyperparameter), and where we are writing the weights of the ANN as a vector (we unraveled the weight matrix into a vector).\n",
    "\n",
    "* If we wrote the weights as a matrix, then we can use double contraction of the weight matrix with itself (leading to the same result as the inner product of the weight vector, of course)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Performing the calculation, the point estimate when assuming a Gaussian prior becomes:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\hat{\\mathbf{z}}_{\\text{map}} &= \\underset{z}{\\mathrm{argmin}}\\left[\\frac{1}{2\\sigma^2}\\sum_{n=1}^{N}\\left[\\mathbf{y}_n-\\mathbf{f}(\\mathbf{x}_n;\\, \\mathbf{z})\\right]^T\\left[\\mathbf{y}_n-\\mathbf{f}(\\mathbf{x}_n;\\, \\mathbf{z})\\right] + \\frac{1}{2\\overset{\\scriptscriptstyle <}{\\sigma}_w^2}\\mathbf{w}^T\\mathbf{w} \\right] \\\\\n",
    "&= \\underset{z}{\\mathrm{argmin}}\\left[\\frac{1}{2}\\sum_{n=1}^{N}\\left[\\mathbf{y}_n-\\mathbf{f}(\\mathbf{x}_n;\\, \\mathbf{z})\\right]^T\\left[\\mathbf{y}_n-\\mathbf{f}(\\mathbf{x}_n;\\, \\mathbf{z})\\right] + \\alpha\\mathbf{w}^T\\mathbf{w} \\right] \\\\\n",
    "&= \\underset{z}{\\mathrm{argmin}}\\left[\\frac{1}{2} ||\\mathbf{Y}-\\mathbf{f}(\\mathbf{X},\\mathbf{z}) ||_2^2 + \\alpha ||\\mathbf{w}||_2^2\\right]\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where $\\alpha = \\frac{\\sigma^2}{\\overset{\\scriptscriptstyle <}{2\\sigma}_w^2}$ is a hyperparameter (proportional to the strength of the prior and the value assumed for the variance of the observation distribution), and\n",
    "\n",
    "$||\\mathbf{w}||_2 = \\sqrt{\\sum |w_m|^2} = \\sqrt{\\mathbf{w}^T\\mathbf{w}} \\quad \\text{where the sum is over all weights (excluding biases) of the ANN.}$\n",
    "\n",
    "This is called the $l_2$ norm of the vector $\\mathbf{w}$. As we saw in Lecture 12 for Ridge regression, this causes a penalty on the weights that become too large in magnitude. In ML literature this is usually called $l_2$ **regularization** or **weight decay**, and is very widely used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Conclusion**: as in the other ML models, finding the point estimate means that we need to minimize the loss function (finding minimum location of the NLL or of the negative log joint likelihood if the prior is not uniform)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* However, there is a **major** problem now: this function is no longer convex because the mean $\\mathbf{f}(\\mathbf{x}_n;\\, \\mathbf{z})$ is a composition of nonlinear activation functions..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "So, what type of optimization algorithm would you choose??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* As we discussed before, first-order optimizers (gradient-based optimizers) are quite efficient.\n",
    "    * But... in order for us to use them **we need to calculate the gradient of the loss function** and then provide both the loss function and its gradient to an optimizer that will find an optimum."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Minimizing the loss function with gradient-based optimizers\n",
    "\n",
    "There are many optimization algorithms that can be used to minimize the loss function.\n",
    "\n",
    "* The most common choices involve the use of first-order optimizers, i.e. optimizers that use the gradient of the loss function.\n",
    "* Among the optimizers we learned on Lectures 24 to 26, ADAM and SGD are particularly common.\n",
    "    * Note: L-BFGS is not so common in ANN literature because the loss landscapes are very high dimensional and, especially, they tend to be noisy (line search does not like that, as we discussed in class)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## How to calculate the gradient of the loss function?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Last Lecture we saw that the mean function of an ANN with $L$ hidden layers is a composition of $L+1$ functions:\n",
    "\n",
    "$\\mathbf{f} = \\mathbf{f}_{L+1} \\circ \\mathbf{f}_L \\circ \\ldots \\circ \\mathbf{f}_{2} \\circ \\mathbf{f}_1 $\n",
    "\n",
    "And if we denote $f_{L+2} = \\mathcal{L}$ as the loss function, then we can write the loss function as follows:\n",
    "\n",
    "$\\mathcal{L} = f_{L+2} \\circ \\mathbf{f}_{L+1} \\circ \\mathbf{f}_L \\circ \\ldots \\circ \\mathbf{f}_{l} \\circ \\ldots \\circ \\mathbf{f}_{2} \\circ \\mathbf{f}_1 $\n",
    "\n",
    "Where $l=1,\\ldots, L$ are the hidden layers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The next slide makes it less abstract, by considering the example of an ANN with **1 hidden layer** and the $\\ell_2$ loss (i.e. regression case)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Neurons of ANN with **1 hidden layer** and using $\\ell_2$ loss (i.e. regression case).\n",
    "\n",
    "* Layer 0 (input layer): $\\mathbf{n}_0 \\equiv \\mathbf{x}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* Layer $L=1$ (first hidden layer):\n",
    "\n",
    "$$\\begin{align}\n",
    "\\mathbf{n}_L &= \\mathbf{n}_1 \\\\\n",
    "             &= \\mathbf{f}_L(\\mathbf{n}_{L-1}, \\mathbf{z}_{L-1}) \\\\\n",
    "             &= \\mathbf{f}_1(\\mathbf{n}_0, \\mathbf{z}_0) \\\\\n",
    "             &= \\mathbf{f}_1(\\mathbf{b}_0 + \\mathbf{W}_0\\mathbf{n}_0)\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* Layer $L+1=2$ (output layer):\n",
    "\n",
    "$$\\begin{align}\n",
    "\\mathbf{n}_{L+1} &= \\mathbf{n}_2 \\equiv \\mathbf{y} \\\\\n",
    "           &= \\mathbf{f}_{L+1}(\\mathbf{n}_{L}, \\mathbf{z}_{L}) = \\mathbf{f}_2(\\mathbf{n}_1, \\mathbf{z}_1) \\\\ \n",
    "           &= \\mathbf{f}_2(\\mathbf{b}_1 + \\mathbf{W}_1\\mathbf{n}_1) & \\\\\n",
    "           &= \\mathbf{b}_1 + \\mathbf{W}_1\\mathbf{n}_1\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Note: the last line results from assuming this is a **regression problem**, and in this case the **output layer** (layer $L+1$) uses a linear activation function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "* Layer $L+2$ can be considered the loss function:\n",
    "\n",
    "$$\\begin{align}\n",
    "\\mathcal{L} &= f_{L+2}(\\mathbf{n}_{L+1}, \\mathbf{y})\\\\\n",
    "            &= f_3(\\mathbf{n}_2, \\mathbf{y}) \\\\\n",
    "            &= \\frac{1}{2}|| \\mathbf{y} - \\mathbf{n}_2  ||_2^2\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where, again, we are considering here the $\\ell_2$ loss because we are focusing on a regression case. If it was a classification problem we would consider the cross-entropy loss (see Lecture 22)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Now we can calculate the derivatives of the loss function wrt to the unknowns $\\mathbf{z}$ very easily!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "In order to calculate the derivatives of the loss function:\n",
    "\n",
    "$\\mathcal{L} = f_{L+2} \\circ \\mathbf{f}_{L+1} \\circ \\mathbf{f}_L \\circ \\ldots \\circ \\mathbf{f}_{2} \\circ \\mathbf{f}_1 $\n",
    "\n",
    "we just need to use the chain rule!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's start by taking the derivative of the loss wrt the unkowns of the last layer, i.e. $\\mathbf{z}_L$:\n",
    "\n",
    "$$\\begin{align}\n",
    "\\nabla_{\\mathbf{z}_{L}} \\mathcal{L} = \\frac{\\partial \\mathcal{L}}{\\partial \\mathbf{z}_L} &= \\frac{\\partial f_{L+2}}{\\partial \\mathbf{z}_{L}} \\\\\n",
    "&= \\frac{\\partial f_{L+2}}{\\partial \\mathbf{f}_{L+1}} \\frac{\\partial \\mathbf{f}_{L+1}}{\\partial \\mathbf{z}_{L}}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Which is usually presented as: $\\frac{\\partial \\mathcal{L}}{\\partial \\mathbf{z}_L} = \\frac{\\partial \\mathcal{L}}{\\partial \\mathbf{n}_{L+1}} \\frac{\\partial \\mathbf{n}_{L+1}}{\\partial \\mathbf{z}_{L}} = \\frac{\\partial \\mathcal{L}}{\\partial \\mathbf{y}} \\frac{\\partial \\mathbf{y}}{\\partial \\mathbf{z}_{L}}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "but we should not stop here! Remember, we want all the derivatives of the loss wrt to all the unknowns $\\mathbf{z}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "So, taking the derivative of the loss wrt the unkowns of the next hidden layer, i.e. wrt $\\mathbf{z}_{L-1}$:\n",
    "\n",
    "$$\\begin{align}\n",
    "\\nabla_{\\mathbf{z}_{L-1}} \\mathcal{L} = \\frac{\\partial \\mathcal{L}}{\\partial \\mathbf{z}_{L-1}} &= \\frac{\\partial f_{L+2}}{\\partial \\mathbf{z}_{L-1}} \\\\\n",
    "&= \\frac{\\partial f_{L+2}}{\\partial \\mathbf{f}_{L+1}} \\frac{\\partial \\mathbf{f}_{L+1}}{\\partial \\mathbf{f}_{L}} \\frac{\\partial \\mathbf{f}_{L}}{\\partial \\mathbf{z}_{L-1}}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Which is usually presented as:\n",
    "\n",
    "$$\n",
    "\\frac{\\partial \\mathcal{L}}{\\partial \\mathbf{z}_{L-1}} = \\frac{\\partial \\mathcal{L}}{\\partial \\mathbf{n}_{L+1}} \\frac{\\partial \\mathbf{n}_{L+1}}{\\partial \\mathbf{n}_{L}} \\frac{\\partial \\mathbf{n}_{L}}{\\partial \\mathbf{z}_{L-1}} = \\frac{\\partial \\mathcal{L}}{\\partial \\mathbf{y}} \\frac{\\partial \\mathbf{y}}{\\partial \\mathbf{n}_{L}} \\frac{\\partial \\mathbf{n}_{L}}{\\partial \\mathbf{z}_{L-1}}\n",
    "$$\n",
    "\n",
    "but we should not stop here!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "So, taking the derivative of the loss wrt $\\mathbf{z}_{l}$:\n",
    "\n",
    "$$\\begin{align}\n",
    "\\nabla_{\\mathbf{z}_{l}} \\mathcal{L} = \\frac{\\partial \\mathcal{L}}{\\partial \\mathbf{z}_{l}} &= \\frac{\\partial f_{L+2}}{\\partial \\mathbf{z}_{l}} \\\\\n",
    "&= \\frac{\\partial f_{L+2}}{\\partial \\mathbf{f}_{L+1}} \\frac{\\partial \\mathbf{f}_{L+1}}{\\partial \\mathbf{f}_{L}} \\cdots \\frac{\\partial \\mathbf{f}_{l+2}}{\\partial \\mathbf{f}_{l+1}}\\frac{\\partial \\mathbf{f}_{l+1}}{\\partial \\mathbf{z}_{l}}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Which is usually presented as:\n",
    "\n",
    "$$\n",
    "\\frac{\\partial \\mathcal{L}}{\\partial \\mathbf{z}_{L-1}} = \\frac{\\partial \\mathcal{L}}{\\partial \\mathbf{n}_{L+1}} \\frac{\\partial \\mathbf{n}_{L+1}}{\\partial \\mathbf{n}_{L}} \\cdots \\frac{\\partial \\mathbf{n}_{l+2}}{\\partial \\mathbf{n}_{l+1}}\\frac{\\partial \\mathbf{n}_{l+1}}{\\partial \\mathbf{z}_{l}} = \\frac{\\partial \\mathcal{L}}{\\partial \\mathbf{y}} \\frac{\\partial \\mathbf{y}}{\\partial \\mathbf{n}_{L}} \\cdots \\frac{\\partial \\mathbf{n}_{l+2}}{\\partial \\mathbf{n}_{l+1}}\\frac{\\partial \\mathbf{n}_{l+1}}{\\partial \\mathbf{z}_{l}}\n",
    "$$\n",
    "\n",
    "but we should not stop here! We need to keep going until we reach the first layer of the network!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "So, taking the derivative of the loss wrt $\\mathbf{z}_{0}$:\n",
    "\n",
    "$$\\begin{align}\n",
    "\\nabla_{\\mathbf{z}_{0}} \\mathcal{L} = \\frac{\\partial \\mathcal{L}}{\\partial \\mathbf{z}_{0}} &= \\frac{\\partial f_{L+2}}{\\partial \\mathbf{z}_{l}} \\\\\n",
    "&= \\frac{\\partial f_{L+2}}{\\partial \\mathbf{f}_{L+1}} \\frac{\\partial \\mathbf{f}_{L+1}}{\\partial \\mathbf{f}_{L}} \\cdots \\frac{\\partial \\mathbf{f}_{l+2}}{\\partial \\mathbf{f}_{l+1}}\\frac{\\partial \\mathbf{f}_{l+1}}{\\partial \\mathbf{z}_{l}} \\cdots \\frac{\\partial \\mathbf{f}_{2}}{\\partial \\mathbf{f}_{1}}\\frac{\\partial \\mathbf{f}_{1}}{\\partial \\mathbf{z}_{0}} \n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Which is usually presented as:\n",
    "\n",
    "$$\n",
    "\\frac{\\partial \\mathcal{L}}{\\partial \\mathbf{z}_{L-1}} = \\frac{\\partial \\mathcal{L}}{\\partial \\mathbf{n}_{L+1}} \\frac{\\partial \\mathbf{n}_{L+1}}{\\partial \\mathbf{n}_{L}} \\cdots \\frac{\\partial \\mathbf{n}_{l+2}}{\\partial \\mathbf{n}_{l+1}}\\frac{\\partial \\mathbf{n}_{l+1}}{\\partial \\mathbf{z}_{l}} \\cdots \\frac{\\partial \\mathbf{n}_{2}}{\\partial \\mathbf{n}_{1}}\\frac{\\partial \\mathbf{n}_{1}}{\\partial \\mathbf{z}_{0}} = \\frac{\\partial \\mathcal{L}}{\\partial \\mathbf{y}} \\frac{\\partial \\mathbf{y}}{\\partial \\mathbf{n}_{L}} \\cdots \\frac{\\partial \\mathbf{n}_{l+2}}{\\partial \\mathbf{n}_{l+1}}\\frac{\\partial \\mathbf{n}_{l+1}}{\\partial \\mathbf{z}_{l}} \\cdots \\frac{\\partial \\mathbf{n}_{2}}{\\partial \\mathbf{n}_{1}}\\frac{\\partial \\mathbf{n}_{1}}{\\partial \\mathbf{z}_{0}} \n",
    "$$\n",
    "\n",
    "* And with this we finally computed all the derivatives!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "This is what we call **Backpropagation**.\n",
    "* The invention of backprop is credited to David Rumelhart. Geoff Hinton, a giant in the field of Artificial Intelligence, later showed that training ANNs using backprop (optimization) led to interesting internal representations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The partial derivatives of the neuron output with respect to the previous neuron outputs are organized in a matrix called the **Jacobian**:\n",
    "\n",
    "$$\n",
    "\\mathbf{J}_{\\mathbf{f}_{l+1}}(\\mathbf{n}_l) = \\frac{\\partial \\mathbf{f}_{l+1}(\\mathbf{n}_{l})}{\\partial \\mathbf{n}_l} = \\begin{bmatrix}\\frac{\\partial f_{1, l+1}}{\\partial n_{1,l}} & \\cdots & \\frac{\\partial f_{D_{l+1}, l+1}}{\\partial n_{D_l,l}}\\\\\n",
    "\\vdots & \\ddots & \\vdots\\\\\n",
    "\\frac{\\partial f_{D_{l+1}, l+1}}{\\partial n_{1,l}} & \\cdots & \\frac{\\partial f_{D_{l+1}, l+1}}{\\partial n_{D_l,l}}\\end{bmatrix} =\n",
    "\\begin{bmatrix}\\nabla f_{1,l+1}(\\mathbf{n}_l)^T\\\\\n",
    "\\vdots\\\\\n",
    "\\nabla f_{D_{l+1}, l+1}(\\mathbf{n}_l)^T\n",
    "\\end{bmatrix}=\n",
    "\\begin{bmatrix}\n",
    "\\frac{\\partial \\mathbf{f}_{l+1}}{\\partial n_{0,l}}, & \\cdots & , \\frac{\\partial \\mathbf{f}_{l+1}}{\\partial n_{D_l,l}}\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "where $D_{l+1}$ and $D_l$ are the number of neurons in layer $l+1$ and $l$, respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "So, doing backpropagation (chain rule) using Jacobians becomes:\n",
    "\n",
    "$$\n",
    "\\mathbf{J}_{\\mathbf{f}_{L+1}}(\\mathbf{n}_0) = \\mathbf{J}_{\\mathbf{y}}(\\mathbf{x}) = \\mathbf{J}_{\\mathbf{f}_{L+1}}(\\mathbf{n}_{L}) \\mathbf{J}_{\\mathbf{f}_{L}}(\\mathbf{n}_{L-1}) \\cdots \\mathbf{J}_{\\mathbf{f}_{l+1}}(\\mathbf{n}_{l}) \\cdots \\mathbf{J}_{\\mathbf{f}_{1}}(\\mathbf{n}_{0})\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Backpropagation\n",
    "\n",
    "* Backpropagation is simply the computation of the gradients (derivatives) of the loss wrt to the unknowns $\\mathbf{z}$.\n",
    "\n",
    "    * Note: remember, \"loss function\" just means NLL (when using a uniform prior) or negative log joint likelihood (when using any other prior)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "* These gradients are then provided to **a gradient-based optimizer** (e.g. in Lecture 27 we used ADAM), that will then provide a point estimate (the weights and biases that minimize the error; or, alternatively, that maximize the joint likelihood)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### What is the consequence of the loss function not being convex?\n",
    "\n",
    "When finding the minimum of a non-convex function, we basically **cross our fingers** and hope that we land in a \"decent\" optimum (which almost never is the \"global\" optimum, contrary to what happened in linear models)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The conundrum we face in machine learning is to:\n",
    "\n",
    "* **\"Make everything as simple as possible, but not simpler\"** (paraphrasing Albert Einstein)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* Linear models are nice because you can find the global optimum and you usually reach it very quickly (convex optimization problems are easy to solve). **BUT** linear models usually are **not sufficiently expressive** to predict complicated maps..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* Nonlinear models with many parameters (and hyperparameters...) are **more expressive** (after all, a nonlinear model includes a linear model!), but the optimization problem to find a point estimate is more difficult to solve and we are not guaranteed to reach the global optimum."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Still, there are countless examples where very expressive nonlinear models (different types of ANNs) outperform simpler models for complex tasks:\n",
    "\n",
    "* For example, classification problems on large datasets (e.g. ImageNet dataset), learning how to play chess or GO, ChatGPT, etc.\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### See you next class\n",
    "\n",
    "Have fun!"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "hide_code_all_hidden": false,
  "kernelspec": {
   "display_name": "Python [conda env:3dasm]",
   "language": "python",
   "name": "conda-env-3dasm-py"
  },
  "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.10.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}