3_stochastic_gradient_descent.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0afd1513",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Notebook 3: Stochastic Gradient Descent\n",
    "\n",
    "Gradient descent is the workhorse of machine learning. The goal of this notebook is to develop the basic algorithms in the context of two common problems: a simple linear regression and logistic regression for binary classification. \n",
    "\n",
    "This notebook was written by [Joe Boyd](https://github.com/jcboyd), [Judith Abecassis](https://github.com/judithabk6) and [Chloé-Agathe Azencott](http://cazencott.info), with inspiration from [the AM207 course at Harvard](http://am207.github.io/2016/) and [Francis Bach](http://www.di.ens.fr/%7Efbach/orsay2018.html).\n",
    "\n",
    "\n",
    "__If you are taking this class for credit, you will need to turn in your completed notebook on Campus.__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a9138c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa2d3802-30c5-477f-90da-0ea0c137c1ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rc('font', **{'size': 12}) # règle la taille de police globalement pour les plots (en pt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a576cfb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Obtaining the colors that matplotlib cycles on\n",
    "# to be able to use them directly\n",
    "# see: https://matplotlib.org/stable/gallery/color/color_cycle_default.html\n",
    "prop_cycle = plt.rcParams['axes.prop_cycle']\n",
    "notebook_colors = prop_cycle.by_key()['color']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38ed9c27",
   "metadata": {},
   "source": [
    "## 1. Linear regression\n",
    "\n",
    "In this section, we will solve a simple linear regression between two variables: the number of features is 1."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "304176f5",
   "metadata": {},
   "source": [
    "### 1.1 Generating data\n",
    "\n",
    "Here we use the `make_regression` method of `sklearn.datasets` ([doc](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_regression.html)) to generate our data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6423d1d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5f62563",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_samples = 100\n",
    "n_features = 1\n",
    "X, y = datasets.make_regression(n_samples=n_samples, n_features=n_features, n_informative=1, noise=20, random_state=37)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "069d79f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "X.shape, y.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2dbfdd8f",
   "metadata": {},
   "source": [
    "Let us visualize the data we have generated:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8bbd7a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create figure\n",
    "fig = plt.figure(figsize=(4, 4))\n",
    "plt.xlabel(\"x (feature)\")\n",
    "plt.ylabel(\"y (output)\")\n",
    "plt.title(\"Synthetic data set\")\n",
    "\n",
    "plt.scatter(X, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2699ec5",
   "metadata": {},
   "source": [
    "### 1.2 Optimization problem"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88b00bfa",
   "metadata": {},
   "source": [
    "Given a data matrix of $n$ samples in $p$ dimensions $X \\in \\mathbb{R}^{nxp}$ (in our simulation, $p=1$) and an output vector $\\boldsymbol{y} \\in \\mathbb{R}^{n}$, the objective of a linear regression is to learn a linear function,\n",
    "\n",
    "$$f: \\boldsymbol{x} \\in \\mathbb{R}^p \\mapsto \\beta_0 + \\boldsymbol{\\beta}^\\top \\boldsymbol{x}$$\n",
    "\n",
    "where $\\beta_0 \\in \\mathbb{R}$ and $\\boldsymbol\\beta \\in \\mathbb{R}^p$.\n",
    "\n",
    "This is achieved by __minimizing the empirical risk__, computed with the quadratic error as loss function: $\\frac1n \\sum_{i=1}^n \\left(y_i - f(\\boldsymbol{x}_i)\\right)^2$. \n",
    "\n",
    "The $\\beta_0$ parameter corresponds to an intercept (or bias) term."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03991cc9",
   "metadata": {},
   "source": [
    "#### Bias trick\n",
    "To simplify the equations, we append a column of 1s to the matrix $X$. \n",
    "\n",
    "$$\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "    x_{11} & x_{12} & \\dots & x_{1p} \\\\\n",
    "    x_{21} & x_{22} & \\dots & x_{2p} \\\\\n",
    "    \\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "    x_{n1} & x_{n2} & \\dots & x_{np}\n",
    "\\end{bmatrix}}_{X} \\to\n",
    "\\underbrace{\n",
    "\\begin{bmatrix}\n",
    "    1 & x_{11} & x_{12} & \\dots & x_{1p} \\\\\n",
    "    1 & x_{21} & x_{22} & \\dots & x_{2p} \\\\\n",
    "    \\vdots & \\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "    1 & x_{n1} & x_{n2} & \\dots & x_{np}\n",
    "\\end{bmatrix}}_{\\text{$X$ with bias variable}}\n",
    "$$\n",
    "\n",
    "Now we can simply use \n",
    "$f: \\boldsymbol{x} \\in \\mathbb{R}^{p+1} \\mapsto \\boldsymbol{\\beta}^\\top \\boldsymbol{x}$ with $\\boldsymbol{\\beta} \\in \\mathbb{R}^{p+1}$.\n",
    "\n",
    "This is called the __bias trick.__\n",
    "\n",
    "Let us now compute the numpy array `X_bt` that corresponds to our data matrix `X`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2d78d2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_bt = np.hstack((np.ones((n_samples, 1)), X))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6cfb3cc",
   "metadata": {},
   "source": [
    "We will now work with matrix `X_bt`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "052fa999",
   "metadata": {},
   "source": [
    "#### Closed-form solution\n",
    "\n",
    "The empirical risk minimization can be rewritten as:\n",
    "$$\\boldsymbol{\\beta}^* \\in \\arg\\min \\frac{1}{n}\\sum_{i=1}^n \\left(y_i - \\boldsymbol{\\beta}^\\top \\boldsymbol{x}_i \\right)^2 = \\frac{1}{n}\\left(\\boldsymbol{y} - X \\boldsymbol\\beta\\right)^\\top(\\boldsymbol{y} - X\\boldsymbol\\beta)$$\n",
    "\n",
    "This is a convex optimization problem, so to solve it we can differentiate the function to minimize and set its gradient to zero.\n",
    "\n",
    "We obtain \n",
    "$(X^\\top X) \\boldsymbol\\beta^* =  X^\\top \\boldsymbol{y},$\n",
    "\n",
    "which, __if $X^\\top X$ can be inverted__, leads to a unique solution\n",
    "$\\boldsymbol\\beta^* = (X^\\top X)^{-1} X^\\top \\boldsymbol{y}.$\n",
    "\n",
    "If $X^\\top X$ cannot be inverted, the above problem has an infinity of solutions. They can also be expressed analytically by replacing $(X^\\top X)^{-1}$ with a [pseudoinverse](https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_inverse) of $X^\\top X.$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7869c8c",
   "metadata": {},
   "source": [
    "##### __Question 1__ \n",
    "Use `np.linalg.inv` ([doc](https://numpy.org/doc/stable/reference/generated/numpy.linalg.inv.html)) to compute the vector `beta_star` that is solution to the ordinary least squares in our data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03a9079f",
   "metadata": {},
   "outputs": [],
   "source": [
    "beta_star = # TODO\n",
    "print(beta_star)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9aa95fb4",
   "metadata": {},
   "source": [
    "#### Quality of the model\n",
    "\n",
    "We can now evaluate the quality of the fit of the closed-form solution model. Note that here we are only concerned with fitting the data, and not about generalization."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2cb5ab36",
   "metadata": {},
   "source": [
    "##### __Question 2__ \n",
    "Compute the vector `y_pred` of outputs predicted by the obtained model on the data "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0113f868",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = # TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a35098d0",
   "metadata": {},
   "source": [
    "We can now visualize the regression line we have just learned, and plot the predicted values of $y$ versus the true values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f888714f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create figure\n",
    "fig = plt.figure(figsize=(8, 4))\n",
    "\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.xlabel(\"x (feature)\")\n",
    "plt.ylabel(\"y (output)\")\n",
    "plt.title(\"Synthetic data set with regression line\")\n",
    "plt.scatter(X, y)\n",
    "plt.plot(X, y_pred, lw=2, color=notebook_colors[1])\n",
    "\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.xlabel(\"true output\")\n",
    "plt.ylabel(\"predicted output\")\n",
    "plt.xlim([-200, 200])\n",
    "plt.ylim([-200, 200])\n",
    "plt.title(\"Quality of predictions\")\n",
    "plt.scatter(y, y_pred)\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9787b383",
   "metadata": {},
   "source": [
    "We can also quantify the error with the RMSE, or Root Mean Squared Error, defined by \n",
    "$$\\sqrt{\\sum_{i=1}^n \\frac1n \\left(y_i - f(\\boldsymbol{x}_i)\\right)^2}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d2a893b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def root_mean_squared_error(y_true, y_predicted):\n",
    "    return np.sqrt(np.sum((y_true - y_predicted)**2)/np.size(y))\n",
    "\n",
    "print(f'RMSE: {root_mean_squared_error(y, y_pred):.2f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "127adb03",
   "metadata": {},
   "source": [
    "#### Relationship between noise and RMSE"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c95b49b",
   "metadata": {},
   "source": [
    "##### __Question 3__ \n",
    "Complete the code below to compute the RMSE of a linear regression fitted to data as a function of the value of the `noise` parameter. Comment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1eeb3d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def play_with_noise(noise_value):\n",
    "    \"\"\"\n",
    "    Return the RMSE of a univariate linear regression on data generated with noise parameter noise_value.\n",
    "    \"\"\"\n",
    "    X, y = datasets.make_regression(n_samples=n_samples, n_features=n_features, n_informative=1, noise=noise_value, random_state=37)\n",
    "    X_bt = np.hstack((np.ones((n_samples, 1)), X))\n",
    "    beta_star = # TODO\n",
    "    y_pred = # TODO\n",
    "    return(root_mean_squared_error(y, y_pred))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d47f64c",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(4, 4))\n",
    "plt.xlabel(\"Noise in the data\")\n",
    "plt.ylabel(\"RMSE of the linear regression\")\n",
    "\n",
    "noise_values = # TODO\n",
    "plt.scatter(noise_values, [play_with_noise(x) for x in noise_values])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2895103c-25a0-4631-98f4-f52f0ed118b3",
   "metadata": {},
   "source": [
    "__Comment:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f41e1382",
   "metadata": {},
   "source": [
    "### 1.3 Gradient descent"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27522e69",
   "metadata": {},
   "source": [
    "It will not always be possible to use a closed-form solution of the empirical risk minimization problem.\n",
    "\n",
    "In the case of linear regression, there are two reasons why we might not want to use the closed-form solution:\n",
    "* if $X^\\top X$ cannot be inverted (which happens if some of the features are correlated, or if the number of features is larger than the number of columns);\n",
    "* if $X^\\top X$ is invertible but of large dimension, ie. when the number of features is large. Indeed, matrix inversion is _cubic_ in the dimension of the matrix.\n",
    "\n",
    "When the empirical risk minimization problem is convex, __gradient descent__ finds a solution by progressing iteratively in the direction opposite the gradient. Specifically, at each iteration $t$, the vector $\\boldsymbol\\beta$ is updated by:\n",
    "$$\\boldsymbol\\beta^{(t+1)} = \\boldsymbol\\beta^{(t)} - \\alpha \\nabla_{\\boldsymbol\\beta} J(\\boldsymbol\\beta^{(t)}),$$\n",
    "where $J: \\mathbb{R}^d \\rightarrow \\mathbb{R} $ is the function to minimize.\n",
    "\n",
    "In the case of the ordinary least squares, $d=p+1$ and $J(\\boldsymbol\\beta) = \\frac{1}{n}\\left(\\boldsymbol{y} - X \\boldsymbol\\beta\\right)^\\top(\\boldsymbol{y} - X\\boldsymbol\\beta).$\n",
    "\n",
    "$\\alpha$ is the step size of the gradient descent, also called the __learning rate.__\n",
    "\n",
    "The gradient descent is stopped either when a preset maximum number of iterations is reached, or when the norm of the gradient is less than a preset error level (called __tolerance__)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e41f1a3",
   "metadata": {},
   "source": [
    "#### Gradient of the empirical risk of the linear regression\n",
    "For ordinary least squares, the gradient of the function to minimize is:\n",
    "\n",
    "$$\\nabla_{\\boldsymbol\\beta} J(\\boldsymbol\\beta) = \\frac1n X^\\top \\left(X \\boldsymbol\\beta - y \\right).$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c422160",
   "metadata": {},
   "source": [
    "##### __Question 4__ \n",
    "Complete the code below to evaluate the gradient of the empirical risk of the linear regression:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fefb163c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def evaluate_gradient(X, y, b_vector):\n",
    "    return # TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf2b0a18",
   "metadata": {},
   "source": [
    "Let us check the value of the gradient at our solution `beta_star`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80dae587",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"In the closed-form solution beta_star, the norm of the gradient is {np.linalg.norm(evaluate_gradient(X_bt, y, beta_star)):.2e}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba46fe35",
   "metadata": {},
   "source": [
    "##### __Question 5__ \n",
    "Does this value match your expectations?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bc1395b-5eb5-4130-974f-1861ea931cd1",
   "metadata": {},
   "source": [
    "__Answer:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59bda0a0",
   "metadata": {},
   "source": [
    "#### Batch gradient descent\n",
    "\n",
    "In batch gradient descent, we use the entire data set at each iteration."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f12d5f1a-b4d1-41b9-9b52-6ff1daf61197",
   "metadata": {},
   "source": [
    "##### __Question 6__ \n",
    "Complete the `batch_gradient_descent` function below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d9ba553",
   "metadata": {},
   "outputs": [],
   "source": [
    "def squared_loss(X, y, b_vector):\n",
    "    \"\"\"\n",
    "    Squared loss of a linear regression.\n",
    "    Parameters:\n",
    "    X: (n_samples, n_features) numpy array\n",
    "        The data matrix.\n",
    "        \n",
    "    y: (n_samples, ) numpy array \n",
    "        The output/target vector.\n",
    "        \n",
    "    b_vector: (n_features, ) numpy array\n",
    "        The weight vector of the linear regression.    \n",
    "    \"\"\"\n",
    "    return(np.sum((y - X.dot(b_vector))**2)/np.size(y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e7aaa6f",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-03-27T17:38:54.314594Z",
     "start_time": "2019-03-27T17:38:54.288955Z"
    }
   },
   "outputs": [],
   "source": [
    "def batch_gradient_descent(X, y, learning_rate=1e-1, max_iters=30, tol=1e-2):\n",
    "    \"\"\"\n",
    "    Batch gradient descent procedure for a linear regression.\n",
    "    Parameters:\n",
    "    X: (n_samples, n_features) numpy array\n",
    "        The data matrix.\n",
    "        \n",
    "    y: (n_samples, ) numpy array \n",
    "        The output/target vector.\n",
    "        \n",
    "    learning_rate: float\n",
    "        The learning rate of the gradient descent.\n",
    "        \n",
    "    max_iters: int\n",
    "        The maximum number of iterations of the gradient descent.\n",
    "        \n",
    "    tol: float\n",
    "        The tolerance of the gradient descent:\n",
    "        the gradient descent will stop when the norm of the gradient is smaller than the tolerance.\n",
    "    \"\"\"\n",
    "    # Random initialisation of beta\n",
    "    beta_current = np.random.rand(# TODO)\n",
    "    \n",
    "    # Keep track of the values of the loss and of beta at each iteration\n",
    "    losses = [squared_loss(X, y, beta_current)]\n",
    "    betas = [beta_current.copy()]\n",
    "\n",
    "    for iter_idx in range(1, max_iters+1):        \n",
    "        # compute the gradient\n",
    "        gradient = # TODO\n",
    "        \n",
    "        # update beta_current \n",
    "        beta_current = # TODO\n",
    "        \n",
    "        # compute the loss and append to losses\n",
    "        losses.append(squared_loss(X, y, beta_current))\n",
    "        \n",
    "        # append beta_current to losses\n",
    "        betas.append(beta_current.copy())\n",
    "\n",
    "        # check whether tolerance is reached\n",
    "        # (you can use np.linalg.norm to compute the norm of a vector)\n",
    "        if # TODO: \n",
    "            print(f\"Gradient descent stopped after {iter_idx} iterations\")\n",
    "            break\n",
    "\n",
    "    return(np.array(losses), np.array(betas))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f49e97f",
   "metadata": {},
   "source": [
    "We can now run the batch gradient descent on the data and observe how beta and the loss evolve:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de52deb1",
   "metadata": {},
   "outputs": [],
   "source": [
    "tolerance = 1e-4\n",
    "iterations = 35\n",
    "lr = 2e-2\n",
    "losses, betas = batch_gradient_descent(X_bt, y, learning_rate=lr, max_iters=iterations, tol=tolerance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa3565ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(np.arange(len(losses)), losses, 'o-')\n",
    "\n",
    "beta_star_loss = squared_loss(X_bt, y, beta_star)\n",
    "plt.plot([0, iterations], [beta_star_loss, beta_star_loss], 'k-', label=\"Optimal loss\")\n",
    "\n",
    "plt.xlabel(\"number of iterations\")\n",
    "plt.ylabel(\"empirical risk\")\n",
    "plt.title(\"Batch gradient descent\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed879ce9",
   "metadata": {},
   "source": [
    "##### __Question 7__ \n",
    "Has your gradient descent stopped because you've reached the maximum number of parameters or the tolerance? If yes, what parameter(s) do you think you could change to actually go below tolerance? Experiment below and comment."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b16f4a8f-56c2-4e57-92cc-aaec0ca45679",
   "metadata": {},
   "source": [
    "__Comment:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e87adf1",
   "metadata": {},
   "source": [
    "#### Quality of the model\n",
    "\n",
    "We can now evaluate the quality of the fit of the model learned by our gradient descent, and compare it to the model we obtained with the closed-form solution."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b709076-bead-450c-87f7-d27cec867422",
   "metadata": {},
   "source": [
    "##### __Question 8__\n",
    "Use the code below to look at the quality of a gradient descent that hasn't converged yet. Comment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fced76f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "tolerance = 1e-4\n",
    "iterations = 35\n",
    "lr = 2e-2\n",
    "losses, betas = batch_gradient_descent(X_bt, y, learning_rate=lr, max_iters=iterations, tol=tolerance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "888827a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred_gd = X_bt.dot(betas[-1])\n",
    "# create figure\n",
    "fig = plt.figure(figsize=(10, 4))\n",
    "\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.xlabel(\"x (feature)\")\n",
    "plt.ylabel(\"y (output)\")\n",
    "plt.title(\"Synthetic data set with regression line\")\n",
    "plt.scatter(X, y, alpha=0.5)\n",
    "plt.plot(X, y_pred, lw=2, color=notebook_colors[1], label='Closed-form solution')\n",
    "plt.plot(X, y_pred_gd, lw=2, color=notebook_colors[2], label='Gradient descent solution')\n",
    "plt.legend()\n",
    "\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.xlabel(\"true output\")\n",
    "plt.ylabel(\"predicted output\")\n",
    "plt.xlim([-200, 200])\n",
    "plt.ylim([-200, 200])\n",
    "plt.title(\"Quality of predictions\")\n",
    "plt.scatter(y, y_pred, label=\"Closed-form solution\", alpha=0.5)\n",
    "plt.scatter(y, y_pred_gd, label=\"Gradient descent solution\", alpha=0.5)\n",
    "plt.legend()\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95770cb1-78a6-4c21-913b-0b83721d62a4",
   "metadata": {},
   "source": [
    "__Comment:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "165aefbd",
   "metadata": {},
   "source": [
    "#### Gradient descent visualization\n",
    "\n",
    "You can use the following code to visualize the gradient descent. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5531b88f",
   "metadata": {},
   "outputs": [],
   "source": [
    "tolerance = 1e-4\n",
    "iterations = 35\n",
    "lr = 5e-1\n",
    "losses, betas = batch_gradient_descent(X_bt, y, learning_rate=lr, max_iters=iterations, tol=tolerance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bfe1c187",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-03-27T17:39:01.369834Z",
     "start_time": "2019-03-27T17:38:56.002975Z"
    }
   },
   "outputs": [],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "# Create mesh grid\n",
    "ms = np.linspace(beta_star[0] - 20 , beta_star[0] + 20, 20)\n",
    "bs = np.linspace(beta_star[1] - 80 , beta_star[1] + 80, 20)\n",
    "M, B = np.meshgrid(ms, bs)\n",
    "\n",
    "# Compute the squared loss in each point of the mesh grid\n",
    "zs = np.array([squared_loss(X_bt, y, theta)\n",
    "               for theta in zip(np.ravel(M), np.ravel(B))])\n",
    "Z = zs.reshape(M.shape)\n",
    "\n",
    "# create 3D axis object\n",
    "fig = plt.figure(figsize=(10, 6))\n",
    "ax = fig.add_subplot(1, 1, 1, projection='3d', xlabel='Intercept',\n",
    "                     ylabel='Slope', zlabel='Loss')\n",
    "\n",
    "# plot loss hypersurface and contours\n",
    "ax.plot_surface(M, B, Z, rstride=1, cstride=1, alpha=0.3)\n",
    "ax.contour(M, B, Z, 20, alpha=0.5, offset=0, stride=30)\n",
    "\n",
    "# plot gradient descent \n",
    "ax.plot(betas[:, 0], betas[:, 1], losses, '.-')\n",
    "\n",
    "# set viewpoint\n",
    "ax.view_init(elev=10, azim=30)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1720a56",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### Newton-Raphson \n",
    "\n",
    "When the function $J$ to optimize is twice differentiable, the Newton-Raphson procedure consists in using $\\| \\nabla^2 J(\\boldsymbol\\beta^{(t)}) \\|$ as learning rate at update $t$ :\n",
    "\n",
    "$$\\boldsymbol\\beta^{(t+1)} = \\boldsymbol\\beta^{(t)} - \\|\\nabla^2 J(\\boldsymbol\\beta^{(t)})\\|.$$\n",
    "\n",
    "In our case, $\\|\\nabla^2 J(\\boldsymbol\\beta)\\| = \\|X^\\top X\\|= \\sigma_{\\max}$, where $\\sigma_{\\max}$ is the largest eigenvalue of $X^\\top X$.\n",
    "\n",
    "The learning rate suggested by this procedure is constant and equal to $\\frac1{\\sigma_{\\max}}.$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba4380e1",
   "metadata": {
    "tags": []
   },
   "source": [
    "##### __Question 9__\n",
    "Define `learning_rates_list` below to compare several learning rates, including $\\frac1{\\sigma_{\\max}}$ and rates both lower and smaller than $\\frac1{\\sigma_{\\max}}$. Comment."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20a79f13",
   "metadata": {},
   "source": [
    "__Remark:__ Note that you can replace `plot` with `semilogy` to use a logarithmic scale on the y-axis and better see differences in the lower values of the loss."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32551df3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the eigenvalues of X'X\n",
    "eigens = np.linalg.eigvals(X_bt.T.dot(X_bt)/X_bt.shape[0])\n",
    "sigma_max = np.max(eigens)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8317993",
   "metadata": {},
   "outputs": [],
   "source": [
    "learning_rates_list = # TODO\n",
    "for lr in learning_rates_list:\n",
    "    # Run the batch gradient descent\n",
    "    losses, betas = batch_gradient_descent(X_bt, y, learning_rate=lr, max_iters=iterations, tol=tolerance)\n",
    "    \n",
    "    # Plot loss\n",
    "    plt.plot(np.arange(len(losses)), losses, 'o-', label=\"Learning rate = %.2e\" % lr)\n",
    "    \n",
    "beta_star_loss = squared_loss(X_bt, y, beta_star)\n",
    "plt.plot([0, iterations], [beta_star_loss, beta_star_loss], 'k-', label=\"Minimum loss\")\n",
    "\n",
    "plt.xlabel(\"number of iterations\")\n",
    "plt.ylabel(\"empirical risk\")\n",
    "plt.title(\"Batch gradient descent\")\n",
    "plt.legend(loc=(1.01, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d2a0778-3fe0-4a10-8e83-f33246e26713",
   "metadata": {},
   "source": [
    "__Comment:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17154009",
   "metadata": {},
   "source": [
    "##### __Question 10__\n",
    "How does this relate to the suggested learning rate for the case where the gradient of $J$ is L-Lipschitz (seen in class)?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2a53251-ec10-4cdd-bcf4-a37e4ba2d911",
   "metadata": {},
   "source": [
    "__Answer:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fdd961e",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 1.4 Minibatch Gradient Descent\n",
    "\n",
    "In larger applications, computing the full gradient can be expensive. __Minibatch gradient descent__ uses a subsample of the data of size $m \\ll n$ (called a __batch__) at each iteration. If the dataset is large, this is often sufficient to make an accurate descent step.\n",
    "\n",
    "Computing the gradient therefore goes from a complexity of $\\mathcal{O}(np)$ to $\\mathcal{O}(mp)$. In the training of deep neural networks, minibatch gradient descent (and its variants) is overwhelmingly the most popular approach.\n",
    "\n",
    "The most straightforward strategy to build batches is to iterate through the data and take the samples slice of $m$ samples by slice of $m$ samples. In our case, at the iteration that uses batch $k$, we compute the gradient as  \n",
    "$$\\nabla_{\\boldsymbol\\beta} J(\\beta) = \\frac1m X_{k:k+m}^\\top \\left(X_{k:k+m} \\boldsymbol\\beta - y_{k:k+m} \\right).$$\n",
    "\n",
    "An __epoch__ corresponds to having seen the entire data set once."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9c21606",
   "metadata": {},
   "source": [
    "##### __Question 11__\n",
    "Complete the code below to implement minibatch gradient descent for the linear regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a750368",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-03-27T17:39:08.231736Z",
     "start_time": "2019-03-27T17:39:08.198510Z"
    }
   },
   "outputs": [],
   "source": [
    "def minibatch_gradient_descent(X, y, batch_size=10, learning_rate=1e-1, max_iters=100, tol=1e-5):    \n",
    "    \"\"\"\n",
    "    Minibatch gradient descent procedure for a linear regression.\n",
    "    Parameters:\n",
    "    X: (n_samples, n_features) numpy array\n",
    "        The data matrix.\n",
    "        \n",
    "    y: (n_samples, ) numpy array \n",
    "        The output/target vector.\n",
    "        \n",
    "    batch_size: int\n",
    "        The number of samples in a batch.\n",
    "        \n",
    "    learning_rate: float\n",
    "        The learning rate of the gradient descent.\n",
    "        \n",
    "    max_iters: int\n",
    "        The maximum number of iterations of the gradient descent.\n",
    "        \n",
    "    tol: float\n",
    "        The tolerance of the gradient descent\n",
    "    \"\"\"\n",
    "    \n",
    "    # Random initialisation of beta\n",
    "    beta_current = np.random.rand(# TODO)\n",
    "    \n",
    "    # Keep track of the values of the loss and of beta at each epoch\n",
    "    losses = [squared_loss(X, y, beta_current)]\n",
    "\n",
    "    # Find total number of samples\n",
    "    num_samples = X.shape[0]\n",
    "    \n",
    "    # Initialize index of start of batch\n",
    "    idx_start_of_batch = 0\n",
    "    for iter_idx in range(1, max_iters+1):\n",
    "        # find index of end of batch\n",
    "        idx_end_of_batch = min(idx_start_of_batch + batch_size, num_samples)\n",
    "        \n",
    "        # construct batch\n",
    "        X_batch = # TODO\n",
    "        y_batch = # TODO\n",
    "        \n",
    "        # find index of start of next batch \n",
    "        idx_start_of_batch = idx_end_of_batch if idx_end_of_batch < num_samples else 0\n",
    "        \n",
    "        # compute the gradient\n",
    "        gradient = # TODO\n",
    "        \n",
    "        # update beta_current\n",
    "        beta_current = # TODO\n",
    "\n",
    "        # check whether tolerance is reached\n",
    "        if # TODO: \n",
    "            print(\"Gradient descent stopped after {iter_idx} iterations\")\n",
    "            break\n",
    "    \n",
    "        # compute the loss and append to losses\n",
    "        losses.append(squared_loss(X, y, beta_current))\n",
    "\n",
    "    return(np.array(losses), beta_current)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2af14a29",
   "metadata": {},
   "source": [
    "#### Impact of the learning rate\n",
    "\n",
    "We will now plot the evolution of the loss as a function of the number of iterations. Larger dots (circles) will mark the epochs."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21c22e31",
   "metadata": {},
   "source": [
    "##### __Question 12__\n",
    "Complete the code below to run the minibatch gradient descent for a few different choices of learning rates. Comment on what you observe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e5dc22f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# List of values of learning rates to test\n",
    "lr_list = # TODO\n",
    "\n",
    "tolerance = 1e-4\n",
    "iterations = 20\n",
    "bs = 20\n",
    "\n",
    "# Loop over the learning rates\n",
    "for lr in lr_list:\n",
    "    # Run the minibatch gradient descent\n",
    "    losses, beta_mgd = minibatch_gradient_descent(X_bt, y, batch_size=bs, learning_rate=lr, max_iters=iterations, tol=tolerance)\n",
    "    \n",
    "    # Plot loss\n",
    "    plt.plot(np.arange(len(losses)), losses, label=\"Learning rate = %.2e\" % lr)\n",
    "    # Mark loss after each epoch\n",
    "    n_iterations_per_epoch = int(np.ceil(n_samples/bs))\n",
    "    n_epochs = int(iterations/n_iterations_per_epoch)\n",
    "    plt.scatter([int(np.ceil(n_samples/bs))*idx for idx in range(n_epochs)], \n",
    "                [losses[int(np.ceil(n_samples/bs))*idx] for idx in range(n_epochs)])\n",
    "    \n",
    "plt.plot([0, iterations], [beta_star_loss, beta_star_loss], '-k', label=\"Minimum loss\")\n",
    "\n",
    "plt.xlabel(\"number of iterations\")\n",
    "plt.ylabel(\"empirical risk\")\n",
    "plt.title(\"Minibatch gradient descent\")\n",
    "plt.legend(loc=(1.01, 0.))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70abf48b-b733-4f76-8385-79c0a025f8f0",
   "metadata": {},
   "source": [
    "__Comment:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0f83e6a",
   "metadata": {},
   "source": [
    "#### Impact of the batch size"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38f8a5c7",
   "metadata": {},
   "source": [
    "##### __Question 13__\n",
    "Now complete the code below to set the learning rate to a fixed value, and run the minibatch gradient descent for a few different choices of batch sizes. Comment on what you observe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55c42388",
   "metadata": {},
   "outputs": [],
   "source": [
    "lr = 1e-1\n",
    "tolerance = 1e-4\n",
    "iterations = 20\n",
    "\n",
    "batch_size_list = # TODO\n",
    "\n",
    "# Loop over the learning rates\n",
    "for bs in batch_size_list:\n",
    "    # Run the minibatch gradient descent\n",
    "    losses, beta_mgd = minibatch_gradient_descent(X_bt, y, batch_size=bs, learning_rate=lr, max_iters=iterations, tol=tolerance)\n",
    "    \n",
    "    # Plot loss\n",
    "    plt.plot(np.arange(len(losses)), losses, label=\"Batch size = %d\" % bs)\n",
    "    # Mark loss after each epoch\n",
    "    n_iterations_per_epoch = int(np.ceil(n_samples/bs))\n",
    "    n_epochs = int(iterations/n_iterations_per_epoch)\n",
    "    plt.scatter([n_iterations_per_epoch*idx for idx in range(n_epochs)], \n",
    "                [losses[int(n_iterations_per_epoch)*idx] for idx in range(n_epochs)])\n",
    "    \n",
    "plt.plot([0, iterations], [beta_star_loss, beta_star_loss], '-k', label=\"Minimum loss\")\n",
    "\n",
    "plt.xlabel(\"number of iterations\")\n",
    "plt.ylabel(\"empirical risk\")\n",
    "plt.title(\"Minibatch gradient descent\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e273cdbb-fc2a-47f0-b1b5-0eaf5c9d66a6",
   "metadata": {},
   "source": [
    "__Comment:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07b4eb2b",
   "metadata": {},
   "source": [
    "### 1.5 Stochastic Gradient Descent"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15e24ab0",
   "metadata": {},
   "source": [
    "Stochastic gradient descent is a special case of minibatch gradient descent:\n",
    "* the batch size is set to 1\n",
    "* one cycles randomly through the samples.\n",
    "\n",
    "In the case of neural networks, which have a non-convex optimization problem, stochasticity also may help to avoid ending in local minima."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20010a4a",
   "metadata": {},
   "source": [
    "##### __Question 14__\n",
    "Complete the implementation of the stochastic gradient descent below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f6bc143",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-03-27T17:39:08.231736Z",
     "start_time": "2019-03-27T17:39:08.198510Z"
    }
   },
   "outputs": [],
   "source": [
    "def stochastic_gradient_descent(X, y, learning_rate=1e-1, max_iters=20):    \n",
    "    \"\"\"\n",
    "    Stochastic gradient descent procedure for a linear regression.\n",
    "    Parameters:\n",
    "    X: (n_samples, n_features) numpy array\n",
    "        The data matrix.\n",
    "        \n",
    "    y: (n_samples, ) numpy array \n",
    "        The output/target vector.\n",
    "        \n",
    "    learning_rate: float\n",
    "        The learning rate of the gradient descent.\n",
    "        \n",
    "    max_iters: int\n",
    "        The maximum number of iters of the gradient descent.\n",
    "        \n",
    "    \"\"\"\n",
    "    \n",
    "    # Random initialisation of beta\n",
    "    beta_current = # TODO\n",
    "    \n",
    "    # Keep track of the values of the loss and of beta at each epoch\n",
    "    losses = [squared_loss(X, y, beta_current)]\n",
    "    \n",
    "    # Vector of sample indices\n",
    "    num_samples = X.shape[0]\n",
    "    samples_indices = np.arange(num_samples)\n",
    "    \n",
    "    # Determine number of epochs\n",
    "    n_epochs = int(ceil(max_iters/num_samples))\n",
    "    \n",
    "    for epoch in range(n_epochs):\n",
    "        # Shuffle the samples\n",
    "        np.random.shuffle(samples_indices)\n",
    "        \n",
    "        # The last epoch may not be full\n",
    "        samples_indices = samples_indices[:min(num_samples, (max_iters-epoch*num_samples))]\n",
    "        \n",
    "        # Loop over the shuffled samples one by one\n",
    "        for idx, sample_idx in enumerate(samples_indices):\n",
    "            # update beta_current\n",
    "            beta_current = # TODO\n",
    "               \n",
    "            # compute the loss and append to losses \n",
    "            losses.append(squared_loss(X, y, beta_current))\n",
    "\n",
    "    return(np.array(losses), beta_current)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88094410",
   "metadata": {
    "tags": []
   },
   "source": [
    "##### __Question 15__\n",
    "Why aren't we using the norm of the gradient as stopping criterion any longer? "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05bbae25-9c8e-4dd3-a3ca-15fc4252350f",
   "metadata": {},
   "source": [
    "__Answer:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "351f41b7",
   "metadata": {},
   "source": [
    "##### __Question 16__\n",
    "Complete the code below by chosing a learning rate `lr` and observe how the loss evolves as a function of the number of iterations. Comment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd306272",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run stochastic gradient descent\n",
    "iterations = 200\n",
    "lr = # TODO\n",
    "losses, betas = stochastic_gradient_descent(X_bt, y, learning_rate=lr, max_iters=iterations)\n",
    "\n",
    "# Plot loss\n",
    "plt.plot(np.arange(len(losses)), losses, label=\"Stochastic gradient descent\")\n",
    "\n",
    "# Mark loss after each epoch\n",
    "n_epochs = int(iterations/n_samples)\n",
    "plt.scatter([n_samples*idx for idx in range(n_epochs+1)], [losses[n_samples*idx] for idx in range(n_epochs+1)])\n",
    "\n",
    "plt.xlabel(\"number of iterations\")\n",
    "plt.ylabel(\"empirical risk\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3547745e-5bce-42a3-9e54-500e5f8f96a3",
   "metadata": {},
   "source": [
    "__Comment:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00c207c8",
   "metadata": {},
   "source": [
    "#### Comparing the three gradient descents"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a10010c",
   "metadata": {},
   "source": [
    "##### __Question 17__\n",
    "Complete and run the following code to compare the three gradient descents by looking at their loss as a function of the number of iterations. Comment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70cc575c",
   "metadata": {},
   "outputs": [],
   "source": [
    "tolerance = 1e-4\n",
    "iterations = # TODO\n",
    "\n",
    "# Run batch gradient descent with Newton-Raphson learning rate\n",
    "lr = 1./sigma_max\n",
    "losses, betas = batch_gradient_descent(X_bt, y, learning_rate=lr, max_iters=iterations, tol=tolerance)\n",
    "# Plot loss\n",
    "plt.plot(np.arange(len(losses)), losses, label=\"Batch gradient descent\")\n",
    "\n",
    "# Run minibatch gradient descent\n",
    "bs = # TODO\n",
    "lr = # TODO\n",
    "losses, beta_mgd = minibatch_gradient_descent(X_bt, y, batch_size=bs, learning_rate=lr, max_iters=iterations, tol=tolerance)\n",
    "# Plot loss\n",
    "plt.plot(np.arange(len(losses)), losses, label=\"Minibatch gradient descent\")\n",
    "    \n",
    "# Run stochastic gradient descent\n",
    "losses, beta_sgd = stochastic_gradient_descent(X_bt, y, learning_rate=lr,  max_iters=iterations)\n",
    "# Plot loss\n",
    "plt.plot(np.arange(len(losses)), losses, label=\"Stochastic gradient descent\")\n",
    "\n",
    "# Plot closed-form loss\n",
    "plt.plot([0, iterations], [beta_star_loss, beta_star_loss], '-k', label=\"Minimum loss\")\n",
    "\n",
    "plt.xlabel(\"number of iterations\")\n",
    "plt.ylabel(\"empirical risk\")\n",
    "plt.title(\"Comparing gradient descents\")\n",
    "plt.legend(loc=(1.01, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "153c7e48-3afa-4936-85e2-eddb2d2c7719",
   "metadata": {},
   "source": [
    "__Comment:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54b6dc59",
   "metadata": {},
   "source": [
    "Another way to look at this is to look at the loss as a function of the *number of epochs*, that is to say the number of times the algorithm looked at the entire data set.\n",
    "\n",
    "##### __Question 18__\n",
    "Complete and run the following code to compare the three gradient descents by looking at their loss as a function of the number of epochs. Comment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6193f103",
   "metadata": {},
   "outputs": [],
   "source": [
    "tolerance = 1e-4\n",
    "epochs = # TODO \n",
    "\n",
    "# Run batch gradient descent\n",
    "lr = 1./sigma_max\n",
    "losses, betas = batch_gradient_descent(X_bt, y, learning_rate=lr, max_iters=epochs, tol=tolerance)\n",
    "# Plot loss\n",
    "plt.semilogy(np.arange(len(losses)), losses, label=\"Batch gradient descent\")\n",
    "\n",
    "# Run minibatch gradient descent\n",
    "bs = # TODO use the same value as for the comparison per number of iterations\n",
    "lr = # TODO use the same value as for the comparison per number of iterations\n",
    "n_iterations_per_epoch = int(np.ceil(n_samples/bs))\n",
    "iterations = epochs * n_iterations_per_epoch\n",
    "losses, beta_mgd = minibatch_gradient_descent(X_bt, y, batch_size=bs, learning_rate=lr, max_iters=iterations, tol=tolerance)\n",
    "# Plot loss\n",
    "plt.plot(np.arange(len(losses))/n_iterations_per_epoch, losses, label=\"Minibatch gradient descent\")\n",
    "    \n",
    "# Run stochastic gradient descent\n",
    "iterations = epochs * n_samples\n",
    "losses, beta_sgd = stochastic_gradient_descent(X_bt, y, learning_rate=lr,  max_iters=iterations)\n",
    "# Plot loss\n",
    "plt.plot(np.arange(len(losses))/n_samples, losses, label=\"Stochastic gradient descent\")\n",
    "\n",
    "# Plot closed-form loss\n",
    "plt.plot([0, epochs], [beta_star_loss, beta_star_loss], '-k', label=\"Minimum loss\")\n",
    "\n",
    "plt.xlabel(\"number of epochs\")\n",
    "plt.ylabel(\"empirical risk\")\n",
    "plt.title(\"Comparing gradient descents\")\n",
    "plt.legend(loc=(1.01, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe740b7d-3e26-4924-b642-f65d784fcdee",
   "metadata": {},
   "source": [
    "__Comment:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ddbcb8ac",
   "metadata": {},
   "source": [
    "### 1.6 Multivariate linear regression\n",
    "\n",
    "The code we have written so far can be applied to data sets with more than one feature."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c62b7ab",
   "metadata": {},
   "source": [
    "##### __Question 19__\n",
    "Complete the following code to compare the three types of gradient descent for a larger dataset, with 30 features. Play with the learning rates and batch sizes. Note and discuss what you observe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df760be8",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_samples = 1000\n",
    "n_features = 30\n",
    "X, y = datasets.make_regression(n_samples=n_samples, n_features=n_features, n_informative=n_features, noise=50, random_state=37)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8544f96",
   "metadata": {},
   "outputs": [],
   "source": [
    "# bias trick\n",
    "X_bt = np.hstack((np.ones((n_samples, 1)), X))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f01d109e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# closed-form solution\n",
    "beta_star = # TODO\n",
    "beta_star_loss = # TODO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f130e589",
   "metadata": {},
   "outputs": [],
   "source": [
    "# optimal learning rate for batch\n",
    "lr_optim = # TODO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d493859",
   "metadata": {},
   "outputs": [],
   "source": [
    "tolerance = 1e-4\n",
    "epochs = 10 # TODO change this value if you wish\n",
    "\n",
    "# Run batch gradient descent\n",
    "lr = 1./sigma_max\n",
    "losses, betas = batch_gradient_descent(X_bt, y, learning_rate=lr, max_iters=epochs, tol=tolerance)\n",
    "# Plot loss\n",
    "plt.plot(np.arange(len(losses)), losses, label=\"Batch gradient descent\")\n",
    "\n",
    "# Run minibatch gradient descent\n",
    "bs = 20 # TODO try changing this value \n",
    "lr = 5e-2 # TODO try changing this value \n",
    "n_iterations_per_epoch = int(np.ceil(n_samples/bs))\n",
    "iterations = epochs * n_iterations_per_epoch\n",
    "losses, beta_mgd = minibatch_gradient_descent(X_bt, y, batch_size=bs, learning_rate=lr, max_iters=iterations, tol=tolerance)\n",
    "# Plot loss\n",
    "plt.plot(np.arange(len(losses))/n_iterations_per_epoch, losses, label=\"Minibatch gradient descent\")\n",
    "    \n",
    "# Run stochastic gradient descent\n",
    "lr = 5e-3 # TODO try changing this value \n",
    "iterations = epochs * n_samples\n",
    "losses, beta_sgd = stochastic_gradient_descent(X_bt, y, learning_rate=lr,  max_iters=iterations)\n",
    "# Plot loss\n",
    "plt.plot(np.arange(len(losses))/n_samples, losses, label=\"Stochastic gradient descent\")\n",
    "\n",
    "# Plot closed-form loss\n",
    "plt.plot([0, epochs], [beta_star_loss, beta_star_loss], '-k', label=\"Minimum loss\")\n",
    "\n",
    "plt.xlabel(\"number of epochs\")\n",
    "plt.ylabel(\"empirical risk\")\n",
    "plt.title(\"Comparing gradient descents\")\n",
    "plt.legend(loc=(1.01, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b71910a-7538-42a6-9b6d-2e97d7f73661",
   "metadata": {},
   "source": [
    "__Comment:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "138b6e2e",
   "metadata": {},
   "source": [
    "### 1.7 SGDRegressor in scikit-learn\n",
    "\n",
    "The `SGDRegressor` class of `scikit-learn` ([doc](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.SGDRegressor.html#sklearn.linear_model.SGDRegressor)) makes it possible to learn linear regression models with various losses and regularizers, using stochastic gradient descent. \n",
    "\n",
    "We can use it here instead of our hand-written procedures.\n",
    "\n",
    "##### __Question 20__\n",
    "Complete the code below to use `SGDRegressor` and compare it to our procedures. Comment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7415b126",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import linear_model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "400ae616",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize an instance of SGDRegressor\n",
    "regressor = linear_model.SGDRegressor(loss='squared_loss', \n",
    "                                      penalty='l2', alpha=1e-5, # you have to add a regularizer, so we use a very small regularization parameter\n",
    "                                      max_iter=20, tol=1e-4, fit_intercept=False, \n",
    "                                      learning_rate='constant', eta0=1e-3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90bb9cd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Learn\n",
    "regressor.fit(X_bt, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5da6c47c",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"SGDRegressor stopped after {regressor.n_iter_} epochs.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa62cc63-aed9-4d60-bac3-6dbfb3031db8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compare the loss of the model learned by SGDRegressor with the minimum loss\n",
    "print(f\"Closed-form loss: {beta_star_loss:.2e}\")\n",
    "print(f\"SGDRegressor loss: {squared_loss(X_bt, y, regressor.coef_):.2e}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d7bf004-6a11-4fb3-a38a-1360a2fcbf17",
   "metadata": {},
   "source": [
    "__Comment:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "293cb689",
   "metadata": {},
   "source": [
    "## 2. Logistic regression\n",
    "\n",
    "We will now work on binary classification with logistic regression."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d28509a3",
   "metadata": {},
   "source": [
    "### 2.1 Generating data "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3da49332",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f858b43",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_samples = 1000\n",
    "n_features = 100\n",
    "X, y = datasets.make_classification(n_samples=n_samples, n_features=n_features, n_informative=int(n_features/10), random_state=59)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d06deb1",
   "metadata": {},
   "outputs": [],
   "source": [
    "X.shape, y.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b82b95c1",
   "metadata": {},
   "source": [
    "Let us standardize the data (each feature will be rescaled to have a mean of 0 and a standard deviation of 1):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8539674d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20d8e730",
   "metadata": {},
   "outputs": [],
   "source": [
    "scaler = preprocessing.StandardScaler()\n",
    "X_scaled = scaler.fit_transform(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75ac034f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bias trick\n",
    "X_bt = np.hstack((np.ones((n_samples, 1)), X_scaled))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f96bd647",
   "metadata": {},
   "source": [
    "### 2.2 Optimization problem\n",
    "\n",
    "In logistic regression, the model is \n",
    "$\\mathbb{P}(Y = 1|X) = \\sigma(\\boldsymbol\\beta^\\top \\boldsymbol{x}),$ with $Y$ being a binary random variable from which the outputs are drawn, and $X$ a $p$-dimensional random vector from which the features are drawn.\n",
    "$\\sigma: u \\mapsto \\frac1{1+e^{-u}}$ is the __sigmoid__ or __logistic__ function.\n",
    "\n",
    "The loss function is the __logistic loss__: $$\\mathcal{L}(y, f(\\boldsymbol{x})) = - y \\log(f(\\boldsymbol{x})- (1-y) \\log(1-f(\\boldsymbol{x})) $$\n",
    "\n",
    "Hence our goal is to minimize\n",
    "$$J(\\boldsymbol\\beta) = \\frac1n \\sum_{i=1}^n - y_i \\log(\\sigma(\\boldsymbol\\beta^\\top \\boldsymbol{x}_i))- (1-y_i) \\log(1-\\sigma(\\boldsymbol\\beta^\\top \\boldsymbol{x}_i))  $$\n",
    "\n",
    "Because $\\sigma^\\prime(u) = \\sigma(u) (1-\\sigma(u))$, we can easily obtain the following gradient for $J$:\n",
    "$$\\nabla J(\\boldsymbol\\beta) = - \\frac1n \\sum_{i=1}^n \\left(y_i - \\sigma(\\boldsymbol\\beta^\\top \\boldsymbol{x}_i) \\right) \\boldsymbol{x}_i$$\n",
    "\n",
    "Because we do not know how to solve \n",
    "$$\\frac1n \\sum_{i=1}^n \\left(y_i - \\sigma(\\boldsymbol\\beta^\\top \\boldsymbol{x}_i) \\right) \\boldsymbol{x}_i = 0,$$\n",
    "we do not have a closed-form solution.\n",
    "\n",
    "However we can use gradient descent!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8193c6aa",
   "metadata": {},
   "source": [
    "Here are an implementation of the evaluation of the logistic loss and of the gradient of the logistic regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de8ba9e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sigmoid(u):\n",
    "    return (1./(1+np.exp(-u)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "932cfbe7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def logistic_loss(X, y, beta):\n",
    "    \"\"\"\n",
    "    Logistic loss of a logistic regression.\n",
    "    Parameters:\n",
    "    X: (n_samples, n_features) numpy array\n",
    "        The data matrix.\n",
    "        \n",
    "    y: (n_samples, ) numpy array \n",
    "        The output/target vector (0/1).\n",
    "        \n",
    "    beta: (n_features, ) numpy array\n",
    "        The weight vector of the logistic regression.    \n",
    "    \"\"\"\n",
    "    where_y_pos = np.where(y)[0]\n",
    "    loss_pos = - np.sum(np.log(sigmoid(X[where_y_pos, :].dot(beta))))\n",
    "    where_y_neg = np.where(y==0)[0]\n",
    "    loss_neg = - np.sum(np.log(1 - sigmoid(X[where_y_neg, :].dot(beta))))\n",
    "    return (loss_pos + loss_neg)/np.size(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50e81389",
   "metadata": {},
   "outputs": [],
   "source": [
    "def evaluate_gradient_logistic(X, y, b_value):\n",
    "    \"\"\" \n",
    "    Gradient of the empirical risk of the logistic regression.\n",
    "    \"\"\"\n",
    "    num_samples = X.shape[0]\n",
    "    diff = sigmoid(X.dot(b_value)) - y \n",
    "    return np.sum(np.multiply(diff, X.T), axis=1)/num_samples"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02347fba",
   "metadata": {},
   "source": [
    "### 2.2 Gradient descent"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ab82429",
   "metadata": {},
   "source": [
    "##### __Question 21__\n",
    "Complete the following code to implement the batch gradient descent procedure for logistic regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ebcf67f4",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-03-27T17:38:54.314594Z",
     "start_time": "2019-03-27T17:38:54.288955Z"
    }
   },
   "outputs": [],
   "source": [
    "def batch_gradient_descent_logistic(X, y, learning_rate=1e-1, max_iters=30, tol=1e-2):\n",
    "    \"\"\"\n",
    "    Batch gradient descent procedure for a logistic regression.\n",
    "    Parameters:\n",
    "    X: (n_samples, n_features) numpy array\n",
    "        The data matrix.\n",
    "        \n",
    "    y: (n_samples, ) numpy array \n",
    "        The output/target vector.\n",
    "        \n",
    "    learning_rate: float\n",
    "        The learning rate of the gradient descent.\n",
    "        \n",
    "    max_iters: int\n",
    "        The maximum number of iterations of the gradient descent.\n",
    "        \n",
    "    tol: float\n",
    "        The tolerance of the gradient descent.\n",
    "        \n",
    "    Returns:\n",
    "    losses: (n_iters, ) numpy array \n",
    "        The loss at each iteration.\n",
    "    \n",
    "    beta_vector: (n_features, ) numpy array\n",
    "        The fitted regression vector (i.e. the solution).\n",
    "    \"\"\"\n",
    "    # TODO\n",
    "    return(losses, beta_vector)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a2f21b2",
   "metadata": {},
   "source": [
    "Now run the code below to run the batch gradient descent on the data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9362c0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "tolerance = 1e-4\n",
    "iterations = 200\n",
    "lr = 5.\n",
    "losses, beta_final = batch_gradient_descent_logistic(X_bt, y, learning_rate=lr, max_iters=iterations, tol=tolerance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b13d575",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(np.arange(len(losses)), losses)\n",
    "\n",
    "plt.xlabel(\"number of iterations\")\n",
    "plt.ylabel(\"value of the loss\")\n",
    "plt.title(\"Batch gradient descent\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4bb7a8c1",
   "metadata": {},
   "source": [
    "#### Quality of the model\n",
    "\n",
    "##### __Question 22__\n",
    "Complete the following code to compute the ROC curve of your model on the training data. Note again that here we're only trying to evaluate how well the model fits the data, and not whether it generalizes well. Comment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0fff21af",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7872e0a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred_gd = # TODO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6a4de38",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(4, 4))\n",
    "\n",
    "fpr, tpr, thresholds = metrics.roc_curve(y, y_pred_gd)\n",
    "plt.plot(fpr, tpr, label=\"Batch gradient descent (AUC=%.3f)\" % metrics.auc(fpr, tpr))\n",
    "\n",
    "plt.xlabel(\"TPR\")\n",
    "plt.ylabel(\"FPR\")\n",
    "plt.title(\"ROC curve (on the train set)\")\n",
    "plt.legend(loc=(1.01, .9))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34bf4b71-ad49-46a4-8f1e-7ceda71fa894",
   "metadata": {},
   "source": [
    "__Comment:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c97f8ec",
   "metadata": {},
   "source": [
    "### 2.3 Comparing batch, minibatch and stochastic gradient descents"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb205db3",
   "metadata": {},
   "source": [
    "#### Minibatch gradient descent"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e49bc6f",
   "metadata": {
    "tags": []
   },
   "source": [
    "##### __Question 23__\n",
    "Complete the code below that implements minibatch gradient descent for the logistic regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "629f58da",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-03-27T17:39:08.231736Z",
     "start_time": "2019-03-27T17:39:08.198510Z"
    }
   },
   "outputs": [],
   "source": [
    "def minibatch_gradient_descent_logistic(X, y, batch_size=10, learning_rate=1e-1, max_iters=100, tol=1e-5):    \n",
    "    \"\"\"\n",
    "    Minibatch gradient descent procedure for a logistic regression.\n",
    "    Parameters:\n",
    "    X: (n_samples, n_features) numpy array\n",
    "        The data matrix.\n",
    "        \n",
    "    y: (n_samples, ) numpy array \n",
    "        The output/target vector.\n",
    "        \n",
    "    batch_size: int\n",
    "        The number of samples in a batch.\n",
    "        \n",
    "    learning_rate: float\n",
    "        The learning rate of the gradient descent.\n",
    "        \n",
    "    max_iters: int\n",
    "        The maximum number of iterations of the gradient descent.\n",
    "        \n",
    "    tol: float\n",
    "        The tolerance of the gradient descent\n",
    "        \n",
    "    Returns:\n",
    "    losses: (n_iters, ) numpy array \n",
    "        The loss at each iteration.\n",
    "    \n",
    "    beta_vector: (n_features, ) numpy array\n",
    "        The fitted regression vector (i.e. the solution).\n",
    "    \"\"\"\n",
    "    # TODO\n",
    "    return(losses, beta_vector)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9011daba",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### Stochastic gradient descent"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90044ba8",
   "metadata": {},
   "source": [
    "##### __Question 24__\n",
    "Complete the code below that implements stochastic gradient descent for the logistic regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3883a34b",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-03-27T17:39:08.231736Z",
     "start_time": "2019-03-27T17:39:08.198510Z"
    }
   },
   "outputs": [],
   "source": [
    "def stochastic_gradient_descent_logistic(X, y, learning_rate=1e-1, max_iters=20):    \n",
    "    \"\"\"\n",
    "    Stochastic gradient descent procedure for a logistic regression.\n",
    "    Parameters:\n",
    "    X: (n_samples, n_features) numpy array\n",
    "        The data matrix.\n",
    "        \n",
    "    y: (n_samples, ) numpy array \n",
    "        The output/target vector.\n",
    "        \n",
    "    learning_rate: float\n",
    "        The learning rate of the gradient descent.\n",
    "        \n",
    "    max_iters: int\n",
    "        The maximum number of iters of the gradient descent.\n",
    "        \n",
    "    Returns:\n",
    "    losses: (n_iters, ) numpy array \n",
    "        The loss at each iteration.\n",
    "    \n",
    "    beta_vector: (n_features, ) numpy array\n",
    "        The fitted regression vector (i.e. the solution).\n",
    "        \n",
    "    \"\"\"\n",
    "    # TODO\n",
    "    return(losses, beta_vector)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "184b243a",
   "metadata": {},
   "source": [
    "#### Comparison per number of iterations"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0165aac",
   "metadata": {},
   "source": [
    "##### __Question 25__\n",
    "Complete the code below to compare the evolution of the loss, for the batch, minibatch, and stochastic gradient descents, as a function of the number of iterations. Comment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83370f4b",
   "metadata": {},
   "outputs": [],
   "source": [
    "tolerance = 1e-4\n",
    "iterations = 50 # TODO change this value if you wish\n",
    "\n",
    "# Run batch gradient descent\n",
    "lr = 1. # TODO try changing this value\n",
    "bgd_losses, bgd_beta = batch_gradient_descent_logistic(X_bt, y, learning_rate=lr, max_iters=iterations, tol=tolerance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55da210d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run minibatch gradient descent\n",
    "bs = 50 # TODO try changing this value\n",
    "lr = 0.1 # TODO try changing this value\n",
    "mgd_losses, mgd_beta = minibatch_gradient_descent_logistic(X_bt, y, batch_size=bs, learning_rate=lr, max_iters=iterations, tol=tolerance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f497486",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run stochastic gradient descent\n",
    "lr = 1e-2 # TODO try changing this value\n",
    "sgd_losses, sgd_beta = stochastic_gradient_descent_logistic(X_bt, y, learning_rate=lr, max_iters=iterations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04b6f6ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot batch gradient loss\n",
    "plt.plot(np.arange(len(bgd_losses)), bgd_losses, label=\"Batch gradient descent\")\n",
    "\n",
    "# Plot loss\n",
    "plt.plot(np.arange(len(mgd_losses)), mgd_losses, label=\"Minibatch gradient descent\")\n",
    "\n",
    "# Plot loss\n",
    "plt.plot(np.arange(len(sgd_losses)), sgd_losses, label=\"Stochastic gradient descent\")\n",
    "\n",
    "plt.xlabel(\"number of iterations\")\n",
    "plt.ylabel(\"value of the loss\")\n",
    "plt.title(\"Comparing gradient descents\")\n",
    "plt.legend(loc=(1.01, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1d5cee1-b99b-4897-b6fc-c316e681f2fe",
   "metadata": {},
   "source": [
    "__Comment:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f829929",
   "metadata": {},
   "source": [
    "#### Comparison per number of epochs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f600f879",
   "metadata": {},
   "source": [
    "##### __Question 26__\n",
    "Complete the code below to compare the evolution of the loss, for the batch, minibatch, and stochastic gradient descents, as a function of the number of epochs. Comment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a544383",
   "metadata": {},
   "outputs": [],
   "source": [
    "tolerance = 1e-4\n",
    "epochs = 5 # TODO change this is you wish\n",
    "\n",
    "# Run batch gradient descent\n",
    "lr = # TODO use the same value as for the comparison per number of iterations\n",
    "bgd_losses, bgd_beta = batch_gradient_descent_logistic(X_bt, y, learning_rate=lr, max_iters=epochs, tol=tolerance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3b060fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run minibatch gradient descent\n",
    "bs = # TODO use the same value as for the comparison per number of iterations\n",
    "lr = # TODO use the same value as for the comparison per number of iterations\n",
    "n_iterations_per_epoch = int(np.ceil(n_samples/bs))\n",
    "iterations = epochs * n_iterations_per_epoch\n",
    "mgd_losses, mgd_beta = minibatch_gradient_descent_logistic(X_bt, y, batch_size=bs, learning_rate=lr, max_iters=iterations, tol=tolerance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8972570f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run stochastic gradient descent\n",
    "lr = # TODO use the same value as for the comparison per number of iterations\n",
    "iterations = epochs * n_samples\n",
    "sgd_losses, sgd_beta = stochastic_gradient_descent_logistic(X_bt, y, learning_rate=lr, max_iters=iterations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a592a01d",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.semilogy(np.arange(len(bgd_losses)), bgd_losses, label=\"Batch gradient descent\")\n",
    "plt.semilogy(np.arange(len(mgd_losses))/n_iterations_per_epoch, mgd_losses, label=\"Minibatch gradient descent\")\n",
    "plt.semilogy(np.arange(len(sgd_losses))/n_samples, sgd_losses, label=\"Stochastic gradient descent\")\n",
    "\n",
    "plt.xlabel(\"number of epochs\")\n",
    "plt.ylabel(\"empirical risk\")\n",
    "plt.title(\"Comparing gradient descents\")\n",
    "plt.legend(loc=(1.01, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ce48dd4-601c-4f89-8fd5-0dd6d2bee5d9",
   "metadata": {
    "tags": []
   },
   "source": [
    "__Comment:__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45e58da1",
   "metadata": {},
   "source": [
    "### 2.4 SGDClassifier in scikit-learn"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f261e6c0-ca79-4694-9897-e6fa836c9b35",
   "metadata": {},
   "source": [
    "##### __Question 27__\n",
    "Complete the code below to use the `SGDClassifier` class of scikit-learn ([doc](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.SGDClassifier.html)) instead of the above hand-written procedures. Comment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac8f0533",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize an instance of SGDClassifier\n",
    "classifier = linear_model.SGDClassifier(# TODO)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5d38db9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Learn\n",
    "classifier.fit(X_bt, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b53b0e71",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"SGDClassifier stopped after {classifier.n_iter_} epochs.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb3362b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compare the loss of the model learned by SGDClassifier with the minimum loss\n",
    "sgd_loss = logistic_loss(X_bt, y, classifier.coef_.reshape((X_bt.shape[1], )))\n",
    "print(f\"SGDClassifier loss: {sgd_loss:.2e}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffd15bef",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.semilogy(np.arange(len(bgd_losses)), bgd_losses, label=\"Batch gradient descent\")\n",
    "plt.semilogy(np.arange(len(mgd_losses))/n_iterations_per_epoch, mgd_losses, label=\"Minibatch gradient descent\")\n",
    "plt.semilogy(np.arange(len(sgd_losses))/n_samples, sgd_losses, label=\"Stochastic gradient descent\")\n",
    "\n",
    "plt.semilogy([0, epochs], [sgd_loss, sgd_loss], 'k--', label=('SGDClassifier after %d epochs' % classifier.n_iter_))\n",
    "\n",
    "plt.xlabel(\"number of epochs\")\n",
    "plt.ylabel(\"empirical risk\")\n",
    "plt.title(\"Comparing gradient descents\")\n",
    "plt.legend(loc=(1.01, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "081c837a-c0b5-4810-b66f-7814a112754d",
   "metadata": {
    "tags": []
   },
   "source": [
    "__Comment:__"
   ]
  }
 ],
 "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}