NN-from-Scratch.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>Neural Networks from Scratch</center>\n",
    "<center> Chia-Hung Yuan & DataLab </center>\n",
    "<center> Fall 2020 </center>\n",
    "\n",
    "In this tutorial, you will learn the fundamentals of how you can build neural networks without the help of the deep learning frameworks, and instead by using NumPy.\n",
    "<img src=\"./figs/deep_nn-1.png\" width=\"800\"/>\n",
    "\n",
    "Creating complex neural networks with different architectures in Python should be a standard practice for any Machine Learning Engineer and Data Scientist. But a genuine understanding of how a neural network works is equally as valuable. This is what we aim to expand on in this article, the very fundamentals on how we can build neural networks, without the help of the frameworks that make it easy for us.\n",
    "\n",
    "## Model architecture\n",
    "We are building a basic deep neural network with 3 layers in total: 1 input layer, 1 hidden layers and 1 output layer. All layers will be fully connected. We implement ReLU and sigmoid activation functions. SGD and Momentum optimizer are available.\n",
    "\n",
    "Let's try to define the layers in an exact way. To be able to classify digits, we must end up with the probabilities of an image belonging to a certain class, after running the neural network, because then we can quantify how well our neural network performed.\n",
    "\n",
    "1. Input layer: In this layer, we input our dataset, consisting of 28x28 images. We flatten these images into one array with $28×28=784$ elements. This means our input layer will have 784 nodes.\n",
    "2. Hidden layer: In this layer, we have decided to reduce the number of nodes from 784 in the input layer to 64 nodes.\n",
    "3. Output layer: In this layer, we are reducing the 64 nodes to a total of 10 nodes, so that we can evaluate the nodes against the label. This label is received in the form of an array with 10 elements, where one of the elements is 1, while the rest is 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sklearn.datasets import fetch_openml\n",
    "import matplotlib.pyplot as plt\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Helper functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Helper function\n",
    "def show_images(image, num_row=2, num_col=5):\n",
    "    # plot images\n",
    "    image_size = int(np.sqrt(image.shape[-1]))\n",
    "    image = np.reshape(image, (image.shape[0], image_size, image_size))\n",
    "    fig, axes = plt.subplots(num_row, num_col, figsize=(1.5*num_col,2*num_row))\n",
    "    for i in range(num_row*num_col):\n",
    "        ax = axes[i//num_col, i%num_col]\n",
    "        ax.imshow(image[i], cmap='gray', vmin=0, vmax=1)\n",
    "        ax.axis('off')\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "    \n",
    "def one_hot(x, k, dtype=np.float32):\n",
    "    \"\"\"Create a one-hot encoding of x of size k.\"\"\"\n",
    "    return np.array(x[:, None] == np.arange(k), dtype)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dataset preparation\n",
    "In this tutorial, we will use [MNIST dataset](http://yann.lecun.com/exdb/mnist/) to demo how to build a neural network. MNIST contains 70,000 images of hand-written digits, 60,000 for training while 10,000 for testing, each $28×28$ pixels, in greyscale with pixel-values from 0 to 255.\n",
    "\n",
    "You might notice we not only import Numpy, but also Scikit-learn and Matplotlib. Matplotlib helps us visualize the training data, while Scikit-learn let us load the dataset more easily, but they are not used for any of the actual neural network. There are many ways to load dataset, you can choose any one you prefer.\n",
    "\n",
    "\n",
    "### Note\n",
    "To load dataset, `fetch_openml` requires `scikit-learn >= 0.22`. For older version, please change `from sklearn.datasets import fetch_openml` to `from sklearn.datasets import fetch_mldata`.\n",
    "\n",
    "Now we have to load the dataset and preprocess it, so that we can use it in NumPy. We do normalization by dividing all images by 255, and make it such that all images have values between 0 and 1, since this removes some of the numerical stability issues with activation functions later on. We choose to go with one-hot encoded labels, since we can more easily subtract these labels from the output of the neural network. We also choose to load our inputs as flattened arrays of 28 * 28 = 784 elements, since that is what the input layer requires."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Load data\n",
    "mnist_data = fetch_openml(\"mnist_784\")\n",
    "x = mnist_data[\"data\"]\n",
    "y = mnist_data[\"target\"]\n",
    "\n",
    "# Normalize\n",
    "x /= 255.0\n",
    "\n",
    "# One-hot encode labels\n",
    "num_labels = 10\n",
    "examples = y.shape[0]\n",
    "y_new = one_hot(y.astype('int32'), num_labels)\n",
    "\n",
    "# Split, reshape, shuffle\n",
    "train_size = 60000\n",
    "test_size = x.shape[0] - train_size\n",
    "x_train, x_test = x[:train_size], x[train_size:]\n",
    "y_train, y_test = y_new[:train_size], y_new[train_size:]\n",
    "shuffle_index = np.random.permutation(train_size)\n",
    "x_train, y_train = x_train[shuffle_index], y_train[shuffle_index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training data: (60000, 784) (60000, 10)\n",
      "Test data: (10000, 784) (10000, 10)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhQAAADqCAYAAAD6fdylAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAfM0lEQVR4nO3de7xPVf7H8XXcIwol9xBTGTOuTc1DyUTEKZSTkqQ6iUwpDdFFTeU2MkkKo+sYmTQRSZR6dJHbJBmVSTEI5RIVp4vr+f3xezw+fdbq7K/tu763fb6v51/vZe3v3iv7fI/VXnutlVNYWGgAAAB8lEh3AwAAQPTRoQAAAN7oUAAAAG90KAAAgDc6FAAAwBsdCgAA4K1UrMqcnBzmlKZRYWFhTpjjuE/pxX2KBu5TNHCfoqGo+8QTCgAA4I0OBQAA8EaHAgAAeKNDAQAAvNGhAAAA3uhQAAAAb3QoAACANzoUAADAGx0KAADgjQ4FAADwRocCAAB4o0MBAAC8xdwcDACOVYMGDazyU089Jblt27ZWXUFBgeSKFSsmtV2IbejQoZJXrVpl1S1atCjVzUEE8YQCAAB4o0MBAAC8Ze2QxznnnCN52bJlVt2RI0cklyxZMmVtyiZ16tQJddy+ffsC6yZNmiS5sLDQqnvhhRckz5079xhbh6OpV6+eVR4yZIjknj17WnUnnHCCZP3dMsaYUqV+/hXUqVMnyQsWLEhEM+EoW7as5HHjxll1vXv3lrxhw4bAz73zzjuSY30/kX14QgEAALzRoQAAAN7oUAAAAG857tizVZmTE1yZAPo9httuu82qW7FiheTx48cn9FrGGPP8889Ldsfz9Thv6dKlva8dr8LCwpwwxyX7PoXlTvu77LLLJHfu3Nmqy8vLkxzrZ/Dzzz+XPHPmTKvunnvuCfycPrZXr16BxyVC1O7Tsahbt67k7t27Sx41apR1XJkyZbyv9eGHH0pu3bq1Vbd//37v8xfn+6Tf9apZs6Zk9zty0UUXSa5du3Zc13r44Ycl33HHHVZdrO9yWMX5PhUnRd0nnlAAAABvdCgAAIC3tA556Gu708lKlPi5r3P55ZdbdS+++OIxX8s9hx7yyMmxn9zodqVz2mgUHv2VL19e8rRp06y6bt26BX5O/52HfUzq3qcPPvhAcosWLQI/17hxY8mfffZZqGsdiyjcp1j0d23ChAlWnZ5KGHYlyyVLlljlV199VfLhw4etujFjxhR5jubNm1vlNWvWhLp2LFG/T8cdd5xkd6ghNzdXcsuWLUOd74cffrDKGzdulOwO9f7qV78q8hzlypWzygcPHgx17Viidp9OPvlkq9y/f3/JU6dOterOPPNMyfn5+ZL1MJUxxqxbt06yHmIyxpj169fH39gEYsgDAAAkBR0KAADgLa0rZephDnfIQ0vEm8Oxrq0f+R6tLbDp2RR6pUPXgAEDrLJ+DB6LHsp46aWXAuti0Y+DkzHkEXV6KMmdDaWHOfbu3St5xowZ1nFz5syRrFdSNMaYAwcOSL7zzjsD27FlyxbJO3bsOFqziz09xGGM/fcadljDtXbtWskPPvigVadXl3VngGzevDmu62UDPXxujL0B3qBBg6y6ChUqSNarxLrDufoc7nfh/vvvD2yLvt6sWbMkf/HFF4GfSSSeUAAAAG90KAAAgDc6FAAAwFta36Fw310IqnPHlxJ9bff8sdqV7dwxwYsvvliy++7JzTffLNmdPhXW1q1bJZ944olWnV5B9brrrrPqCgoKijwOv6Sncvbt29eq0ytl6qmbmzZtCn3+Jk2aSNY/E66lS5dKztZ3KE4//XTJI0aMsOrCvjehx8v/+c9/WnV6CuLXX38deA49vu86dOhQqHYUJ+7yARdccIHkVq1aBX7u3XfftcqzZ8+WXL9+/cDPDR48WLL+/rhteeCBB6w6/Y6SfqetQ4cOgddKJP7lBAAA3uhQAAAAb0wbNUwb9aH/rtz7pDd7SgR3A6oLL7ww8NrxDrFku9WrV8csx0NPc6tevbpV991330m+8cYbva8VdXpVV725nuvbb7+1ynrK57hx4yRv2LAhrnbE2nhPr26aiJUxo+C8886zygsXLpTsDpk/99xzkt0hxJ9++inU9Ro2bCjZncr95z//WfKwYcMCz5GIDfuOFU8oAACANzoUAADAGx0KAADgLeXvUOjxIKaNFm96h8Lly5d7n0/vbGqMMbVq1Qo8NlN25MtG7rjuRRddFHisXv5ZT/XNVnqaoTvFVk/bdad86t0p46WngPfp0yfwOP2+RnFWqVIlybHeyXJ3173pppskh31nwrVq1SrJV111lVV31113BX5u27Ztkm+44Ya4ru2DfzkBAIA3OhQAAMBbyoc8brvtNsmpnDaqr+tej2mj4T355JNWWU9Rcz3++OOSd+/ebdXNnz//mK89adKkwDp3F1F3B0DEp3Tp0pJjDT3eeuutkt1HsuXKlZPsPh7WPyOwvyeTJ09O6rX0KqjGGPPII49Idu/12LFjJSdieCUKqlWrJllP43S5U0rj4e78mpeXJznW904PcRhjzDXXXCM5HcO+PKEAAADe6FAAAABvdCgAAIC3lL9DsWLFCsk9evSQ7L7HoMeN3GlKus59vyKozh2HilXHtNFg+/bts8p6eWC9k54x9jRPvVysMfb0OPec2r333is5NzfXqtP38Prrr7fq9JLOiE2P39aoUcOq69Kli+SqVavGdf5XXnlFsn7XwhhjDhw4ENc54S8/P98q690v9bRFY+xdLLNRvO/xNW3a1CrrqdFDhgyR7C47r6/nXlvv9uoukf7WW2/F1c5E4V9OAADgjQ4FAADwlhPrUU5OTk7it/lUDh8+LNmdqqmHHRJRF2tqaKw6PW0u1QoLC0MtEZrs+xSLHtb4+9//btVdeumlgZ9bunRpkce1adPGOm7atGlFXssYY95++23JevjMmF9OU02mKNynWD755BPJZ5xxRsLP37x5c8l6tcdUi/p9SoRrr71W8kMPPWTVffPNN5LbtWtn1W3ZsiWp7dIy5T7Vq1dPsv6OGGNPhXZ3dN24caPk888/36rTK2dWrFhRcqwhedejjz4qedCgQYHHJVtR94knFAAAwBsdCgAA4C2tQx76sffZZ5/tXlty2JkcseqOZZaHritZsmTwf0CSZcqjv7D0IzxjjHn66aclu8Mf+u9cr3LpzjI4/vjjJeufF/ecqRzicEXtPrl69eoluUqVKnGdY8yYMZL142BjjKldu7bkr776Kq7zJ0LU71NYepj297//vVW3cOFCye5Qr97U6plnnklS644uE+/TsmXLrPLvfvc73Q6rTg9rvPnmm1bde++9J1n/znI3H9P/Br3xxhtWXffu3SWnc0M9hjwAAEBS0KEAAADe6FAAAABvKV8pU9NT/dx3KPTqmOmcNorw3BUv+/TpI7ls2bJWnV71slGjRoHn1Cv2uStlxlphE+E999xzx/yZU0891SqPGjUqUc0p9mrVqiVZr05pjDGnnXaa5E6dOgWeY9asWVZZv5vSuHFjybF2LHV373311Vclu++O6Sn+2ehvf/ubVdbvBenp68YYM2HCBMn//e9/rbrvv/9ecqwdl7/++mvJU6ZMserS+d7E0fCEAgAAeKNDAQAAvKV12mgqXX755Vb5+eefl8y00eTTm3wZY28WFutn8Mknn5Tcr1+/hLcrEYrTfYqlTJkykvVjXWPszY3cjeBGjhwpOZ3Diem8T/3795d89913S65Zs2aiLxWafvxujDF79uyR/O2331p1Qfdt8eLFVnn16tWS3cf9y5cvD9Wu4vx9atKkiWS9aqz7b1BeXp5kd3grUzBtFAAAJAUdCgAA4I0OBQAA8JbWaaPpFHbaqLub2/jx45PbsAhzl95+8MEHJefn51t1+r2JWO9Q6M+tXLnSqnviiSfiaifio3eg1O9MuN5//32rnI3TsBs2bGiVJ06cKFn/vpk+fbp13Nq1awPPqZelv+WWW0K1Y9euXVb5wIEDgcfq72+dOnVCnb9p06aBdYcOHbLK7tTxbPT4449L1r/33N1MY00pzWQ8oQAAAN7oUAAAAG9ZO+ShHzu6U3Z0HUMcsZUvX16yu0Nht27dAj+nj9WP99wV6fTul/fcc49VV1yHPEqVsr+WeifPK664wqrr2bOnZL2TYaI0a9ZMsp5q7Vq/fr1kd+XAbHTJJZdYZf07RU8hjfUz7A6buN8NbdOmTUWe0/3MN998E3gOvfppixYtAo/TWrdubZX1z4ueFpmt9PfTGGNatWolWa/0606J1zuWRglPKAAAgDc6FAAAwBsdCgAA4C1rl96eMWOG5FjTRkuXLp3chsUQhSVo582bJznW7ogDBgywylOnTi3yuCVLllhlvQvttm3brDp3x8t0SfR9cqfi6R1Xt2/fbtU1b95c8s6dO8OcPqZKlSpZZf3eRMeOHQM/17JlS8l6+eVMksrvk94t0hhjKleuLFnvKLp3717ruKFDh0q+5pprrLpq1apJnjt3rlWnp7dv3rw5jhZnjij83ovlpJNOkuy+R3LKKadI1rv8uvc6Clh6GwAAJAUdCgAA4I1poyb2tFHYzj//fKucm5sref/+/VadnlKmH9sbY0+N1FPU6tWrZx2n7417n7JR1apVrXKPHj0ku9N2g+gVF40xpmvXrpI7dOhg1bVv316y3o1y+PDh1nEfffRRqGtnC/c+6WHUf/3rX5Ld4S29u7E7HP3UU09Jvv/++606dzgQ6aOnhrqrB2sbN25MRXNSin85AQCANzoUAADAW9YOebA5WHzcFfT0Y1m98psx9hvr7pDHueeeKznW7BB9/uK6MqZry5YtVlnPfHFXJpwwYUKROVH0UMbMmTMlT5kyJeHXKk4ee+wxq6xnOcVahVJvgDd27FirbtasWQlqHZJJ/z6rUKGCVffZZ59JjrXyaVTxhAIAAHijQwEAALzRoQAAAN6y9h0Kpo0mnt4Z1BhjunfvXmQ2xv471+9JbN261Tpu0qRJkt0x5eJKT880xp7W6e42etddd0muVatWXNf7+OOPJb/00ktWnf47/+GHH+I6fzb605/+ZJXDvt/y5ZdfSubvOxpOPPFEq3zllVdKdqf+jhw5UrK+18UF/3ICAABvdCgAAIC3rNkcrHbt2lZZbzr1wgsvWHVsDhbMXflt3LhxkvPz8wM/504bLSgokDx//nzJzz77rHXc7t2742lmSmXifcIvcZ+iIWr3SU+PN8aYUaNGSd60aZNV16ZNG8lRX92UzcEAAEBS0KEAAADe6FAAAABvWTNt1J2OqMt6hz/E5i6v3a9fvyIzAGSDRo0aBdZNmzbNKkf9vYmj4QkFAADwRocCAAB4y5ohDwAAEu3OO++0yvXq1ZO8Zs2aFLcmvXhCAQAAvNGhAAAA3uhQAAAAb1mz9HYURW0J2mzFfYoG7lM0cJ+igaW3AQBAUtChAAAA3mIOeQAAAITBEwoAAOCNDgUAAPBGhwIAAHijQwEAALzRoQAAAN7oUAAAAG90KAAAgDc6FAAAwBsdCgAA4I0OBQAA8EaHAgAAeKNDAQAAvNGhAAAA3uhQAAAAb3QoAACANzoUAADAGx0KAADgjQ4FAADwRocCAAB4o0MBAAC80aEAAADe6FAAAABvdCgAAIA3OhQAAMAbHQoAAOCNDgUAAPBGhwIAAHijQwEAALzRoQAAAN5KxarMyckpTFVD8EuFhYU5YY7jPqUX9ykauE/RwH2KhqLuE08oAACANzoUAADAGx0KAADgjQ4FAADwRocCAAB4o0MBAAC80aEAAADe6FAAAABvdCgAAIA3OhQAAMAbHQoAAOCNDgUAAPBGhwIAAHijQwEAALzRoQAAAN5KpbsBQCq1bdvWKg8dOlRyiRJ2//rBBx+U/N577yW1XbF07NhRcuvWrRN+/h07dkhet26dVbd48eJQ5zj99NMl5+XlBR43cOBAq1ypUiXJy5Ytk/yHP/zBOu7AgQOh2pEt8vPzrfKpp54q+e6775bs/kwfOXIk1Pnvueceq7x161bJ//jHP0K3E7aSJUtKHjNmjOTBgwdbxz311FOSb7jhhuQ3LEF4QgEAALzRoQAAAN7oUAAAAG/F7h0KPWZYt25dye6YbPXq1SVffvnlVl2zZs0k79mzx6rTY4tTpkzxayyOSU5OjlUuXbq05MaNG1t17dq1k9ypUyfJ5513XuA5XMcdd5zkNm3aHFtjE2ju3LmS3fYePHgwsC4K9LsR77zzTpF/nq1+/etfW+Vrr71Wsvsuih6bLywslOy+M6HrYtHvDxljv2ezceNGycuXL7eOO3ToUKjzZyv9b83tt98uedWqVdZxAwYMSFmbEoknFAAAwBsdCgAA4C0n1iOwnJyccM/H0qhGjRpW+b777pPct2/fwM8VFBRI1tPVjDFmyZIlkt1H3RdccIFkPVVu/fr1IVscXmFhYc7Rj4rGferTp49VPvnkkyW/8sorkj/99FPruNzcXMk9evSw6nr37u3dru3bt0vevXu3Vbd69epQ10r2fTrnnHMkN2/e3KrTj59nz54tuWzZstZx+uf98OHDVt0JJ5wQT7NC+/bbbyV/8MEHVt3o0aMlv/XWW0ltR9S+Ty+//LJV7ty5c+Cx+p7qe12mTBnrOD2Mlwi1atWyynpoJF5Ru0+xVKhQwSp/8sknkuvUqSP5N7/5jXXc2rVrk9uwBCjqPvGEAgAAeKNDAQAAvEVyyKN9+/aSp0+fbtVVq1ZN8vfffy/ZXfntxRdflLxt27bAa+lhDWOMWblypeSrr75asn4TP1Gi/uhv+PDhkm+66SarTs+y0W+GuytStmzZUnLFihVDX1s/Lt68ebPkOXPmWMfpoaotW7aEPr+WzvukZzXp1fWuueYa6zg9Q2D+/PlW3RVXXBF4fv1YVs+QmTFjRug26mO/++670J9LtCh8nxo2bCjZ/Z3i/i7Sxo0bJ3nYsGGSzz77bOs4/T3Uv79cr732mlXWq7Vq7s+B+3MXjyjcp7CqVKlilXft2lXkcXqWTlQw5AEAAJKCDgUAAPBGhwIAAHiLxEqZt9xyi1X+61//Ktkde9K7I1533XWS//e//8V1bXf3xcmTJ0tesGBBXOcsToKmfxpjv//g7nqolSr184+hnvJmjDGTJk2SvGLFCqvu/fffDzznV199JTnsDotRVLlyZcl6/NqdAqunRrurv+qfaaSWfmfCGGPmzZsnuVGjRoGfGzFihFXW02819zujp0K77yR17dpV8qhRo6w6PT20SZMmge1CduMJBQAA8EaHAgAAeMvYIQ89vWnChAlWnd4QSQ9/GGPMHXfccczXch/h6Y1xpk2b5n3+4qxDhw6SzzrrrLjOoR/z6g1zjDFmw4YN8TUsS3Tp0qXIPy9XrpxVPuWUUyQnY1VXhKeH//Lz8626WMMcesVRPU3UGGP2798f6tr6uFjT5d3fiQxzIAyeUAAAAG90KAAAgDc6FAAAwFvGvEPhTmEaMGCAZHd5cF2nlxuO93o33nijVafHpZcuXRrX+Ysrd5put27dAo/96aefJP/lL3+x6vSue//+978lf/HFF75NzCq//e1vi/xzd5dDvWx2vMtf62mpb775plUXazz+888/lxx2rL840797hgwZEvpzY8eOlexOr46He+2tW7dKTucS6cWJO2Vd//y7OwIXBzyhAAAA3uhQAAAAbxmz2+gjjzxilQcOHCjZfVx+3333ST5w4ECo89erV88q65Xl3N0W9S6l7m59a9euDXW9RMjEXfdat25tlRctWiTZnaqoH9HqHRCLm3Tep/Hjx0vW35lM8sYbb0hesmSJVffss89KTvZwV6Z8n/RuyQsXLgz9uR9//FFyzZo1rbp9+/b5N0zJzc21ykG7Kbvt2Llzp/e1M+U+JcPEiRMl66F7Pf3emF8OKWYidhsFAABJQYcCAAB4o0MBAAC8Zcy00eOPPz6w7tRTT7XKnTp1krx9+3arrm7dupL/+Mc/StY7XxpjTPny5QOvp3exTOU7E1FwzjnnWGX3vQlNLxWM7KXfGdDZGHtH4Ouvv16yu+R61KcT62m8t956a6jPDB061CovX75cst4Z1Bhjpk+f7tG6X3KXwA+ip4bj6FatWlXkn+fl5VnlKLxDURSeUAAAAG90KAAAgLeMGfJwp5NdffXVkq+88kqrzi1rehrs5s2bJU+ePNk6bseOHZLdnfv0yo2wPf/881ZZrzLq7pQ4YsQIyZdeeqlV9/DDD0ueOXNmIpuYVfTQwJdffim5atWq1nErVqyI6/wzZsyQvHHjxlCfcVee1UOU7lCjHqLU00v1dY0xZvjw4ZI3bdoUqh2ZpFKlSpJPOumkUJ9xh1v170j392WiNWjQILBuzpw5khnyODaLFy8u8s+rV6+e4pYkB08oAACANzoUAADAW8YMeTzzzDNWWb8N27Nnz9Dn0Zt5vfzyy5LdR63/+c9/JMdaLTRbXXzxxZL1G/Zr1qyxjmvXrp1kdzW9pk2bSj7rrLOsumnTpkmeMGGC5KlTpwa2yV1Ndc+ePYHHZovHHntMsl51sUqVKtZxqRzG00MXxhjTqlUryfPmzbPqqlWrVuQ5rrrqKqt87rnnSq5fv75vE9OqRImf/z8uJyfUopAp57ZLl/UsuLArFeP/6U309N+p/o4YY0yZMmUkR+nvmCcUAADAGx0KAADgjQ4FAADwljG7jSabu0vmu+++K/mHH36w6vR41rp165LbsBjSueve559/Llm/f9K5c2frOP0uiqtLly6S3d0L9dTfihUrhmrTxx9/bJXvv/9+ybNmzQp1jmQozrsjJtrTTz9tlfXPQdmyZUOdo2TJknFdO533Sb8D8vbbb4f6zCWXXGKVFyxYkMgm/UKdOnUku9OMDx8+LLlFixaSd+3alfB2FOfv03HHHSe5oKBAsvvOiv6dm6lTc9ltFAAAJAUdCgAA4C1jpo0m27XXXhtY99JLL1nldA5zZIqDBw9KrlGjhuTevXtbx8Ua8tDTdnU2xpiHHnpIsl7VsW/fvtZxetikSZMmVt3o0aMlu4+Rd+/eHdgupI/eAMwYe+hxypQpkkuXLp2yNuH/9erVS7I7nVdPHU/GMAeKB55QAAAAb3QoAACANzoUAADAW7F+h6J27dqSr7jiCqtOT9N56623UtamqNDLaJ9xxhmSb775Zus4vcyyXpLXmNjvMaxfv77IrHeINcaePuXuMqvHcr/55pvAayFz7d+/X7JekhrJV65cOavcpk2bNLUke+gpoK+//rrkjh07pqM5Ccc3GAAAeKNDAQAAvBW7IY9SpX7+Txo8eLDkChUqWMctXrxYsrvTKYwZNWqU5OrVq0t2h45effVVyV999ZVVN3nyZMl6J0xjjFm5cmWR192+fbtV1qsDurvO6hU23fu7b9++Is+PzHLiiSdKjrUC5kcffZSK5iSN3hlX7zhZq1atdDTHGGNP3TbGmA4dOkjesWOHVdevX7+UtKm40ytT6x2A3SEPvdvziy++mPyGJQhPKAAAgDc6FAAAwBsdCgAA4K3YvUPRoEEDybfcckvgcVOnTk1FcyJLv4Nw3XXXSd66dat13GWXXSb5zDPPtOoeeOAByXfddZdVt2HDhlDt0GPs7s64enlwnZG53Km/Q4cODfU5PY05itauXStZ/7cMGDAg8DNVqlTxvq67k++4ceMkd+rUKfBzekqjMcYsWrTIuy2wudPsNf3eWpTwhAIAAHijQwEAALzluI+RrcqcnODKDPXCCy9I7t69u2Q9RccYY1q3bi35yJEjyW9YHAoLC3OOflR671P79u0ld+3a1aq76qqrJFeuXNn7WnrHQ2OMefLJJyWPGDHC+/zxisJ90vT0QGOMOe200xJ6/saNG1tl/T086aSTrLpYU0U1/bl4V0XNlPt0wgknSI61muzevXutst6NVQ8nGmOvwKjpadfGGHPhhReGamP9+vWt8pYtW0J9LhEy5T6lkvtv0Keffiq5RYsWVl3QvU61ou4TTygAAIA3OhQAAMBb5Ic8ypYta5X14yD93+ausjhz5szkNiwBov7or2HDhpL1BmPG2I/Ba9SoIdl9m3zdunWSX3vtNasuU2Z2pPM+tW3bVnKzZs0k65X23Dp3VdEyZcokulneHn30Uat8++23S471OyuWTPk+6Y0J7733Xqtu+PDhoc7hrgQb9Heih1fc45YsWWLVPfHEE5JnzJhh1aVyWDhT7lMq6d9zxti/O91hwkzZCJEhDwAAkBR0KAAAgDc6FAAAwFvk36GYNGmSVe7fv7/kzZs3S3ansv3444/JbVgCZONYYhSl8z6NHz9e8sCBAxN9+qT78MMPJY8cOVLy/PnzreMOHDjgfa1M/D7Vrl3bKufm5koeM2aMVeeuehmGfl/DGGMKCgoku6uWulNM0yUT71Oy6Z99Y4wZNmyYZN6hAAAAWYUOBQAA8Bb5IY+dO3daZf14aPbs2ZLz8vJS1qZEycZHf1GUzvukV5osUSL4/w+6desmuUmTJoluRmgTJ060yvrx7eHDh5N67ah9n9zNu/Rmh+5qp0EWL15slR955BHJmbrhWtTuUyLoqfPGGLN06VLJ7vf1+++/T0mbjoYhDwAAkBR0KAAAgDc6FAAAwFvk36EYNGiQVR49erTkLl26SH799ddT1qZEycaxxCjiPkUD9ykauE/RwDsUAAAgKehQAAAAb5Ef8ijOePQXDdynaOA+RQP3KRoY8gAAAElBhwIAAHijQwEAALzRoQAAAN7oUAAAAG90KAAAgLeY00YBAADC4AkFAADwRocCAAB4o0MBAAC80aEAAADe6FAAAABvdCgAAIC3/wMAzmeZcE4upwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 540x288 with 10 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(\"Training data: {} {}\".format(x_train.shape, y_train.shape))\n",
    "print(\"Test data: {} {}\".format(x_test.shape, y_test.shape))\n",
    "show_images(x_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initialization\n",
    "The specific problem that arises, when trying to implement the feedforward neural network, is that we are trying to transform from 784 nodes all the way down to 10 nodes. When instantiating the `DeepNeuralNetwork` class, we pass in an array of sizes that defines the number of activations for each layer.\n",
    "\n",
    "```python\n",
    "dnn = DeepNeuralNetwork(sizes=[784, 64, 10])\n",
    "```\n",
    "\n",
    "This initializes the `DeepNeuralNetwork` class by the init function.\n",
    "\n",
    "```python\n",
    "def __init__(self, sizes, activation='sigmoid'):\n",
    "    self.sizes = sizes\n",
    "\n",
    "    # Choose activation function\n",
    "    if activation == 'relu':\n",
    "        self.activation = self.relu\n",
    "    elif activation == 'sigmoid':\n",
    "        self.activation = self.sigmoid\n",
    "\n",
    "    # Save all weights\n",
    "    self.params = self.initialize()\n",
    "    # Save all intermediate values, i.e. activations\n",
    "    self.cache = {}\n",
    "```\n",
    "\n",
    "To smarten up our initialization, we shrink the variance of the weights in each layer. Following [this nice video](https://www.coursera.org/lecture/deep-neural-network/weight-initialization-for-deep-networks-RwqYe) by Andrew Ng, we’ll set the variance for each layer to 1/n, where n is the number of inputs feeding into that layer. We use `np.random.randn()` function to get our initial weights, which draws from the standard normal distribution. So to adjust the variance to 1/n, we just divide by √n. \n",
    "\n",
    "The initialization of weights in the neural network is kind of hard to think about. To really understand how and why the following approach works, you need a grasp of linear algebra, specifically dimensionality when using the dot product operation, which is beyond the scope of this class.\n",
    "\n",
    "```python\n",
    "def initialization(self):\n",
    "    # Number of nodes in each layer\n",
    "    input_layer=self.sizes[0]\n",
    "    hidden_1=self.sizes[1]\n",
    "    hidden_2=self.sizes[2]\n",
    "    output_layer=self.sizes[3]\n",
    "\n",
    "    params = {\n",
    "        'W1':np.random.randn(hidden_1, input_layer) * np.sqrt(1./hidden_1),\n",
    "        'W2':np.random.randn(hidden_2, hidden_1) * np.sqrt(1./hidden_2),\n",
    "        'W3':np.random.randn(output_layer, hidden_2) * np.sqrt(1./output_layer)\n",
    "    }\n",
    "\n",
    "    return params\n",
    "```\n",
    "\n",
    "## Feedforward\n",
    "The forward pass consists of the dot operation in NumPy, which turns out to be just matrix multiplication. As described in the [introduction to neural networks](https://mlfromscratch.com/neural-networks-explained/#/) article, we have to multiply the weights by the activations of the previous layer. Then we have to apply the activation function to the outcome. \n",
    "\n",
    "To get through each layer, we sequentially apply the dot operation, followed by the sigmoid/relu activation function. In the last layer we use the softmax activation function, since we wish to have probabilities of each class, so that we can measure how well our current forward pass performs.\n",
    "\n",
    "<img src=\"./figs/backprop_algo_forward.png\" width=\"300\"/>\n",
    "\n",
    "```python\n",
    "def feed_forward(self, x):\n",
    "    self.cache[\"X\"] = x\n",
    "    self.cache[\"Z1\"] = np.matmul(self.params[\"W1\"], self.cache[\"X\"].T) +\\\n",
    "                                                        self.params[\"b1\"]\n",
    "    self.cache[\"A1\"] = self.activation(self.cache[\"Z1\"])\n",
    "    self.cache[\"Z2\"] = np.matmul(self.params[\"W2\"], self.cache[\"A1\"]) +\\\n",
    "                                                        self.params[\"b2\"]\n",
    "    self.cache[\"A2\"] = self.softmax(self.cache[\"Z2\"])\n",
    "    return self.cache[\"A2\"]\n",
    "```\n",
    "\n",
    "The following are the activation functions used for this article. As you can see, we provide a derivative version of the relu and sigmoid, since we will need that later on when backpropagating through the neural network.\n",
    "\n",
    "```python\n",
    "def relu(self, x, derivative=False):\n",
    "    if derivative:\n",
    "        x = np.where(x < 0, 0, x)\n",
    "        x = np.where(x >= 0, 1, x)\n",
    "        return x\n",
    "    return np.maximum(0, x)\n",
    "\n",
    "def sigmoid(self, x, derivative=False):\n",
    "    if derivative:\n",
    "        return (np.exp(-x))/((np.exp(-x)+1)**2)\n",
    "    return 1/(1 + np.exp(-x))\n",
    "\n",
    "def softmax(self, x):\n",
    "    # Numerically stable with large exponentials\n",
    "    exps = np.exp(x - x.max())\n",
    "    return exps / np.sum(exps, axis=0)\n",
    "```\n",
    "\n",
    "Note: A numerical stable version of the softmax function was chosen, you can read more from [Stanford CS231n](https://cs231n.github.io/linear-classify/#softmax) course.\n",
    "\n",
    "## Backpropagation\n",
    "The backward pass is hard to get right, because there are so many sizes and operations that have to align, for all the operations to be successful. Here is the full function for the backward pass; we will go through each weight update below.\n",
    "\n",
    "<img src=\"./figs/backprop_algo_backward.png\" width=\"400\"/>\n",
    "\n",
    "For people who are interested in the magic of backpropagation, please refer to [this nice article](https://mlfromscratch.com/neural-networks-explained/). For people who want to deep dive into methmetics and understand thoroughly, please refer to [NTHU CS565600 Deep Learning](https://nthu-datalab.github.io/ml/index.html). Here is [slide](https://nthu-datalab.github.io/ml/slides/10_NN_Design.pdf) and [video](https://www.youtube.com/watch?v=uYRUbvyKXAo&list=PLlPcwHqLqJDk3A0qFgFUDlyALzaF44NTL&index=4&ab_channel=Shan-HungWu) for backpropagation.\n",
    "\n",
    "```python\n",
    "def back_propagate(self, y, output):\n",
    "    current_batch_size = y.shape[0]\n",
    "\n",
    "    dZ2 = output - y.T\n",
    "    dW2 = (1./current_batch_size) * np.matmul(dZ2, self.cache[\"A1\"].T)\n",
    "    db2 = (1./current_batch_size) * np.sum(dZ2, axis=1, keepdims=True)\n",
    "\n",
    "    dA1 = np.matmul(self.params[\"W2\"].T, dZ2)\n",
    "    dZ1 = dA1 * self.activation(self.cache[\"Z1\"], derivative=True)\n",
    "    dW1 = (1./current_batch_size) * np.matmul(dZ1, self.cache[\"X\"])\n",
    "    db1 = (1./current_batch_size) * np.sum(dZ1, axis=1, keepdims=True)\n",
    "\n",
    "    self.grads = {\"W1\": dW1, \"b1\": db1, \"W2\": dW2, \"b2\": db2}\n",
    "    return self.grads\n",
    "```\n",
    "\n",
    "## Training (Stochastic Gradient Descent)\n",
    "We have defined a forward and backward pass, but how can we start using them? We have to make a training loop and choose to use Stochastic Gradient Descent (SGD) as the optimizer to update the parameters of the neural network.\n",
    "\n",
    "There are two main loops in the training function. One loop for the number of epochs, which is the number of times we run through the whole dataset, and a second loop for running through each batch one by one.\n",
    "\n",
    "For each batch, we do a forward pass by calling `self.feedforward()` with `x`, which is one batch in an array with the length 784, as explained earlier. The `output` of the forward pass is used along with `y`, which are the one-hot encoded labels (the ground truth), in the backward pass. `self.back_propagate()` returns the gradients of weights.\n",
    "\n",
    "```python\n",
    "def train(self, x_train, y_train, x_test, y_test):\n",
    "    for i in range(self.epochs):\n",
    "        # Shuffle\n",
    "        permutation = np.random.permutation(x_train.shape[0])\n",
    "        x_train_shuffled = x_train[permutation]\n",
    "        y_train_shuffled = y_train[permutation]\n",
    "\n",
    "        for j in range(num_batches):\n",
    "            # Batch\n",
    "            begin = j * self.batch_size\n",
    "            end = min(begin + self.batch_size, x_train.shape[0]-1)\n",
    "            x = x_train_shuffled[begin:end]\n",
    "            y = y_train_shuffled[begin:end]\n",
    "\n",
    "            # Forward\n",
    "            output = self.feed_forward(x)\n",
    "            # Backprop\n",
    "            grad = self.back_propagate(y, output)\n",
    "            # Optimize\n",
    "            self.optimize(l_rate=l_rate, beta=beta)\n",
    "```\n",
    "\n",
    "The `optimize()` function has the code for both SGD and momentum update rules. SGD algorithm is relatively straightforward, updating the networks by calculated gradient directly, while momentum make the same movement in the last iteration, corrected by negative gradient.\n",
    "\n",
    "```python\n",
    "def optimize(self, l_rate=0.1, beta=.9):\n",
    "    '''\n",
    "        Stochatic Gradient Descent (SGD):\n",
    "        θ^(t+1) <- θ^t - η∇L(y, ŷ)\n",
    "\n",
    "        Momentum:\n",
    "        v^(t+1) <- βv^t + (1-β)∇L(y, ŷ)^t\n",
    "        θ^(t+1) <- θ^t - ηv^(t+1)\n",
    "    '''\n",
    "    if self.optimizer == \"sgd\":\n",
    "        for key in self.params:\n",
    "            self.params[key] = self.params[key] -\\\n",
    "                                        l_rate*self.grads[key]\n",
    "    elif self.optimizer == \"momentum\":\n",
    "        for key in self.params:\n",
    "            self.momemtum_opt[key] = (beta*self.momemtum_opt[key] +\\\n",
    "                                      (1.-beta)*self.grads[key])\n",
    "            self.params[key] = self.params[key] -\\\n",
    "                                        l_rate * self.momemtum_opt[key]\n",
    "\n",
    "```\n",
    "\n",
    "After having updated the parameters of the neural network, we can measure the accuracy on a validation set that we conveniently prepared earlier, to validate how well our network performs after each iteration over the whole dataset.\n",
    "\n",
    "```python\n",
    "def accuracy(self, y, output):\n",
    "    return np.mean(np.argmax(y, axis=-1) == np.argmax(output.T, axis=-1))\n",
    "```\n",
    "\n",
    "Finally, we can call the training function, after knowing what will happen. We use the training and validation data as input to the training function, and then we wait.\n",
    "\n",
    "```python\n",
    "dnn.train(x_train, y_train, x_val, y_val)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3-layer neural network with Numpy\n",
    "Here is the full code, for an easy copy-paste and overview of what's happening."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "class DeepNeuralNetwork():\n",
    "    def __init__(self, sizes, activation='sigmoid'):\n",
    "        self.sizes = sizes\n",
    "        \n",
    "        # Choose activation function\n",
    "        if activation == 'relu':\n",
    "            self.activation = self.relu\n",
    "        elif activation == 'sigmoid':\n",
    "            self.activation = self.sigmoid\n",
    "        else:\n",
    "            raise ValueError(\"Activation function is currently not support, please use 'relu' or 'sigmoid' instead.\")\n",
    "        \n",
    "        # Save all weights\n",
    "        self.params = self.initialize()\n",
    "        # Save all intermediate values, i.e. activations\n",
    "        self.cache = {}\n",
    "        \n",
    "    def relu(self, x, derivative=False):\n",
    "        '''\n",
    "            Derivative of ReLU is a bit more complicated since it is not differentiable at x = 0\n",
    "        \n",
    "            Forward path:\n",
    "            relu(x) = max(0, x)\n",
    "            In other word,\n",
    "            relu(x) = 0, if x < 0\n",
    "                    = x, if x >= 0\n",
    "\n",
    "            Backward path:\n",
    "            ∇relu(x) = 0, if x < 0\n",
    "                     = 1, if x >=0\n",
    "        '''\n",
    "        if derivative:\n",
    "            x = np.where(x < 0, 0, x)\n",
    "            x = np.where(x >= 0, 1, x)\n",
    "            return x\n",
    "        return np.maximum(0, x)\n",
    "\n",
    "    def sigmoid(self, x, derivative=False):\n",
    "        '''\n",
    "            Forward path:\n",
    "            σ(x) = 1 / 1+exp(-z)\n",
    "            \n",
    "            Backward path:\n",
    "            ∇σ(x) = exp(-z) / (1+exp(-z))^2\n",
    "        '''\n",
    "        if derivative:\n",
    "            return (np.exp(-x))/((np.exp(-x)+1)**2)\n",
    "        return 1/(1 + np.exp(-x))\n",
    "\n",
    "    def softmax(self, x):\n",
    "        '''\n",
    "            softmax(x) = exp(x) / ∑exp(x)\n",
    "        '''\n",
    "        # Numerically stable with large exponentials\n",
    "        exps = np.exp(x - x.max())\n",
    "        return exps / np.sum(exps, axis=0)\n",
    "\n",
    "    def initialize(self):\n",
    "        # number of nodes in each layer\n",
    "        input_layer=self.sizes[0]\n",
    "        hidden_layer=self.sizes[1]\n",
    "        output_layer=self.sizes[2]\n",
    "        \n",
    "        params = {\n",
    "            \"W1\": np.random.randn(hidden_layer, input_layer) * np.sqrt(1./input_layer),\n",
    "            \"b1\": np.zeros((hidden_layer, 1)) * np.sqrt(1./input_layer),\n",
    "            \"W2\": np.random.randn(output_layer, hidden_layer) * np.sqrt(1./hidden_layer),\n",
    "            \"b2\": np.zeros((output_layer, 1)) * np.sqrt(1./hidden_layer)\n",
    "        }\n",
    "        return params\n",
    "    \n",
    "    def initialize_momemtum_optimizer(self):\n",
    "        momemtum_opt = {\n",
    "            \"W1\": np.zeros(self.params[\"W1\"].shape),\n",
    "            \"b1\": np.zeros(self.params[\"b1\"].shape),\n",
    "            \"W2\": np.zeros(self.params[\"W2\"].shape),\n",
    "            \"b2\": np.zeros(self.params[\"b2\"].shape),\n",
    "        }\n",
    "        return momemtum_opt\n",
    "\n",
    "    def feed_forward(self, x):\n",
    "        '''\n",
    "            y = σ(wX + b)\n",
    "        '''\n",
    "        self.cache[\"X\"] = x\n",
    "        self.cache[\"Z1\"] = np.matmul(self.params[\"W1\"], self.cache[\"X\"].T) + self.params[\"b1\"]\n",
    "        self.cache[\"A1\"] = self.activation(self.cache[\"Z1\"])\n",
    "        self.cache[\"Z2\"] = np.matmul(self.params[\"W2\"], self.cache[\"A1\"]) + self.params[\"b2\"]\n",
    "        self.cache[\"A2\"] = self.softmax(self.cache[\"Z2\"])\n",
    "        return self.cache[\"A2\"]\n",
    "    \n",
    "    def back_propagate(self, y, output):\n",
    "        '''\n",
    "            This is the backpropagation algorithm, for calculating the updates\n",
    "            of the neural network's parameters.\n",
    "\n",
    "            Note: There is a stability issue that causes warnings. This is \n",
    "                  caused  by the dot and multiply operations on the huge arrays.\n",
    "                  \n",
    "                  RuntimeWarning: invalid value encountered in true_divide\n",
    "                  RuntimeWarning: overflow encountered in exp\n",
    "                  RuntimeWarning: overflow encountered in square\n",
    "        '''\n",
    "        current_batch_size = y.shape[0]\n",
    "        \n",
    "        dZ2 = output - y.T\n",
    "        dW2 = (1./current_batch_size) * np.matmul(dZ2, self.cache[\"A1\"].T)\n",
    "        db2 = (1./current_batch_size) * np.sum(dZ2, axis=1, keepdims=True)\n",
    "\n",
    "        dA1 = np.matmul(self.params[\"W2\"].T, dZ2)\n",
    "        dZ1 = dA1 * self.activation(self.cache[\"Z1\"], derivative=True)\n",
    "        dW1 = (1./current_batch_size) * np.matmul(dZ1, self.cache[\"X\"])\n",
    "        db1 = (1./current_batch_size) * np.sum(dZ1, axis=1, keepdims=True)\n",
    "\n",
    "        self.grads = {\"W1\": dW1, \"b1\": db1, \"W2\": dW2, \"b2\": db2}\n",
    "        return self.grads\n",
    "    \n",
    "    def cross_entropy_loss(self, y, output):\n",
    "        '''\n",
    "            L(y, ŷ) = −∑ylog(ŷ).\n",
    "        '''\n",
    "        l_sum = np.sum(np.multiply(y.T, np.log(output)))\n",
    "        m = y.shape[0]\n",
    "        l = -(1./m) * l_sum\n",
    "        return l\n",
    "                \n",
    "    def optimize(self, l_rate=0.1, beta=.9):\n",
    "        '''\n",
    "            Stochatic Gradient Descent (SGD):\n",
    "            θ^(t+1) <- θ^t - η∇L(y, ŷ)\n",
    "            \n",
    "            Momentum:\n",
    "            v^(t+1) <- βv^t + (1-β)∇L(y, ŷ)^t\n",
    "            θ^(t+1) <- θ^t - ηv^(t+1)\n",
    "        '''\n",
    "        if self.optimizer == \"sgd\":\n",
    "            for key in self.params:\n",
    "                self.params[key] = self.params[key] - l_rate * self.grads[key]\n",
    "        elif self.optimizer == \"momentum\":\n",
    "            for key in self.params:\n",
    "                self.momemtum_opt[key] = (beta * self.momemtum_opt[key] + (1. - beta) * self.grads[key])\n",
    "                self.params[key] = self.params[key] - l_rate * self.momemtum_opt[key]\n",
    "        else:\n",
    "            raise ValueError(\"Optimizer is currently not support, please use 'sgd' or 'momentum' instead.\")\n",
    "\n",
    "    def accuracy(self, y, output):\n",
    "        return np.mean(np.argmax(y, axis=-1) == np.argmax(output.T, axis=-1))\n",
    "\n",
    "    def train(self, x_train, y_train, x_test, y_test, epochs=10, \n",
    "              batch_size=64, optimizer='momentum', l_rate=0.1, beta=.9):\n",
    "        # Hyperparameters\n",
    "        self.epochs = epochs\n",
    "        self.batch_size = batch_size\n",
    "        num_batches = -(-x_train.shape[0] // self.batch_size)\n",
    "        \n",
    "        # Initialize optimizer\n",
    "        self.optimizer = optimizer\n",
    "        if self.optimizer == 'momentum':\n",
    "            self.momemtum_opt = self.initialize_momemtum_optimizer()\n",
    "        \n",
    "        start_time = time.time()\n",
    "        template = \"Epoch {}: {:.2f}s, train acc={:.2f}, train loss={:.2f}, test acc={:.2f}, test loss={:.2f}\"\n",
    "        \n",
    "        # Train\n",
    "        for i in range(self.epochs):\n",
    "            # Shuffle\n",
    "            permutation = np.random.permutation(x_train.shape[0])\n",
    "            x_train_shuffled = x_train[permutation]\n",
    "            y_train_shuffled = y_train[permutation]\n",
    "\n",
    "            for j in range(num_batches):\n",
    "                # Batch\n",
    "                begin = j * self.batch_size\n",
    "                end = min(begin + self.batch_size, x_train.shape[0]-1)\n",
    "                x = x_train_shuffled[begin:end]\n",
    "                y = y_train_shuffled[begin:end]\n",
    "                \n",
    "                # Forward\n",
    "                output = self.feed_forward(x)\n",
    "                # Backprop\n",
    "                grad = self.back_propagate(y, output)\n",
    "                # Optimize\n",
    "                self.optimize(l_rate=l_rate, beta=beta)\n",
    "\n",
    "            # Evaluate performance\n",
    "            # Training data\n",
    "            output = self.feed_forward(x_train)\n",
    "            train_acc = self.accuracy(y_train, output)\n",
    "            train_loss = self.cross_entropy_loss(y_train, output)\n",
    "            # Test data\n",
    "            output = self.feed_forward(x_test)\n",
    "            test_acc = self.accuracy(y_test, output)\n",
    "            test_loss = self.cross_entropy_loss(y_test, output)\n",
    "            print(template.format(i+1, time.time()-start_time, train_acc, train_loss, test_acc, test_loss))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Results\n",
    "The results completely dependent on how the weights are initialized and the activation function we use. Experimentally, due to non-bounded behavior of `relu()`, the learning rate should be set much smaller than the one for `sigmoid()` (bounded). In addition, training with SGD optimizer with momentum should have better result since it avoids from getting stuck in local minima or saddle points.\n",
    "\n",
    "The reason behind this phenomenon is complicated and beyond the scope of this class. In short, the training results will be more stable and consistent as the batch size increases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 1: 0.90s, train acc=0.95, train loss=0.16, test acc=0.95, test loss=0.17\n",
      "Epoch 2: 1.75s, train acc=0.97, train loss=0.10, test acc=0.96, test loss=0.12\n",
      "Epoch 3: 2.58s, train acc=0.98, train loss=0.08, test acc=0.97, test loss=0.10\n",
      "Epoch 4: 3.42s, train acc=0.98, train loss=0.07, test acc=0.97, test loss=0.09\n",
      "Epoch 5: 4.27s, train acc=0.98, train loss=0.05, test acc=0.97, test loss=0.08\n",
      "Epoch 6: 5.13s, train acc=0.99, train loss=0.04, test acc=0.98, test loss=0.08\n",
      "Epoch 7: 5.97s, train acc=0.99, train loss=0.04, test acc=0.97, test loss=0.08\n",
      "Epoch 8: 6.82s, train acc=0.99, train loss=0.03, test acc=0.97, test loss=0.08\n",
      "Epoch 9: 7.69s, train acc=0.99, train loss=0.03, test acc=0.98, test loss=0.08\n",
      "Epoch 10: 8.55s, train acc=0.99, train loss=0.02, test acc=0.98, test loss=0.08\n"
     ]
    }
   ],
   "source": [
    "# Sigmoid + Momentum\n",
    "dnn = DeepNeuralNetwork(sizes=[784, 64, 10], activation='sigmoid')\n",
    "dnn.train(x_train, y_train, x_test, y_test, batch_size=128, optimizer='momentum', l_rate=4, beta=.9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 1: 0.70s, train acc=0.89, train loss=0.41, test acc=0.89, test loss=0.39\n",
      "Epoch 2: 1.27s, train acc=0.90, train loss=0.34, test acc=0.91, test loss=0.32\n",
      "Epoch 3: 1.82s, train acc=0.91, train loss=0.31, test acc=0.91, test loss=0.30\n",
      "Epoch 4: 2.34s, train acc=0.92, train loss=0.29, test acc=0.92, test loss=0.28\n",
      "Epoch 5: 2.88s, train acc=0.92, train loss=0.28, test acc=0.92, test loss=0.27\n",
      "Epoch 6: 3.42s, train acc=0.92, train loss=0.27, test acc=0.92, test loss=0.27\n",
      "Epoch 7: 3.94s, train acc=0.92, train loss=0.27, test acc=0.92, test loss=0.26\n",
      "Epoch 8: 4.48s, train acc=0.93, train loss=0.26, test acc=0.92, test loss=0.26\n",
      "Epoch 9: 5.00s, train acc=0.93, train loss=0.25, test acc=0.93, test loss=0.25\n",
      "Epoch 10: 5.56s, train acc=0.93, train loss=0.25, test acc=0.93, test loss=0.25\n"
     ]
    }
   ],
   "source": [
    "# ReLU + SGD\n",
    "dnn = DeepNeuralNetwork(sizes=[784, 64, 10], activation='relu')\n",
    "dnn.train(x_train, y_train, x_test, y_test, batch_size=128, optimizer='sgd', l_rate=0.05)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Good exercises in NumPy\n",
    "You might have noticed that the code is very readable, but takes up a lot of space and could be optimized to run in loops. Here is a chance to optimize and improve the code. For example, you can optimize the forward and backward pass, such that they run in a for loop in each function. This makes the code easier to modify and possibly easier to maintain. \n",
    "\n",
    "More challenging exercises including implement any other activation function from this overview of activation functions, and remember to implement the derivatives as well. Different optimizers, e.g. Adam, RMSProp, etc, are also worth to try.\n",
    "\n",
    "## Reference \n",
    "1. [CS565600 Deep Learning](https://nthu-datalab.github.io/ml/index.html), National Tsing Hua University\n",
    "2. [Building a Neural Network from Scratch: Part 1](https://jonathanweisberg.org/post/A%20Neural%20Network%20from%20Scratch%20-%20Part%201/)\n",
    "3. [Building a Neural Network from Scratch: Part 2](https://jonathanweisberg.org/post/A%20Neural%20Network%20from%20Scratch%20-%20Part%202/)\n",
    "4. [Neural networks from scratch](https://developer.ibm.com/technologies/artificial-intelligence/articles/neural-networks-from-scratch), IBM Developer\n",
    "5. [The Softmax Function Derivative (Part 1)](https://aimatters.wordpress.com/2019/06/17/the-softmax-function-derivative/)"
   ]
  }
 ],
 "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}