notebooks_en/8_More_about_optimization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bb6b004f",
   "metadata": {},
   "source": [
    "# More about optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96a88e50",
   "metadata": {},
   "source": [
    "In previous lessons, we use *optimization*, *learning*, and *training* interchangeably. They all mean finding the best set of model parameters. However, optimization has a broader meaning -- any applications in which we want to find parameters that best meet our goals belong to optimization. It is itself a topic of study and research.\n",
    "\n",
    "Gradient descent is an entry-level optimization method. The first part of this lesson will introduce stochastic gradient descent, an optimization method derived from vanilla gradient descent. In the second half of this lesson, we will talk about using the optimizers from SciPy. SciPy is a well-known Python package for scientific\n",
    "computing. It has several advanced optimizers that may be too complicated for us to write from scratch. Using third-party packages like SciPy can save us a lot of time.\n",
    "\n",
    "Our goal is to give you a taste of different optimizations and how to use third-party packages. We will refrain from digging too deep into the math.\n",
    "\n",
    "##### Warm-up exercise\n",
    "\n",
    "Consider the following problem: *a to-be-built rectangular fence will enclose $120$ square feet. The material cost for the long and short sides will be $\\$5$ and $\\$6$ per foot, respectively. What are the dimensions of the fence that will minimize cost?*\n",
    "\n",
    "Let the lengths of the long and short sides be $l$ and $s$. This is an optimization problem as we want to find a set of $l$ and $s$ that gives us the best cost $c = 2\\left(5l+6s\\right)$. Try to solve this problem with gradient descent. The solutions are $l=12$ and $s=10$. (Hint: substitute $s = 120l^{-1}$ into the cost first.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2f2c099",
   "metadata": {},
   "source": [
    "### Before we start"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb69a275",
   "metadata": {},
   "source": [
    "We again use the application of identifying defective metal-casting parts. Let's import the required packages and normalized data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b6ca0b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from autograd import numpy\n",
    "from autograd import grad\n",
    "from matplotlib import pyplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c61563d",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"The following lines imports variables and functions from lesson 5 and 7. The\n",
    "imported variables and functions include : images_train, images_val, images_test,\n",
    "labels_train, labels_val, labels_test, neural_network_model, model_loss,\n",
    "regularized_loss, and performance.\"\"\"\n",
    "import sys\n",
    "sys.path.insert(0, \"../scripts\")\n",
    "from load_casting_data import *\n",
    "from lesson_7_functions import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a22e81b",
   "metadata": {},
   "source": [
    "In addition, we define a function to return a new set of parameters with Xavier initialization. Every time we train the model with a new optimization method, we want to re-initialize the parameters to the same values to have a fair comparison between different methods. Note that we use a length-one NumPy array for $b^1$ (which is supposed to be a scalar) because then we can use many NumPy functions on it later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de891498",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_new_parameters(n0, n1):\n",
    "    \"\"\"Get a new set of parameters using a fixed seed and the Xavier initialization.\n",
    "    \n",
    "    Arguments\n",
    "    ---------\n",
    "    n0, n1 : int\n",
    "        The number of elements in vector x and the first intermediate vector z1.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    A list of 4 elementsb\n",
    "        W0 : a 2D array of shape (n0, n1).\n",
    "        b0 : a 1D array of length n1.\n",
    "        W1 : a 1D array of length n1.\n",
    "        b1 : a 1D array of only one element (i.e., to mimic a scalar).\n",
    "    \"\"\"\n",
    "    numpy.random.seed(0)  # fix the random seed to avoid \"true randomness\"\n",
    "    \n",
    "    scale = numpy.sqrt(6/(n0+n1))\n",
    "    W0 = numpy.random.uniform(-scale, scale, (n0, n1))\n",
    "    b0 = numpy.zeros(n1)\n",
    "\n",
    "    scale = numpy.sqrt(6/(n1+1))\n",
    "    W1 = numpy.random.uniform(-scale, scale, n1)\n",
    "    b1 = numpy.zeros(1)  # using a length-one vector as a scalar\n",
    "    \n",
    "    return [W0, b0, W1, b1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "636db73b",
   "metadata": {},
   "source": [
    "The lengths of each input vector $\\mathbf{x}$ and the intermediate vector $\\mathbf{z}^1$ are fixed throughout this lesson."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e6002dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "n0 = images_train.shape[1]  # i.e., total num. of pixels per image\n",
    "n1 = 64  # the length of z1 vector"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2926e745",
   "metadata": {},
   "source": [
    "Don't forget the function that calculates the gradients for us: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6fb9e003",
   "metadata": {},
   "outputs": [],
   "source": [
    "gradients = grad(regularized_loss, argnum=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb6671ef",
   "metadata": {},
   "source": [
    "## Revisiting gradient descent"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f96c0067",
   "metadata": {},
   "source": [
    "In gradient descent, the gradients are the derivatives of a total loss $L$ with respect to parameters. The total loss can be the sum or the mean of the losses from all training samples. For example, we use the mean when doing logistic regression over $N$ training samples:\n",
    "\n",
    "$$\n",
    "L\\left(\\hat{\\mathbf{y}}, \\mathbf{y}\\right)\n",
    "=\n",
    "\\frac{1}{N} \\sum_{i=1}^{N} l\\left(\\hat{y}_i, y_i\\right)\n",
    "=\n",
    "\\frac{1}{N}\\sum_{i=1}^{N} −y_i\\log(\\hat{y}_i)-(1-y_i)\\log(1-\\hat{y}_i)\n",
    "$$\n",
    "\n",
    "The gradients, in this case, are also the mean of the gradients calculated from the single-sample loss, $l\\left(\\hat{y}_i, y_i\\right)$:\n",
    "\n",
    "$$\n",
    "\\frac{\\mathrm{d}~L}{\\mathrm{d}~W}\\left(\\hat{\\mathbf{y}}, \\mathbf{y}\\right) = \\frac{1}{N}\\sum_{i=1}^{N}\\frac{\\mathrm{d}~l}{\\mathrm{d}~W}\\left(\\hat{y}_i, y_i\\right)\n",
    "$$\n",
    "\n",
    "And the gradient descent algorithm is roughly\n",
    "\n",
    "```\n",
    "loop:\n",
    "    Y_hat 🠘 model(X, W)\n",
    "    L 🠘 loss(Y_hat, Y)\n",
    "    W 🠘 W - step_size x (dL/dW)\n",
    "    exit optimization if model performance is satisfying\n",
    "```\n",
    "\n",
    "(Note we never implemented code to calculate `dL/dW` (i.e., $\\frac{1}{N}\\sum_{i=1}^{N}\\frac{\\mathrm{d}~l}{\\mathrm{d}~W}\\left(y_i, \\hat{y}_i\\right)$) because `autograd` did that for us.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec92537a",
   "metadata": {},
   "source": [
    "## Batched and stochastic gradient descent"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e1e7510",
   "metadata": {},
   "source": [
    "When we have more samples in the training dataset, an iteration of gradient descent becomes more computationally expensive, i.e., it needs more time and computer memory. This is when ***batched gradient descent*** comes into play."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8419ea89",
   "metadata": {},
   "source": [
    "Batched gradient descent divides an entire training dataset into several batches. Next, it assumes the gradients calculated from the loss of each batch are about the same as the overall gradients:\n",
    "\n",
    "$$\n",
    "\\frac{\\mathrm{d}~L}{\\mathrm{d}~W}\\left(\\hat{\\mathbf{y}}, \\mathbf{y}\\right)\n",
    "= \\frac{1}{N}\\sum_{i=1}^{N}\\frac{\\mathrm{d}~l}{\\mathrm{d}~W}\\left(\\hat{y}_i, y_i\\right)\n",
    "\\approx \\frac{1}{N_b}\\sum_{i=1}^{N_b}\\frac{\\mathrm{d}~l}{\\mathrm{d}~W}\\left(\\hat{y}_i, y_i\\right)\n",
    "\\approx \\frac{1}{N_b}\\sum_{i=N_b+1}^{2N_b}\\frac{\\mathrm{d}~l}{\\mathrm{d}~W}\\left(\\hat{y}_i, y_i\n",
    "\\right)\n",
    "\\approx \\cdots\n",
    "$$\n",
    "\n",
    "$N_b$ denotes the number of samples in each batch. Batches don't need to have the same number of samples. For example, the last batch can be smaller because $N$ may not be divisible by $N_b$.\n",
    "\n",
    "Some people use batched gradient descent and ***stochastic gradient descent*** interchangeably. Some people, on the other hand, specifically use stochastic gradient descent for the case $N_b=1$. \n",
    "\n",
    "The data distribution should be similar in each batch to make batch gradients close enough to overall gradients. For example, let's say we have 800 good casting parts and 200 defective casting parts in our training dataset. In that case, each batch should also have approximately 80%/20% of good/defective casting parts. We can achieve this by ***shuffling*** the training data before dividing them into batches."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d39f6e42",
   "metadata": {},
   "outputs": [],
   "source": [
    "# fix the random seed to avoid \"true randomness\"\n",
    "numpy.random.seed(12345)\n",
    "\n",
    "# create an array for original indices\n",
    "order = numpy.arange(len(images_train))\n",
    "\n",
    "# shuffle the indices (`order` is re-ordered in-place)\n",
    "numpy.random.shuffle(order)\n",
    "\n",
    "# get re-ordered training data\n",
    "images_train = images_train[order]\n",
    "labels_train = labels_train[order]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3171d76d",
   "metadata": {},
   "source": [
    "Note that in the above code cell, the function `numpy.random.shuffle` can shuffle `images_train` directly. However, we first created `new_order` instead because we need to re-order `labels_train` with the same order as `images_train`.\n",
    "\n",
    "We use $N_b=100$ here. The total number of samples in the training dataset is $782$, so we have $8$ batches (with the last batch only having $82$ samples)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b42d46d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the batch size Nb\n",
    "Nb = 100\n",
    "\n",
    "# define the start and end index of each batch\n",
    "starts = [0, 100, 200, 300, 400, 500, 600, 700]  # python starts from 0\n",
    "ends = [100, 200, 300, 400, 500, 600, 700, 782]  # python uses end+1 for slicing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bacc31f0",
   "metadata": {},
   "source": [
    "Let's examine if each batch has a similar ratio between good and defective casting parts:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35161eda",
   "metadata": {},
   "outputs": [],
   "source": [
    "# the overall good-defective ratio in the training dataset\n",
    "Ndef = numpy.count_nonzero(labels_train)\n",
    "print(\"The overall good-defective ratio in images_train: {}\".format(\n",
    "    (len(images_train)-Ndef)/Ndef))\n",
    "\n",
    "for i, (bg, ed) in enumerate(zip(starts, ends)):\n",
    "    # calculate and print the ratio\n",
    "    Ndef = numpy.count_nonzero(labels_train[bg:ed])\n",
    "    print(\"Batch {} (samples {} to {}): the ratio is {}\".format(\n",
    "        i, bg, ed-1, (ed-bg-Ndef)/Ndef))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67affb21",
   "metadata": {},
   "source": [
    "There's no golden rule how close these numbers should be. We assume the good-defective ratios from the batches are close enough to the overall one for the sake of convenience."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf295565",
   "metadata": {},
   "source": [
    "When updating the parameters, we loop through batches, calculate the loss of the current batch, get batch gradients, and then update the parameters using the batch gradients:\n",
    "\n",
    "```\n",
    "loop:\n",
    "    loop i = [0, 1, ..., number_of_batches-1]:\n",
    "        X_b 🠘 training data in the i-th batch\n",
    "        Y_b 🠘 labels in the i-th batch\n",
    "        Y_hat_b 🠘 model(X_b, W)\n",
    "        L_b 🠘 loss(Y_hat_b, Y_b)\n",
    "        W 🠘 W - step_size x (dL_b/dW)\n",
    "        exit optimization if model performance is satisfying\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bfc5b329",
   "metadata": {},
   "source": [
    "The `L_b` in the pseudo-code denotes the loss of `i`-th batch, $L_b = \\frac{1}{N_b}\\sum_{j=iN_b}^{(i+1)N_b}l\\left(\\hat{y}_j, y_j\\right)$, and `dL_b/dW` denotes its gradients."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e840bc88",
   "metadata": {},
   "source": [
    "Let's try the new method to optimize our neural network model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "038a8169",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "# get a new set of parameters (params=[W0, b0, W1, b1])\n",
    "params = get_new_parameters(n0, n1)\n",
    "\n",
    "# step size\n",
    "lr = 1e-1\n",
    "\n",
    "# variables to store parameters at when the validation loss is at its lowest\n",
    "best_val_loss = numpy.inf\n",
    "best_iter = 0\n",
    "best_params = None\n",
    "\n",
    "# variables to store the history of training and validation loss\n",
    "hist_train_loss = []\n",
    "hist_val_loss = []\n",
    "\n",
    "# a counter for the numbero of total iterations\n",
    "i = 0\n",
    "\n",
    "# limit the total number of iterations to 2000\n",
    "while i < 2000:\n",
    "    \n",
    "    # loop through each batch\n",
    "    for j, (bg, ed) in enumerate(zip(starts, ends)):\n",
    "        \n",
    "        # alias for data in this batch to shorten the code\n",
    "        x_b = images_train[bg:ed]\n",
    "        labels_b = labels_train[bg:ed]\n",
    "        \n",
    "        # calculate gradients and use gradient descents\n",
    "        grads = gradients(x_b, labels_b, params)\n",
    "        params = [p - lr * g for p, g in zip(params, grads)]\n",
    "    \n",
    "        # save losses for future investigation\n",
    "        hist_train_loss.append(model_loss(x_b, labels_b, params))\n",
    "        hist_val_loss.append(model_loss(images_val, labels_val, params))\n",
    "    \n",
    "        # if the validation loss is lower than the current best, save it\n",
    "        if hist_val_loss[-1] < best_val_loss:\n",
    "            best_val_loss = hist_val_loss[-1]\n",
    "            best_iter = i\n",
    "            best_params = params\n",
    "    \n",
    "        # update iteration counter\n",
    "        i = i + 1\n",
    "        \n",
    "        # print something every 100 iterations so we know it's not dead\n",
    "        if i % 100 == 0:\n",
    "            print(\"{}...\".format(i), end=\"\")\n",
    "\n",
    "# just to add a new line\n",
    "print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e68a769",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"The best validation loss:\", best_val_loss)\n",
    "print(\"It happened at iteration {}\".format(best_iter))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e729593e",
   "metadata": {},
   "source": [
    "The optimization took about $100$ seconds to finish $2000$ iterations on a personal desktop with 6 physical cores. Compared to the $2.5$ minutes in lesson 7, it is about $1.5\\texttt{x}$ faster. Moreover, we now get the best parameters earlier at the $390$th iteration. If we monitor and stop the optimization in real-time, it means the time-to-solution is $19.5$ seconds. ($100 \\times 390 \\div 2000 = 19.5$). The time-to-solution of the vanilla gradient descent in lesson 7 is $34.2$ seconds ($2.5 \\times 60 \\times 456 \\div 2000 = 34.2$). So using batched gradient descent improves the time-to-solution with a $1.75\\texttt{x}$ speedup!\n",
    "\n",
    "Just so that you know: it's also possible to further improve the performance by dividing validation data into batches."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e87a110",
   "metadata": {},
   "source": [
    "We examine the loss history in the following figures:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50a27396",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt, ax = pyplot.subplots(1, 2, figsize=(10, 5))\n",
    "ax[0].plot(range(len(hist_train_loss)), hist_train_loss, label=\"Training loss\")\n",
    "ax[0].plot(range(len(hist_val_loss)), hist_val_loss, label=\"Validation loss\")\n",
    "ax[0].set_xlabel(\"Iteration\")\n",
    "ax[0].set_ylabel(\"Loss\")\n",
    "ax[0].set_title(\"Iteration 0 to 1999\")\n",
    "ax[0].legend(loc=0)\n",
    "ax[1].plot(range(100), hist_train_loss[:100], label=\"Training loss\")\n",
    "ax[1].plot(range(100), hist_val_loss[:100], label=\"Validation loss\")\n",
    "ax[1].set_xlabel(\"Iteration\")\n",
    "ax[1].set_ylabel(\"Loss\")\n",
    "ax[1].set_title(\"Iteration 0 to 100\")\n",
    "ax[1].legend(loc=0)\n",
    "plt.suptitle(\"Training and validation loss history\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "358af9c8",
   "metadata": {},
   "source": [
    "We have the complete loss history on the left. Oscillation happened from the beginning to the end. To have a clearer view of the oscillation, we plotted the first 100 iterations on the right. This oscillation is normal for batched gradient descent. In a nutshell, though we assume the characteristics of a batch (distributions, means, gradients, losses, etc.) are about the same as those of the whole training dataset, they are still different more or less. So using one batch at a time usually gives a bumpier loss history than using the entire dataset.\n",
    "\n",
    "The smaller the $N_b$ is, the more oscillation we get. The strictly stochastic gradient descent ($N_b=1$) hence usually shows a considerable fluctuation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0fa1f24",
   "metadata": {},
   "source": [
    "The model's final performance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a8d74cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# final accuracy\n",
    "pred_labels_test = classify(images_test, best_params, neural_network_model)\n",
    "perf = performance(pred_labels_test, labels_test)\n",
    "print(\"Final precision: {:.1f}%\".format(perf[0]*100))\n",
    "print(\"Final recall: {:.1f}%\".format(perf[1]*100))\n",
    "print(\"Final F-score: {:.1f}%\".format(perf[2]*100))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fc01a53",
   "metadata": {},
   "source": [
    "The performance is slightly different from what we got in lesson 7. Many factors contribute to the difference. The most obvious reason is that we were not using the best hyperparameters in both lessons."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c979e70",
   "metadata": {},
   "source": [
    "Though there are still many other gradient-descent-based optimization methods, we stop the discussion of gradient descent here. Among so many different optimization methods, gradient descent and its variants are the most important ones in modern deep learning. They require fewer computing resources, so optimizing models against millions or even billions of data becomes feasible. We encourage those who want to go deeper with deep learning to check reference [1]."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "419d711a",
   "metadata": {},
   "source": [
    "## Using advanced optimizers from SciPy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1496c17",
   "metadata": {},
   "source": [
    "We have been implementing gradient descent from scratch in all lessons so far. Gradient descent is easy to implement and powerful in many use cases.  However, sometimes we may want to use more advanced methods, but implementing advanced methods may be challenging. For example, [optimizers based on Newton's method](https://en.wikipedia.org/wiki/Quasi-Newton_method) are commonly used for small-size problems but are complicated. It may be more feasible to use third-party packages in such a situation instead of writing the code from scratch.\n",
    "\n",
    "We want to show you how to use the optimizers from a well-known Python package called SciPy. We will use the [L-BFGS-B optimizer](https://en.wikipedia.org/wiki/Limited-memory_BFGS) to optimize our neural network. L-BFGS-B is based on Newton's method. The primary goal is to show how to use SciPy, and L-BFGS-B serves as an example, so we'll skip the math.\n",
    "\n",
    "SciPy provides a uniform interface for all optimization methods through [`scipy.optimize.minimize`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html). In this lesson, we will use this function with\n",
    "\n",
    "```\n",
    "results = scipy.optimize.minimize(\n",
    "    fun      = <a function returning loss and gradients>,\n",
    "    x0       = <initial values for parameters>,\n",
    "    args     = <arguments to `fun`>,\n",
    "    method   = <optimization method>,\n",
    "    jac      = <whether `fun` returns gradients,\n",
    "    callback = <a function to call after every iteration>\n",
    ")\n",
    "```\n",
    "\n",
    "Let's take a look at the arguments:\n",
    "\n",
    "##### 1. $\\texttt{fun}$\n",
    "   \n",
    "`fun` is a function that returns both the loss and gradients. This function must have the model parameters as its first argument. All other arguments should be after model parameters. For example, we usually use `params` to represent model parameters in our lessons and use `x` and `y` to represent input data and labels. So if the function name is `loss_and_grad`, then the definition of the function should be something like:\n",
    "```python\n",
    "def loss_and_grad(params, x, y):\n",
    "    ...\n",
    "    return loss, grads\n",
    "```\n",
    "\n",
    "And then `fun=loss_and_grad`.\n",
    "\n",
    "SciPy expects all parameters to be packed into a flattened 1D array. Previously, our `params` is a `list` of several `numpy.ndarray`. Now we need to flatten each `numpy.ndarray` and concatenate them into a single vector.\n",
    "\n",
    "##### 2. $\\texttt{x0}$\n",
    "\n",
    "SciPy uses different notation. It uses `x` to denote model parameters while we use `params`. `x0` represents the initial values of model parameters.\n",
    "\n",
    "##### 3. $\\texttt{args}$\n",
    "\n",
    "All arguments except for the model parameters should be provided through `args`. For example, `loss_and_grad` requires three arguments: model parameters, input data, and labels. During training, input data and labels are `images_train` and `labels_train`, respectively. So `args=(images_train, labels_train)`.\n",
    "\n",
    "##### 4. $\\texttt{method}$\n",
    "\n",
    "This is a string indicating what optimization method we want to use. [This SciPy tutorial](https://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html#unconstrained-minimization-of-multivariate-scalar-functions-minimize) lists the methods allowed in SciPy. Here we use `mothod='L-BFGS-B'`.\n",
    "   \n",
    "##### 5. $\\texttt{jac}$\n",
    "\n",
    "We use `jac=True` to let SciPy know the function provided in `fun` also returns gradients. (`fun` may return the loss only in some situations, and SciPy will then use other approaches to calculate gradients.) Since we already use `autograd` to get gradients, we don't need SciPy to calculate the gradients for us.\n",
    "\n",
    "##### 6. $\\texttt{callback}$\n",
    "\n",
    "If we provide a function to `callback`, SciPy will call this function after each optimization iteration. The `callback` function takes only one argument: the model parameters at the current optimization iteration. We can then do things like calculating the current validation loss or printing iteration information."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfdaaa8b",
   "metadata": {},
   "source": [
    "Now let's see how to use SciPy. First, all optimizers are under submodule `scipy.optimize`. We import it first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98941e50",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.optimize"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bbeed1e",
   "metadata": {},
   "source": [
    "As mentioned, SciPy needs all model parameters to be flattened and packed into a 1D array. We define `join_parameters` to achieve it and `split_parameters` to recover the 1D array back to a list of `numpy.ndarray`s. `split_parameters` takes in three arguments: model parameters, input data sample size, and the length of the intermediate vector. We need the sizes to know the shapes of $W^0$, $\\mathbf{b}^0$, $W^1$, and $\\mathbf{b}^1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f14a5e2d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def join_parameters(params):\n",
    "    \"\"\"Flatten and concatenate several ndarrays into one single vector.\n",
    "    \n",
    "    Arguments\n",
    "    ---------\n",
    "    params : a list of numpy.ndarray\n",
    "        For example, in our neural network, params=[W0, b0, W1, b1].\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    params : a 1D numpy.ndarray\n",
    "    \"\"\"\n",
    "    params = numpy.concatenate([p.flatten() for p in params])\n",
    "    return params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dec20312",
   "metadata": {},
   "outputs": [],
   "source": [
    "def split_parameters(params, n0, n1):\n",
    "    \"\"\"Split an 1D array or parameters into W0, b0, W1, b1.\n",
    "    \n",
    "    Arguments\n",
    "    ---------\n",
    "    params : 1D numpy.ndarray\n",
    "        The parameters of a neural network packed in an 1D array.\n",
    "    n0 : integer\n",
    "        The number of element in the input vector x.\n",
    "    n1 : integer\n",
    "        The number of element in the intermediate vector z1.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    A list of 4 elements\n",
    "        W0 : a 2D array of shape (n0, n1).\n",
    "        b0 : a 1D array of length n1.\n",
    "        W1 : a 1D array of length n1.\n",
    "        b1 : a scalar.\n",
    "    \"\"\"\n",
    "    \n",
    "    W0 = params[:n0*n1].reshape((n0, n1))\n",
    "    b0 = params[n0*n1:(n0+1)*n1].reshape((n1,))\n",
    "    W1 = params[(n0+1)*n1:(n0+2)*n1].reshape((n1,))\n",
    "    b1 = params[-1]\n",
    "    return [W0, b0, W1, b1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f8e68f6",
   "metadata": {},
   "source": [
    "Next is what we will provide to `fun`. It is a function taking five arguments: model parameters `params`, input data `x`, labels `y`, size of each input data sample `n0`, and the length of intermediate vector `n1`. The function returns both the loss and gradients. Gradients also need to be flattened and packed into a 1D array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8cabaa27",
   "metadata": {},
   "outputs": [],
   "source": [
    "def loss_and_grad(params, x, y, n0, n1):\n",
    "    \"\"\"A function that returns loss and gradients for using SciPy's optimizers.\n",
    "    \n",
    "    Arguments\n",
    "    ---------\n",
    "    params : 1D numpy.ndarray\n",
    "        The parameters of a neural network packed in an 1D array.\n",
    "    x : numpy.ndarray\n",
    "        The input data. If a 1D array, it represents the features of one sample. If a\n",
    "        2D array, the first dimension represents the number of samples, and the second\n",
    "        dimension represents the features.\n",
    "    y : numpy.ndarray\n",
    "        The output of the model. If x is a 1D array, y is a scalar. If x is a 2D array,\n",
    "        y is a 1D array.\n",
    "    n0 : integer\n",
    "        The number of element in the input vector x.\n",
    "    n1 : integer\n",
    "        The number of element in the intermediate vector z1.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    loss : a scalar\n",
    "        The total loss (regularized).\n",
    "    grads : a 1D numpy.ndarray\n",
    "        The gradients of w.r.t. parameters packed into a single vector.\n",
    "    \"\"\"\n",
    "    # split and recover the shape of parameters\n",
    "    params = split_parameters(params, n0, n1)\n",
    "    \n",
    "    # get the regularized loss (params now is [W0, b0, W1, b1])\n",
    "    loss = regularized_loss(x, y, params)\n",
    "    \n",
    "    # get the gradients (grads is now [dW0, db0, dW1, db1])\n",
    "    grads = gradients(x, y, params)\n",
    "    \n",
    "    # concatenate all gradients into a 1D array\n",
    "    grads = join_parameters(grads)\n",
    "    \n",
    "    return loss, grads"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a3dacd7",
   "metadata": {},
   "source": [
    "We want to calculate the validation loss at the end of every iteration and record the best model parameters (the ones that give the lowest validation loss). We define the function `best_param_recorder` for this purpose:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "113e00e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def best_param_recorder(params):\n",
    "    \"\"\"A callback function to record the parameters that give the best validation loss.\n",
    "    \n",
    "    Arguments\n",
    "    ---------\n",
    "    params : 1D numpy.ndarray\n",
    "        The parameters of a neural network packed in an 1D array.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    A boolean. This function always return False to let SciPy optimizers know they can\n",
    "    continue the optimization.\n",
    "    \n",
    "    Notes\n",
    "    -----\n",
    "    Variables `best_val_loss` and `best_params` should already exist outside this\n",
    "    function when this function is called.\n",
    "    \"\"\"\n",
    "    \n",
    "    # `global` indicates these variables will be the same as those outside the function\n",
    "    global best_val_loss, best_params\n",
    "    \n",
    "    # split and recover the shape of parameters\n",
    "    params = split_parameters(params, n0, n1)\n",
    "\n",
    "    # calculate the validation loss\n",
    "    val_loss = model_loss(images_val, labels_val, params)\n",
    "\n",
    "    # if the new validation loss is smaller than the current best, save it\n",
    "    if val_loss <= best_val_loss:\n",
    "        best_val_loss = val_loss\n",
    "        best_params = params.copy()\n",
    "\n",
    "    return False"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6d56f39",
   "metadata": {},
   "source": [
    "`best_param_recorder` takes in only model parameters and can only return a boolean. This limitation comes from SciPy. So in order to record the best parameters, we use a keyword `global` in the function. `global` notifies the function that the subsequent variables will be obtained and synchronized with the same variables outside the function, i.e., in the global scope. In this way, we don't have to return the best parameters using `return` but are still able to access them from outside the function. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d0c3d00",
   "metadata": {},
   "source": [
    "In other words, we need to decalre and initialize those global variables first before using the function `best_param_recorder`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df951612",
   "metadata": {},
   "outputs": [],
   "source": [
    "# we need these variables to be decalred outside the function `best_param_recorder`\n",
    "best_val_loss = numpy.inf\n",
    "best_params = None"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38f06fd9",
   "metadata": {},
   "source": [
    "Finally, we can use `scipy.optimize.minimize` to optimize the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce58f241",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "# get a new set of parameters\n",
    "params = get_new_parameters(n0, n1)  # params=[W0, b0, W1, b1]\n",
    "params = join_parameters(params)  # flattened & concatenated\n",
    "\n",
    "# optimize the model using SciPy\n",
    "results = scipy.optimize.minimize(\n",
    "    loss_and_grad, params, args=(images_train, labels_train, n0, n1),\n",
    "    method='L-BFGS-B', jac=True, tol=1e-9, callback=best_param_recorder)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "251d826b",
   "metadata": {},
   "source": [
    "The returned object `results` contains several information. We can access them through a dot `.`. For example, the number of iterations is `results.nit`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d381a3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Number of iterations: {}\".format(results.nit))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eed5afa7",
   "metadata": {},
   "source": [
    "And the final model parameters are in `results.x` (recall that SciPy uses `x` to denote model parameters). [This webpage](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.OptimizeResult.html#scipy.optimize.OptimizeResult) shows other information carried by `results`. Let's see the performance against the test dataset using the model parameters from the last iteration:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "346a3d36",
   "metadata": {},
   "outputs": [],
   "source": [
    "# final accuracy using the parameters from the last iteration\n",
    "params = split_parameters(results.x, n0, n1)\n",
    "pred_labels_test = classify(images_test, params, neural_network_model)\n",
    "perf = performance(pred_labels_test, labels_test)\n",
    "print(\"Final precision: {:.1f}%\".format(perf[0]*100))\n",
    "print(\"Final recall: {:.1f}%\".format(perf[1]*100))\n",
    "print(\"Final F-score: {:.1f}%\".format(perf[2]*100))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f877bbfd",
   "metadata": {},
   "source": [
    "It doesn't look as nice as what we got previously using gradient descent and the batched version. However, remember these numbers are obtained using the model parameters from the last iteration. They are not the parameters that give the lowest validation. Let's calculate the test performance again using the best parameters (i.e., `best_params`):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16617ec6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# final accuracy using the best parameters\n",
    "pred_labels_test = classify(images_test, best_params, neural_network_model)\n",
    "perf = performance(pred_labels_test, labels_test)\n",
    "print(\"Final precision: {:.1f}%\".format(perf[0]*100))\n",
    "print(\"Final recall: {:.1f}%\".format(perf[1]*100))\n",
    "print(\"Final F-score: {:.1f}%\".format(perf[2]*100))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae8ea7c3",
   "metadata": {},
   "source": [
    "These numbers are different from what we got from gradient descent. It's possible to tune the L-BFGS-B optimizer to get closer numbers. We can find all tunable options for L-BFGS-B optimizer [here](https://docs.scipy.org/doc/scipy/reference/optimize.minimize-lbfgsb.html). We will stop here as our primary goal is to show you how to use SciPy's optimizers."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b500410e",
   "metadata": {},
   "source": [
    "## Challenge"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "335c1c7c",
   "metadata": {},
   "source": [
    "Recall the 120-square-feet fence we wanted to build at the beginning of this notebook. Can you write code to solve the same problem using SciPy's optimizers?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "084352c0",
   "metadata": {},
   "source": [
    "## What we've learned"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f95e6918",
   "metadata": {},
   "source": [
    "1. Using batched gradient descent to accelerate the optimization.\n",
    "2. Using third-party libaraies, like SciPy, to optimize a model with advanced optimizers."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "772c6efd",
   "metadata": {},
   "source": [
    "## Reference\n",
    "1. Goodfellow, I., Bengio, Y., & Courville, A. (2016). Chapter 8. Optimization for Training Deep Models. In Deep Learning (pp. 271–325). MIT Press. http://www.deeplearningbook.org"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0be9b20d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute this cell to load the notebook's style sheet, then ignore it\n",
    "from IPython.core.display import HTML\n",
    "css_file = '../style/custom.css'\n",
    "HTML(open(css_file, \"r\").read())"
   ]
  }
 ],
 "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.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}