Deep.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 7: Deep Learning"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ % Latex macros\n",
    "\\newcommand{\\mat}[1]{\\begin{pmatrix} #1 \\end{pmatrix}}\n",
    "\\newcommand{\\p}[2]{\\frac{\\partial #1}{\\partial #2}}\n",
    "\\renewcommand{\\b}[1]{\\boldsymbol{#1}}\n",
    "\\newcommand{\\w}{\\boldsymbol{w}}\n",
    "\\newcommand{\\x}{\\boldsymbol{x}}\n",
    "\\newcommand{\\y}{\\boldsymbol{y}}\n",
    "\\newcommand{\\z}{\\boldsymbol{z}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "* 7.1 Multi-Layer Neural Networks (Bishop, Chater 5)\n",
    "    * Universality of two-layer networks\n",
    "    * Recurrent neural networks\n",
    "\n",
    "\n",
    "* 7.2 Back-Propagation Learning\n",
    "    * Function approximation\n",
    "    * Classification\n",
    "\n",
    "\n",
    "* 7.3 Deep Generative Models\n",
    "    \n",
    "* 7.4 Restricted Boltzmann Machines (RBM)\n",
    "\n",
    "* 7.5 Variational Auto-Encoders (VAE)\n",
    "\n",
    "\n",
    "* 7.6 Deep Learning Tools\n",
    "    * Appendix: Deep learning by PyTorch\n",
    "\n",
    "\n",
    "* 7.6 Visual Cortex and Convolutional Networks\n",
    "    * Visual cortex\n",
    "    * Neocognitron\n",
    "    * Deep neural networks for neuroscience"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.1 Multi-Layer Neural Networks\n",
    "\n",
    "In the classic perceptrons and standard linear regression models, the input features $\\b{\\phi}(\\x)$ were fixed and only the output connection weights $\\w$ were changed by learning to compute the output\n",
    "\n",
    "$$\n",
    "    y = \\w^T \\b{\\phi}(\\x).\n",
    "$$\n",
    "\n",
    "The capability of such networks is dependent on what features we prepare, either by many randomly connected units in perceptrons or hand-crafted features in conventional patter classification.\n",
    "\n",
    "An alternative approach is to learn features that suite the required input-output mapping. \n",
    "Let us consider a two-layer network\n",
    "\n",
    "$$\n",
    "y = w^o_{0} + \\sum_{i=1}^M w^o_{i} h_i\n",
    "$$\n",
    "\n",
    "$$\n",
    "    h_i = g(w^h_{i0} + \\sum_{j=1}^D w^h_{ij} x_j),\n",
    "$$\n",
    "\n",
    "where $(h_1,...,h_M)$ are the outputs of *hidden units*.\n",
    "\n",
    "The *activation function* $g(\\ )$ is usually the logistic sigmoid function\n",
    "\n",
    "$$\n",
    "    g(u)=\\frac{1}{1+e^{-u}}\n",
    "$$\n",
    "\n",
    "or the rectified linear unit (ReLU)\n",
    "\n",
    "$$\n",
    "    g(u)=\\max(u, 0).\n",
    "$$\n",
    "\n",
    "Learning of the hidden unit weights $w^h_{ij}$ is possible by error gradient descent when the nonlinear function $g(\\ )$ is differentiable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_net(M=[3,4,2], label=\"\", color=None, offs=0):\n",
    "    \"\"\"Draw a multi-layer network with M units\"\"\"\n",
    "    L = len(M)    # number of layers\n",
    "    Y = [ np.linspace(0,1,M[l]+2)[1:-1] for l in range(L)]  # cell positions\n",
    "    for l in range(L):  # from bottom to top layers\n",
    "        if l<L-1:  # cross lines\n",
    "            y = np.row_stack((np.tile(Y[l],M[l+1]), np.repeat(Y[l+1],M[l])))\n",
    "            plt.plot([l, l+1], y, c=color)  # connections\n",
    "        plt.plot(l, Y[l].reshape((1,-1)), 'ow', ms=20, mec='k')  # cells\n",
    "        if type(label) is list:  # each string\n",
    "            plt.text(l, Y[l][0]/2, label[l], fontsize=16, ha='center')\n",
    "        else:    # format with l\n",
    "            plt.text(l, Y[l][0]/2, label.format(l), fontsize=16, ha='center')\n",
    "    plt.axis('off');\n",
    "    return Y  # cell positions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "draw_net([2,5,1], \"y{}\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Universality of multi-layer neural networks\n",
    "\n",
    "It has been shown that the two-layer network can approximate any nonlinear function to any desired accuracy using sufficiently large number of hidden units, called *universality* of multi-layer neural networks.\n",
    "\n",
    "While just two layers are theoretically enough, it has been shown in applications of multi-layer neural networks that three or more deeper neural networks are more capable in learning complex function approximation or classification tasks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Recurrent neural networks\n",
    "\n",
    "A discrete-time recurrent neural networks can be considered, by spatial unrolling, as a deep neural network with the same weights shared across all layers. \n",
    "\n",
    "Thus back-propagation algorithm below can be used for supervised learning of recurrent neural networks.\n",
    "It is known as *back-propatation through time*.\n",
    "\n",
    "As a corollary to the universality of two-layer neural networks, recurrent neural networks can approximate arbitrary dynamical systems.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.2 Back-Propagation Learning\n",
    "\n",
    "Here we consider a general $L$-layer network\n",
    "\n",
    "$$\n",
    "    y^l_i = g^l(w^l_{i0} + \\sum_{j=1}^{M^l} w^l_{ij} y^{l-1}_j)\n",
    "$$\n",
    "for $l=(1,...,L)$.  \n",
    "$\\y^0=\\x$ is the input vector and $\\y=\\y^L$ is the output vector.\n",
    "\n",
    "For function approximation, the output function of the last layer is usually linear, i.e. $g^L(u)=u$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we consider a *stochastic gradient* learning, in which we change the weight parameters toward the descending gradient of the error for each input data:\n",
    "\n",
    "The basic way of online learning is to minimize the output error for each input \n",
    "\n",
    "$$\n",
    "E = \\frac{1}{2}||\\y - \\y^*||^2 = \\frac{1}{2}\\sum_{i=1}^{M^L}(y^L_i - y^*_i)^2.\n",
    "$$\n",
    "\n",
    "The error gradient for the output unit weights are computes as in the linear regression\n",
    "\n",
    "$$\n",
    "\\p{E}{w^L_{ij}} = \\p{E}{y^L_i}\\p{y^L_i}{w^L_{ij}}\n",
    "  = (y^L_i - y^*_i)y^{L-1}_j. \n",
    "$$\n",
    "\n",
    "By further applying the chain rule of derivatives, the gradient for the hidden unit weights can be computed by following the network top to bottom:\n",
    "\n",
    "$$\n",
    "\\p{E}{w^l_{ij}} = \\p{E}{y^l_i}\\p{y^l_i}{w^l_{ij}}\n",
    "  = \\p{E}{y^l_i} g'^l_i y^{l-1}_j, \n",
    "$$\n",
    "where $g'^l_i$ is the derivative of the output function. For the logistic sigmoid, it is given as\n",
    "\n",
    "$$\n",
    "g'^l_i = y^l_i(1-y^l_i). \n",
    "$$\n",
    "\n",
    "In the above, $\\p{E}{y^l_i}$ takes the role of an effective error for the $i$-th unit in layer $l$, and computed iteratively\n",
    "\n",
    "$$\n",
    "\\p{E}{y^l_i} = \\sum_{k=1}^{M^{l+1}}\\p{E}{y^{l+1}_k}g'^{l+1}_k w^{l+1}_{ki}. \n",
    "$$\n",
    "\n",
    "In the example below, the red lines show the propagation of errors from the two output units to a unit in layer 2.  \n",
    "The blue lines show the propagation of errors through layer 2 to one unit in layer 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y = draw_net([2, 4, 3, 2], label=\"y{}\", color=\"gray\")\n",
    "# to layer 2\n",
    "plt.plot([3,2], np.row_stack((Y[3],np.repeat(Y[2][1],2)))+0.01, \"r\")\n",
    "# to layer 1\n",
    "plt.plot([3,2], np.row_stack((np.repeat(Y[3],3),np.tile(Y[2],2))), \"b\")\n",
    "plt.plot([2,1], np.row_stack((Y[2],np.repeat(Y[1][1],3))), \"b\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using these error gradients, all the weights can be updated by\n",
    "\n",
    "$$\n",
    "    \\Delta w^l_{ij} = - \\alpha \\p{E}{w^l_{ij}} \n",
    "$$\n",
    "where $\\alpha>0$ is a learning rate paramter.\n",
    "\n",
    "This is called *error back-propagation* learning algorithm.\n",
    "\n",
    "Below is a sample implementation of back-propagation with sigmoid hidden units and linear output units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class DNN:\n",
    "    \"\"\"Simple class for a deep neural network\"\"\"\n",
    "\n",
    "    def __init__(self, units=[2,2,1], winit=0.1):\n",
    "        \"\"\"Create a new netowrk: units:list of numbers of units\"\"\"\n",
    "        self.L = len(units)-1   # number of layers (excluding input)\n",
    "        self.M = units  # numbers of units in layers\n",
    "        # output and error vectors: layer 0 to L\n",
    "        self.Y = [np.zeros(self.M[l]) for l in range(self.L+1)]\n",
    "        self.err = [np.zeros(self.M[l]) for l in range(self.L+1)]\n",
    "        # initialize small random weights and bias\n",
    "        self.W = [np.random.normal(scale=winit, size=(self.M[l+1],self.M[l]+1)) for l in range(self.L)]\n",
    "    \n",
    "    def act(self, u):\n",
    "        \"\"\"activation function y=g(u)\"\"\"\n",
    "        return 1/(1+np.exp(-u))  # logistic sigmoid\n",
    "    def dact(self, y):\n",
    "        \"\"\"derivative of activation function dy/du|y\"\"\"\n",
    "        return y*(1-y)     # dy/du\n",
    "\n",
    "    def output(self, u):\n",
    "        \"\"\"activation function for the output unit\"\"\"\n",
    "        return u  # identity\n",
    "    def outerror(self, y, yt):\n",
    "        \"\"\"output error based on negative log likelihood\"\"\"\n",
    "        return y-yt  # difference\n",
    "    \n",
    "    def forward(self, input):\n",
    "        \"\"\"Compute the output\"\"\"\n",
    "        self.Y[0] = np.array(input).reshape((-1)) # input vector\n",
    "        for l in range(self.L):\n",
    "            z = self.W[l][:,0] + self.W[l][:,1:]@self.Y[l]  # include bias\n",
    "            self.Y[l+1] = (self.output(z) if l+1==self.L else self.act(z))  # linear output\n",
    "        return self.Y[-1].squeeze()  # last layer, as scalar if 1D\n",
    "    \n",
    "    def backprop(self, target, alpha=0.01):\n",
    "        \"\"\"Error backpropagation learning\"\"\"\n",
    "        for l in range(self.L, 0, -1):  # from top to bottom\n",
    "            if l==self.L:  # output layer\n",
    "                self.err[l] = self.outerror(self.Y[l], target)  # output error\n",
    "            else:\n",
    "                # error propapation from the layer above\n",
    "                self.err[l] = self.err[l+1]@self.W[l][:,1:] # exclude bias\n",
    "                self.err[l] *= self.dact(self.Y[l])  # error before the gain\n",
    "            # error gradient\n",
    "            dW = np.outer(self.err[l], np.concatenate(([1], self.Y[l-1])))\n",
    "            # update weights by the error gradient\n",
    "            self.W[l-1] -= alpha*dW\n",
    "        return np.dot(self.err[-l],self.err[-l])  # sum of squared error\n",
    "    \n",
    "    def train(self, inputs, targets, alpha=0.01, repeat=1):\n",
    "        \"\"\"train by a dataset\"\"\"\n",
    "        N = inputs.shape[0]  # data size\n",
    "        if repeat>1:\n",
    "            mse = np.zeros(repeat)  # record of mean square errors\n",
    "            for t in range(repeat):\n",
    "                sse = self.train(inputs, targets, alpha, repeat=1)\n",
    "                mse[t] = np.mean(sse)\n",
    "            return mse\n",
    "        else:\n",
    "            sse = np.zeros(N)  # record of sum square errors\n",
    "            for n in np.random.permutation(N):\n",
    "                y = self.forward(inputs[n])\n",
    "                sse[n] = self.backprop(targets[n], alpha)\n",
    "            return sse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "source": [
    "### Example: Sine wave"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sine wave dataset\n",
    "N = 100\n",
    "sigma = 0.05  # noise\n",
    "xr = 6  # range of x\n",
    "X = np.random.uniform(-xr, xr, N) #.reshape((N,1))\n",
    "t = np.sin(X) + np.random.normal(scale=sigma, size=N) #(N,1))\n",
    "Np = 100  # data for test/plot\n",
    "Xp = np.linspace(-xr, xr, Np) #.reshape((N,1))\n",
    "fp = np.sin(Xp)\n",
    "plt.plot(Xp, fp, \"g\")  # target\n",
    "plt.plot(X, t, \"ro\")   # training data\n",
    "plt.xlabel(\"x\"); plt.ylabel(\"y\"); plt.title(\"training data\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a network\n",
    "M = 10\n",
    "dn = DNN([1, M, 1], winit=0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize hidden and output units\n",
    "yp = np.zeros(Np)\n",
    "hp = np.zeros((Np, dn.M[1]))\n",
    "for n in range(Np):\n",
    "    yp[n] = dn.forward(Xp[n])\n",
    "    hp[n] = dn.Y[1]  # hidden units\n",
    "plt.plot(Xp, yp, \"b\")  # output unit\n",
    "plt.plot(Xp, hp)  # hidden units\n",
    "plt.xlabel(\"x\"); plt.ylabel(\"h, y\"); plt.title(\"hidden & output units\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Train with the data\n",
    "mse = dn.train(X, t, alpha=0.01, repeat=1000)\n",
    "#print(dn.W)\n",
    "print(\"mse =\", mse[-1])  # final mse\n",
    "plt.plot(mse); \n",
    "plt.xlabel(\"iteration\"); plt.ylabel(\"mse\"); plt.title(\"learning curve\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# see fitting\n",
    "yp = [dn.forward(Xp[n]) for n in range(Np)]\n",
    "plt.plot(Xp, fp, \"g\")  # target\n",
    "plt.plot(X, t, \"r.\")   # training data\n",
    "plt.plot(Xp, yp, \"b\")  # testing\n",
    "plt.xlabel(\"x\"); plt.ylabel(\"y\"); plt.title(\"testing\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Visualize hidden and output units\n",
    "for n in range(Np):\n",
    "    yp[n] = dn.forward(Xp[n])\n",
    "    hp[n] = dn.Y[1]  # hidden units\n",
    "plt.plot(Xp, yp, \"b\")  # output unit\n",
    "plt.plot(Xp, hp)  # hidden units\n",
    "plt.xlabel(\"x\"); plt.ylabel(\"h, y\"); plt.title(\"hidden & output units\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Deep NN for classification\n",
    "\n",
    "For binary classification, the target output is $y^*\\in\\{0,1\\}$ and we use a network with sigmoid output to predict the probablity \n",
    "\n",
    "$$\n",
    "p(y^*=1)=y=\\sigma(u)=\\frac{1}{1+e^{-u}}.\n",
    "$$\n",
    "\n",
    "The output error calculated as negative log likelihood\n",
    "\n",
    "$$\n",
    "E = -y^*\\log y -(1-y^*)\\log(1-y) \n",
    "$$\n",
    "is called *cross entropy error*.\n",
    "\n",
    "Its gradient with respect to the output $y$ is\n",
    "\n",
    "$$\n",
    "\\p{E}{y} = -\\frac{y^*}{y} + \\frac{1-y^*}{1-y}. \n",
    "$$\n",
    "\n",
    "The gradient with respect to the input sum $z$ is\n",
    "\n",
    "$$\n",
    "\\p{E}{z} = \\p{E}{y}\\p{y}{z}\n",
    "   = (-\\frac{y^*}{y} + \\frac{1-y^*}{1-y}) y(1-y) \n",
    "$$\n",
    "\n",
    "$$\n",
    "= -y^*(1-y) + (1-y^*)y = y - y^*. \n",
    "$$\n",
    "\n",
    "Below is a modified class for binary classification."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class DNNb(DNN):\n",
    "    \"\"\"A deep neural network for classification\"\"\"\n",
    "    # override the output function and error\n",
    "    def output(self, u):\n",
    "        \"\"\"activation function for the output unit\"\"\"\n",
    "        return self.act(u)  # sigmoid\n",
    "    # the output error dE/dz stays the same"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example: Exclusive OR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ExOr data\n",
    "N = 4\n",
    "#sigma = 0.1  # input spread\n",
    "X = np.array([[0,0], [0,1], [1,0], [1,1]])  # four corners\n",
    "yt = np.array([0, 1, 1, 0])  # ExOr\n",
    "plt.scatter(X[:,0], X[:,1], c=yt)\n",
    "plt.xlabel(\"x1\"); plt.ylabel(\"x2\"); plt.axis('square'); plt.title(\"training data\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a network\n",
    "M = 5  # hidden units\n",
    "dn = DNNb([2, M, 1], winit=0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data for test/plot\n",
    "Np = 10  \n",
    "xp = np.linspace(0, 1, Np)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Output\n",
    "yp = np.array([[dn.forward(np.array([xp1,xp2])) for xp2 in xp] for xp1 in xp])\n",
    "plt.pcolormesh(xp, xp, yp.squeeze())\n",
    "plt.colorbar()\n",
    "plt.xlabel(\"x1\"); plt.ylabel(\"x2\"); plt.axis('square');\n",
    "# Hidden unit activation boundaries\n",
    "# w0+w1*x1+w2*x2=0 -> x2=-(w0+w1*x1)/w2\n",
    "h0 = -dn.W[0][:,0]/dn.W[0][:,2]\n",
    "h1 = -(dn.W[0][:,0]+dn.W[0][:,1])/dn.W[0][:,2]\n",
    "plt.scatter(X[:,0], X[:,1], c=yt)\n",
    "plt.plot([0,1], np.row_stack((h0,h1)));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Train with the data\n",
    "mse = dn.train(X, yt, alpha=0.1, repeat=3000)\n",
    "#print(dn.W)\n",
    "print(\"mse =\", mse[-1])  # final mse\n",
    "plt.plot(mse); \n",
    "plt.xlabel(\"iteration\"); plt.ylabel(\"mse\"); plt.title(\"learning curve\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Output\n",
    "yp = np.array([[dn.forward(np.array([xp1,xp2])) for xp2 in xp] for xp1 in xp])\n",
    "plt.pcolormesh(xp, xp, yp.squeeze())\n",
    "plt.colorbar()\n",
    "plt.xlabel(\"x1\"); plt.ylabel(\"x2\"); plt.axis('square');\n",
    "# Hidden unit activation boundaries\n",
    "# w0+w1*x1+w2*x2=0 -> x2=-(w0+w1*x1)/w2\n",
    "h0 = -dn.W[0][:,0]/dn.W[0][:,2]\n",
    "h1 = -(dn.W[0][:,0]+dn.W[0][:,1])/dn.W[0][:,2]\n",
    "plt.scatter(X[:,0], X[:,1], c=yt)\n",
    "plt.plot([0,1], np.row_stack((h0,h1)));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example: Donut"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 200\n",
    "xr = 4\n",
    "X = np.random.uniform(-xr, xr, (N,2))\n",
    "y = (X[:,0]**2+X[:,1]**2)\n",
    "yt = (y>4) * (y<9)\n",
    "plt.scatter(X[:,0], X[:,1], c=yt)\n",
    "plt.xlabel(\"x1\"); plt.ylabel(\"x2\"); plt.axis('square'); plt.title(\"training data\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a network\n",
    "M = 10  # hidden units\n",
    "dn = DNNb([2, M, 1], winit=0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data for test/plot\n",
    "Np = 20  \n",
    "xp = np.linspace(-xr, xr, Np)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Output\n",
    "yp = np.array([[dn.forward(np.array([xp1,xp2])) for xp2 in xp] for xp1 in xp])\n",
    "plt.pcolormesh(xp, xp, yp.squeeze())\n",
    "plt.xlabel(\"x1\"); plt.ylabel(\"x2\"); plt.axis('square');\n",
    "plt.colorbar()\n",
    "# Hidden unit activation boundary\n",
    "# w0+w1*x1+w2*x2=0 -> x2=-(w0+w1*x1)/w2\n",
    "h0 = -(dn.W[0][:,0]-xr*dn.W[0][:,1])/dn.W[0][:,2]\n",
    "h1 = -(dn.W[0][:,0]+xr*dn.W[0][:,1])/dn.W[0][:,2]\n",
    "plt.plot([-xr,xr], np.row_stack((h0,h1)))\n",
    "# training data\n",
    "plt.scatter(X[:,0], X[:,1], c=yt, edgecolors='w');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Train with the data\n",
    "mse = dn.train(X, yt, alpha=0.01, repeat=3000)\n",
    "#print(dn.W)\n",
    "print(\"mse =\", mse[-1])  # final mse\n",
    "plt.plot(mse); \n",
    "plt.xlabel(\"iteration\"); plt.ylabel(\"mse\"); plt.title(\"learning curve\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Output\n",
    "yp = np.array([[dn.forward(np.array([xp1,xp2])) for xp2 in xp] for xp1 in xp])\n",
    "plt.pcolormesh(xp, xp, yp.squeeze())\n",
    "plt.xlabel(\"x1\"); plt.ylabel(\"x2\"); plt.axis('square');\n",
    "plt.colorbar()\n",
    "# Hidden unit activation boundary\n",
    "# w0+w1*x1+w2*x2=0 -> x2=-(w0+w1*x1)/w2\n",
    "h0 = -(dn.W[0][:,0]-xr*dn.W[0][:,1])/dn.W[0][:,2]\n",
    "h1 = -(dn.W[0][:,0]+xr*dn.W[0][:,1])/dn.W[0][:,2]\n",
    "plt.plot([-xr,xr], np.row_stack((h0,h1)))\n",
    "# training data\n",
    "plt.scatter(X[:,0], X[:,1], c=yt, edgecolors='w');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "source": [
    "### Multi-class classification\n",
    "\n",
    "For multiple classes $K$, the target output is $\\y^*=(y^*_1,...,y^*_K)$ where only one of the component $y^*_k=1$ and others are zero.\n",
    "\n",
    "In this case the standard activation function is the *softmax* function\n",
    "\n",
    "$$ y_k = p(y^*_k=1) = \\mbox{softmax}_k(\\b{u}) = \\frac{e^{u_k}}{\\sum_{j=1}^K e^{u_j}} $$\n",
    "and the cross entropy error is given as\n",
    "\n",
    "$$ E = - \\sum_{k=1}^K y^*_k \\log y_k $$\n",
    "\n",
    "Its gradient with respect to the output $y_k$ is\n",
    "\n",
    "$$ \\p{E}{y_k} = -\\frac{y^*_k}{y_k} $$\n",
    "\n",
    "The derivative of the softmax function is\n",
    "\n",
    "$$ \\p{y_k}{u_i} = \\delta_{ki}\\frac{e^{u_k}}{\\sum_{j=1}^K e^{u_j}}\n",
    "- \\frac{e^{u_k}e^{u_i}}{(\\sum_{j=1}^K e^{u_j})^2} \n",
    " = \\delta_{ki}y_k + y_k y_i $$\n",
    "\n",
    "Thus the error gradient with respect to the input sum $u_i$ is\n",
    "\n",
    "$$ \\p{E}{u_i} = \\sum_{k=1}^K \\p{E}{y_k}\\p{y_k}{u_i}\n",
    "   = \\sum_{k=1}^K -\\frac{y^*_k}{y_k}(\\delta_{ki}y_k - y_k y_i) $$\n",
    "\n",
    "$$ = \\sum_{k=1}^K -y^*_k(\\delta_{ki} + y_i) \n",
    "   = y_i - y^*_i, $$\n",
    "same as in the case of linear output and sigmoid output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.3 Deep Generative Models\n",
    "\n",
    "In back-propagation, we considered input-output mapping\n",
    "\n",
    "$$\n",
    "    \\y = f(\\x)\n",
    "$$\n",
    "to approximate the target output $\\y^*$.\n",
    "\n",
    "An opposite approach is *generative models*\n",
    "\n",
    "$$\n",
    "    \\x \\sim g(\\x|\\z)\n",
    "$$\n",
    "which assumes that the data $\\x$ are produced by a hierarchical probabilistic model with the higher level *latent variable* $\\z$, such as the class labels or low-dimensional parameter space.\n",
    "For a given input $\\x$, we consider what is the latent variable behind the data, such as the MAP estimate\n",
    "\n",
    "$$\n",
    "    \\y = \\arg\\max_\\z g(\\x|\\z)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Boltzmann machine\n",
    "\n",
    "The Boltzmann machine is a stochastic binary recurrent network consisting of *visible* and *hidden* units. The joint distribution over the visible and hidden units is given by the *energy function*. \n",
    "\n",
    "With sufficient number of hidden units and setting of the connection weights, Boltzmann machines can represent arbitrary distributions over visible units. However, learning of Boltzmann machine is *intractable*, requiring computation of posterior probabilities for exponentially large number of states.\n",
    "\n",
    "The *restricted Bolzmann machine (RBM)* is a Boltzmann machine having a layered structure, with connections only between subsequent layers and no connections within each layer.\n",
    "As described below, RBM can be trained efficiently by an algorithm called *contrastive divergence*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Deep generative networks\n",
    "\n",
    "Early deep generative models were produced by stacking RBMs and train them from the bottom to top. Examples are *deep belief networks (DBN)* and *deep Boltzmann machine (DBM)*.\n",
    "They were used for pre-training of deep networks for pattern classification, although recently simple back-propagation without pre-training has become popular (Salakhutdinov & Hinton 2012).\n",
    "\n",
    "Currently *variational auto encoders (VAE)* and *generative adversarial networks (GAN)* are the most popular and successful architectures in image generation.\n",
    "\n",
    "VAE is composed of two networks, a differentiable *generator* network that converts hidden variable $\\b{z}$ to data $\\x$, and an *encoder* network that learns to approximate the posterior probability of the hidden variables $p(\\b{z}|\\x)$.\n",
    "\n",
    "GAN also uses a differentiable generator network but is coupled with a *discriminator* network that learns to distinguish the real data and the data produced by the generator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.4 Restricted Boltzmann Machines\n",
    "\n",
    "A restricted Boltzmann machine (RBM) is a two-layer network with undirected connections between visible and hidden layers, with no connections among visible or hidden layers. This produces conditional independence of distributions of hidden units given visible units, and visible units given hidden units, which allows efficient learning (Hinton 2010)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "draw_net([4,3], label=[\"visible\", \"hidden\"])\n",
    "plt.title(\"Restricted Boltzmann machine (RBM)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we represent the state of visible units by $\\b{v}\\in\\{0,1\\}^{M_v}$ and hidden units by $\\b{h}\\in\\{0,1\\}^{M_h}$. \n",
    "For the connection weights $W\\in R^{M_v\\times M_h}$ and the bias inputs $\\b{b}$ and $\\b{c}$ for visible and hidden units, respectively, \n",
    "the *energy* function of the RBM is defined as\n",
    "\n",
    "$$ E(\\b{v},\\b{h}) = −\\b{b}^T\\b{v} − \\b{c}^T\\b{h} − \\b{v}^T W\\b{h}. $$\n",
    "\n",
    "The joint distribution of the states of the visible and and hidden units are given by\n",
    "\n",
    "$$ p(\\b{v},\\b{h}) = \\frac{1}{Z} e^{−E(\\b{v},\\b{h})}, $$\n",
    "where $Z$ is the normalizing constant called *partition function*, given by\n",
    "\n",
    "$$ Z = \\sum_\\b{v}\\sum_\\b{h} e^{−E(\\b{v},\\b{h})}. $$\n",
    "\n",
    "The marginal distribution for the visibile units is given by summing over all states of hidden units\n",
    "\n",
    "$$ p(\\b{v}) = \\frac{1}{Z}\\sum_\\b{h} e^{−E(\\b{v},\\b{h})}. $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the analytic form of maximul likelihood estimate is not available, we apply the stochastic gradient with respect to the parameters $\\theta=(W,\\b{b},\\b{c})$ . \n",
    "The gradient of the log probability of the energy-based model is given by\n",
    "\n",
    "$$ \\p{\\log p(\\b{v},\\b{h}|\\theta)}{\\theta} \n",
    " = -\\p{E(\\b{v},\\b{h};\\theta)}{\\theta} \n",
    " + \\sum_{\\b{v}}\\sum_{\\b{h}}p(\\b{v},\\b{h})\\p{E(\\b{v},\\b{h};\\theta)}{\\theta}. $$\n",
    "\n",
    "Thus the log likelihood for the visible unit state $\\b{v}_n$ $(n=1,...,N)$ is given as\n",
    "\n",
    "$$ \\p{\\log p(\\b{v}_n|\\theta)}{\\theta} \n",
    " = -\\sum_\\b{h}p(\\b{h}|\\b{v}_n)\\p{E(\\b{v}_n,\\b{h};\\theta)}{\\theta} \n",
    " + \\sum_{\\b{v}}\\sum_{\\b{h}}p(\\b{v},\\b{h})\\p{E(\\b{v},\\b{h};\\theta)}{\\theta} $$\n",
    "The first term is data-dependent statistics and the second term is data-independent statistics coming from the partition function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Learning by contrastive divergence (CD)\n",
    "\n",
    "Because RBM has no connections between hidden units, the probability of the hidden unit states is conditionally independent given the visible unit state\n",
    "\n",
    "$$ p(\\b{h}|\\b{v}) = \\prod_{j=1}^{M_h} p(h_j|\\b{v}). $$\n",
    "The probability for each hidden unit to be active is given by the sigmoid function of its input\n",
    "\n",
    "$$ p(h_j=1|\\b{v}) = \\sigma(c_j + \\sum_{i=1}^{M_v}v_i w_{ij}).$$\n",
    "This allows evaluation of the first term of the log likelihood gradient by simple sampling.\n",
    "\n",
    "Similarly, the probability for the visible units is also conditionally independent given the hidden units and each visible unit satate is given by\n",
    "\n",
    "$$ p(v_i=1|\\b{h}) = \\sigma(b_i + \\sum_{j=1}^{N_h}w_{ij}h_j). $$\n",
    "\n",
    "In order to evaluate the second term, we have to take samples from the unconstrained model dynamics. The *contrastive divergence* method takes a crude approximation of the model distribution but has been shown to perform surprizingly well.\n",
    "\n",
    "First we set the initial state of the visible units by a sample $\\b{v}^0=\\b{v}_n$. Then we repeatedly sample hidden units\n",
    "\n",
    "$$ \\b{h}^k \\sim p(\\b{h}|\\b{v}^k) $$\n",
    "and the visible units\n",
    "\n",
    "$$ \\b{v}^{k+1} \\sim p(\\b{v}|\\b{h}^k) $$\n",
    "for $K$ times. \n",
    "\n",
    "The unconstrained joint distribution is approximated by\n",
    "\n",
    "$$ p_{CD_K}(\\b{v},\\b{h})=\\delta(\\b{v}-\\b{v}_n^K)p(\\b{h}|\\b{v}_n^K), $$\n",
    "\n",
    "where $\\delta(\\x)$ is a single point distribution at $\\x=0$.\n",
    "\n",
    "The parameters are updated with a learning rate parameter $\\alpha$ as follows.\n",
    "\n",
    "$$ \\Delta W = \\alpha [\\b{v}_n p(\\b{h}|\\b{v}_n) - \\b{v}_n^K p(\\b{h}|\\b{v}_n^K)]  $$\n",
    "\n",
    "$$ \\Delta\\b{b} = \\alpha [\\b{v}_n - \\b{v}_n^K]  $$\n",
    "\n",
    "$$ \\Delta\\b{c} = \\alpha [p(\\b{h}|\\b{v}_n) - p(\\b{h}|\\b{v}_n^K)]. $$\n",
    "\n",
    "This contrastive divergence method has been shown to work well even for small $K=1$.\n",
    "Instead of this stochastic gradient for each data, *mini-batch* method to average gradients for tens of data points are often used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# RBM\n",
    "class RBM:\n",
    "    \"\"\"Restricted Boltzmann machine [Mv,Mh]\"\"\"\n",
    "\n",
    "    def __init__(self, units=[4,3], winit=0.1):\n",
    "        \"\"\"Create a new RBM: units:list of numbers of units\"\"\"\n",
    "        self.Mv, self.Mh = units  # number of visible/hidden units\n",
    "        # visible and hidden units\n",
    "        self.v = np.zeros(self.Mv)\n",
    "        self.h = np.zeros(self.Mh)\n",
    "        # initialize small random weights and bias\n",
    "        self.W = np.random.normal(scale=winit, size=(self.Mv,self.Mh))\n",
    "        self.b = np.random.normal(scale=winit, size=self.Mv)\n",
    "        self.c = np.random.normal(scale=winit, size=self.Mh)\n",
    "\n",
    "    def sigmoid(self, x):\n",
    "        return 1/(1+np.exp(-x))\n",
    "\n",
    "    def phv(self, v):\n",
    "        \"\"\"p(h|v)\"\"\"\n",
    "        return self.sigmoid(self.c + v@self.W)\n",
    "        \n",
    "    def pvh(self, h):\n",
    "        \"\"\"p(v|h)\"\"\"\n",
    "        return self.sigmoid(self.b + (self.W@h).T)\n",
    "\n",
    "    def sample(self, vin, K=1):\n",
    "        \"\"\"sample hidden and visible units K time\"\"\"\n",
    "        pv = np.zeros((K+1,self.Mv))  # probabilities\n",
    "        ph = np.zeros((K+1,self.Mh))\n",
    "        v = np.zeros((K+1,self.Mv))  # binary samples\n",
    "        h = np.zeros((K+1,self.Mh))\n",
    "        v[0] = pv[0] = vin\n",
    "        for k in range(K):\n",
    "            ph[k] = self.phv(v[k])\n",
    "            h[k] = np.random.random(self.Mh)<ph[k]\n",
    "            pv[k+1] = self.pvh(h[k])\n",
    "            v[k+1] = np.random.random(self.Mv)<pv[k+1]\n",
    "        ph[K] = self.phv(v[K])\n",
    "        return ph, h, pv, v\n",
    "    \n",
    "    def train(self, V, K=1, alpha=0.01, T=100):\n",
    "        \"\"\"train with data V=[v1,...,vN]\"\"\"\n",
    "        N = V.shape[0]  # data size\n",
    "        mse = np.zeros(T)  # record of mean square errors\n",
    "        for t in range(T): # repeat\n",
    "            for n in np.random.permutation(N):  # simple SGD without minibatch\n",
    "                ph, h, pv, v = self.sample(V[n], K)\n",
    "                self.W += alpha*(np.outer(v[0],ph[0]) - np.outer(v[-1],ph[-1]))\n",
    "                self.b += alpha*(v[0] - v[-1])\n",
    "                self.c += alpha*(ph[0] - ph[-1])\n",
    "                mse[t] += np.dot(v[0]-v[-1], v[0]-v[-1])\n",
    "            mse[t] /= N\n",
    "        return mse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Digits\n",
    "Let's try RBM with the `digits` dataset in `sklearn`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "digits = datasets.load_digits()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def imgrid(x, Nh=None):\n",
    "    \"\"\"x: N*Pv*Ph image array\"\"\"\n",
    "    N, Pv, Ph = x.shape\n",
    "    if Nh == None:\n",
    "        Nh = (int)(np.ceil(np.sqrt(N)/10)*10) if N>10 else N\n",
    "    Nv = (int)(np.ceil(N/Nh))  # rows\n",
    "    if N < Nv*Nh:   # pad by zeros\n",
    "        x = np.vstack((x, np.zeros((Nv*Nh-N,Pv,Ph))))\n",
    "    x = x.reshape((Nv,Nh,Pv,Ph))\n",
    "    x = np.transpose(x, (0,2,1,3))\n",
    "    x = x.reshape((Nv*Pv, Nh*Ph))\n",
    "    plt.imshow(x, extent=(0,Nh,Nv,0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# select only 0-3\n",
    "yt = digits.target[digits.target<4]\n",
    "x = digits.data[digits.target<4]\n",
    "x = x/np.max(x)  # 0-1\n",
    "x.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imgrid(x.reshape((-1,8,8)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a RBM for 8*8 images\n",
    "Mv = 8*8\n",
    "Mh = 4\n",
    "rb = RBM([Mv, Mh], winit=0.1)\n",
    "imgrid(rb.W.T.reshape((-1,8,8)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 500   # training set size\n",
    "mse = rb.train(x[:N], K=1, alpha=0.01, T=100)\n",
    "plt.plot(mse)\n",
    "plt.xlabel(\"t\"); plt.ylabel(\"mse\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Learned weights\n",
    "imgrid(rb.W.T.reshape((-1,8,8)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test by new data\n",
    "M = 20\n",
    "v = x[N:N+M]\n",
    "plt.subplot(3,1,1)\n",
    "imgrid(v.reshape((-1,8,8)))\n",
    "plt.subplot(3,1,2)\n",
    "h = rb.phv(v)\n",
    "plt.imshow(h.T)\n",
    "plt.subplot(3,1,3)\n",
    "vr = rb.pvh(h.T)  # reconstructed\n",
    "imgrid(vr.reshape((-1,8,8)));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.5 Variational Autoencoders (VAE)\n",
    "\n",
    "The objective of learning in a generative model is to minimize the discrepancy between the observed data distribution $p(\\x)$ and the generated data distribution\n",
    "\n",
    "$$\n",
    "    \\int_\\z p_{\\theta}(\\x|\\z) p(\\z) d\\z\n",
    "$$\n",
    "while assuming a simple low-dimensional prior distribution $p(\\z)$ of the latent varible, such as independent Gaussians.\n",
    "\n",
    "A standard way to do this is to maximize the log likelihood for data $(\\x_1,...,\\x_N)$:\n",
    "\n",
    "$$\n",
    "    \\sum_{n=1}^N\\log p_{\\theta}(\\x_n) = \\sum_n\\log\\int_\\z p_\\theta(\\x_n|\\z) p(\\z) d\\z \n",
    "$$\n",
    "\n",
    "A problem with this approach is that computing $p_{\\theta}(x)$ is often intractable: it is hard to compute the integration over all possible ranges of $\\z$ with a high dimension.\n",
    "\n",
    "In the *autoencoder* framework, we consider two networks:\n",
    "\n",
    "* The *encoder* network that maps the data $\\x$ to the latent variable$\\z$:\n",
    "\n",
    "$$\n",
    "    q_\\phi(\\z|\\x)\n",
    "$$\n",
    "\n",
    "* The *decoder* network to regenerate the data $\\z$ from the latent variable $\\z$:\n",
    "\n",
    "$$\n",
    "    p_\\theta(\\x|\\z).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "draw_net([3,5,2,5,3], [\"x\",r\"$q_\\phi(z|x)$\",\"z\",r\"$p_\\theta(x|z)$\",\"x\"]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal of learning of autoencoder is to reduce the reconstruction error, or to maximize the likelihood of regenerated data\n",
    "\n",
    "$$\n",
    "    \\sum_{n=1}^N\\log p_\\theta(\\x_n) = \n",
    "    \\sum_{n=1}^N\\log\\int_\\z p_\\theta(\\x_n|\\z) q_\\phi(\\z|\\x_n) d\\z \n",
    "$$\n",
    "while keeping the data-constrained latent variable distribution $q_\\phi(\\z|\\x_n)$ close to a desired prior distribution $p(\\z)$, such as normal gaussian.\n",
    "\n",
    "This can be achieved by maximizing the *expected variational lower bound (ELBO):\n",
    "\n",
    "$$\n",
    "    \\mathcal{L} := \n",
    "    \\sum_{n=1}^N\\log\\int_\\z p_\\theta(\\x_n|\\z) q_\\phi(\\z|\\x_n) d\\z\n",
    "    - \\mbox{KL}[q_\\phi(\\z|\\x_n) || p(\\z)]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reparametrization trick\n",
    "\n",
    "For optimizing the parameters $\\phi$ and $\\theta$ during training the generative model, we want to apply the gradient descent algorithm. We consider a model consisting of two neural networks; one encoder and one decoder network. Since our latent representation $z$ is a stochastic variable, it is mathematically not possible to backpropagate the error through the stochastic network nodes. Therefore, we need to apply a \"reparametrization trick\" (Kingma and Welling 2014). \n",
    "\n",
    "Instead of encoding the stochastic variable z, it will generate its' mean and standard deviation plus adding some Gaussian noise (see Figure 4). In that case, $\\mu$ and $\\sigma$ are deterministic, so that we can use backpropagation to update our network parameters during training. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "draw_net([3,1], [r\"$f_\\phi(x)_i$\",r\"$z_i$\"])\n",
    "plt.text(-0.1, 0.75, r\"$\\mu_i$\")\n",
    "plt.text(-0.1, 0.5, r\"$\\sigma_i$\")\n",
    "plt.text(-0.2, 0.25, r\"$\\epsilon\\sim N(0,1)$\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.6 Deep learning tools and pre-trained weights\n",
    "\n",
    "An important factor in today's boom in deep learning is the availability of open-source software tools that are optimized for general-purpose graphics processing units (GPUs). *Tensorflow* and *Pytorch* are examples of most popular tools today.\n",
    "\n",
    "In addition, some of the network weights that were trained by large data sets and achieved high performance in competitions are publicly available. Training of deep neural networks require heavy computing, but running a deep neural network for classification or prediction is much less demanding.\n",
    "By downloading such pre-trained weights, anybody can reproduce state-of-the-art performance in deep neural networks. You can further customize or improve the model by adding your own higher layers on top of those pre-trained networks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.7 Visual cortex and covolutional networks\n",
    "\n",
    "### Simple and complex receptive fields\n",
    "Soon after the discovery of orientation selective tuning of cat visual cortex neurons, Hubel & Wiesel further found that there are neurons that show *complex* receptive fields. Those neurons respond to presentation of bars in particular orientations in different positions in their large receptive fields. \n",
    "\n",
    "They suggested a networks architecutre that simple cells sum together inputs from on- or off-center cells in the thalamus, while complex cells sum together inputs from simple cells with similar orientation selectivity.\n",
    "\n",
    "% ![HW62](figures/HW62.png)\n",
    "% Simple and complex receptive fields (Hubel, Wiese 1962).\n",
    "\n",
    "### Neocognitron\n",
    "\n",
    "Inspired by the simple and complex receptive fields found by Hubel & Wiesel, Fukushima proposed a pattern classification neural network model called *Cognitron* and then its extended version *Neocognitron* (Fukushima 1980).\n",
    "\n",
    "% ![Fukushima80](figures/Fukushima80.png)\n",
    "% The neocognitron architecture (Fukushima 1980)\n",
    "\n",
    "### Analyzing neural coding by trained deep networks\n",
    "\n",
    "Inspired by the success of deep neural networks in visual object recognition, Yamins and colleagues explored how such trained deep neural networks can be used to characterize the response properties of higher visual cortex neurons, which are not easy to express by mathematical formulas.\n",
    "\n",
    "They first trained a deep neural network to perform visual object recognition task and then used the responses of the higher layers of the deep neural network to predict the activity of visual cortex neurons when the same stimulus was presented. \n",
    "\n",
    "They found that the response propereties of the highest visual cortical area, inferior temporal (IT) cortex, was well predicted using the activities of the highest layer of the trained deep network, while those in the middle level of the visual cortex, area V4, were better predicted by the responses of the intermediate layers of the deep network (Yamins et al. 2014, 2016)\n",
    "\n",
    "<img src=\"https://www.pnas.org/cms/10.1073/pnas.1403112111/asset/7a11a5c7-9ea1-44fc-8e46-b050d1abda83/assets/graphic/pnas.1403112111fig02.jpeg\" width=\"500px\">\n",
    "Figure 7.3: Prediction of visual cortical neural activities by trained deep neural networks (Yamins et al. 2014)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "source": [
    "## References\n",
    "\n",
    "* Bishop CM (2006) Pattern Recognition and Machine Learning. Springer. \n",
    "https://www.microsoft.com/en-us/research/people/cmbishop/prml-book/\n",
    "    * Chapter 5: Neural networks\n",
    "    \n",
    "    \n",
    "* Goodfellow I, Bengio Y, Courville A (2016) Deep Learning. MIT Press. (http://www.deeplearningbook.org)\n",
    "    * 3.13 Information Theory\n",
    "    * 19.4 Variational Inference and Learning\n",
    "    * Chapter 20 Deep Generative Models\n",
    "    \n",
    "    \n",
    "* Backpropagation\n",
    "    * Rumelhart DE, Hinton GE, Williams RJ (1986). Learning Representations by Back-Propagating Errors. Nature, 323, 533-536. https://doi.org/10.1038/323533a0\n",
    "\n",
    "\n",
    "* Boltzmann machines\n",
    "    * Ackley DH, Hinton GE, Sejnowski TJ (2010). A Learning Algorithm for Boltzmann Machines*. Cogn Sci, 9, 147-169. https://doi.org/10.1207/s15516709cog0901_7\n",
    "    * Salakhutdinov R, Hinton G (2012) An efficient learning procedure for deep Boltzmann machines. Neural Computation 24:1967-2006. https://doi.org/10.1162/NECO_a_00311\n",
    "\n",
    "\n",
    "* Variational autoencoders\n",
    "    * Kingma DP, Welling M (2014). Auto-encoding variational Bayes. International Conference on Learning Representations (ICLR). https://doi.org/10.48550/arXiv.1312.6114  \n",
    "    * Doersch, C. (2016). Tutorial on variational autoencoders. arXiv:1606.05908. https://doi.org/10.48550/arXiv.1606.05908  \n",
    "\n",
    "\n",
    "* Visual cortex and convolutional neural networks\n",
    "    * Hubel DH, Wiesel TN (1962). Receptive fields, binocular interaction and functional architecture in the cat's visual cortex. Journal of Physiology, 160, 106-154. https://doi.org/10.1113/jphysiol.1962.sp006837\n",
    "    * Fukushima K (1980). Neocognitron: a self organizing neural network model for a mechanism of pattern recognition unaffected by shift in position. Biol Cybern, 36, 193-202. https://doi.org/10.1007/BF00344251\n",
    "    * Yamins DL, Hong H, Cadieu CF, Solomon EA, Seibert D, DiCarlo JJ (2014). Performance-optimized hierarchical models predict neural responses in higher visual cortex. Proceedings of the National Academy of Sciences USA, 111, 8619-24. https://doi.org/10.1073/pnas.1403112111\n",
    "    * Horikawa T, Kamitani Y (2017). Generic decoding of seen and imagined objects using hierarchical visual features. Nature Communications, 8. https://doi.org/10.1038/ncomms15037\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}