SelfImplementedNeuralNetwork.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initiliase\n",
    "Initialises the weights and biases of the network as all zero ~~with uniform distributions from -1 to 1~~. The network comprises 4 layers having 28² = 784, 128, 128 and 10 neurons respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def init_network():\n",
    "    # The layer sizes\n",
    "    n0 = 28*28\n",
    "    n1 = 128\n",
    "    n2 = 128\n",
    "    n3 = 10\n",
    "    \n",
    "    init_weights = lambda m, n: np.random.uniform(-1,1,(n,m))\n",
    "    init_biases = lambda m: np.random.uniform(-1,1,(m,))\n",
    "    # init_weights = lambda m, n: np.zeros((n,m)) # <- Attention: This function flips n and m for convenient input\n",
    "    #  init_biases = lambda m: np.zeros((m,))\n",
    "    \n",
    "    layer0_1_weights = init_weights(n0,n1)\n",
    "    layer0_1_biases = init_biases(n1)\n",
    "    layer0_1 = [layer0_1_weights, layer0_1_biases]\n",
    "    \n",
    "    layer1_2_weights = init_weights(n1,n2)\n",
    "    layer1_2_biases = init_biases(n2)\n",
    "    layer1_2 = [layer1_2_weights, layer1_2_biases]\n",
    "    \n",
    "    layer2_3_weights = init_weights(n2,n3)\n",
    "    layer2_3_biases = init_biases(n3)\n",
    "    layer2_3 = [layer2_3_weights, layer2_3_biases]\n",
    "    \n",
    "    return [layer0_1, layer1_2, layer2_3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Activation functions\n",
    "The *ReLU* (Rectified Linear Unit) activation function is defined as $\\text{ReLU}(x)=\\text{max}(0,x)$ and is the go-to activation function for hidden layers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def relu(x):\n",
    "    return np.maximum(0,x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The *Sotfmax* acitvation function is a generalization of the *Sigmoid* function. It is applied to a 1D-Array of values which are squished to the interval $[0,1)$ such that the sum of all values is 1. Thereby, we get a true probability distribution which is quite convenient. It is defined as\n",
    "\n",
    "$$\\sigma(x)_i=\\frac{e^{x_i}}{\\sum_{j}e^{x_j}}\\quad\\forall\\;i.$$\n",
    "\n",
    "To be compatible with forward propagation on a 2D-Array of multiple training examples we extend the function to also take 2D-Arrays as a parameter and perform the softmax row-wise. (How this works in detail is annotated in the comments.)\n",
    "\n",
    "Moreover, if the arrays contain big values we quickly get an overflow for `np.exp`. Therefore, we subtract the maximum value of each row (or for the 1D case: simply the maximum) from all values which doesn't alter the endresult of our computation, but rids us of some possible (pseudo-)infinities along the way. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def softmax(x):\n",
    "    # If x is a 2D-array of multiple training exmaples, we transpose so that axis=0 is the axis of a single \n",
    "    # training example. (We can't simply use axis=1 because that would break if we use the function for a \n",
    "    # 1D-array, i.e. only one training example. Transposing first and then using axis=0 works in both cases.)\n",
    "    x = x.T\n",
    "    # To avoid overflow of np.exp (doesn't alter the value)\n",
    "    x = x - np.max(x, axis=0)\n",
    "    # The main calculation\n",
    "    result = np.exp(x) / np.sum(np.exp(x), axis=0)\n",
    "    # Tranpsose back to normal form where each row is a training example (only makes a difference if array is 2D;\n",
    "    # doesn't change anything if array is 1D.)\n",
    "    return result.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Forward propagation\n",
    "Here we implement forward propagation. This function can either take a 1D-Array of a single training example or a 2D-Array of multiple training examples where each row is a training example. The transformations (`np.array.T`) are necessary so that the dot product and the vector addition work in the 2D case. They can be ignored in the 1D case.\n",
    "\n",
    "For backpropagation we not only need the endresult, but all intermediate results. If the option `intermediate_results` is set to `True` (default is `False`) we return a tuple $(y^{(1)},o^{(1)},y^{(2)},o^{(2)},y^{(3)},o^{(3)})$. Otherwise we simply return $o^{(3)}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def feed_forward(data, network, intermediate_results=False):\n",
    "    # You can either pass only one training example,\n",
    "    # then we add a new axis so that axis=0 is the axis\n",
    "    # of training examples.\n",
    "    if data.ndim == 1:\n",
    "        x = data[np.newaxis, :]\n",
    "    elif data.ndim == 2:\n",
    "        x = data\n",
    "    else:\n",
    "        raise Exception(\"Too many dimensions\")\n",
    "        \n",
    "    y1 = np.einsum(\"jk,ik\", network[0][0], x) + network[0][1]\n",
    "    o1 = relu(y1)\n",
    "    \n",
    "    y2 = np.einsum(\"jk,ik\", network[1][0], o1) + network[1][1]\n",
    "    o2 = relu(y2)\n",
    "    \n",
    "    y3 = np.einsum(\"jk,ik\", network[2][0], o2) + network[2][1]\n",
    "    o3 = softmax(y3)\n",
    "    \n",
    "    if not intermediate_results:\n",
    "        return o3\n",
    "    else:\n",
    "        return y1, o1, y2, o2, y3, o3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def feed_forward_old(data, network, intermediate_results=False): \n",
    "    y1s, o1s, y2s, o2s, y3s, o3s = [], [], [], [], [], []\n",
    "    \n",
    "    for datapoint in data:\n",
    "        x = datapoint\n",
    "        \n",
    "        y1 = np.dot(network[0][0], x) + network[0][1]\n",
    "        o1 = relu(y1)\n",
    "\n",
    "        y2 = np.dot(network[1][0], o1) + network[1][1]\n",
    "        o2 = relu(y2)\n",
    "\n",
    "        y3 = np.dot(network[2][0], o2) + network[2][1]\n",
    "        o3 = softmax(y3)\n",
    "        \n",
    "        y1s.append(y1)\n",
    "        y2s.append(y2)\n",
    "        y3s.append(y3)\n",
    "        o1s.append(o1)\n",
    "        o2s.append(o2)\n",
    "        o3s.append(o3)\n",
    "        \n",
    "    if not intermediate_results:\n",
    "        return np.array(o3s)\n",
    "    else:\n",
    "        return (np.array(y1s), \n",
    "               np.array(o1s), \n",
    "               np.array(y2s), \n",
    "               np.array(o2s),\n",
    "               np.array(y3s), \n",
    "               np.array(o3s))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get data\n",
    "We use the mnist handwritten digit classification data which is shipped with keras."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "network = init_network()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tensorflow as tf\n",
    "(x_train_raw, y_train), (x_test_raw, y_test) = tf.keras.datasets.mnist.load_data()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Flatten the data to make it a well-behaved input for our neural network:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_train = x_train_raw.reshape(60000, 28*28)\n",
    "x_test = x_test_raw.reshape(10000, 28*28)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparing old and new implementation of feed forward"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 11.9 s\n"
     ]
    }
   ],
   "source": [
    "%time a1, a2, a3, a4, a5, a6 = feed_forward(x_train, network, intermediate_results=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 22.9 s\n"
     ]
    }
   ],
   "source": [
    "%time b1, b2, b3, b4, b5, b6 = feed_forward_old(x_train, network, intermediate_results=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 6.21 s\n"
     ]
    }
   ],
   "source": [
    "%time a = feed_forward(x_train, network)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 4.35 s\n"
     ]
    }
   ],
   "source": [
    "%time b = feed_forward_old(x_train, network)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "New implementation in NumPy seems to be a little better than the python-loop implementation, but the difference is smaller than I expected. Let's check whether the `einsum`s are really correct:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-24895.13949769, -43819.82993956,  -8045.26730106,  21888.33303494,\n",
       "       -32186.98355666,  20470.37117075,  10200.70150049, -66351.75941768,\n",
       "        12147.13547264, -21356.35795113])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a5[5000]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-24895.13949769, -43819.82993956,  -8045.26730106,  21888.33303494,\n",
       "       -32186.98355666,  20470.37117075,  10200.70150049, -66351.75941768,\n",
       "        12147.13547264, -21356.35795113])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b5[5000]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(a5, b5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(a6,b6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also check whether the softmax function works correctly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "a = np.array([softmax(elem) for elem in b5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = softmax(b5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(a,b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The feed forward function can now either take a single example or an array of examples!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cost function\n",
    "First, the loss function, i.e. the error for a single training example. There are a multitude of loss functions, though there is an \"industry-standard\" for each type of problem. While the simpler and better-known MSE (*mean squared error*) is used for regression problems, one uses *cross entropy loss* for categorization problems (which fits our case precisely).\n",
    "\n",
    "Given a prediction $\\hat{y_i}$ and the correct solution $y_i$ (in our case for $0\\leq i<10$), we can define *cross entropy loss* as \n",
    "$$H(y,\\hat{y})=-\\sum_{i=0}^9 y_i\\cdot \\text{log}(\\hat{y}_i).$$\n",
    "Since we know that $y_i$ is zero for all but one $i$ (for which it has value 1), we can simplify this to\n",
    "$$H(y,\\hat{y})=-\\text{log}(\\hat{y}_k)$$\n",
    "where $k$ is the index of the correct solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cross_entropy_loss(y_hat, y):\n",
    "    y_val = y_hat[y]\n",
    "    # To avoid division by zero\n",
    "    if y_val == 0:\n",
    "        # #The numpy epsilon (np.finfo(float).eps) is apparently not really \n",
    "        # # the smallest number possible; by trying out I found the limit for \n",
    "        # # np.log to not throw an error was at about eps*10^(-307).\n",
    "        # y_val = np.finfo(\"float64\").eps*10**-307\n",
    "        # Update: \n",
    "        # np.spacing(0) gives the smallest number bigger than 0\n",
    "        y_val = np.spacing(0)\n",
    "    return -np.log(y_val)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = feed_forward(x_train, network)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "744.4400719213812"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cross_entropy_loss(result[10], y_train[10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we define the cost function, i.e. the error for a whole training set (or subset). In other words\n",
    "$$E(Y,\\hat{Y})=\\sum_{\\forall\\;y,\\,\\hat{y}}H(y,\\hat{y}).$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cross_entropy_cost(y_hats, ys):\n",
    "    return np.sum(cross_entropy_loss(y_hat, y) for y_hat, y in zip(y_hats, ys))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "39289585.71953629"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cross_entropy_cost(result, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Backpropagation\n",
    "Now, finally the training will be implemented.\n",
    "\n",
    "### Definitions\n",
    "\n",
    "- Let the layers be labeled as $n=0,\\dots,3$.\n",
    "- Let $x_i$ be the $i$-th $x$-value, that is the $i$-th input of the network ($0\\leq i<28^2$).\n",
    "- Let $y_i^{(n)}$ be the values of the $i$-th neuron of layer $n$ *before* the activation function ($0\\leq i<128$ for $n\\in\\{1,2\\}$, $0\\leq i<10$ for $n=3$).\n",
    "- Let $o_i^{(n)}$ be the values of the $i$-th neuron of layer $n$ *after* the activation function.\n",
    "    - $o_i^{(n)}=\\text{ReLU}(y_i^{n})$ for $n\\in\\{1,2\\}$\n",
    "    - $o_i^{(3)}=\\sigma(y^{(3)})_i$\n",
    "- Let $w_{ij}^{(n)}$ represent the weight on the connection of the $i$-th neuron in layer n to the $j$-th neuron in layer n+1. (Can be found in `network[n][0][j][i]`.)  \n",
    "- Let $b_i^{(n)}$ represent the bias of the $i$-th neuron in layer n+1. (Can be found in `network[n][1][i]`.)\n",
    "- Within the context of the cost/loss function we use the following definitions:\n",
    "    - $x$ and $y$ are the input and expected output of a single training example;\n",
    "    - Within the subsequent text, I will also use $k$ (as in **k**orrect solution) synonymously for $y$\n",
    "    - $X$ and $Y$ are the set (or a subset) of training examples; \n",
    "    - $\\hat{y}$ is the actual output of our network, in other words: $\\hat{y}=o^{(3)}$; \n",
    "    - $\\hat{Y}$ is the set of actual outputs for a set of inputs $X$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stochastic Gradient Descent\n",
    "\n",
    "We split the computation of the gradient into six parts:\n",
    "- Gradient of the layer 2 to 3 weights ($w^{(2)}$)\n",
    "- Gradient of the layer 3 biases ($b^{(2)}$)\n",
    "- Gradient of the layer 1 to 2 weights ($w^{(1)}$)\n",
    "- Gradient of the layer 2 biases ($b^{(1)}$)\n",
    "- Gradient of the layer 0 to 1 weights ($w^{(0)}$)\n",
    "- Gradient of the layer 1 biases ($b^{(0)}$)\n",
    "\n",
    "#### Layer 2 to 3 weights\n",
    "\n",
    "##### Overview\n",
    "Our goal is to determine\n",
    "$$\\frac{\\partial}{\\partial w_{ij}^{(2)}} \\;E(Y,\\hat{Y})$$\n",
    "$$=\\quad\\sum_{\\forall y,\\hat{y}} \\frac{\\partial}{\\partial w_{ij}^{(2)}} \\;H(y,\\hat{y})$$\n",
    "$$=\\quad\\sum_{\\forall y,\\hat{y}} \\frac{\\partial}{\\partial w_{ij}^{(2)}} \\;\\left(-\\text{log}(\\hat{y}_k)\\right),$$\n",
    "where $k$ is the correct solution for a given sample."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Cost function derivative with respect to a single weight\n",
    "Hence, let us consider\n",
    "$$\\frac{\\partial}{\\partial w_{ij}^{(2)}} \\;\\left(-\\text{log}(\\hat{y}_k)\\right)$$\n",
    "$$=\\quad \\frac{\\partial}{\\partial w_{ij}^{(2)}} \\;\\left(-\\text{log}(o^{(3)}_k)\\right)$$\n",
    "$$=\\quad \\frac{\\partial}{\\partial o^{(3)}_k}\\;\\left(-\\text{log}(o^{(3)}_k)\\right)\\cdot\\frac{\\partial o^{(3)}_k}{\\partial w_{ij}^{(2)}}$$\n",
    "$$=\\quad -\\frac{1}{o^{(3)}_k}\\cdot \\frac{\\partial}{\\partial w_{ij}^{(2)}} \\;\\sigma(y^{(3)})_k$$\n",
    "$$=\\quad -\\frac{1}{o^{(3)}_k}\\cdot \\frac{\\partial}{\\partial w_{ij}^{(2)}} \\;\\frac{e^{y^{(3)}_k}}{\\sum_{\\forall r}e^{y^{(3)}_r}}.$$\n",
    "$$=\\quad -\\frac{1}{o^{(3)}_k}\\cdot \\left(\\sum_{\\forall s}\\frac{\\partial}{\\partial y^{(3)}_s} \\;\\frac{e^{y^{(3)}_k}}{\\sum_{\\forall r}e^{y^{(3)}_r}} \\cdot \\frac{\\partial y^{(3)}_s}{\\partial w_{ij}^{(2)}} \\right).$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To simplify it more, let us focus on\n",
    "$$\\frac{\\partial}{\\partial y^{(3)}_s} \\;\\frac{e^{y^{(3)}_k}}{\\sum_{\\forall r}e^{y^{(3)}_r}}$$\n",
    "$$=\\quad \\frac{\\partial}{\\partial e^{y^{(3)}_s}} \\;\\frac{e^{y^{(3)}_k}}{\\sum_{\\forall r}e^{y^{(3)}_r}} \\cdot \\frac{\\partial e^{y^{(3)}_s}}{\\partial y^{(3)}_s}$$\n",
    "$$=\\quad e^{y^{(3)}_s}\\cdot\\frac{\\partial}{\\partial e^{y^{(3)}_s}} \\;\\frac{e^{y^{(3)}_k}}{\\sum_{\\forall r}e^{y^{(3)}_r}}$$\n",
    "$$=\\quad e^{y^{(3)}_s} \\cdot \n",
    "\\frac{ \n",
    "    \\frac{\\partial}{ \\partial e^{y^{(3)}_s} } \\; e^{y^{(3)}_k} \\cdot \\left(\\sum_{\\forall r}e^{y^{(3)}_r}\\right) \n",
    "    -  e^{y^{(3)}_k} \\cdot \\frac{\\partial}{ \\partial e^{y^{(3)}_s} } \\; \\left(\\sum_{\\forall r}e^{y^{(3)}_r}\\right)\n",
    "}{\n",
    "    \\left(\\sum_{\\forall r}e^{y^{(3)}_r}\\right)^2\n",
    "}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we need to make a distinction between $s=k$ and $s\\neq k$. If $s=k$, we get\n",
    "$$e^{y^{(3)}_s} \\cdot\n",
    "\\frac{ \n",
    "    \\sum_{\\forall r}e^{y^{(3)}_r} \n",
    "    -  e^{y^{(3)}_k}\n",
    "}{\n",
    "    \\left(\\sum_{\\forall r}e^{y^{(3)}_r}\\right)^2\n",
    "}$$\n",
    "$$=\\quad \n",
    "\\frac{\n",
    "    e^{y^{(3)}_s}\n",
    "}{\n",
    "    \\sum_{\\forall r}e^{y^{(3)}_r}\n",
    "}\n",
    "\\cdot\n",
    "\\frac{\n",
    "    \\sum_{\\forall r}e^{y^{(3)}_r} \n",
    "    -  e^{y^{(3)}_k}\n",
    "}{\n",
    "    \\sum_{\\forall r}e^{y^{(3)}_r}\n",
    "}\n",
    "$$\n",
    "$$=\\quad \\sigma(y^{(3)})_s \\cdot \\left( 1 - \\sigma(y^{(3)})_k \\right)$$\n",
    "$$=\\quad o^{(3)}_s \\cdot \\left( 1 - o^{(3)}_k \\right).$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If $s\\neq k$, we get\n",
    "$$e^{y^{(3)}_s} \\cdot\n",
    "\\frac{ \n",
    "    -  e^{y^{(3)}_k}\n",
    "}{\n",
    "    \\left(\\sum_{\\forall r}e^{y^{(3)}_r}\\right)^2\n",
    "}$$\n",
    "$$=\\quad - \\sigma(y^{(3)})_s \\cdot \\sigma(y^{(3)})_k$$\n",
    "$$=\\quad - o^{(3)}_s \\cdot o^{(3)}_k.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Furthermore, we need to evaluate\n",
    "$$\\frac{\\partial y^{(3)}_s}{\\partial w_{ij}^{(2)}}.$$\n",
    "We know that\n",
    "$$y^{(3)}_s = \\sum_{\\forall t} w_{ts}^{(2)}o^{(2)}_t + b_s^{(2)}.$$\n",
    "If $j=s$, we have\n",
    "$$\\frac{\\partial y^{(3)}_s}{\\partial w_{ij}^{(2)}} = o^{(2)}_i.$$\n",
    "Else, if $j\\neq s$, we have\n",
    "$$\\frac{\\partial y^{(3)}_s}{\\partial w_{ij}^{(2)}} = 0.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The equation we considered before was\n",
    "$$-\\frac{1}{o^{(3)}_k}\\cdot \\left(\\sum_{\\forall s}\\frac{\\partial}{\\partial y^{(3)}_s} \\;\\frac{e^{y^{(3)}_k}}{\\sum_{\\forall r}e^{y^{(3)}_r}} \\cdot \\frac{\\partial y^{(3)}_s}{\\partial w_{ij}^{(2)}} \\right),$$\n",
    "where -- as we see from our last result -- all terms but one turn out to be zero. What remains is\n",
    "$$-\\frac{1}{o^{(3)}_k}\\cdot \\left(\\frac{\\partial}{\\partial y^{(3)}_j} \\;\\frac{e^{y^{(3)}_k}}{\\sum_{\\forall r}e^{y^{(3)}_r}}\\cdot o_i^{(2)}\\right).$$\n",
    "Now, we have to distinguish between $j=k$ and $j\\neq k$, which yields us for $j=k$\n",
    "$$\\frac{\\partial}{\\partial w_{ij}^{(2)}} \\;H(y,\\hat{y}) = \n",
    "    -\\frac{1}{o^{(3)}_k}\n",
    "    \\cdot \n",
    "    o^{(3)}_j \\cdot \\left( 1 - o^{(3)}_k \\right)\n",
    "    \\cdot \n",
    "    o_i^{(2)},$$\n",
    "and for $j\\neq k$\n",
    "$$\\frac{\\partial}{\\partial w_{ij}^{(2)}} \\;H(y,\\hat{y}) = \n",
    "    \\frac{1}{o^{(3)}_k}\n",
    "    \\cdot\n",
    "    o^{(3)}_j \\cdot o^{(3)}_k\n",
    "    \\cdot \n",
    "    o_i^{(2)}$$\n",
    "$$=\\quad o^{(3)}_j\\cdot o_i^{(2)}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Unifying it into a single matrix\n",
    "Our weights matrix $W^{(2)}$ has dimensions $n^{(3)} \\times n^{(2)}$. (This means the indices $i$ and $j$ in $w_{ij}^{(2)}$ are actually the wrong way around.) Let us call the gradient of $H(y,\\hat y)$ w.r.t. the weight matrix \n",
    "$$\\nabla_{W^{(2)}} \\,H(y,\\hat y).$$ \n",
    "This gradient shall be a $n^{(3)} \\times n^{(2)}$ matrix as well. (Hence, contrary to popular usage, we extend the notion of gradient from simply being a vector to also being a matrix; this is useful because we also treat the weights as matrices which simplifies both updating the weights according to the gradient and computing the gradient non-iteratively easier.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is a distinction to be made between $j=k$ and $j\\neq k$, i.e. the $k$-th row in the gradient matrix needs a special treating. Luckily, considering the case $j=k$ as derived above we see\n",
    "$$\\frac{\\partial}{\\partial w_{ij}^{(2)}} \\;H(y,\\hat{y}) = \n",
    "    -\\frac{1}{o^{(3)}_k}\n",
    "    \\cdot \n",
    "    o^{(3)}_j \\cdot \\left( 1 - o^{(3)}_k \\right)\n",
    "    \\cdot \n",
    "    o_i^{(2)}$$\n",
    "$$= \\quad\n",
    "    -\\frac{1}{o^{(3)}_k} \\cdot o^{(3)}_j \\cdot o_i^{(2)}\n",
    "    \\; + \\;\n",
    "    \\frac{1}{o^{(3)}_k} \\cdot o^{(3)}_j \\cdot o^{(3)}_k \\cdot o_i^{(2)}$$\n",
    "$$= \\quad\n",
    "    - o_i^{(2)}\n",
    "    \\; + \\;\n",
    "    o^{(3)}_j \\cdot o_i^{(2)},$$\n",
    "where the second summand coincides with the case $j\\neq k$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus, we are able to split the unification of the many partial derivatives into a single gradient matrix into two steps\n",
    "1. Compute the matrix where the element in the $i$-th column and $j$-th row is $o^{(3)}_j \\cdot o_i^{(2)}$.\n",
    "2. Subtract from the $k$-th row a (row) vector where the $i$-th column is $o_i^{(2)}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The computation in step 1 is\n",
    "$$\\quad o^{(3)} \\otimes o^{(2)}$$\n",
    "where $a\\otimes b = a\\cdot b^T$ denotes the *outer product*, *matrix product* or *tensor product* of two vectors (called `np.outer`). In the second step, we simply subtract $o^{(2)}$ (as a row vector) from the $k$-th column of the result of step 1. (Expressing this operation mathematically is not necessary, here we just have to consider implementing it efficiently.) \n",
    "\n",
    "Lastly, we can sum this matrix for multiple training examples\n",
    "$$\\nabla_{W^{(2)}} \\;E(Y,\\hat{Y}) = \\sum_{\\forall\\;y,\\hat y} \\nabla_{W^{(2)}} \\;H(y,\\hat y).$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Layer 3 biases\n",
    "Here, we examine\n",
    "$$\\frac{\\partial}{\\partial b_{j}^{(2)}} \\;E(Y,\\hat{Y})$$\n",
    "$$=\\quad\\sum_{\\forall y,\\hat{y}} \\frac{\\partial}{\\partial b_{j}^{(2)}} \\;H(y,\\hat{y}),$$\n",
    "which can be derived in a very analogous way. The key difference arises when determining\n",
    "$$\\frac{\\partial y^{(3)}_s}{\\partial b_{j}^{(2)}}.$$\n",
    "If $j=s$, we have\n",
    "$$\\frac{\\partial y^{(3)}_s}{\\partial b_{j}^{(2)}} = 1,$$\n",
    "instead of $o_i^{(2)}$ which we would have got for for $w_{ij}^{(2)}$. Else, if $j\\neq s$, we still have\n",
    "$$\\frac{\\partial y^{(3)}_s}{\\partial b_{j}^{(2)}} = 0.$$\n",
    "Replacing every occasion of $o_i^{(2)}$ with $1$, we see that $\\nabla_{b^{(2)}} \\;H(y,\\hat y)$ is a $n^{(3)}$-dimensional vector where the $j$-th row is $o_j^{(3)}$ if $j\\neq k$ and $o_j^{(3)}-1$ if $j=k$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Layer 1 to 2 weights\n",
    "Now, let's go a layer deeper! For determining \n",
    "$$\\frac{\\partial}{\\partial w_{ij}^{(1)}} \\;E(Y,\\hat{Y})$$\n",
    "we can follow the same steps as for $w_{ij}^{(2)}$ up until\n",
    "$$\\quad -\\frac{1}{o^{(3)}_k}\\cdot \\left(\\sum_{\\forall s}\\frac{\\partial}{\\partial y^{(3)}_s} \\;\\frac{e^{y^{(3)}_k}}{\\sum_{\\forall r}e^{y^{(3)}_r}} \\cdot \\frac{\\partial y^{(3)}_s}{\\partial w_{ij}^{(1)}} \\right)$$\n",
    "$$= \\quad -\\frac{1}{o^{(3)}_k}\\cdot \n",
    "\\left(\n",
    "    \\sum_{\\forall s / k} \\left(- o^{(3)}_s \\cdot o^{(3)}_k\\right)\n",
    "    \\cdot \\frac{\\partial y^{(3)}_s}{\\partial w_{ij}^{(1)}}\n",
    "    \\;+\\;\n",
    "    o^{(3)}_k \\cdot \\left( 1 - o^{(3)}_k \\right)\n",
    "    \\cdot \\frac{\\partial y^{(3)}_k}{\\partial w_{ij}^{(1)}}\n",
    "\\right)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$= \\quad -\\frac{1}{o^{(3)}_k}\\cdot \n",
    "\\left(\n",
    "    \\sum_{\\forall s} \\left(- o^{(3)}_s \\cdot o^{(3)}_k\\right)\n",
    "    \\cdot \\frac{\\partial y^{(3)}_s}{\\partial w_{ij}^{(1)}}\n",
    "    \\;+\\;\n",
    "    o^{(3)}_k\n",
    "    \\cdot \\frac{\\partial y^{(3)}_k}{\\partial w_{ij}^{(1)}}\n",
    "\\right)$$\n",
    "$$= \\quad \n",
    "\\sum_{\\forall s} o^{(3)}_s\n",
    "\\cdot \\frac{\\partial y^{(3)}_s}{\\partial w_{ij}^{(1)}}\n",
    "\\;-\\;\n",
    "\\frac{\\partial y^{(3)}_k}{\\partial w_{ij}^{(1)}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us examine\n",
    "$$\\frac{\\partial}{\\partial w_{ij}^{(1)}} \\;y^{(3)}_s$$\n",
    "$$=\\quad \\frac{\\partial}{\\partial w_{ij}^{(1)}} \\; \\left(\\sum_{\\forall t} w_{ts}^{(2)}o^{(2)}_t + b_s^{(2)}\\right)$$$$=\\quad \\frac{\\partial}{\\partial w_{ij}^{(1)}} \\; \\left(\\sum_{\\forall t} w_{ts}^{(2)}o^{(2)}_t + b_s^{(2)}\\right)$$\n",
    "$$=\\quad \\sum_{\\forall t} \\frac{\\partial}{\\partial w_{ij}^{(1)}} \\; w_{ts}^{(2)}o^{(2)}_t $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Applying the chain rule)\n",
    "$$=\\quad \\sum_{\\forall t} \\frac{\\partial}{\\partial o^{(2)}_t} \\; \\left(w_{ts}^{(2)}o^{(2)}_t\\right) \\cdot \\frac{\\partial o^{(2)}_t}{\\partial w_{ij}^{(1)}}$$\n",
    "$$=\\quad \\sum_{\\forall t} w_{ts}^{(2)} \\cdot \\frac{\\partial o^{(2)}_t}{\\partial w_{ij}^{(1)}}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us, furthermore, concentrate on\n",
    "$$\\frac{\\partial}{\\partial w_{ij}^{(1)}} \\;o^{(2)}_t$$\n",
    "$$=\\quad \\frac{\\partial}{\\partial w_{ij}^{(1)}} \\; \\text{ReLU}(y^{(2)}_t)$$\n",
    "$$=\\quad \\frac{\\partial}{\\partial y^{(2)}_t} \\; \\text{ReLU}(y^{(2)}_t) \\cdot \\frac{\\partial y^{(2)}_t}{\\partial w_{ij}^{(1)}}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the derivative of the ReLU function we have to consider that it is not everywhere differentiable. We (arbitrarily) choose that the derivative at point 0 is 1. If $y_t^{(2)}< 0$, we then get\n",
    "$$\\frac{\\partial}{\\partial w_{ij}^{(1)}} \\;o^{(2)}_t$$\n",
    "$$=\\quad 0.$$\n",
    "Else, if $y_t^{(2)} \\geq 0$, we get\n",
    "$$\\frac{\\partial}{\\partial w_{ij}^{(1)}} \\;o^{(2)}_t$$\n",
    "$$=\\quad 1 \\cdot \\frac{\\partial y^{(2)}_t}{\\partial w_{ij}^{(1)}}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, lastly, let us have a look at\n",
    "$$\\frac{\\partial}{\\partial w_{ij}^{(1)}} y^{(2)}_t$$\n",
    "$$=\\quad \\frac{\\partial}{\\partial w_{ij}^{(1)}} \n",
    "\\left( \\sum_{\\forall r} w_{rt}^{(1)}o_{r}^{(1)} + b_{t}^{(1)}\\right)$$\n",
    "$$=\\quad \\frac{\\partial}{\\partial w_{ij}^{(1)}} w_{it}^{(1)}o_{i}^{(1)}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus, if $j\\neq t$\n",
    "$$\\frac{\\partial}{\\partial w_{ij}^{(1)}} y^{(2)}_t = 0$$\n",
    "and if $j=t$\n",
    "$$\\frac{\\partial}{\\partial w_{ij}^{(1)}} y^{(2)}_t = o_{i}^{(1)}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Going back up the rabbit hole\n",
    "We evaluated all the partial derivatives. Now, we have to put them together to get back to where we started off. As an intermediate step, we had\n",
    "$$\\frac{\\partial}{\\partial w_{ij}^{(1)}} \\;y^{(3)}_s \n",
    "\\;=\\; \n",
    "\\sum_{\\forall t} w_{ts}^{(2)} \\cdot \\frac{\\partial o^{(2)}_t}{\\partial w_{ij}^{(1)}}.$$\n",
    "We found out that $\\frac{\\partial o^{(2)}_t}{\\partial w_{ij}^{(1)}}$ is definitely zero if $j\\neq t$, which yields\n",
    "$$\\frac{\\partial}{\\partial w_{ij}^{(1)}} \\;y^{(3)}_s \n",
    "\\;=\\; \n",
    "w_{js}^{(2)} \\cdot \\frac{\\partial o^{(2)}_j}{\\partial w_{ij}^{(1)}}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therefore, if $y^{(2)}_j<0$ we get\n",
    "$$\\frac{\\partial}{\\partial w_{ij}^{(1)}} \\;y^{(3)}_s \n",
    "\\;=\\; 0.$$\n",
    "Else, if $y^{(2)}_j\\geq 0$ we have\n",
    "$$\\frac{\\partial}{\\partial w_{ij}^{(1)}} \\;y^{(3)}_s \n",
    "\\;=\\; \n",
    "w_{js}^{(2)} \\cdot o^{(1)}_i.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we are able to take a look at the uppermost level, where we had\n",
    "$$\\frac{\\partial}{\\partial w_{ij}^{(1)}} \\;E(Y,\\hat{Y})\n",
    "\\;=\\;\n",
    "\\sum_{\\forall s} o^{(3)}_s\n",
    "\\cdot \\frac{\\partial y^{(3)}_s}{\\partial w_{ij}^{(1)}}\n",
    "\\;-\\;\n",
    "\\frac{\\partial y^{(3)}_k}{\\partial w_{ij}^{(1)}}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If $y^{(2)}_j<0$ we get\n",
    "$$\\frac{\\partial}{\\partial w_{ij}^{(1)}} \\;E(Y,\\hat{Y})\n",
    "\\;=\\;0$$\n",
    "Else, if $y^{(2)}_j\\geq 0$ we get\n",
    "$$\\frac{\\partial}{\\partial w_{ij}^{(1)}} \\;E(Y,\\hat{Y})\n",
    "\\;=\\;\n",
    "\\sum_{\\forall s} o^{(3)}_s\n",
    "\\cdot \\left(w_{js}^{(2)} o^{(1)}_i\\right)\n",
    "\\;-\\;\n",
    "w_{jk}^{(2)} \\cdot o^{(1)}_i.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Unifying it into a single matrix\n",
    "Our weights matrix $W^{(1)}$ has dimensions $n^{(2)}\\times n^{(1)}$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "def stochastic_gradient_descent(x_train, y_train, network, batch_size=1000, epochs=3, learning_rate=0.0001):\n",
    "    \n",
    "    num_examples = len(x_train)\n",
    "    iterations_per_epoch = num_examples//batch_size\n",
    "    \n",
    "    ###########################\n",
    "    # Iterating through epochs\n",
    "    for _ in range(epochs):\n",
    "        ###############################\n",
    "        # Iterate through the dataset within an epoch\n",
    "        # where i is the current position/offset\n",
    "        # within the dataset\n",
    "        for i in range(iterations_per_epoch):\n",
    "            X = x_train[i:i+batch_size]\n",
    "            Y = y_train[i:i+batch_size]\n",
    "\n",
    "            y1, o1, y2, o2, y3, o3 = feed_forward(X, network, intermediate_results=True)\n",
    "\n",
    "            print(cross_entropy_cost(o3,Y))\n",
    "\n",
    "            ## Layer 2to3 weights\n",
    "            # Step 1\n",
    "            W2_grad = np.einsum(\"ij,ik\",o3, o2)\n",
    "            # Step 2\n",
    "            for i, y in enumerate(Y):\n",
    "                W2_grad[y] -= o2[i]\n",
    "            # Update weights\n",
    "            network[2][0] -= W2_grad*learning_rate\n",
    "\n",
    "            ## Layer 3 biases\n",
    "            # Step 1\n",
    "            b2_grad = np.sum(o3, axis=0)\n",
    "            # Step 2\n",
    "            for i, y in enumerate(Y):\n",
    "                b2_grad[y] -= 1\n",
    "            # Update biases\n",
    "            network[2][1] -= b2_grad*learning_rate\n",
    "    \n",
    "    # Layer 1to2 weights\n",
    "    \n",
    "    # Layer 2 biases\n",
    "    \n",
    "    # Layer 0to1 weights\n",
    "    \n",
    "    # Layer 1 biases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "o3 = feed_forward(x_train, network)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "network = init_network()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "646084.1600977055\n",
      "657340.5835065747\n",
      "529296.891136106\n",
      "620118.5799105082\n",
      "630540.7409174069\n",
      "599274.257896711\n",
      "542696.8124306899\n",
      "585129.8965302057\n",
      "611929.7391193736\n",
      "628307.4207016429\n",
      "483141.6066769803\n",
      "517385.8499853641\n",
      "426564.1612109548\n",
      "571729.9752356218\n",
      "358820.11466610845\n",
      "393808.7980464137\n",
      "368497.8356010865\n",
      "416142.00020405534\n",
      "379664.4366799073\n",
      "313409.2702789038\n",
      "337231.3525803882\n",
      "355842.3543784229\n",
      "293309.3883370263\n",
      "267998.4258916991\n",
      "263531.78546017077\n",
      "256087.38474095688\n",
      "273953.9464670702\n",
      "256831.82481287827\n",
      "227201.2429156116\n",
      "215887.62085720192\n",
      "212165.42049759498\n",
      "213654.30064143776\n",
      "235243.06272715802\n",
      "221843.14143257303\n",
      "210737.57879067073\n",
      "189832.21833995334\n",
      "162287.93567886198\n",
      "167499.0161823117\n",
      "171965.65661384002\n",
      "198021.0591310886\n",
      "188343.33819611056\n",
      "192809.9786276389\n",
      "182387.81762073946\n",
      "187598.89812418917\n",
      "209187.66020990943\n",
      "228543.10207986552\n",
      "222587.58150449442\n",
      "221843.14143257303\n",
      "224820.90172025858\n",
      "210519.08587806852\n",
      "205465.4598503025\n",
      "207868.05049641547\n",
      "189832.21833995334\n",
      "169732.33639807586\n",
      "164521.25589462614\n",
      "145165.81402467005\n",
      "166754.5761103903\n",
      "162287.93567886198\n",
      "180154.4974049753\n",
      "175687.85697344696\n",
      "176432.29704536835\n",
      "181643.37754881807\n",
      "191321.09848379612\n",
      "213654.30064143776\n",
      "233009.74251139385\n",
      "239709.70315868634\n",
      "221098.70136065164\n",
      "241943.0233744505\n",
      "212909.86056951637\n",
      "226309.78186410136\n",
      "262042.905316328\n",
      "234498.62265523663\n",
      "235243.06272715802\n",
      "252365.18438134994\n",
      "217376.5010010447\n",
      "238965.26308676496\n",
      "224076.4616483372\n",
      "198021.0591310886\n",
      "147399.1342404342\n",
      "118365.97143550008\n",
      "114643.77107589314\n",
      "116132.65121973591\n",
      "111666.01078820758\n",
      "112410.45086012897\n",
      "102732.72992515092\n",
      "102732.72992515092\n",
      "101243.98578213555\n",
      "102732.72992515092\n",
      "102732.72992515092\n",
      "104221.6100689937\n",
      "101988.28985322954\n",
      "107943.81042860064\n",
      "110177.1306443648\n",
      "125065.93208279257\n",
      "151121.33460004115\n",
      "203043.90152427956\n",
      "217376.5010010447\n",
      "206209.89992222388\n",
      "230031.9822237083\n",
      "260554.0251724852\n",
      "253854.06452519272\n",
      "270231.74610746326\n",
      "250876.30423750717\n",
      "199509.93927493138\n",
      "199509.93927493138\n",
      "153354.65481580532\n",
      "130277.01258624229\n",
      "128043.69237047812\n",
      "129532.5725143209\n",
      "121343.73172318563\n",
      "130277.01258624229\n",
      "131021.45265816367\n",
      "141443.6136650631\n",
      "135488.093089692\n",
      "153354.65481580532\n",
      "151865.77467196254\n",
      "195043.29884340306\n",
      "190576.65841187473\n",
      "193554.41869956028\n",
      "186854.4580522678\n",
      "198765.49920301\n",
      "195043.29884340306\n",
      "206209.89992222388\n",
      "205465.4598503025\n",
      "236731.9428710008\n",
      "238965.26308676496\n",
      "227054.22193602275\n",
      "192065.5385557175\n",
      "192065.5385557175\n",
      "159504.57848654286\n",
      "145910.25409659144\n",
      "128043.69237047812\n",
      "116877.0912916573\n",
      "93055.00899017288\n",
      "92310.56891825149\n",
      "80399.52776750928\n",
      "83377.28805519483\n",
      "79655.08769558789\n",
      "84121.72812711622\n",
      "82632.84798327344\n",
      "84121.72812711622\n",
      "88588.36855864455\n",
      "87843.92848672316\n",
      "90821.68877440871\n",
      "90077.24870248733\n",
      "96032.76927785843\n",
      "99010.52956554398\n",
      "109432.69057244342\n",
      "125065.93208279257\n",
      "148888.014384277\n",
      "181643.37754881807\n",
      "200998.81941877416\n",
      "200998.81941877416\n",
      "203976.5797064597\n",
      "213654.30064143776\n",
      "203976.5797064597\n",
      "211420.9804256736\n",
      "206209.89992222388\n",
      "225565.34179217997\n",
      "221843.14143257303\n",
      "209932.10028183082\n",
      "211420.9804256736\n",
      "203232.13963453833\n",
      "203976.5797064597\n",
      "193554.41869956028\n",
      "185365.577908425\n",
      "160799.0555350192\n",
      "160054.6154630978\n",
      "143676.93388082727\n",
      "126554.81222663535\n",
      "116132.65121973591\n",
      "101988.28985322954\n",
      "104221.6100689937\n",
      "93799.44906209427\n",
      "93799.44906209427\n",
      "90077.24870248733\n",
      "87099.48841480177\n",
      "88588.36855864455\n",
      "80399.52776750928\n",
      "81143.96783943067\n"
     ]
    }
   ],
   "source": [
    "stochastic_gradient_descent(x_train, y_train, network)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "5//3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "y1, o1, y2, o2, y3, o3 = feed_forward(x_train[0:50], network, intermediate_results=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(50, 128)"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "o2.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "50"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(o3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10, 128)"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.einsum(\"ij,ik\",o3,o2).shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10, 128)"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "network[2][0].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "y1, o1, y2, o2, y3, o3 = feed_forward(x_train, network, intermediate_results=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.einsum(\"ij,ik\",o3[0:1],o2[0:1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = np.outer(o3[0],o2[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(a,b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}