task_6_regularisation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Task 6: Regularisation\n",
    "\n",
    "_All credit for the code examples of this notebook goes to the book \"Hands-On Machine Learning with Scikit-Learn & TensorFlow\" by A. Geron. Modifications were made and text was added by K. Zoch in preparation for the hands-on sessions._"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, import a few common modules, ensure MatplotLib plots figures inline and prepare a function to save the figures:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Common imports\n",
    "import numpy as np\n",
    "import os\n",
    "\n",
    "# to make this notebook's output stable across runs\n",
    "np.random.seed(42)\n",
    "\n",
    "# To plot pretty figures\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "mpl.rc('axes', labelsize=14)\n",
    "mpl.rc('xtick', labelsize=12)\n",
    "mpl.rc('ytick', labelsize=12)\n",
    "\n",
    "# Function to save a figure. This also decides that all output files \n",
    "# should stored in the subdirectorz 'classification'.\n",
    "PROJECT_ROOT_DIR = \".\"\n",
    "EXERCISE = \"regularisation\"\n",
    "\n",
    "def save_fig(fig_id, tight_layout=True):\n",
    "    path = os.path.join(PROJECT_ROOT_DIR, \"output\", EXERCISE, fig_id + \".png\")\n",
    "    print(\"Saving figure\", fig_id)\n",
    "    if tight_layout:\n",
    "        plt.tight_layout()\n",
    "    plt.savefig(path, format='png', dpi=300)\n",
    "\n",
    "# Ignore useless warnings (see SciPy issue #5998)\n",
    "import warnings\n",
    "warnings.filterwarnings(action=\"ignore\", message=\"^internal gelsd\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Polynomial regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even for non-linear data we can make use of linear regression – we simply need to add higher degrees of features to the set of features. Using those 'new' extended features, linear regression can still give us good results. Let's get started by generating some random data, with a maximum degree of 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.random as rnd\n",
    "\n",
    "np.random.seed(42)\n",
    "\n",
    "m = 100\n",
    "X = 6 * np.random.rand(m, 1) - 3\n",
    "y = 0.5 * X**2 + X + 2 + np.random.randn(m, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's have a look at how the data is distributed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(X, y, \"b.\")\n",
    "plt.xlabel(\"$x_1$\", fontsize=18)\n",
    "plt.ylabel(\"$y$\", rotation=0, fontsize=18)\n",
    "plt.axis([-3, 3, 0, 10])\n",
    "save_fig(\"quadratic_data_plot\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now comes your first task. We would like to create a feature set which also includes all features to the power of 2. Can you create this new feature set and then perform a linear regression as we already did in the previous task? One helpful class in Scikit-Learn is the [PolynominalFeatures](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.PolynomialFeatures.html) class. You will have to create an instance of that class and somehow get a 1-dimensional array of the new features. Then, you can use the [LinearRegression](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html) to fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "\n",
    "X_poly # To be implemented: polynomial (degree=2) version of our data\n",
    "\n",
    "# Let's test if that actually worked. If we look at instance 0,\n",
    "# do we have the value to the power of 2 in there?\n",
    "print(X[0])\n",
    "print(X_poly[0])\n",
    "\n",
    "# Now perform the linear regression.\n",
    "from sklearn.linear_model import LinearRegression\n",
    "lin_reg # To be implemented: linear regression instance, can you fit the data X_poly?\n",
    "\n",
    "# And print the fit results.\n",
    "print(\"Fitted values: intercept = %s, coefficient = %s\" % (lin_reg.intercept_[0], lin_reg.coef_[0][0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will see if your implementation worked with the following piece of code. It will again print the dataset created above, and then also plot the fitted model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a set of 100 X values in the interval [-3, 3], which\n",
    "# is the area we want to plot, and create a 100x1 array from them.\n",
    "X_new=np.linspace(-3, 3, 100).reshape(100, 1)\n",
    "\n",
    "# Now use the PolynomialFeatures class to create a feature matrix.\n",
    "X_new_poly = poly_features.transform(X_new)\n",
    "\n",
    "# Make predictions based on this new feature matrix.\n",
    "y_new = lin_reg.predict(X_new_poly)\n",
    "\n",
    "plt.plot(X, y, \"b.\")\n",
    "plt.plot(X_new, y_new, \"r-\", linewidth=2, label=\"Predictions\")\n",
    "plt.xlabel(\"$x_1$\", fontsize=18)\n",
    "plt.ylabel(\"$y$\", rotation=0, fontsize=18)\n",
    "plt.legend(loc=\"upper left\", fontsize=14)\n",
    "plt.axis([-3, 3, 0, 10])\n",
    "save_fig(\"quadratic_predictions_plot\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pre-step: learning curves\n",
    "\n",
    "Before we can go into detail about the regularisation technqiues, we need to come back to another basic performance measure: the learning curve. Learning curves plot the model's performance on data with increasing size of the training set. They are particularly interesting to check convergence behaviour of the model, but also to compare the model's performance on training and validation set. Reminder: our typical cost function on the y axis is the mean squared error (MSE), or the square-root of that (RSME).\n",
    "\n",
    "Can you try to implement the needed computation of the RSME as a function of the training set size? The first step will be to split the dataset into training and test, a good function to use for that is [train_test_split](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.train_test_split.html). Maybe having a training:test ratio of 80:20 would be a good start. Then, create a for loop which iterates over the number of training cycles, and – for each number of cycles – make a fit with the model. Then, predict the model performance on the training set (only looking at the instances it has already seen!) and the complete validation set. You can then store the RSME values of the prediction using the [mean_squared_error](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.mean_squared_error.html) function. We will be using the [LinearRegression](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html) class for the defined function, so make sure to check out the `fit` and `predict` functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import mean_squared_error\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "# Define a function which takes a training model, and \n",
    "# our sets of X and y values.\n",
    "def plot_learning_curves(model, X, y):\n",
    "    # First, split into training and validation sets.\n",
    "    # Implemnt instances such as X_train, X_val, y_train and y_val\n",
    "    \n",
    "    # Make sure to store the RSME values for training\n",
    "    # and validation – we want to plot them later.\n",
    "    train_errors, val_errors = [], []\n",
    "    \n",
    "    # Then start the loop, where we want to go from 1 to the\n",
    "    # total number of training cycles, but for each iteration,\n",
    "    # we want to evaluate the instances _up to_ that point!\n",
    "    for m in range(1, len(X_train)):\n",
    "        # implement here\n",
    "        \n",
    "        # Store the values into the lists.\n",
    "        train_errors.append(mean_squared_error())  # to be implemented\n",
    "        val_errors.append(mean_squared_error())    # to be implemented\n",
    "\n",
    "    # And the actual plotting commands. Plot the RSME for the\n",
    "    # training set in red plus signs, and the RSME for the \n",
    "    # validation set in blue.\n",
    "    plt.plot(np.sqrt(train_errors), \"r-+\", linewidth=2, label=\"train\")\n",
    "    plt.plot(np.sqrt(val_errors), \"b-\", linewidth=3, label=\"val\")\n",
    "    plt.legend(loc=\"upper right\", fontsize=14)\n",
    "    plt.xlabel(\"Training set size\", fontsize=14)\n",
    "    plt.ylabel(\"RMSE\", fontsize=14)\n",
    "\n",
    "# Now create an instance of the LinearRegression class, and \n",
    "# create a plot with the function we just defined.\n",
    "lin_reg = LinearRegression()\n",
    "plot_learning_curves(lin_reg, X, y)\n",
    "plt.axis([0, 80, 0, 3])\n",
    "save_fig(\"underfitting_learning_curves_plot\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just for the fun of it, the following code creates a polynomial set of feature up to the power of 10, which we can then fit using the LinearRegression class. Is the fitting behaviour different?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.pipeline import Pipeline\n",
    "\n",
    "polynomial_regression = Pipeline([\n",
    "        (\"poly_features\", PolynomialFeatures(degree=10, include_bias=False)),\n",
    "        (\"lin_reg\", LinearRegression()),\n",
    "    ])\n",
    "\n",
    "plot_learning_curves(polynomial_regression, X, y)\n",
    "plt.axis([0, 80, 0, 3])\n",
    "save_fig(\"learning_curves_plot\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regularized models\n",
    "\n",
    "There are two basic types of regularisation, both of which are based on mathematical norms. The first one, the Ridge regression, uses the Euclidian L2 norm. This enters the cost function for the training as an additional 'penalty' term, scaled with a parameter alpha. With small alphas, we essentially 'turn off' regularisation. Do you have an idea what happens with large values for alpha? How does the regularisation kick in?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's start by generating some random data.\n",
    "np.random.seed(42)\n",
    "m = 40\n",
    "X = 3 * np.random.rand(m, 1)\n",
    "y = 1 + 0.5 * X + np.random.randn(m, 1) / 2.5\n",
    "\n",
    "# We will need these data points later for our predictions.\n",
    "X_new = np.linspace(0, 3, 100).reshape(100, 1)\n",
    "\n",
    "# And plot it.\n",
    "plt.plot(X, y, \"b.\", linewidth=3)\n",
    "plt.xlabel(\"$x_1$\", fontsize=18)\n",
    "plt.axis([0, 3, 0, 4])\n",
    "save_fig(\"ridge_regression_plot_data\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This looks like some good test data to try Ridge regularisation with. Below you can find two functions, which still have some functionality missing. Can you complete them? The first one is meant to plot linear regressions for different values of alpha, the latter the same for polynomial regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import Ridge\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "# Function to take multiple values of alpha and then\n",
    "# create instances of Ridge models for each of them.\n",
    "#   - model_class is a flexible parameter and could take\n",
    "#     various models from sklearn.linear_model\n",
    "#   - alphas should be a tuple of alpha values\n",
    "def plot_model_lin(model_class, alphas):\n",
    "    # Let's combine the alpha values with different plotting styles.\n",
    "    alpha_styles = zip(alphas, (\"b-\", \"g--\", \"r:\"))\n",
    "    \n",
    "    # Now we can start our loop over the zipped alphas. What\n",
    "    # we will need here is to instantiate a model object of\n",
    "    # the desired class (maybe alo give it a fixed random_state\n",
    "    # value), and then perform the fit on our X, y data. Then,\n",
    "    # use the X_new to make a prediction, which we then want to\n",
    "    # plot together with the data. Can you implement that?\n",
    "    for alpha, style in alpha_styles:\n",
    "        # Implement a model with parameter alpha and a fixed\n",
    "        # random state. Fit the data X, y, then make a prediction\n",
    "        # on the X_new data points created above (which will\n",
    "        # give us something like y_new_regul).\n",
    "        \n",
    "        # Plot the results.\n",
    "        plt.plot(X_new, y_new_regul, style, linewidth=2, label=r\"$\\alpha = {}$\".format(alpha))\n",
    "    \n",
    "    # This will also plot the data, create a legend etc.\n",
    "    plt.plot(X, y, \"b.\", linewidth=3)\n",
    "    plt.legend(loc=\"upper left\", fontsize=15)\n",
    "    plt.xlabel(\"$x_1$\", fontsize=18)\n",
    "    plt.axis([0, 3, 0, 4])\n",
    "\n",
    "# Function to take multiple values of alpha and then\n",
    "# create instances of Ridge models for each of them. This\n",
    "# will automatically expand the feature set with polynomial\n",
    "# features up to the power of 10.\n",
    "#   - model_class is a flexible parameter and could take\n",
    "#     various models from sklearn.linear_model\n",
    "#   - alphas should be a tuple of alpha values\n",
    "#   - **model_args in case we want to forward arguments to the\n",
    "#     instance of the model class (ignore this for now).\n",
    "def plot_model_poly(model_class, alphas, **model_kargs):\n",
    "    # Let's combine the alpha values with different plotting styles.\n",
    "    alpha_styles = zip(alphas, (\"b-\", \"g--\", \"r:\"))\n",
    "    \n",
    "    # Now we can start our loop over the zipped alphas.\n",
    "    for alpha, style in alpha_styles:\n",
    "        # Implement a model with parameter alpha, a fixed random\n",
    "        # state, and also forward the **model_kargs arguments. \n",
    "        \n",
    "        # This is just to make our life easier. We make a few\n",
    "        # transformations of our model, by adding polynomial\n",
    "        # features up the power of 10, implement a standard \n",
    "        # scaler, and eventually get back the updated model.\n",
    "        # All of these could be done in individual steps, but\n",
    "        # the sklearn.pipeline.Pipeline class makes this a lot\n",
    "        # easier.\n",
    "        model = Pipeline([\n",
    "                (\"poly_features\", PolynomialFeatures(degree=10, include_bias=False)),\n",
    "                (\"std_scaler\", StandardScaler()),\n",
    "                (\"regul_reg\", model),\n",
    "            ])\n",
    "\n",
    "        # Implement: perform a fit on data X, y and make a\n",
    "        # prediction on X_new. The result should be something like\n",
    "        # y_new_regul.\n",
    "\n",
    "        # Plot the results.\n",
    "        plt.plot(X_new, y_new_regul, style, linewidth=2, label=r\"$\\alpha = {}$\".format(alpha))\n",
    "    # This will also plot the data, create a legend etc.\n",
    "    plt.plot(X, y, \"b.\", linewidth=3)\n",
    "    plt.legend(loc=\"upper left\", fontsize=15)\n",
    "    plt.xlabel(\"$x_1$\", fontsize=18)\n",
    "    plt.axis([0, 3, 0, 4])\n",
    "\n",
    "# Now call the above functions and make two comparison plots for\n",
    "# exemplary values of alpha.\n",
    "plt.figure(figsize=(8,4))\n",
    "plt.subplot(121)\n",
    "plot_model_lin(Ridge, alphas=(0, 10, 100))\n",
    "plt.ylabel(\"$y$\", rotation=0, fontsize=18)\n",
    "plt.subplot(122)\n",
    "plot_model_poly(Ridge, alphas=(0, 10**-5, 1))\n",
    "\n",
    "save_fig(\"ridge_regression_plot\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks good! To get a better idea of the predicted values, let's look at one example point at 1.5. In principle, the Ridge model can be implemented with a closed-form solution (remember the lecture), or a gradient descent method. Can you compare these two approaches? The two classes for this are [linear_model.Ridge](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Ridge.html) and [linear_model.SGDRegressor](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.SGDRegressor.html), which we have both used before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import Ridge\n",
    "from sklearn.linear_model import SGDRegressor\n",
    "\n",
    "# Create an instance of Ridge regression following the\n",
    "# closed-form solution.\n",
    "ridge_reg # To be implemented\n",
    "print(\"Closed form predicts: %s\" % ridge_reg.predict([[1.5]])[0][0])\n",
    "\n",
    "# Create an SGD regressor and implement Ridge regularisation.\n",
    "sgd_reg # To be implemented\n",
    "print(\"SGD regressor predicts: %s\" % sgd_reg.predict([[1.5]])[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A second regularisation technique is the Lasso regression, which uses the L2 norm (i.e. \"least absolute deviations\" instead of \"least squares\"). In case you're not, please make yourself familiar with these two norms, it is important to understand their impact on the cost function. Again, Lasso regression is implemented as a 'penalty' term with a parameter alpha to control the impact. What would be the typical behaviour of a regressor with Lasso regression implemented? Do the weights tend to be small/large? Are they equally distributed or far apart?\n",
    "\n",
    "The following piece of code uses our previously defined functions for linear and polynomial regression, but plots Lasso regression for different values of alpha instead. Play around with those values to see what happens. We also finally get to use what we implemented earlier: we can forward additional arguments to the model. Here `tol=1` is set to one. Can you check what it does? Class reference: [linear_model.Lasso](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Lasso.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import Lasso\n",
    "\n",
    "plt.figure(figsize=(8,4))\n",
    "plt.subplot(121)\n",
    "plot_model_lin(Lasso, alphas=(0, 0.1, 1))\n",
    "plt.ylabel(\"$y$\", rotation=0, fontsize=18)\n",
    "plt.subplot(122)\n",
    "plot_model_poly(Lasso, alphas=(0, 10**-7, 1), tol=1)\n",
    "\n",
    "save_fig(\"lasso_regression_plot\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, as a last step before jumping into the next topic, it would be nice to compare the above predictions for Ridge regression with those of Lasso regression. A third model is [linear_model.ElasticNet](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.ElasticNet.html) which essentially implements the best of both worlds, i.e. it is a mixed version of L1 and L2 regularisation. Can you try to implement Lasso and ElasticNet and make a prediction for the value 1.5?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import Lasso\n",
    "from sklearn.linear_model import ElasticNet\n",
    "\n",
    "# Make a prediction with Lasso regression.\n",
    "lasso_reg # To be implemented\n",
    "print(\"Lasso predicts: %s\" % lasso_reg.predict([[1.5]])[0])\n",
    "\n",
    "# Make a prediction with ElasticNet.\n",
    "elastic_net # To be implemented\n",
    "print(\"ElasticNet predicts: %s\" % elastic_net.predict([[1.5]])[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Early stopping\n",
    "\n",
    "Another important, but completely different aspect of regularisation is 'early stopping'. The idea is to stop the training as soon as the cost function of the validation set reaches its minimum. This is basically the turning point, where the validation cost function will start to increase again and the model starts to overfit the training data. This first bit of code only visualises the idea of early stopping, but doesn't actually implement it. We will do that later. But can you add the missing parts to the code below? You might have to look into the definitions of [train_test_split](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.train_test_split.html) and [linear_model.SGDRegressor](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.SGDRegressor.html) again for this. For the latter, please look up the `penalty` argument. What does it do? In our case, we probably want to set it to `None`. Other important arguments to look up and include are `max_iter` (how many do we want?), `eta0`, `warm_start` (do we need this?), and `learning_rate` (we want the learning rate to be constant here!)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Generate random data, both for training and validation.\n",
    "np.random.seed(42)\n",
    "m = 100\n",
    "X = 6 * np.random.rand(m, 1) - 3\n",
    "y = 2 + X + 0.5 * X**2 + np.random.randn(m, 1)\n",
    "\n",
    "# As before, we will need to split into training and validation\n",
    "# set. Can you implement that? This time, maybe a ratio of 50:50\n",
    "# would be desirable. After splitting we should have arrays\n",
    "# such as X_train, X_val, y_train and y_val.\n",
    "\n",
    "# Create a small pipeline which adds polynomial features to our\n",
    "# feature set and applies a standard scaler (which removes the\n",
    "# mean and scales to unit variance).\n",
    "poly_scaler = Pipeline([\n",
    "        (\"poly_features\", PolynomialFeatures(degree=90, include_bias=False)),\n",
    "        (\"std_scaler\", StandardScaler()),\n",
    "    ])\n",
    "\n",
    "# Apply the poly_scaler pipeline to our training and validation sets.\n",
    "X_train_poly_scaled = poly_scaler.fit_transform(X_train)\n",
    "X_val_poly_scaled = poly_scaler.transform(X_val)\n",
    "\n",
    "# Instantiate an SGDRegressor. Make sure to set the iteration\n",
    "# to 1, because we want to control the number of epochs by hand\n",
    "# (see below).\n",
    "sgd_reg # To be implemented\n",
    "\n",
    "# Create arrays to store the mean squared errors on training\n",
    "# and validation data.\n",
    "train_errors, val_errors = [], []\n",
    "\n",
    "# Now let's loop over 500 epochs, fit the training data once\n",
    "# per epoch, and make predictions on the training and the\n",
    "# validation dataset. Remember: what is one epoch for the \n",
    "# stochastic GD method? How many instances does the model see\n",
    "# per epoch?\n",
    "for epoch in range(500):\n",
    "    # Implement fit and predictions here.\n",
    "    \n",
    "    # Store the mean squared errors in the arrays.\n",
    "    train_errors.append(mean_squared_error())  # To be implemented\n",
    "    val_errors.append(mean_squared_error())    # To be implemented\n",
    "\n",
    "# Now let's get some info about the 'best' epoch. We want to\n",
    "# know the epoch which had the smallest mean sqared error on\n",
    "# the validation set. numpy has a function for that ... Then,\n",
    "# with the best epoch extracted, we also want to know which\n",
    "# value was the best.\n",
    "best_epoch = np.argmin(val_errors)\n",
    "best_val_rmse = np.sqrt(val_errors[best_epoch])\n",
    "\n",
    "# Now this is just plotting. Make an annotation where the best\n",
    "# model was actually located, and point it out in the plot.\n",
    "plt.annotate('Best model',\n",
    "             xy=(best_epoch, best_val_rmse),\n",
    "             xytext=(best_epoch, best_val_rmse + 1),\n",
    "             ha=\"center\",\n",
    "             arrowprops=dict(facecolor='black', shrink=0.05),\n",
    "             fontsize=16,\n",
    "            )\n",
    "\n",
    "# And do the plotting.\n",
    "best_val_rmse -= 0.03  # just to make the graph look better\n",
    "plt.plot([0, 500], [best_val_rmse, best_val_rmse], \"k:\", linewidth=2)\n",
    "plt.plot(np.sqrt(val_errors), \"b-\", linewidth=3, label=\"Validation set\")\n",
    "plt.plot(np.sqrt(train_errors), \"r--\", linewidth=2, label=\"Training set\")\n",
    "plt.legend(loc=\"upper right\", fontsize=14)\n",
    "plt.xlabel(\"Epoch\", fontsize=14)\n",
    "plt.ylabel(\"RMSE\", fontsize=14)\n",
    "save_fig(\"early_stopping_plot\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now of couse it would be interesting to actually make the model _stop_ after the best epoch. Of course, Scikit-Learn provides functionality for that, but it's not really difficult to implement this ourselves. Let's try this in the following piece of code. So be able to save the model state after each epoch, we need to import the [base.clone](https://scikit-learn.org/stable/modules/generated/sklearn.base.clone.html) method, which can create copies of our model instance. Can you implement the rest yourself? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.base import clone\n",
    "\n",
    "# Instantiate an SGDRegressor. Make sure to set the iteration\n",
    "# to 1, because we want to control the number of epochs by hand\n",
    "# (see below). As above, also check that you set values for\n",
    "# 'warm_start', 'penalty', 'learning_rate' and 'eta'.\n",
    "sgd_reg # To be implemented\n",
    "\n",
    "# As before, we will need to split into training and validation\n",
    "# set. Can you implement that? This time, maybe a ratio of 50:50\n",
    "# would be desirable.\n",
    "\n",
    "# Create a small pipeline which adds polynomial features to our\n",
    "# feature set and applies a standard scaler (which removes the\n",
    "# mean and scales to unit variance).\n",
    "poly_scaler = Pipeline([\n",
    "        (\"poly_features\", PolynomialFeatures(degree=90, include_bias=False)),\n",
    "        (\"std_scaler\", StandardScaler()),\n",
    "    ])\n",
    "\n",
    "# Apply the poly_scaler pipeline to our training and validation sets.\n",
    "X_train_poly_scaled = poly_scaler.fit_transform(X_train)\n",
    "X_val_poly_scaled = poly_scaler.transform(X_val)\n",
    "\n",
    "# Our reference value for the validation error, which we will update\n",
    "# for every epoch, starting with 'inf'.\n",
    "minimum_val_error = float(\"inf\")\n",
    "\n",
    "# Store our best epochs and model instances in these variables.\n",
    "best_epoch = None\n",
    "best_model = None\n",
    "\n",
    "# Now perform the loop. Let's go through 1000 epochs. Make sure\n",
    "# that the regressor continues where it left off in the epoch\n",
    "# before (check the class documentation). Then, perform a fit on\n",
    "# the poly-scaled version of our training data. Afterwards, make\n",
    "# a prediction on our poly-scaled validation data, and store the\n",
    "# mean_squared_error value of that. Is it smaller than the\n",
    "# minimum value observed so far? What should we do in that case?\n",
    "for epoch in range(1000):\n",
    "    # Implement fit, prediction and calculation of error.\n",
    "\n",
    "    # Take appropriate actions if the observed error is smaller\n",
    "    # than the one previously observed.\n",
    "    if val_error < minimum_val_error:\n",
    "        # Implement code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Did it work? What is the current best_epoch and best_model?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "best_epoch, best_model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following (large) bit of code is directly taken from Geron's book. You don't need to understand the details of the implementation, but rather take the plots as a nice inspiration. They nicely show the different characteristics of Ridge and Lasso regression. Explanations are found below the plots. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t1a, t1b, t2a, t2b = -1, 3, -1.5, 1.5\n",
    "\n",
    "# ignoring bias term\n",
    "t1s = np.linspace(t1a, t1b, 500)\n",
    "t2s = np.linspace(t2a, t2b, 500)\n",
    "t1, t2 = np.meshgrid(t1s, t2s)\n",
    "T = np.c_[t1.ravel(), t2.ravel()]\n",
    "Xr = np.array([[-1, 1], [-0.3, -1], [1, 0.1]])\n",
    "yr = 2 * Xr[:, :1] + 0.5 * Xr[:, 1:]\n",
    "\n",
    "J = (1/len(Xr) * np.sum((T.dot(Xr.T) - yr.T)**2, axis=1)).reshape(t1.shape)\n",
    "\n",
    "N1 = np.linalg.norm(T, ord=1, axis=1).reshape(t1.shape)\n",
    "N2 = np.linalg.norm(T, ord=2, axis=1).reshape(t1.shape)\n",
    "\n",
    "t_min_idx = np.unravel_index(np.argmin(J), J.shape)\n",
    "t1_min, t2_min = t1[t_min_idx], t2[t_min_idx]\n",
    "\n",
    "t_init = np.array([[0.25], [-1]])\n",
    "\n",
    "def bgd_path(theta, X, y, l1, l2, core = 1, eta = 0.1, n_iterations = 50):\n",
    "    path = [theta]\n",
    "    for iteration in range(n_iterations):\n",
    "        gradients = core * 2/len(X) * X.T.dot(X.dot(theta) - y) + l1 * np.sign(theta) + 2 * l2 * theta\n",
    "\n",
    "        theta = theta - eta * gradients\n",
    "        path.append(theta)\n",
    "    return np.array(path)\n",
    "\n",
    "plt.figure(figsize=(12, 8))\n",
    "for i, N, l1, l2, title in ((0, N1, 0.5, 0, \"Lasso\"), (1, N2, 0,  0.1, \"Ridge\")):\n",
    "    JR = J + l1 * N1 + l2 * N2**2\n",
    "    \n",
    "    tr_min_idx = np.unravel_index(np.argmin(JR), JR.shape)\n",
    "    t1r_min, t2r_min = t1[tr_min_idx], t2[tr_min_idx]\n",
    "\n",
    "    levelsJ=(np.exp(np.linspace(0, 1, 20)) - 1) * (np.max(J) - np.min(J)) + np.min(J)\n",
    "    levelsJR=(np.exp(np.linspace(0, 1, 20)) - 1) * (np.max(JR) - np.min(JR)) + np.min(JR)\n",
    "    levelsN=np.linspace(0, np.max(N), 10)\n",
    "    \n",
    "    path_J = bgd_path(t_init, Xr, yr, l1=0, l2=0)\n",
    "    path_JR = bgd_path(t_init, Xr, yr, l1, l2)\n",
    "    path_N = bgd_path(t_init, Xr, yr, np.sign(l1)/3, np.sign(l2), core=0)\n",
    "\n",
    "    plt.subplot(221 + i * 2)\n",
    "    plt.grid(True)\n",
    "    plt.axhline(y=0, color='k')\n",
    "    plt.axvline(x=0, color='k')\n",
    "    plt.contourf(t1, t2, J, levels=levelsJ, alpha=0.9)\n",
    "    plt.contour(t1, t2, N, levels=levelsN)\n",
    "    plt.plot(path_J[:, 0], path_J[:, 1], \"w-o\")\n",
    "    plt.plot(path_N[:, 0], path_N[:, 1], \"y-^\")\n",
    "    plt.plot(t1_min, t2_min, \"rs\")\n",
    "    plt.title(r\"$\\ell_{}$ penalty\".format(i + 1), fontsize=16)\n",
    "    plt.axis([t1a, t1b, t2a, t2b])\n",
    "    if i == 1:\n",
    "        plt.xlabel(r\"$\\theta_1$\", fontsize=20)\n",
    "    plt.ylabel(r\"$\\theta_2$\", fontsize=20, rotation=0)\n",
    "\n",
    "    plt.subplot(222 + i * 2)\n",
    "    plt.grid(True)\n",
    "    plt.axhline(y=0, color='k')\n",
    "    plt.axvline(x=0, color='k')\n",
    "    plt.contourf(t1, t2, JR, levels=levelsJR, alpha=0.9)\n",
    "    plt.plot(path_JR[:, 0], path_JR[:, 1], \"w-o\")\n",
    "    plt.plot(t1r_min, t2r_min, \"rs\")\n",
    "    plt.title(title, fontsize=16)\n",
    "    plt.axis([t1a, t1b, t2a, t2b])\n",
    "    if i == 1:\n",
    "        plt.xlabel(r\"$\\theta_1$\", fontsize=20)\n",
    "\n",
    "save_fig(\"lasso_vs_ridge_plot\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start with the top left plot, on the x and y axis are two model parameters. The background ellipses (contours) show the behaviour of the MSE cost function _without_ regularisation. The foreground diamond-shaped contours represent the penalty that the L1 norm applies. The white circles represent the path that an _unregularised_ regressor would go to find the minimum (i.e. with `alpha = 0`). The yellow triangles show the path that a _purely_ penalty-based regressor would go (i.e. with `alpha -> infty`). What's interesting to see is that the regressor first 'walks back' to `theta_1 = 0`, and then walks along the y axis to reach the minimum of the diamond-shaped contours.\n",
    "\n",
    "The ellipse contours in the top right plot show the cost function with an L1 regularisation `alpha = 0.5`. Compared to the left-hand plot, you will notice, how the global minimum is shifted to `theta_2 = 0`, and overall is closer to smaller values of `theta`. The white circles are the path the regressor would go with this regularisation term.\n",
    "\n",
    "The two plots on the bottom show the same for Ridge regression, i.e. using L2 regularisation. First of all, notice how the contours for `alpha -> infty` on the bottom left are now circular. This also changes the path the regressor would take for a _purely_ penalty-based model: the yellow triangles show a path essentially orthogonal to the circular contour lines. \n",
    "\n",
    "The plot on the bottom right again shows the regression behaviour for L2 regularisation with `alpha = 0.5`. Notice how the global minimum now it _not_ at `theta_2 = 0` is it is for L1 regularisation. But still, overall it is closer to the smaller values of `theta` than in the unregularised case."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "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": 1
}