{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Summary\n\nThis note contains an **explicit** routine of a neural network for MNIST digit recognizer implemented in numpy from scratch.\n\nThis is a simple demonstration mainly for pedagogical purposes, which shows the basic workflow of a machine learning algorithm using a simple feedforward neural network. The derivative at the backpropagation stage is computed explicitly through the chain rule.\n\n* The model is a 3-layer feedforward neural network (FNN), in which the input layer has 784 units, and the 256-unit hidden layer is activated by ReLU, while the output layer is activated by softmax function to produce a discrete probability distribution for each input. \n\n* The loss function, model hypothesis function, and the gradient of the loss function are all implemented from ground-up in numpy in a highly vectorized fashion (no FOR loops).\n\n* The training is through a standard gradient descent with step size adaptively changing by Root Mean Square prop (RMSprop), and there is no cross-validation set reserved nor model averaging for simplicity.\n\n\nThe code is vectorized and is adapted from the Softmax regression and Neural Network lectures used in [UCI Math 10](https://github.com/scaomath/UCI-Math10). \n\nCaveat lector: fluency in linear algebra and multivariate calculus.\n\n\n#### References:\n* [Stanford Deep Learning tutorial in MATLAB](http://ufldl.stanford.edu/tutorial/supervised/MultiLayerNeuralNetworks/).\n* [Learning PyTorch with examples](https://pytorch.org/tutorials/beginner/pytorch_with_examples.html).\n* [An overview of gradient descent optimization algorithms](http://ruder.io/optimizing-gradient-descent/).","metadata":{}},{"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib.pyplot as plt\n%matplotlib inline","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-04-24T19:58:45.797136Z","iopub.execute_input":"2023-04-24T19:58:45.797981Z","iopub.status.idle":"2023-04-24T19:58:45.809795Z","shell.execute_reply.started":"2023-04-24T19:58:45.797616Z","shell.execute_reply":"2023-04-24T19:58:45.808709Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Data input","metadata":{}},{"cell_type":"code","source":"import os\npath = os.listdir(\"../input\")\nprint(path)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-04-24T19:58:45.812076Z","iopub.execute_input":"2023-04-24T19:58:45.812420Z","iopub.status.idle":"2023-04-24T19:58:45.822382Z","shell.execute_reply.started":"2023-04-24T19:58:45.812368Z","shell.execute_reply":"2023-04-24T19:58:45.821366Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Read the data\ntrain_data = pd.read_csv('../input/train.csv')\ntest_data = pd.read_csv(\"../input/test.csv\")","metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","execution":{"iopub.status.busy":"2023-04-24T19:58:45.823906Z","iopub.execute_input":"2023-04-24T19:58:45.824362Z","iopub.status.idle":"2023-04-24T19:58:54.178722Z","shell.execute_reply.started":"2023-04-24T19:58:45.824291Z","shell.execute_reply":"2023-04-24T19:58:54.177402Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Set up the data\ny_train = train_data['label'].values\nX_train = train_data.drop(columns=['label']).values/255\nX_test = test_data.values/255","metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","execution":{"iopub.status.busy":"2023-04-24T19:58:54.180611Z","iopub.execute_input":"2023-04-24T19:58:54.181026Z","iopub.status.idle":"2023-04-24T19:58:55.127624Z","shell.execute_reply.started":"2023-04-24T19:58:54.180949Z","shell.execute_reply":"2023-04-24T19:58:55.126463Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Some examples of the input data\nWe randomly choose 10 samples from the training set and visualize it.","metadata":{}},{"cell_type":"code","source":"fig, axes = plt.subplots(2,5, figsize=(12,5))\naxes = axes.flatten()\nidx = np.random.randint(0,42000,size=10)\nfor i in range(10):\n axes[i].imshow(X_train[idx[i],:].reshape(28,28), cmap='gray')\n axes[i].axis('off') # hide the axes ticks\n axes[i].set_title(str(int(y_train[idx[i]])), color= 'black', fontsize=25)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-04-24T19:58:55.133680Z","iopub.execute_input":"2023-04-24T19:58:55.134055Z","iopub.status.idle":"2023-04-24T19:58:56.039453Z","shell.execute_reply.started":"2023-04-24T19:58:55.134002Z","shell.execute_reply":"2023-04-24T19:58:56.038409Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Network structures\n\n\"drawing\"\n\nThe figure above is a simplication of the neural network used in this example. The circles labeled \"+1\" are the bias units. Layer 1 is the input layer, and Layer 3 is the output layer. The middle layer, Layer 2, is the hidden layer.\n\nThe neural network in the figure above has 2 input units (not counting the bias unit), 3 hidden units, and 1 output unit. In this actual computation below, the input layer has 784 units, the hidden layer has 256 units, and the output layers has 10 units ($K =10$ classes).\n\nThe weight matrix $W^{(0)}$ mapping input $\\mathbf{x}$ from the input layer (Layer 1) to the hidden layer (Layer 2) is of shape `(784,256)` together with a `(256,)` bias. Then $\\mathbf{a}$ is the activation from the hidden layer (Layer 2) can be written as:\n$$\n\\mathbf{a} = \\mathrm{ReLU}\\big((W^{(0)})^{\\top}\\mathbf{x} + \\mathbf{b}\\big),\n$$\nwhere the ReLU activation function is $\\mathrm{ReLU}(z) = \\max(z,0)$ and can be implemented in a vectorized fashion as follows.","metadata":{}},{"cell_type":"code","source":"# relu activation function\n# THE fastest vectorized implementation for ReLU\ndef relu(x):\n x[x<0]=0\n return x","metadata":{"execution":{"iopub.status.busy":"2023-04-24T19:58:56.041843Z","iopub.execute_input":"2023-04-24T19:58:56.042616Z","iopub.status.idle":"2023-04-24T19:58:56.049138Z","shell.execute_reply.started":"2023-04-24T19:58:56.042537Z","shell.execute_reply":"2023-04-24T19:58:56.047723Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Softmax activation, prediction, and the loss function\n\nFrom the hidden layer (Layer 2) to the output layer (layer 3), the weight matrix $W^{(1)}$ is of shape `(256,10)`, the form of which is as follows:\n$$\nW^{(1)} =\n\\begin{pmatrix}\n| & | & | & | \\\\\n\\boldsymbol{\\theta}_1 & \\boldsymbol{\\theta}_2 & \\cdots & \\boldsymbol{\\theta}_K \\\\\n| & | & | & |\n\\end{pmatrix},\n$$\nwhich maps the activation from Layer 2 to Layer 3 (output layer), and there is no bias because a constant can be freely added to the activation without changing the final output. \n\nAt the last layer, a softmax activation is used, which can be written as follows combining the weights matrix $W^{(1)}$ that maps the activation $\\mathbf{a}$ from the hidden layer to output layer:\n$$\nP\\big(y = k \\;| \\;\\mathbf{a}; W^{(1)}\\big) = \\sigma_k(\\mathbf{a}; W^{(1)}) := \\frac{\\exp\\big(\\boldsymbol{\\theta}^{\\top}_k \\mathbf{a} \\big)}\n{\\sum_{j=1}^K \\exp\\big(\\boldsymbol{\\theta}^{\\top}_j \\mathbf{a} \\big)}.\n$$\n$\\{P\\big(y = k \\;| \\;\\mathbf{a}; W^{(1)}\\big)\\}_{k=1}^K$ is the probability distribution of our model, which estimates the probability of the input $\\mathbf{x}$'s label $y$ is of class $k$. We denote this distribution by a vector \n$$\\boldsymbol{\\sigma}:= (\\sigma_1,\\dots, \\sigma_K)^{\\top}.$$\nWe hope that this estimate is as close as possible to the true probability: $1_{\\{y=k\\}}$, that is $1$ if the sample $\\mathbf{x}$ is in the $k$-th class and 0 otherwise. \n\nLastly, our prediction $\\hat{y}$ for sample $\\mathbf{x}$ can be made by choosing the class with the highest probability:\n$$\n\\hat{y} = \\operatorname{argmax}_{k=1,\\dots,K} P\\big(y = k \\;| \\;\\mathbf{a}; W^{(1)}\\big). \\tag{$\\ast$}\n$$\n\nDenote the label of the $i$-th input as $y^{(i)}$, and then the sample-wise loss function is the cross entropy measuring the difference of the distribution of this model function above with the true one $1_{\\{y^{(i)}=k\\}}$: denote $W = (W^{(0)}, W^{(1)})$, $b = (\\mathbf{b})$, let $\\mathbf{a}^{(i)}$ be the activation for the $i$-th sample in the hidden layer (Layer 2),\n$$\nJ_i:= J(W,b;\\mathbf{x}^{(i)},y^{(i)}) := - \\sum_{k=1}^{K} \\left\\{ 1_{\\left\\{y^{(i)} = k\\right\\} }\n\\log P\\big(y^{(i)} = k \\;| \\;\\mathbf{a}^{(i)}; W^{(1)}\\big)\\right\\}. \\tag{1}\n$$\n\nDenote the data sample matrix $X := (\\mathbf{x}^{(1)}, \\dots, \\mathbf{x}^{(N)})^{\\top}$, its label vector as $\\mathbf{y} := (y^{(1)}, \\dots, y^{(N)})$, and then the final loss has an extra $L^2$-regularization term for the weight matrices (not for bias): \n$$\nL(W,b; X, \\mathbf{y}) := \\frac{1}{N}\\sum_{i=1}^{N} J_i + \\frac{\\alpha}{2} \\Big(\\|W^{(0)}\\|^2 + \\|W^{(1)}\\|^2\\Big),\n\\tag{2}\n$$\nwhere $\\alpha>0$ is a hyper-parameter determining the strength of the regularization, the bigger the $\\alpha$ is, the smaller the magnitudes of the weights will be after training.","metadata":{}},{"cell_type":"code","source":"def h(X,W,b):\n '''\n Hypothesis function: simple FNN with 1 hidden layer\n Layer 1: input\n Layer 2: hidden layer, with a size implied by the arguments W[0], b\n Layer 3: output layer, with a size implied by the arguments W[1]\n '''\n # layer 1 = input layer\n a1 = X\n # layer 1 (input layer) -> layer 2 (hidden layer)\n z1 = np.matmul(X, W[0]) + b[0]\n \n # add one more layer if applicable\n \n # layer 2 activation\n a2 = relu(z1)\n # layer 2 (hidden layer) -> layer 3 (output layer)\n z2 = np.matmul(a2, W[1])\n s = np.exp(z2)\n total = np.sum(s, axis=1).reshape(-1,1)\n sigma = s/total\n # the output is a probability for each sample\n return sigma","metadata":{"execution":{"iopub.status.busy":"2023-04-24T19:58:56.050998Z","iopub.execute_input":"2023-04-24T19:58:56.051837Z","iopub.status.idle":"2023-04-24T19:58:56.063507Z","shell.execute_reply.started":"2023-04-24T19:58:56.051418Z","shell.execute_reply":"2023-04-24T19:58:56.062164Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def softmax(X_in,weights):\n '''\n Un-used cell for demo\n activation function for the last FC layer: softmax function \n Output: K probabilities represent an estimate of P(y=k|X_in;weights) for k=1,...,K\n the weights has shape (n, K)\n n: the number of features X_in has\n n = X_in.shape[1]\n K: the number of classes\n K = 10\n '''\n \n s = np.exp(np.matmul(X_in,weights))\n total = np.sum(s, axis=1).reshape(-1,1)\n return s / total","metadata":{"execution":{"iopub.status.busy":"2023-04-24T19:58:56.065622Z","iopub.execute_input":"2023-04-24T19:58:56.066580Z","iopub.status.idle":"2023-04-24T19:58:56.078499Z","shell.execute_reply.started":"2023-04-24T19:58:56.066191Z","shell.execute_reply":"2023-04-24T19:58:56.076889Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def loss(y_pred,y_true):\n '''\n Loss function: cross entropy with an L^2 regularization\n y_true: ground truth, of shape (N, )\n y_pred: prediction made by the model, of shape (N, K) \n N: number of samples in the batch\n K: global variable, number of classes\n '''\n global K \n K = 10\n N = len(y_true)\n # loss_sample stores the cross entropy for each sample in X\n # convert y_true from labels to one-hot-vector encoding\n y_true_one_hot_vec = (y_true[:,np.newaxis] == np.arange(K))\n loss_sample = (np.log(y_pred) * y_true_one_hot_vec).sum(axis=1)\n # loss_sample is a dimension (N,) array\n # for the final loss, we need take the average\n return -np.mean(loss_sample)","metadata":{"execution":{"iopub.status.busy":"2023-04-24T19:58:56.080760Z","iopub.execute_input":"2023-04-24T19:58:56.081238Z","iopub.status.idle":"2023-04-24T19:58:56.093814Z","shell.execute_reply.started":"2023-04-24T19:58:56.081156Z","shell.execute_reply":"2023-04-24T19:58:56.092959Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Backpropagation (Chain rule)\n\nThe derivative of the cross entropy $J$ in (1), for a single sample and its label $(\\mathbf{x}, y)$ , with respect to the weights and the bias is computed using the following procedure:\n> **Step 1**: Forward pass: computing the activations $\\mathbf{a} = (a_1,\\dots, a_{n_2})$ from the hidden layer (Layer 2), and $\\boldsymbol{\\sigma} = (\\sigma_1,\\dots, \\sigma_K)$ from the output layer (Layer 3). \n>\n> **Step 2**: Derivatives for $W^{(1)}$: recall that $W^{(1)} = (\\boldsymbol{\\theta}_1 ,\\cdots, \\boldsymbol{\\theta}_K)$ and denote \n$$\\mathbf{z}^{(2)} = \\big(z^{(2)}_1, \\dots, z^{(2)}_K\\big) = (W^{(1)})^{\\top}\\mathbf{a} =\n\\big(\\boldsymbol{\\theta}^{\\top}_1 \\mathbf{a} ,\\cdots, \\boldsymbol{\\theta}^{\\top}_K \\mathbf{a}\\big),$$ \nfor the $k$-th output unit in the output layer (Layer 3), then\n$$\n\\delta^{(2)}_k\n:= \\frac{\\partial J}{\\partial z_k^{(2)}} = \\Big\\{ P\\big(y = k \\;| \\;\\mathbf{a}; W^{(1)}\\big)- 1_{\\{ y = k\\}} \\Big\\} = \\sigma_k - 1_{\\{ y = k\\}}\n$$\nand \n$$\n\\frac{\\partial J}{\\partial \\boldsymbol{\\theta}_k}= \\frac{\\partial J}{\\partial z_k^{(2)}}\\frac{\\partial z_k^{(2)}}{\\partial \\boldsymbol{\\theta}_k} = \\delta^{(2)}_k \\mathbf{a}.\n$$\n>\n> **Step 3**: Derivatives for $W^{(0)}$, $\\mathbf{b}$: recall that $W^{(0)} = (\\boldsymbol{w}_1 ,\\cdots, \\boldsymbol{w}_{n_2})$, $\\mathbf{b} = (b_1,\\dots, b_{n_2})$, where $n_2$ is the number of units in the hidden layer (Layer 2), and denote \n$$\\mathbf{z}^{(1)} = (z_1^{(1)}, \\dots, z_{n_2}^{(1)}) = (W^{(0)})^{\\top}\\mathbf{x} + \\mathbf{b} =\n\\big(\\mathbf{w}^{\\top}_1 \\mathbf{x} +b_1 ,\\cdots, \\mathbf{w}^{\\top}_{n_2} \\mathbf{x} + b_{n_2}\\big),$$ \nfor each node $i$ in the hidden layer (Layer $2$), $i=1,\\dots, n_2$, then\n$$\\delta^{(1)}_i : = \\frac{\\partial J}{\\partial z^{(1)}_i} =\n\\frac{\\partial J}{\\partial a_i} \n\\frac{\\partial a_i}{\\partial z^{(1)}_i}=\n\\frac{\\partial J}{\\partial \\mathbf{z}^{(2)}}\n\\cdot\\left(\\frac{\\partial \\mathbf{z}^{(2)}}{\\partial a_i} \n\\frac{\\partial a_i}{\\partial z^{(1)}_i}\\right)\n\\\\\n=\\left( \\sum_{k=1}^{K} \\frac{\\partial J}{\\partial {z}^{(2)}_k}\n\\frac{\\partial {z}^{(2)}_k}{\\partial a_i} \\right) f'(z^{(1)}_i) = \\left( \\sum_{k=1}^{K} w_{ki} \\delta^{(2)}_k \\right) 1_{\\{z^{(1)}_i\\; > 0\\}},\n$$\nwhere $1_{\\{z^{(1)}_i\\; > 0\\}}$ is ReLU activation $f$'s (weak) derivative, and the partial derivative of the $k$-th component (before activated by the softmax) in the output layer ${z}^{(2)}_k$ with respect to the $i$-th activation $a_i$ from the hidden layer is the weight $w^{(1)}_{ki}$. Thus\n>\n$$\n\\frac{\\partial J}{\\partial w_{ji}} = x_j \\delta_i^{(1)} ,\\;\n\\frac{\\partial J}{\\partial b_{i}} = \\delta_i^{(1)}, \\;\\text{ and }\\;\n\\frac{\\partial J}{\\partial \\mathbf{w}_{i}} = \\delta_i^{(1)}\\mathbf{x} ,\\;\n\\frac{\\partial J}{\\partial \\mathbf{b}} = \\boldsymbol{\\delta}^{(1)}.\n$$","metadata":{}},{"cell_type":"code","source":"def backprop(W,b,X,y,alpha=1e-4):\n '''\n Step 1: explicit forward pass h(X;W,b)\n Step 2: backpropagation for dW and db\n '''\n K = 10\n N = X.shape[0]\n \n ### Step 1:\n # layer 1 = input layer\n a1 = X\n # layer 1 (input layer) -> layer 2 (hidden layer)\n z1 = np.matmul(X, W[0]) + b[0]\n # layer 2 activation\n a2 = relu(z1)\n \n # one more layer\n \n # layer 2 (hidden layer) -> layer 3 (output layer)\n z2 = np.matmul(a2, W[1])\n s = np.exp(z2)\n total = np.sum(s, axis=1).reshape(-1,1)\n sigma = s/total\n \n ### Step 2:\n \n # layer 2->layer 3 weights' derivative\n # delta2 is \\partial L/partial z2, of shape (N,K)\n y_one_hot_vec = (y[:,np.newaxis] == np.arange(K))\n delta2 = (sigma - y_one_hot_vec)\n grad_W1 = np.matmul(a2.T, delta2)\n \n # layer 1->layer 2 weights' derivative\n # delta1 is \\partial a2/partial z1\n # layer 2 activation's (weak) derivative is 1*(z1>0)\n delta1 = np.matmul(delta2, W[1].T)*(z1>0)\n grad_W0 = np.matmul(X.T, delta1)\n \n # Possible student project: extra layer of derivative\n \n # no derivative for layer 1\n \n # the alpha part is the derivative for the regularization\n # regularization = 0.5*alpha*(np.sum(W[1]**2) + np.sum(W[0]**2))\n \n \n dW = [grad_W0/N + alpha*W[0], grad_W1/N + alpha*W[1]]\n db = [np.mean(delta1, axis=0)]\n # dW[0] is W[0]'s derivative, and dW[1] is W[1]'s derivative; similar for db\n return dW, db","metadata":{"execution":{"iopub.status.busy":"2023-04-24T19:58:56.095631Z","iopub.execute_input":"2023-04-24T19:58:56.096320Z","iopub.status.idle":"2023-04-24T19:58:56.109995Z","shell.execute_reply.started":"2023-04-24T19:58:56.096254Z","shell.execute_reply":"2023-04-24T19:58:56.108919Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Hyper-parameters and network initialization","metadata":{}},{"cell_type":"code","source":"eta = 5e-1\nalpha = 1e-6 # regularization\ngamma = 0.99 # RMSprop\neps = 1e-3 # RMSprop\nnum_iter = 1000 # number of iterations of gradient descent\nn_H = 256 # number of neurons in the hidden layer\nn = X_train.shape[1] # number of pixels in an image\nK = 10","metadata":{"execution":{"iopub.status.busy":"2023-04-24T19:58:56.111809Z","iopub.execute_input":"2023-04-24T19:58:56.112210Z","iopub.status.idle":"2023-04-24T19:58:56.124823Z","shell.execute_reply.started":"2023-04-24T19:58:56.112131Z","shell.execute_reply":"2023-04-24T19:58:56.123696Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# initialization\nnp.random.seed(1127825)\nW = [1e-1*np.random.randn(n, n_H), 1e-1*np.random.randn(n_H, K)]\nb = [np.random.randn(n_H)]","metadata":{"execution":{"iopub.status.busy":"2023-04-24T19:58:56.126882Z","iopub.execute_input":"2023-04-24T19:58:56.128419Z","iopub.status.idle":"2023-04-24T19:58:56.150733Z","shell.execute_reply.started":"2023-04-24T19:58:56.127237Z","shell.execute_reply":"2023-04-24T19:58:56.149531Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Gradient Descent: training of the network\n\nIn the training, we use a GD-variant of the RMSprop: for $\\mathbf{w}$ which stands for the parameter vector in our model\n> Choose $\\mathbf{w}_0$, $\\eta$, $\\gamma$, $\\epsilon$, and let $g_{-1} = 1$

\n> For $k=0,1,2, \\cdots, M$

\n>      $g_{k} = \\gamma g_{k-1} + (1 - \\gamma)\\, \\left|\\partial_{\\mathbf{w}} L (\\mathbf{w}_k)\\right|^2$

\n>      $\\displaystyle\\mathbf{w}_{k+1} = \\mathbf{w}_k - \\frac{\\eta} {\\sqrt{g_{k}+ \\epsilon}} \\partial_{\\mathbf{w}} L(\\mathbf{w}_k)$ \n\n### Remark: \nThe training takes a while since we use the gradient descent for all samples.","metadata":{}},{"cell_type":"code","source":"%%time\ngW0 = gW1 = gb0 = 1\n\nfor i in range(num_iter):\n dW, db = backprop(W,b,X_train,y_train,alpha)\n \n gW0 = gamma*gW0 + (1-gamma)*np.sum(dW[0]**2)\n etaW0 = eta/np.sqrt(gW0 + eps)\n W[0] -= etaW0 * dW[0]\n \n gW1 = gamma*gW1 + (1-gamma)*np.sum(dW[1]**2)\n etaW1 = eta/np.sqrt(gW1 + eps)\n W[1] -= etaW1 * dW[1]\n \n gb0 = gamma*gb0 + (1-gamma)*np.sum(db[0]**2)\n etab0 = eta/np.sqrt(gb0 + eps)\n b[0] -= etab0 * db[0]\n \n if i % 500 == 0:\n # sanity check 1\n y_pred = h(X_train,W,b)\n print(\"Cross-entropy loss after\", i+1, \"iterations is {:.8}\".format(\n loss(y_pred,y_train)))\n print(\"Training accuracy after\", i+1, \"iterations is {:.4%}\".format( \n np.mean(np.argmax(y_pred, axis=1)== y_train)))\n \n # sanity check 2\n print(\"gW0={:.4f} gW1={:.4f} gb0={:.4f}\\netaW0={:.4f} etaW1={:.4f} etab0={:.4f}\"\n .format(gW0, gW1, gb0, etaW0, etaW1, etab0))\n \n # sanity check 3\n print(\"|dW0|={:.5f} |dW1|={:.5f} |db0|={:.5f}\"\n .format(np.linalg.norm(dW[0]), np.linalg.norm(dW[1]), np.linalg.norm(db[0])), \"\\n\")\n \n # reset RMSprop\n gW0 = gW1 = gb0 = 1\n\ny_pred_final = h(X_train,W,b)\nprint(\"Final cross-entropy loss is {:.8}\".format(loss(y_pred_final,y_train)))\nprint(\"Final training accuracy is {:.4%}\".format(np.mean(np.argmax(y_pred_final, axis=1)== y_train)))","metadata":{"_kg_hide-output":false,"execution":{"iopub.status.busy":"2023-04-24T19:58:56.152706Z","iopub.execute_input":"2023-04-24T19:58:56.153301Z","iopub.status.idle":"2023-04-24T20:18:17.441060Z","shell.execute_reply.started":"2023-04-24T19:58:56.153050Z","shell.execute_reply":"2023-04-24T20:18:17.439656Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Predictions for testing data\nThe prediction labels are generated by $(\\ast)$.","metadata":{}},{"cell_type":"code","source":"# predictions\ny_pred_test = np.argmax(h(X_test,W,b), axis=1)","metadata":{"execution":{"iopub.status.busy":"2023-04-24T20:18:17.443309Z","iopub.execute_input":"2023-04-24T20:18:17.444066Z","iopub.status.idle":"2023-04-24T20:18:17.891213Z","shell.execute_reply.started":"2023-04-24T20:18:17.443981Z","shell.execute_reply":"2023-04-24T20:18:17.889355Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Generating submission using pandas for grading\nsubmission = pd.DataFrame({'ImageId': range(1,len(X_test)+1) ,'Label': y_pred_test })\nsubmission.to_csv(\"simplemnist_result.csv\",index=False)","metadata":{"execution":{"iopub.status.busy":"2023-04-24T20:18:17.899726Z","iopub.execute_input":"2023-04-24T20:18:17.900945Z","iopub.status.idle":"2023-04-24T20:18:18.159706Z","shell.execute_reply.started":"2023-04-24T20:18:17.900483Z","shell.execute_reply":"2023-04-24T20:18:18.158292Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Group members: Alicia Chuang, Vinachau Pham, Tyler Brown","metadata":{}}]}