02_basics_algorithms.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Basics, algorithms\n",
    "\n",
    "Welcome to the first jupyter notebook! In this session, we won't go into too many minute details yet, but will cover the basics of machine learning. We've tried to keep it as simple as possible, also because many might only have little experience with python and/or programming a 'learning machine'. If it looks like you're gonna be through the content of this notebook in ten minutes or so, because you're already familiar with all of its concepts, then feel free to challenge yourself a little more by making your own adjustments to the code, diving deeper into advanced topics or doing some last-minute preparations for your group presentations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup\n",
    "\n",
    "To allow the next code blocks to run smoothly, this section sets a couple of settings."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some imports that we will be using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the random seed to a fixed number. This will guarantee that the notebook output is generated the same way for every run, otherwise the seed would be – random, as the name suggests."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the numpy output precision for some prettier numbers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.set_printoptions(precision=3, suppress=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some figure plotting settings: increase the axis labels of our figures a bit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpl.rc('axes', labelsize=14)\n",
    "mpl.rc('xtick', labelsize=12)\n",
    "mpl.rc('ytick', labelsize=12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear regression with the normal equation\n",
    "\n",
    "Let's start with something simple: linear regression. As we've learnt before, we can use the _normal equation_ to calculate a prediction. In statistics, we usually label this as $\\hat{\\theta}$, because it is an estimator for the parameter vector $\\theta$ of the model. The hat indicates that it's an estimator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " First we need to generate some random data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 100   # number of data points\n",
    "X = 2 * np.random.rand(m, 1)\n",
    "y = 3 + 4 * X + np.random.randn(m, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And let's plot it to get an idea what we're looking at. Of course, we do things the proper way, and put labels on our axes!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(X, y, \".\")\n",
    "plt.xlabel(\"$x_1$\", fontsize=18)\n",
    "plt.ylabel(\"$y$\", rotation=0, fontsize=18)\n",
    "plt.axis([0, 2, 0, 15])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also print a few of the generated data values, just to get an idea what we're talking about. The following command shows us the first three entries of the object `X` that we just generated:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X[0:3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And this shows us the first three corresponding entries in the object `y`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y[0:3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks good. Reminder: what is the normal equation? And what are `X` and `y`?\n",
    "\n",
    "$$ \\hat{\\theta} = \\left( \\mathbf{X}^T \\cdot \\mathbf{X} \\right)^{-1} \\cdot \\mathbf{X}^T \\cdot \\mathbf{y} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quick refresher if your mathematics is a littley rusty:\n",
    "* $\\hat{\\theta}$ is our estimator for the vector of parameter $\\theta$. This is what we want to calculate!\n",
    "* $\\mathbf{x}^{(i)}$, beware that it's lower-case, is a vector which contains all features of the training instance with index i. In the data generated above, we only have one 'feature', which is called $x_1$. So, $\\mathbf{x}^{(0)}$, the feature vector of instance number zero, includes one single value: the value of $x_1$ for that instance, as printed out with the commands above.\n",
    "* $\\mathbf{X}$, now it's upper-case, is a vector of _all feature vectors_ $\\mathbf{x}^{(i)}$. To make things more confusing, the entries of $\\mathbf{X}$ are actually not $\\mathbf{x}^{(i)}$, but the transposed vectors $(\\mathbf{x}^{(i)})^T$.\n",
    "* $\\mathbf{y}$ are the _true_ target values of the instances. So, this is also a vector with the same dimension as $\\mathbf{X}$. But even in more complicated input data structures, every entry of the target vector will just be the one target value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, what can we do with the normal equation? And what _is_ actually this $\\theta$? It's our vector of model parameters. The above case is very simple: we would like to create a model that represents the data as well as possible. With just looking at the plot, and without too much thinking, it's obvious that there is some sort of linear dependence between $x_1$ and $y$.\n",
    "\n",
    "How many parameters do we need to describe this model? Probably two: one for the linear dependence, and one _bias term_ to shift the entire model along the y axis. So, our $\\theta$ in this case is just a vector of two entries, and the goal of 'linear regression' is to find the optimal values of the two. And actually, if you scroll back up for a bit, you can already see the linear dependence, when we generated the data:\n",
    "\n",
    "$$ y(x_1) = 3 + 4 \\cdot x_1 $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But let's pretend we don't know that ... Without using any machine learning yet, we can just use the above normal equation to get estimators for the two values. For that, we can make use of numpy's `linalg.inv()` function to invert matrices. Essentially, we then just need to 'type' the above formula into python and let our computer do the rest. Easy, right?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One more step is necessary: we need to append an additional feature $x_0 = 1$ to all instances, because otherwise we would ignore the bias parameter in our calculation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_b = np.c_[np.ones((m, 1)), X]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cool. Now, here's the typed-out formula for calculating the normal equation. It only uses numpy functions, such as the matrix inversion, or the calculation of dot products:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do we know it worked? One easy thing is to check what the values of the two parameters are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"theta_0 = %.3f\" % theta[0][0])\n",
    "print(\"theta_1 = %.3f\" % theta[1][0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, that's not perfect, but close. Where does the difference come from? Probably from the fact that our data 'only' consists of 100 data points. the more data points we use, the closer we get to the values that we used to generate the data with."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maybe it is also useful to plot the prediction as a line into the plot. For that, we should first calculate the predictions for the value of `y` for all our instances:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_predict = X_b.dot(theta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now we can do the plotting:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(X, y, \".\")\n",
    "plt.plot(X, y_predict, \"-\", linewidth=2, label=\"Predictions\")\n",
    "plt.xlabel(\"$x_1$\", fontsize=18)\n",
    "plt.ylabel(\"$y$\", rotation=0, fontsize=18)\n",
    "plt.legend(loc=\"upper right\", fontsize=14)\n",
    "plt.axis([0, 2, 0, 15])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Great! As a quick summary: we generated random data in python, determined an appropriate model to represent that data (by eye-balling very carefully), and used the normal equation to get estimators for our model parameters. Now, let's start with some actual machine learning."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear regression using batch gradient descent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try and implement the first machine-learning algorithm to solve our linear-regression problem: batch gradient descent. Quick reminder: gradient descent is an iterative approach to find $\\hat{\\theta}$. Using the learning rate $\\eta$, we adjust our estimates for $\\theta$ in each learning step iteratively. The \"direction\" of adjustment is determined by the _gradient_ of the mean square error."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maybe we should have a quick revision of the formula:\n",
    "\n",
    "$$ \\mathit{MSE}(\\theta) = \\frac{1}{m} \\sum_{i=1}^{m} \\left( \\theta^T \\cdot \\mathbf{x}^{(i)} - y^{(i)} \\right)^2 $$\n",
    "\n",
    "Now, most of you will probably know that the gradient of this function just means taking the derivative of it with respect to $\\theta_1,\\dots,\\theta_n$. To refresh your memory, let's write down the formula for the partial derivative as well:\n",
    "\n",
    "$$ \\frac{\\partial}{\\partial \\theta_j} \\mathit{MSE}(\\theta) = \\frac{2}{m} \\sum_{i=1}^{m} \\left( \\theta^T \\cdot \\mathbf{x}^{(i)} - y^{(i)} \\right) x_j^{(i)}$$ \n",
    "\n",
    "Then, the entire gradient is:\n",
    "\n",
    "$$ \\nabla_\\theta \\mathit{MSE}(\\theta) = \\frac{2}{m} \\mathbf{X}^T \\cdot \\left( \\mathbf{X} \\cdot \\theta - \\mathbf{y} \\right) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now it's really just one last step missing: we need to calculate our predictions for $\\theta$. For the very first step, we start with random values of $\\theta$. Then, after calculating the gradient above for a step, we update the value of $\\theta$ according to:\n",
    "\n",
    "$$ \\theta \\rightarrow \\theta - \\eta \\nabla_\\theta \\mathit{MSE}(\\theta) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That wasn't too hard, was it? Writing this out with python is even easier. Let's start with setting a learning rate $\\eta$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eta = 0.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we also need to decide how many steps of calculations we would like to perform:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_iterations = 1000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And initialise our $\\theta$ with random values: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta = np.random.randn(2,1)   # let's start with random values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, it's really just creating a small loop and implementing the calculation of the gradients, and then updating $\\theta$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for iteration in range(n_iterations):\n",
    "    gradients = 2/m * X_b.T.dot(X_b.dot(theta) - y)\n",
    "    theta = theta - eta * gradients"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cool, but did it do anything? We should probably check the values for $\\theta$ once again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"theta_0 = %.3f\" % theta[0][0])\n",
    "print(\"theta_1 = %.3f\" % theta[1][0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You might be surprised to see that these values basically are _exactly_ the same as those obtained with the normal equation. That's because our estimate, as much as before, completely depends on the data points we fed into the model. You can go to the earlier cells, change some of the parameters, and run the code again. Does anything change, for example when adjusting to use a larger/smaller dataset (the `m` parameter)?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, the implementation of batch gradient descent looks rather simple, but it's really not that obvious what happens in each step of the iteration. Remember: we look at the data points one thousand times, calculate some gradient one thousand times, update our estimate for $\\theta$ one thousand times, and only see the final result of that final step.\n",
    "\n",
    "Ignore the details of the following code (unless you're really interested of course). But here's what it does: we repeat the batch gradient descent method on the same dataset as before, but with different learning rates. When you execute the code, you'll see that the model with a very low learning rate only very slowly 'approaches' the dataset. The second one seems to be a lot faster in comparison. The third one, however, 'overshooots' the data with its very large learning rate. The model still converges, but it jumps around quite a bit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Store the path of BGD, we'll need that later.\n",
    "theta_path_bgd = []\n",
    "\n",
    "def plot_gradient_descent(theta, eta, theta_path=None):\n",
    "    plt.plot(X, y, \".\")\n",
    "    n_iterations = 1000\n",
    "    thetas = []\n",
    "    \n",
    "    for iteration in range(n_iterations):\n",
    "        if iteration < 15:\n",
    "            y_predict = X_b.dot(theta)\n",
    "            plt.plot(X, y_predict, color='#ff7f0e', alpha=0.2+0.06*iteration)\n",
    "        gradients = 2/m * X_b.T.dot(X_b.dot(theta) - y)\n",
    "        theta = theta - eta * gradients\n",
    "        thetas.append(theta)\n",
    "        # Again, this is just for later. If a path is provided,\n",
    "        # store it to visualise it in comparison to SGD and MBGD.\n",
    "        if theta_path is not None: theta_path.append(theta)\n",
    "\n",
    "    plt.xlabel(\"$x_1$\", fontsize=18)\n",
    "    plt.axis([0, 2, 0, 15])\n",
    "    plt.title(r\"$\\eta = {}$\".format(eta), fontsize=16)\n",
    "    return thetas\n",
    "\n",
    "np.random.seed(42)\n",
    "theta = np.random.randn(2,1)\n",
    "\n",
    "plt.figure(figsize=(12,8))\n",
    "plt.subplot(231)\n",
    "plt.ylabel(\"$y$\", rotation=0, fontsize=18)\n",
    "predictions_002 = plot_gradient_descent(theta, eta=0.02)\n",
    "plt.subplot(232)\n",
    "predictions_010 = plot_gradient_descent(theta, eta=0.1, theta_path=theta_path_bgd)\n",
    "plt.subplot(233)\n",
    "predictions_040 = plot_gradient_descent(theta, eta=0.4)\n",
    "plt.subplot(234)\n",
    "plt.axis([0, 30, 0, 6])\n",
    "plt.xlabel(\"$n_{iterations}$\", fontsize=18)\n",
    "plt.ylabel(\"$\\\\theta_0$, $\\\\theta_1$\", fontsize=18)\n",
    "plt.plot([i[0] for i in predictions_002])\n",
    "plt.plot([i[1] for i in predictions_002])\n",
    "plt.subplot(235)\n",
    "plt.axis([0, 30, 0, 6])\n",
    "plt.xlabel(\"$n_{iterations}$\", fontsize=18)\n",
    "plt.plot([i[0] for i in predictions_010])\n",
    "plt.plot([i[1] for i in predictions_010])\n",
    "plt.subplot(236)\n",
    "plt.axis([0, 30, 0, 6])\n",
    "plt.xlabel(\"$n_{iterations}$\", fontsize=18)\n",
    "plt.plot([i[0] for i in predictions_040])\n",
    "plt.plot([i[1] for i in predictions_040])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try out different learning rates and see what happens. You'll notice that – at very large rates – the model won't converge."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stochastic Gradient Descent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point, you'll probably know about the most problematic aspect of BGD already: we always evaluate the _entire_ dataset at every step of the training process. That's perfect for smaller datasets, because the model will usually find the 'best' estimate for $\\theta$ possible (similar to what the normal equation does). However, this is not very feasible for huge datasets.\n",
    "\n",
    "An alternative in that case is the _stochastic gradient descent_ (SGD) method. As opposed to batch gradient descent, the stochastic technique picks one instance from the dataset _randomly_ and adjusts the estimate for $\\theta$ according to that instance. This means _a lot_ of jumping around, because we follow the randomness of individual instances. On the other hand, it's computationally very inexpensive. And if it's done right, it can still find the 'best' estimate possible. However, it will _never_ converge to that 'best' estimate per se. One common technique to overcome this problem is to adjust the learning rate according to a _schedule_ during the training process. That is, start with a high learning rate and decrease it constantly to help the model to 'settle' in the global minimum."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Back to hands-on experience. Let's first decide for how many _epochs_ we would like to evaluate. The idea for this is that, despite an instance being picked randomly at every step, we would still want to evaluate them in sets of a certain size (which is usually just the size of the training data). Once we did, we call that an epoch and jump to the next one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_epochs = 50"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we should implement this learning-rate scheduling, because otherwise our model would not converge. Let's write a simple function to return a learning rate with two parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t0, t1 = 5, 50  \n",
    "def learning_schedule(t):\n",
    "    return t0 / (t + t1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once again, we'll have to initialise our $\\theta$ randomly before the first calculation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta = np.random.randn(2,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we can do the actual implementation of the loops:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a loop over the number of epochs.\n",
    "for epoch in range(n_epochs):\n",
    "    # And a loop over the size of the training data.\n",
    "    for i in range(m):\n",
    "        # Pick a random instance: this function takes\n",
    "        # a random value in the range from 0 to m.\n",
    "        random_index = np.random.randint(m)\n",
    "        \n",
    "        # Now, store this random instance and its target\n",
    "        # value into the variables xi and yi ...\n",
    "        xi = X_b[random_index:random_index+1]\n",
    "        yi = y[random_index:random_index+1]\n",
    "\n",
    "        # ... and calculate gradients and theta as before.\n",
    "        # Remember to calculate the value of the learning\n",
    "        # rate from the learning-rate schedule defined\n",
    "        # above (epoch * m + i) is just the number of the\n",
    "        # current training step over all epochs.\n",
    "        gradients = 2 * xi.T.dot(xi.dot(theta) - yi)\n",
    "        eta = learning_schedule(epoch * m + i)\n",
    "        theta = theta - eta * gradients"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We should probably check the values for $\\theta$ once again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"theta_0 = %.3f\" % theta[0][0])\n",
    "print(\"theta_1 = %.3f\" % theta[1][0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Great, it worked! The result seems to be similar, but not the same as before. Try changing some of the parameters. For instance, increase the number of epochs. Will we eventually reach the values from before? Remember: once we reach one thousand (!) epochs, we will have seen the same number of instances as if we had done one thousand iteration steps in the BGD method. Another thing you could try: change the learning-rate schedule (for example: turn it off). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an alternative, we can also use the _Scikit Learn_ library for stochastic gradient descent. Scikit Learn provides a plethora of functions and implementation of machine-learning models, and we will keep using it in more examples later. Of course, it might sometimes be helpful to implement algorithms yourself, but taking them from a library with dozens of experts working on them every day, doesn't sound too bad of an idea either.\n",
    "\n",
    "For stochastic gradient descent, Scikit Learn provides a class called [SGDRegressor](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.SGDRegressor.html). The good thing about this class is that we can just give it a bunch of parameters (many others will also be set by default), and then the class is ready to do exactly what we need."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we need to _import_ the class from Scikit Learn, then we can _instantiate_ it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import SGDRegressor\n",
    "sgd_reg = SGDRegressor(max_iter=50, eta0=0.1, penalty=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have already set a few parameters up there: we would like fifty epochs, and start with a learning rate $\\eta$ of 0.1. The Scikit Learn implementation uses learning-rate scheduling by default. We also need to set `penalty` explicitly, because otherwise Scikit Learn would otherwise use $\\ell_2$ regularisation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, the only thing left is to perform the actual training, by giving it our training data (this will also show you all parameters of our model):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sgd_reg.fit(X, y.ravel())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And what does it predict?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"theta_0: %.3f\" % sgd_reg.intercept_[0])\n",
    "print(\"theta_1: %.3f\" % sgd_reg.coef_[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To visualise the randomness of SGD, let's create some plot again. Again, you may just ignore the details of the following code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta_path_sgd = []\n",
    "theta = np.random.randn(2,1)\n",
    "n_epochs = 50\n",
    "\n",
    "t0, t1 = 10, 50  \n",
    "def learning_schedule(t):\n",
    "    return t0 / (t + t1)\n",
    "\n",
    "thetas = []\n",
    "\n",
    "plt.figure(figsize=(12,6))\n",
    "plt.subplot(121)\n",
    "plt.plot(X, y, \".\")\n",
    "\n",
    "for epoch in range(n_epochs):\n",
    "    for i in range(m):\n",
    "        if epoch == 0 and i < 10:\n",
    "            y_predict = X_b.dot(theta)\n",
    "            plt.plot(X, y_predict, color='#ff7f0e', alpha=0.2+0.06*i)\n",
    "            \n",
    "        random_index = np.random.randint(m)     \n",
    "        xi = X_b[random_index:random_index+1]\n",
    "        yi = y[random_index:random_index+1]\n",
    "        gradients = 2 * xi.T.dot(xi.dot(theta) - yi)\n",
    "        eta = learning_schedule(epoch * m + i)\n",
    "        theta = theta - eta * gradients\n",
    "        thetas.append(theta)\n",
    "        if epoch == 0 and i < 200: theta_path_sgd.append(theta)\n",
    "\n",
    "plt.xlabel(\"$x_1$\", fontsize=18)\n",
    "plt.ylabel(\"$y$\", rotation=0, fontsize=18)\n",
    "plt.axis([0, 2, 0, 15])\n",
    "plt.subplot(122)\n",
    "plt.plot([i[0] for i in thetas])\n",
    "plt.plot([i[1] for i in thetas])\n",
    "plt.axis([0, 500, 1, 7])\n",
    "plt.xlabel(\"$n_{iterations}$\", fontsize=18)\n",
    "plt.ylabel(\"$\\\\theta_0$, $\\\\theta_1$\", fontsize=18)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mini-batch gradient descent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By 'taking the best of both worlds', the _mini-batch gradient descent_ might be the perfect compromise between BGD and SGD. Instead of taking all or just one instance, the gradient is evaluated on a _mini-batch_ of instances. This makes it a little more stable than SGD, especially with large mini-batches. And it allows for vectorisation optimisations in terms of computing. It has the same 'issue' as SGD, however, that it never stops at the optimal values for the estimators, but keeps jumping around the global minimum. Therefore, a good learning schedule is pivotal to implement this technique successfully."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's first fix some parameters: the number of epochs, and the size of the mini batches:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_epochs = 50\n",
    "minibatch_size = 20"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, implement learning-rate scheduling as before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t0, t1 = 200, 1000\n",
    "def learning_schedule(t):\n",
    "    return t0 / (t + t1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make sure to initialise the values of $\\theta$ randomly in the beginning:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta = np.random.randn(2,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A small helper function to give us randomly shuffled instances. We shuffle the instances for every new epoch to make sure we will always look at 'new' mini batches for every epoch. Otherwise, we would end up looking at all instances in the same order in all epochs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def shuffled_instances(m, X_b, y):\n",
    "    shuffled_indices = np.random.permutation(m)\n",
    "    return X_b[shuffled_indices], y[shuffled_indices]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can do the actual looping:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Count the number of steps for the learning-rate schedule.\n",
    "step = 0\n",
    "\n",
    "# Loop over all epochs.\n",
    "for epoch in range(n_epochs):\n",
    "    # Start off every epoch by shuffling around the instances.\n",
    "    X_b_shuffled, y_shuffled = shuffled_instances(m, X_b, y)\n",
    "    \n",
    "    # Loop over all instances in mini batches.\n",
    "    for i in range(0, m, minibatch_size):\n",
    "        step += 1\n",
    "        \n",
    "        # Before we can calculate the gradients, make sure we\n",
    "        # get the mini-batch subset of instances.\n",
    "        xi = X_b_shuffled[step:step+minibatch_size]\n",
    "        yi = y_shuffled[step:step+minibatch_size]\n",
    "        \n",
    "        # Then, calculate the gradient and update theta as before.\n",
    "        gradients = 2/minibatch_size * xi.T.dot(xi.dot(theta) - yi)\n",
    "        eta = learning_schedule(step)\n",
    "        theta = theta - eta * gradients"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And what are the predicted values of the parameters?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"theta_0: %.3f\" % theta[0][0])\n",
    "print(\"theta_1: %.3f\" % theta[1][0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code won't produce any output, but please run it quickly for the last part of this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta_path_mgd = []\n",
    "n_iterations = 50\n",
    "minibatch_size = 20\n",
    "np.random.seed(42)\n",
    "theta = np.random.randn(2,1)\n",
    "\n",
    "t0, t1 = 200, 1000\n",
    "def learning_schedule(t):\n",
    "    return t0 / (t + t1)\n",
    "\n",
    "step = 0\n",
    "\n",
    "for epoch in range(n_iterations):\n",
    "    shuffled_indices = np.random.permutation(m)\n",
    "    X_b_shuffled = X_b[shuffled_indices]\n",
    "    y_shuffled = y[shuffled_indices]\n",
    "    \n",
    "    for i in range(0, m, minibatch_size):\n",
    "        step += 1\n",
    "        xi = X_b_shuffled[step:step+minibatch_size]\n",
    "        yi = y_shuffled[step:step+minibatch_size]\n",
    "        gradients = 2/minibatch_size * xi.T.dot(xi.dot(theta) - yi)\n",
    "        eta = learning_schedule(step)\n",
    "        theta = theta - eta * gradients\n",
    "        theta_path_mgd.append(theta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparing BGD, SGD and mini-batch GD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're almost done, but let's quickly have a look at what the algorithms actually do in the parameter space. The following code produces a plot that points out every step the batch gradient descent, stochastic gradient descent, and mini-batch gradient descent methods took. Do you spot any differences?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta_path_bgd = np.array(theta_path_bgd)\n",
    "theta_path_sgd = np.array(theta_path_sgd)\n",
    "theta_path_mgd = np.array(theta_path_mgd)\n",
    "\n",
    "plt.figure(figsize=(7,4))\n",
    "plt.plot(theta_path_sgd[:, 0], theta_path_sgd[:, 1], \"-s\", linewidth=1, label=\"Stochastic\")\n",
    "plt.plot(theta_path_mgd[:, 0], theta_path_mgd[:, 1], \"-+\", linewidth=2, label=\"Mini-batch\")\n",
    "plt.plot(theta_path_bgd[:, 0], theta_path_bgd[:, 1], \"-o\", linewidth=3, label=\"Batch\")\n",
    "plt.legend(loc=\"upper left\", fontsize=16)\n",
    "plt.xlabel(r\"$\\theta_0$\", fontsize=20)\n",
    "plt.ylabel(r\"$\\theta_1$   \", fontsize=20, rotation=0)\n",
    "plt.axis([1.0, 4.0, 2.5, 6.0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.11"
  },
  "nav_menu": {},
  "toc": {
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 6,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}