#1 MLP for MNIST.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Assignment 1, programming part\n",
    "#### IFT 6135\n",
    "#### February 17, 2019 23:59\n",
    "#### Authors: Leo Boisvert(20032421), Andrew Williams(20125276)\n",
    "#### Link to github repository: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import dataset "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import utils.mnist_reader as mnist_reader\n",
    "import pickle\n",
    "import gzip\n",
    "\n",
    "f = gzip.open(r'C:\\Users\\User\\Desktop\\UDEM\\IFT6135H19_assignment-master\\mnist.pkl.gz', 'rb')\n",
    "train_set, valid_set, test_set = pickle.load(f, encoding='latin1')\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, y_train = train_set\n",
    "X_test, y_test = test_set\n",
    "X_valid, Y_valid = valid_set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Building the Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import random\n",
    "import numpy as np\n",
    "import sys\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "class MLPClassifier(object):\n",
    "    \n",
    "    def __init__(self, d, d_h,d_h2, m, hidden_dims=(512,512),n_hidden=2,mode='train', datapath=None,model_path=None):\n",
    "        self.m = m\n",
    "        self.d_h = d_h\n",
    "        self.hidden_dims = hidden_dims\n",
    "        self.n_hidden = n_hidden\n",
    "        self.mode = mode\n",
    "        self.datapath = datapath\n",
    "        self.model_path = model_path\n",
    "        self.binary_y = 0\n",
    "        \n",
    "        #parameters\n",
    "        self.W1 = np.random.uniform(-1/np.sqrt(d_h), 1/np.sqrt(d_h), (d_h, d))\n",
    "        self.W2 = np.random.uniform(-1/np.sqrt(d_h2), 1/np.sqrt(d_h2), (d_h2, d_h))\n",
    "        self.W3 = np.random.uniform(-1/np.sqrt(m), 1/np.sqrt(m), (m, d_h2))\n",
    "        \n",
    "        self.b1 = np.zeros((d_h, 1))\n",
    "        self.b2 = np.zeros((d_h2, 1))\n",
    "        self.b3 = np.zeros((m, 1))\n",
    "        \n",
    "        #activations\n",
    "        self.os = 0   # (m, K)\n",
    "        self.oa = 0   # (m, K)\n",
    "        \n",
    "        self.hs1 = 0   # (d_h, K)\n",
    "        self.ha1 = 0   # (d_h, K)\n",
    "        \n",
    "        self.hs2 = 0 \n",
    "        self.ha2 = 0\n",
    "        \n",
    "        #gradients\n",
    "        self.dW1 = np.zeros((d_h, d))\n",
    "        self.dW2 = np.zeros((d_h2, d_h))\n",
    "        self.dW3 = np.zeros((m, d_h2))\n",
    "        \n",
    "        self.db1 = np.zeros((d_h, 1))\n",
    "        self.db2 = np.zeros((d_h2, 1))  \n",
    "        self.db3 = np.zeros((m, 1))\n",
    "        \n",
    "        #losses and missclassification rates\n",
    "        self.loss_train =[]\n",
    "        self.missclass_train = []\n",
    "        \n",
    "        self.loss_valid = []\n",
    "        self.missclass_valid = []\n",
    "        \n",
    "        self.loss_test = []\n",
    "        self.missclass_test = []\n",
    "        \n",
    "        \n",
    "    def initialize_weights(self, n_hidden=2, dims=(784,10), initialization = None):\n",
    "        if initialization == None:\n",
    "            pass\n",
    "            \n",
    "        if initialization == 'normal':\n",
    "            self.W1 = np.random.normal(0,1 , (self.hidden_dims[0],dims[0]))\n",
    "            self.W2 = np.random.normal(0,1 , (self.hidden_dims[1],self.hidden_dims[0]))\n",
    "            self.W3 = np.random.normal(0,1 , (dims[1],self.hidden_dims[1]))\n",
    "            \n",
    "        if initialization == 'glorot':\n",
    "            self.W1 = np.random.uniform((-1*np.sqrt(6/(784+512))),(np.sqrt(6/(784+512))), \\\n",
    "                                        (self.hidden_dims[0],dims[0]))\n",
    "            self.W2 = np.random.uniform((-1*np.sqrt(6/(512+512))),(np.sqrt(6/(512+512))), \\\n",
    "                                        (self.hidden_dims[1],self.hidden_dims[0]))\n",
    "            self.W3 = np.random.uniform((-1*np.sqrt(6/(10+512))),(np.sqrt(6/(10+512))), \\\n",
    "                                        (dims[1],self.hidden_dims[1]))\n",
    "    \n",
    "    \n",
    "    def forward(self, X):        \n",
    "        #first layer prop\n",
    "        self.ha = np.dot(self.W1, X.T) + self.b1\n",
    "        self.hs = self.activation(self.ha)\n",
    "        \n",
    "        #second layer prop\n",
    "        self.ha2 = np.dot(self.W2, self.hs) + self.b2\n",
    "        self.hs2 = self.activation(self.ha2)\n",
    "        \n",
    "        #third layer prop\n",
    "        self.oa = np.dot(self.W3, self.hs2) + self.b3\n",
    "        self.os = self.softmax()\n",
    "     \n",
    "    \n",
    "    def activation(self, fct):\n",
    "        activation = [[max(0, i) for i in row] for row in fct]\n",
    "        return np.array(activation)\n",
    "    \n",
    "    \n",
    "    def loss(self, X, Y):\n",
    "        t = self.getbinary(Y)\n",
    "        self.forward(X)\n",
    "        K = self.os.shape[1]\n",
    "        if K == t.shape[1]:\n",
    "            return -1 * self.m * np.mean(t * np.log(self.os))\n",
    "        else :\n",
    "          print(\"The subset of example is not consistent between X and Y\")             \n",
    "    \n",
    "    \n",
    "    def softmax(self):\n",
    "        shift = np.max(self.oa, axis = 0)\n",
    "        shifted_output = self.oa - shift\n",
    "        normal = np.sum(np.exp(shifted_output), axis = 0)\n",
    "        return np.exp(shifted_output) / normal\n",
    "    \n",
    "    \n",
    "    #it will automatically update the parameters and the gradients\n",
    "    def backward(self, X, Y, lamb = np.zeros((3, 2))):\n",
    "        t = self.getbinary(Y)\n",
    "        K = X.shape[0]\n",
    "        \n",
    "        if K == t.shape[1]:\n",
    "            #backprop through activation\n",
    "            doa = self.os - t     # (m, K)\n",
    "            self.dW3 = 1 / K * np.dot(doa, self.hs2.T)   # (m, d_h)\n",
    "            self.db3 = np.mean(doa, axis = 1, keepdims = True)  # (m, 1)\n",
    "            \n",
    "            #backprop through second hidden layer\n",
    "            dha2 = np.dot(self.W3.T, doa) * self.positive(self.ha2) \n",
    "            self.dW2 = 1 / K * np.dot(dha2, self.hs.T)\n",
    "            self.db2 = np.mean(dha2, axis = 1, keepdims = True)\n",
    "            \n",
    "            #backprop through first hidden layer\n",
    "            dha = np.dot(self.W2.T, dha2) * self.positive(self.ha)  # (d_h, K)\n",
    "            self.dW1 = 1 / K * np.dot(dha, X)\n",
    "            self.db1 = np.mean(dha, axis = 1, keepdims = True)  #(d_h, 1)\n",
    "            \n",
    "            #regularization\n",
    "            self.dW3 += lamb[0][0]*np.sign(self.W3) + 2 * lamb[0][1] * self.W3\n",
    "            self.dW2 += lamb[1][0]*np.sign(self.W2) + 2 * lamb[1][1] * self.W2\n",
    "            self.dW1 += lamb[2][0]*np.sign(self.W1) + 2 * lamb[2][1] * self.W1\n",
    "              \n",
    "        else:\n",
    "            print(\"The subset of example is not consistent between X and Y\")\n",
    "    \n",
    "    \n",
    "    def update(self, eta):\n",
    "        self.W3 -= eta * self.dW3\n",
    "        self.b3 -= eta * self.db3\n",
    "        \n",
    "        self.W2 -= eta * self.dW2\n",
    "        self.b2 -= eta * self.db2\n",
    "        \n",
    "        self.W1 -= eta * self.dW1\n",
    "        self.b1 -= eta * self.db1\n",
    "                \n",
    "                \n",
    "    #train with training set X and targets t. The gradients steps are done with batches of K examples \n",
    "    def train(self, X, Y, K, eta=0.07, epsilon=1e-3, max_ite=50000, X_valid=np.array([]), Y_valid=np.array([]), \\\n",
    "              n_rep_early_stopping=2, print_loss = False, lamb = np.zeros((3, 2)), X_test=[], \\\n",
    "              Y_test=[], print_graphs= False, initialization = 'glorot'):\n",
    "        \n",
    "        #bool that decides if early stopping is used\n",
    "        early_stopping = (X_valid.size != 0)\n",
    "        if early_stopping:\n",
    "            valid_error = self.compute_missclass(X_valid, Y_valid)\n",
    "        \n",
    "        #weight and iteration initialization\n",
    "        self.initialize_weights(initialization = initialization)\n",
    "        cond = True\n",
    "        ite = 0\n",
    "        ite_overfit = 0\n",
    "        \n",
    "        while cond:   \n",
    "            ## generate subset of training data\n",
    "            X_minibatch = X[[(i % X.shape[0]) for i in range(ite, ite + K)],:]\n",
    "            Y_minibatch = Y[[(i % X.shape[0]) for i in range(ite, ite + K)]]\n",
    "            \n",
    "            #training iteration\n",
    "            self.forward(X_minibatch)\n",
    "            self.backward(X_minibatch, Y_minibatch, lamb)\n",
    "            self.update(eta = eta)\n",
    "            \n",
    "            #if start of training or epoch finishes\n",
    "            if (ite)%5000 ==0:\n",
    "                if print_loss:\n",
    "                    print('Missclassification rate after epoch ', ite/5000, ' : ' , self.compute_missclass(X, Y))\n",
    "                    print('Average loss after epoch ', ite/5000, ' : ', self.loss(X,Y))\n",
    "                    \n",
    "                if print_graphs:\n",
    "                    self.loss_train.append(self.loss(X, Y))\n",
    "                    self.missclass_train.append(self.compute_missclass(X, Y))\n",
    "                    \n",
    "                    self.loss_valid.append(self.loss(X_valid, Y_valid))\n",
    "                    self.missclass_valid.append(self.compute_missclass(X_valid, Y_valid)) \n",
    "                    \n",
    "                    self.loss_test.append(self.loss(X_test, Y_test))\n",
    "                    self.missclass_test.append(self.compute_missclass(X_test, Y_test)) \n",
    "                \n",
    "                #smaller stepsize\n",
    "                if ite==35000:\n",
    "                    eta = 0.01\n",
    "                \n",
    "                #even smaller stepsize\n",
    "                if ite==40000:\n",
    "                    eta = 0.003\n",
    "                \n",
    "                if early_stopping:\n",
    "                    valid_error_after = self.compute_missclass(X_valid, Y_valid)\n",
    "                    #print(valid_error_after)\n",
    "                    if valid_error_after > valid_error:\n",
    "                        ite_overfit += 1\n",
    "                    else:\n",
    "                        ite_overfit = 0\n",
    "                    valid_error = valid_error_after\n",
    "                \n",
    "                if ite_overfit > n_rep_early_stopping:\n",
    "                    print(\"Exit to avoid overfit\")\n",
    "                    break \n",
    "            \n",
    "            ite += 1\n",
    "            \n",
    "            #training either converges or we run out of iterations\n",
    "            cond = self.norm_grad() > epsilon and ite < max_ite\n",
    "            \n",
    "        print(\"Missclassification ratio on training set\"+ str(self.compute_missclass(X, Y)*100)+ \" %\")\n",
    "        \n",
    "        if early_stopping:\n",
    "            print(\"Missclassification ratio on validation set\"+ str(self.compute_missclass(X_valid, Y_valid)*100)+ \" %\")\n",
    "    \n",
    "    \n",
    "    def test(self, X, Y):\n",
    "        print('test loss: ', self.loss(X,Y))\n",
    "        print('test missclassification rate: ', self.compute_missclass(X,Y))\n",
    "        \n",
    "##############################################################################\n",
    "######################### additional functionalities #########################\n",
    "\n",
    "    def predict(self, X):\n",
    "        self.forward(X)\n",
    "        return np.argmax(self.os, axis = 0)\n",
    "    \n",
    "    \n",
    "    def positive(self, fct):\n",
    "        positive = [[ int(i > 0) for i in row] for row in fct]\n",
    "        return np.array(positive)\n",
    "    \n",
    "    \n",
    "    #convert Y (K X 1) to a target matrix (m X K)\n",
    "    def getbinary(self, Y):\n",
    "        try: \n",
    "            return np.array([[int(Y[j] == i) for j in range(Y.shape[0])] for i in range(self.m)])\n",
    "        except IndexError:\n",
    "            return np.transpose(np.array([[int(Y == i)] for i in range(self.m)]))\n",
    "    \n",
    "    \n",
    "    #computes missclassification ratio      \n",
    "    def compute_missclass(self, X, Y):\n",
    "        t = self.getbinary(Y)\n",
    "        self.forward(X)\n",
    "        K = self.os.shape[1]\n",
    "        if K == t.shape[1]:\n",
    "            predictions = np.argmax(self.os, axis = 0)\n",
    "            return np.mean(predictions!=Y)\n",
    "        else :\n",
    "          print(\"The subset of example is not consistent between X and Y\") \n",
    "    \n",
    "    \n",
    "    #returns the norm of the gradient      \n",
    "    def norm_grad(self):\n",
    "        return np.sqrt(np.sum(self.dW3**2)+np.sum(self.db3**2)+np.sum(self.dW2**2)+np.sum(self.dW1**2)+\\\n",
    "                           np.sum(self.db2**2)+np.sum(self.db1**2))\n",
    "    \n",
    "    \n",
    "    #returns the approximative gradient using mini batches of size K\n",
    "    def finite_differences(self, X, Y, epsilon, print_grads=False):\n",
    "        differences_plus = np.empty(10)\n",
    "        differences_minus = np.empty(10)\n",
    "        finite_grad = np.empty(10)\n",
    "        \n",
    "        for i in range(10):\n",
    "            self.W2[0][i] += epsilon\n",
    "            differences_plus[i] = self.loss(X.reshape(-1,1).T,Y.reshape(-1,1).T)\n",
    "            self.W2[0][i] -= 2*epsilon\n",
    "            differences_minus[i] = self.loss(X.reshape(-1,1).T,Y.reshape(-1,1).T)\n",
    "            self.W2[0][i] += epsilon\n",
    "            finite_grad[i] = ((differences_plus[i]) - (differences_minus[i]))/(2*epsilon)\n",
    "        return finite_grad\n",
    "        \n",
    "        if print_grads:\n",
    "                print('dW3 = ', self.dW3)\n",
    "        else: return (self.dW3)\n",
    "###############################################################################    \n",
    "##example for debugging\n",
    "        \n",
    "# X = np.random.uniform(-1, 1, (100, 2))\n",
    "# Y = np.array([int(i) for i in X[:, 0]**2 + X[:, 1]**2 > 1/2])\n",
    "\n",
    "# Xvalid = np.random.uniform(-1, 1, (100, 2))\n",
    "# Yvalid = np.array([int(i) for i in Xvalid[:, 0]**2 + Xvalid[:, 1]**2 > 1/2])\n",
    "\n",
    "# NN = MLPClassifier(2, 3, 2)\n",
    "# NN.forward(X[0:50, :])\n",
    "# ##standard back prop prints dW2\n",
    "# NN.backward(X[0:50,:], Y[0:50], 0.01, print_grads=True, update=False)\n",
    "# NN.backward_loop(X[0:50,:], Y[0:50], 0.01, print_grads=True, update = False)\n",
    "# ##finite diff returns dW2\n",
    "# NN.finite_diff(X[0:50], Y[0:50], 1e-5)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# #%%\n",
    "# ##train on training set without \n",
    "# NN.train(X, Y, 50, epsilon=0.01, print_loss=False, X_valid=Xvalid, Y_valid=Yvalid)\n",
    "\n",
    "# plt.figure(1)\n",
    "# plt.scatter(X[:,0], X[:,1], c = Y)\n",
    "\n",
    "# plt.figure(2)\n",
    "# plt.scatter(X[:,0], X[:,1], c= NN.predict(X))\n",
    "\n",
    "# plt.figure(3)\n",
    "# plt.scatter(Xvalid[:,0], Xvalid[:,1], c=Yvalid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initializations\n",
    "    Aside from the weight initializations, the models were trained using the optimal hyperparameters foudn in the subsequent question."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Missclassification rate after epoch  0.0  :  0.8066\n",
      "Average loss after epoch  0.0  :  2.294620167995422\n",
      "Missclassification rate after epoch  1.0  :  0.09494\n",
      "Average loss after epoch  1.0  :  0.3544980922889262\n",
      "Missclassification rate after epoch  2.0  :  0.06168\n",
      "Average loss after epoch  2.0  :  0.2206068463090487\n",
      "Missclassification rate after epoch  3.0  :  0.05738\n",
      "Average loss after epoch  3.0  :  0.2114955154374316\n",
      "Missclassification rate after epoch  4.0  :  0.06218\n",
      "Average loss after epoch  4.0  :  0.23073012184585912\n",
      "Missclassification rate after epoch  5.0  :  0.04386\n",
      "Average loss after epoch  5.0  :  0.14928821987452173\n",
      "Missclassification rate after epoch  6.0  :  0.04896\n",
      "Average loss after epoch  6.0  :  0.18284297253315124\n",
      "Missclassification rate after epoch  7.0  :  0.03168\n",
      "Average loss after epoch  7.0  :  0.1090657316474467\n",
      "Missclassification rate after epoch  8.0  :  0.02412\n",
      "Average loss after epoch  8.0  :  0.08193156883655245\n",
      "Missclassification rate after epoch  9.0  :  0.02316\n",
      "Average loss after epoch  9.0  :  0.07447427755931804\n",
      "Missclassification ratio on training set2.168 %\n",
      "Missclassification ratio on validation set2.85 %\n",
      "Missclassification rate after epoch  0.0  :  0.88644\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:100: RuntimeWarning: divide by zero encountered in log\n",
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:100: RuntimeWarning: invalid value encountered in multiply\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average loss after epoch  0.0  :  nan\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:267: RuntimeWarning: overflow encountered in square\n",
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:129: RuntimeWarning: invalid value encountered in multiply\n",
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\numpy\\core\\fromnumeric.py:83: RuntimeWarning: invalid value encountered in reduce\n",
      "  return ufunc.reduce(obj, axis, dtype, out, **passkwargs)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Missclassification ratio on training set90.13600000000001 %\n",
      "Missclassification ratio on validation set90.09 %\n",
      "Missclassification rate after epoch  0.0  :  0.8255\n"
     ]
    }
   ],
   "source": [
    "#Zero weights \n",
    "NN_Mnist_zero = MLPClassifier(784,512,512,10)\n",
    "NN_Mnist_zero.train(X_train, y_train, 10, eta=0.05, max_ite = 50000 , epsilon=1e-6, print_loss=True,\n",
    "                    X_valid = X_valid, Y_valid= Y_valid, X_test = X_test, Y_test = y_test, \n",
    "                    lamb = ((0,0.00005),(0,0.00005),(0,0.00005)), print_graphs=True, initialization = None)\n",
    "\n",
    "#normal weights\n",
    "NN_Mnist_normal = MLPClassifier(784,512,512,10)\n",
    "NN_Mnist_normal.train(X_train, y_train, 10, eta=0.05, max_ite = 50000 , epsilon=1e-6, print_loss=True,\n",
    "                      X_valid = X_valid, Y_valid= Y_valid, X_test = X_test, Y_test = y_test, \n",
    "                      lamb = ((0,0.00005),(0,0.00005),(0,0.00005)), print_graphs=True, initialization = 'normal')\n",
    "\n",
    "#glorot weights\n",
    "NN_Mnist_glorot = MLPClassifier(784,512,512,10)\n",
    "NN_Mnist_glorot.train(X_train, y_train, 10, eta=0.05, max_ite = 50000 , epsilon=1e-6, print_loss=True,\n",
    "                      X_valid = X_valid, Y_valid= Y_valid, X_test = X_test, Y_test = y_test, \n",
    "                      lamb = ((0,0.00005),(0,0.00005),(0,0.00005)), print_graphs=True, initialization = 'glorot')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#average losses at the end of each epoch\n",
    "#and graphs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hyperparameter search"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1. Find a set of hyperparameters with validation accuracy over 97%\n",
    "\n",
    "    -512 neurons on each hidden layer\n",
    "    -learning rate schedule is 0.05 for first 7 epochs, 0.01 for the 8th, 0.003 for the last 2. \n",
    "    -ReLu activation function\n",
    "    -Gradient norm cutoff is 1e-6\n",
    "    -Number of epochs where valid_error increases before early stopping occurs is 3 in a row\n",
    "    -Minibatch of size 10\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Missclassification rate after epoch  0.0  :  0.87054\n",
      "Average loss after epoch  0.0  :  2.2750267300636295\n",
      "Missclassification rate after epoch  1.0  :  0.10668\n",
      "Average loss after epoch  1.0  :  0.4188884878623532\n",
      "Missclassification rate after epoch  2.0  :  0.07028\n",
      "Average loss after epoch  2.0  :  0.2539725050956212\n",
      "Missclassification rate after epoch  3.0  :  0.06082\n",
      "Average loss after epoch  3.0  :  0.22218040479558918\n",
      "Missclassification rate after epoch  4.0  :  0.05648\n",
      "Average loss after epoch  4.0  :  0.21306434965561363\n",
      "Missclassification rate after epoch  5.0  :  0.04566\n",
      "Average loss after epoch  5.0  :  0.1575525836962566\n",
      "Missclassification rate after epoch  6.0  :  0.0369\n",
      "Average loss after epoch  6.0  :  0.12876049213659402\n",
      "Missclassification rate after epoch  7.0  :  0.0289\n",
      "Average loss after epoch  7.0  :  0.09969934357079317\n",
      "Missclassification rate after epoch  8.0  :  0.02972\n",
      "Average loss after epoch  8.0  :  0.09742512583975697\n",
      "Missclassification rate after epoch  9.0  :  0.0282\n",
      "Average loss after epoch  9.0  :  0.08927742681770881\n",
      "Missclassificaiton ratio on training set2.606 %\n",
      "Missclassificaiton ratio on validation set3.06 %\n"
     ]
    }
   ],
   "source": [
    "#successful MLP: valid error is _____ missclassification and test error is ______\n",
    "np.random.seed(123)\n",
    "NN_Mnist = MLPClassifier(784,512,512,10)\n",
    "\n",
    "first = time.time()\n",
    "NN_Mnist.train(X_train, y_train, 10, eta=0.05, max_ite = 50000 , epsilon=1e-12, print_loss=True, \\\n",
    "               X_valid = X_valid, Y_valid= Y_valid, X_test = X_test, Y_test = y_test, \\\n",
    "               lamb = ((0,0.00005),(0,0.00005),(0,0.00005)), print_graphs=True, initialization = 'glorot')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Greedy search\n",
    "First we find the optimal stepsize with no regularization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stepsizes = [1,1e-1,1e-2,1e-3]\n",
    "r_valids = []\n",
    "for i in range(len(stepsizes)):\n",
    "    NN_Mnist_tune = MLPClassifier(784,512,512,10)\n",
    "    NN_Mnist_tune.train(X_train, y_train, 10, eta=stepsizes[i], max_ite = 50000 , epsilon=1e-6, \\print_loss=False,\\\n",
    "                        X_valid = X_valid, Y_valid= Y_valid, X_test = X_test, Y_test = y_test,\\\n",
    "                        lamb = lambdas[k], print_graphs=True, initialization = 'glorot')\n",
    "\n",
    "    r_valids.append((NN_Mnist_tune.compute_missclass(X_valid,Y_valid),stepsizes[i]))\n",
    "\n",
    "print(r_valids)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Best stepsize is 0.01. Test for values around this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stepsizes2 = [0.05,0.01,0.005]\n",
    "r_valids_2 = []\n",
    "for i in range(len(stepsizes)):    \n",
    "    NN_Mnist_tune = MLPClassifier(784,512,512,10)\n",
    "    NN_Mnist_tune.train(X_train, y_train, 10, eta=stepsizes[i], max_ite = 50000 , epsilon=1e-6, print_loss=False,\\\n",
    "                        X_valid = X_valid, Y_valid= Y_valid, X_test = X_test, Y_test = y_test, \\\n",
    "                        print_graphs=True, initialization = 'glorot')\n",
    "\n",
    "    r_valids_2.append((NN_Mnist_tune.compute_missclass(X_valid,Y_valid),stepsizes2[i]))\n",
    "\n",
    "print(r_valids_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Best stepsize is 0.05. Add decreasing regime of stepsizes for later iterations(hardcoded into class).\n",
    "\n",
    "Find best lambda for l2 regularization. l1 regularization is coded but not used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lambdas = [((0,0.1),(0,0.1),(0,0.1)),((0,0.001),(0,0.001),(0,0.001)),((0,0.0001),(0,0.0001),(0,0.0001))]\n",
    "r_valids_3 = []\n",
    "for i in range(len(lambdas)):\n",
    "    NN_Mnist_tune = MLPClassifier(784,512,512,10)\n",
    "    NN_Mnist_tune.train(X_train, y_train, 10, eta=0.05, max_ite = 50000 , epsilon=1e-6, print_loss=False, \\\n",
    "                        X_valid = X_valid, Y_valid= Y_valid, X_test = X_test, Y_test = y_test, \\\n",
    "                        lamb = lambdas[i], print_graphs=True, initialization = 'glorot')\n",
    "\n",
    "        r_valids_3.append((NN_Mnist_tune.compute_missclass(X_valid,Y_valid),lambdas[i]))\n",
    "\n",
    "print(r_valids_3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Best lambda is 0.0001. With additional tuning, the ideal lambda is found to be 0.00005."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Validate Gradients using Finite Differences"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5 different values of N = 1/epsilon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "      <th>5</th>\n",
       "      <th>6</th>\n",
       "      <th>7</th>\n",
       "      <th>8</th>\n",
       "      <th>9</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>N = 5</th>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.059965</td>\n",
       "      <td>-0.030483</td>\n",
       "      <td>-0.031991</td>\n",
       "      <td>-0.01183</td>\n",
       "      <td>-0.030912</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.044868</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.035811</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>N = 50</th>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.089845</td>\n",
       "      <td>-0.030483</td>\n",
       "      <td>-0.032988</td>\n",
       "      <td>-0.01183</td>\n",
       "      <td>-0.030912</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.059055</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.040703</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>N = 500</th>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.089845</td>\n",
       "      <td>-0.030483</td>\n",
       "      <td>-0.032988</td>\n",
       "      <td>-0.01183</td>\n",
       "      <td>-0.030912</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.059055</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.040703</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>N = 5000</th>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.089845</td>\n",
       "      <td>-0.030483</td>\n",
       "      <td>-0.032988</td>\n",
       "      <td>-0.01183</td>\n",
       "      <td>-0.030912</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.059055</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.040703</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>N = 50000</th>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.089845</td>\n",
       "      <td>-0.030483</td>\n",
       "      <td>-0.032988</td>\n",
       "      <td>-0.01183</td>\n",
       "      <td>-0.030912</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.059055</td>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.040703</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             0         1         2         3        4         5    6  \\\n",
       "N = 5      0.0 -0.059965 -0.030483 -0.031991 -0.01183 -0.030912  0.0   \n",
       "N = 50     0.0 -0.089845 -0.030483 -0.032988 -0.01183 -0.030912  0.0   \n",
       "N = 500    0.0 -0.089845 -0.030483 -0.032988 -0.01183 -0.030912  0.0   \n",
       "N = 5000   0.0 -0.089845 -0.030483 -0.032988 -0.01183 -0.030912  0.0   \n",
       "N = 50000  0.0 -0.089845 -0.030483 -0.032988 -0.01183 -0.030912  0.0   \n",
       "\n",
       "                  7    8         9  \n",
       "N = 5     -0.044868  0.0 -0.035811  \n",
       "N = 50    -0.059055  0.0 -0.040703  \n",
       "N = 500   -0.059055  0.0 -0.040703  \n",
       "N = 5000  -0.059055  0.0 -0.040703  \n",
       "N = 50000 -0.059055  0.0 -0.040703  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.seed(111)\n",
    "\n",
    "NN_Mnist = MLPClassifier(784,512,512,10)\n",
    "NN_Mnist.forward(X_train[150].reshape(-1,1).T)\n",
    "NN_Mnist.backward(X_train[150].reshape(-1,1).T, y_train[150].reshape(-1,1).T)\n",
    "\n",
    "finite_differences = []\n",
    "for i in range(5):\n",
    "    finite_differences.append(NN_Mnist.finite_differences(X_train[150],y_train[150], epsilon = 1/((5)*10**(i))))\n",
    "finite_differences = np.array(finite_differences)\n",
    "\n",
    "import pandas \n",
    "row_labels = ['N = 5','N = 50','N = 500','N = 5000','N = 50000']\n",
    "column_labels = []\n",
    "pandas.DataFrame(finite_differences, index = row_labels, \\\n",
    "                 title = 'Values of finite difference approximations for first 10 parameters')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot max differences as a function of N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_differences = np.max(np.abs(finite_differences - NN_Mnist.dW2[0][0:10]), axis = 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3 align=\"center\"> Plot of $ \\left(\\max_{1\\leq i \\leq p} | \\nabla_i^N - \\frac{\\partial L}{\\partial \\theta_i}| \\right)$ as a function of N <h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "testing = pd.DataFrame(max_differences).transpose()\n",
    "testing.columns = [5,50,500,5000,50000]\n",
    "testing.index.name = 'param'\n",
    "testing.columns.name = 'Value of N'\n",
    "testing.plot(style = '.', grid=True, linestyle = 'dotted')"
   ]
  }
 ],
 "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": 2
}