{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lab 1 : First order methods on regression models\n",
    "\n",
    "#### Authors: A. Gramfort, R. Gower, P. Ablin, M. Massias\n",
    "\n",
    "## Aim\n",
    "\n",
    "The aim of this material is to code \n",
    "- proximal gradient descent (PGD)\n",
    "- accelerated gradient descent (APGD) \n",
    "\n",
    "for \n",
    "- linear regression\n",
    "- logistic regression \n",
    "\n",
    "models.\n",
    "\n",
    "The proximal operators we will use are the \n",
    "- L2 penalization\n",
    "- L1 penalization\n",
    "\n",
    "## VERY IMPORTANT\n",
    "\n",
    "- This work **must be done by pairs of students**.\n",
    "- **Each** student must send their work using the **moodle platform**. Please check the deadline on moodle!\n",
    "- This means that **each student in the pair sends the same file**\n",
    "- The **name of the file must be** constructed as in the next cell\n",
    "\n",
    "# Gentle reminder: no evaluation if you don't respect this EXACTLY\n",
    "\n",
    "### How to construct the name of your file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Change here using YOUR first and last names\n",
    "fn1 = \"pierre\"\n",
    "ln1 = \"ablin\"\n",
    "fn2 = \"alexandre\"\n",
    "ln2 = \"gramfort\"\n",
    "\n",
    "filename = \"_\".join(map(lambda s: s.strip().lower(), \n",
    "                        [\"lab1\", ln1, fn1, \"and\", ln2, fn2])) + \".ipynb\"\n",
    "print(filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## to embed figures in the notebook\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 0 : Introduction\n",
    "\n",
    "We'll start by generating sparse vectors and simulating data\n",
    "\n",
    "### Getting sparse coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "np.set_printoptions(precision=2)  # to have simpler print outputs with numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_features = 50\n",
    "n_samples = 1000\n",
    "idx = np.arange(n_features)\n",
    "coefs = ((-1) ** idx) * np.exp(-idx / 10.)\n",
    "coefs[20:] = 0.\n",
    "plt.stem(coefs)\n",
    "plt.title(\"Parameters / Coefficients\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Functions for the simulation of the models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy.random import multivariate_normal\n",
    "from scipy.linalg.special_matrices import toeplitz\n",
    "\n",
    "\n",
    "def simu_linreg(coefs, n_samples=1000, corr=0.5):\n",
    "    \"\"\"Simulation of a linear regression model\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    coefs : `numpy.array`, shape (n_features,)\n",
    "        Coefficients of the model\n",
    "    \n",
    "    n_samples : `int`, default=1000\n",
    "        Number of samples to simulate\n",
    "    \n",
    "    corr : `float`, default=0.5\n",
    "         Correlation between features i and j is corr^|i - j|\n",
    "\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    A : `numpy.ndarray`, shape (n_samples, n_features)\n",
    "        Simulated features matrix. It samples of a centered Gaussian \n",
    "        vector with covariance given by the Toeplitz matrix\n",
    "    \n",
    "    b : `numpy.array`, shape (n_samples,)\n",
    "        Simulated labels\n",
    "    \"\"\"\n",
    "    # Construction of a covariance matrix\n",
    "    cov = toeplitz(corr ** np.arange(0, n_features))\n",
    "    # Simulation of features\n",
    "    A = multivariate_normal(np.zeros(n_features), cov, size=n_samples)\n",
    "    # Simulation of the labels\n",
    "    b = A.dot(coefs) + np.random.randn(n_samples)\n",
    "    return A, b\n",
    "\n",
    "\n",
    "def sigmoid(t):\n",
    "    \"\"\"Sigmoid function\"\"\"\n",
    "    return 1. / (1. + np.exp(-t))\n",
    "\n",
    "\n",
    "def simu_logreg(coefs, n_samples=1000, corr=0.5):\n",
    "    \"\"\"Simulation of a logistic regression model\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    coefs : `numpy.array`, shape (n_features,)\n",
    "        Coefficients of the model\n",
    "    \n",
    "    n_samples : `int`, default=1000\n",
    "        Number of samples to simulate\n",
    "    \n",
    "    corr : `float`, default=0.5\n",
    "        Correlation between features i and j is corr^|i - j|\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    A : `numpy.ndarray`, shape (n_samples, n_features)\n",
    "        Simulated features matrix. It samples of a centered Gaussian \n",
    "        vector with covariance given by the Toeplitz matrix\n",
    "    \n",
    "    b : `numpy.array`, shape (n_samples,)\n",
    "        Simulated labels\n",
    "    \"\"\"\n",
    "    cov = toeplitz(corr ** np.arange(0, n_features))\n",
    "    A = multivariate_normal(np.zeros(n_features), cov, size=n_samples)\n",
    "    p = sigmoid(A.dot(coefs))\n",
    "    b = np.random.binomial(1, p, size=n_samples)\n",
    "    b = 2 * b - 1\n",
    "    return A, b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A, b = simu_linreg(coefs)\n",
    "# A, b = simu_logreg(coefs)  # uncomment when you want to test logistic regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1 : Proximal operators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We remind that the proximal operator of a function $R$ is given by:\n",
    "\n",
    "$$\n",
    "\\text{prox}_{(\\lambda \\,R)}(x) = \\arg\\min_z \\left\\{ \\frac{1}{2} \\Vert x - z\\Vert_2^2 + \\lambda \\,R(z) \\right\\}.\n",
    "$$\n",
    "\n",
    "\n",
    "We have in mind to use the following cases\n",
    "\n",
    "- L2 penalization, where $R(z) = \\frac{1}{2} \\|z\\|_2^2$\n",
    "- L1 penalization, where $R(z) =  \\|z|\\|_1$\n",
    "\n",
    "where $\\lambda > 0$ is a regularization parameter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>Questions</b>:\n",
    "     <ul>\n",
    "      <li>Code a function that computes $R(x)$ in both cases and $\\text{prox}_{\\lambda\\, R}(x)$ for L2 and  L1 penalization (use the slides of the first course to get the formulas), using the prototypes given below</li>\n",
    "      <li>Visualize the functions applied element wise by the proximity operators of the L2 and L1 \n",
    "    </li>\n",
    "    </ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def prox_L1(x, lmbd):\n",
    "    \"\"\"Proximal operator for the L1 at x\"\"\"    \n",
    "    return x # TODO\n",
    "    \n",
    "def L1_norm(x, lmbd):\n",
    "    \"\"\"Value of the L1 penalization at x\"\"\"\n",
    "    return 0. # TODO\n",
    "\n",
    "def prox_L2(x, lmbd):\n",
    "    \"\"\"Proximal operator for the L2 at x\"\"\"    \n",
    "    return x # TODO\n",
    "    \n",
    "def L2_norm(x, lmbd):\n",
    "    \"\"\"Value of the L2 penalization at x\"\"\"\n",
    "    return 0. # TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualization\n",
    "\n",
    "We are now going to visualize the effect of the proximity operators on coefficients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.random.randn(50)\n",
    "l_l1 = 1.\n",
    "l_l2 = 0.5\n",
    "\n",
    "plt.figure(figsize=(15.0, 4.0))\n",
    "plt.subplot(1, 3, 1)\n",
    "plt.stem(x)\n",
    "plt.title(\"Original parameter\", fontsize=16)\n",
    "plt.ylim([-2, 2])\n",
    "plt.subplot(1, 3, 2)\n",
    "plt.stem(prox_L1(x, lmbd=l_l1))\n",
    "plt.title(\"Proximal L1\", fontsize=16)\n",
    "plt.ylim([-2, 2])\n",
    "plt.subplot(1, 3, 3)\n",
    "plt.stem(prox_L2(x, lmbd=l_l2))\n",
    "plt.title(\"Proximal L2\", fontsize=16)\n",
    "plt.ylim([-2, 2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>Question</b>:\n",
    "     <ul>\n",
    "      <li>Comment what you observe (1 or 2 sentences).</li>\n",
    "    </ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2: Gradients\n",
    "\n",
    "The problems we want to minimize take the form:\n",
    "$$\n",
    "\\arg\\min_x f(x) + \\lambda \\,R(x)\n",
    "$$\n",
    "where $f$ is $L$-smooth and $R$ has a simple proximal operator.\n",
    "\n",
    "Consider the following cases:\n",
    "\n",
    "**Linear regression**, where \n",
    "$$\n",
    "f(x) = \\frac{1}{2n} \\sum_{i=1}^n (b_i - a_i^\\top x)^2 = \\frac{1}{2 n} \\| b - A x \\|_2^2,\n",
    "$$\n",
    "where $n$ is the sample size, $b = [b_1 \\cdots b_n]$ is the vector of labels and $A = [a_1,\\ldots, a_n]^{\\top} \\in \\mathbb{R}^{n \\times d}$ is the matrix of samples.\n",
    "\n",
    "**Logistic regression**, where\n",
    "$$\n",
    "f(x) = \\frac{1}{n} \\sum_{i=1}^n \\log(1 + \\exp(-b_i a_i^\\top x)),\n",
    "$$\n",
    "where $n$ is the sample size, and where labels $b_i \\in \\{ -1, 1 \\}$ for all $i$.\n",
    "\n",
    "We need to be able to compute $f$ and its gradient"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>Questions</b>:\n",
    "     <ul>\n",
    "      <li>Compute on paper the gradient $\\nabla f$ of $f$ for both cases (linear and logistic regression)</li>\n",
    "      <li>Code a function that computes $f$ and its gradient $\\nabla f$ in both cases, using the prototypes below.</li>\n",
    "      <li>Check that these functions are correct by numerically checking the gradient, using the function ``<a href=\"https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.check_grad.html\">check_grad</a>`` from ``scipy.optimize``. Remark: use the functions `simu_linreg` and `simu_logreg` to simulate data according to the right model</li>\n",
    "    </ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def loss_linreg(x):\n",
    "    \"\"\"Least-squares loss\"\"\"\n",
    "    # TODO\n",
    "    pass\n",
    "\n",
    "def grad_linreg(x):\n",
    "    \"\"\"Leas-squares gradient\"\"\"\n",
    "    # TODO\n",
    "    pass\n",
    "\n",
    "def loss_logreg(x):\n",
    "    \"\"\"Logistic loss\"\"\"\n",
    "    # TODO\n",
    "    pass\n",
    "\n",
    "def grad_logreg(x):\n",
    "    \"\"\"Logistic gradient\"\"\"\n",
    "    # TODO\n",
    "    pass\n",
    "\n",
    "# TO BE COMPLETED"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 3: Solvers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We know have a function to compute $f$, $\\nabla f$ and $R$ and $\\text{prox}_R$. \n",
    "\n",
    "We want now to code the PGD and APGD solvers to minimize\n",
    "\n",
    "$$\n",
    "\\arg\\min_x f(x) + \\lambda \\, R(x)\n",
    "$$\n",
    "\n",
    "where $\\lambda >0$ is the regularization parameter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>Questions</b>:\n",
    "     <ul>\n",
    "      <li>Implement functions that compute the Lipschitz constants for linear and \n",
    "  logistic regression losses. Note that the operator norm of a matrix can \n",
    "  be computed using the function <a href=\"https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.norm.html\">numpy.linalg.norm</a> (read the documentation\n",
    "  of the function)</li>\n",
    "      <li>Finish the functions `PGD` and `APGD` below that implements the \n",
    "  PGD(Proximal Gradient Descent) and APGD (Accelerated Proximal \n",
    "  Gradient Descent) algorithms</li>\n",
    "    </ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TO BE COMPLETED\n",
    "\n",
    "# Calculate the smoothness constant for L2\n",
    "def smoothness_const_linreg(A):\n",
    "    \"\"\"Lipschitz constant for linear squares loss\"\"\"    \n",
    "    # TODO\n",
    "    pass\n",
    "# Calculate the smoothness constant for Logistic loss   \n",
    "def smoothness_const_loreg(A):\n",
    "    \"\"\"Lipschitz constant for logistic loss\"\"\"    \n",
    "    # TODO\n",
    "    pass\n",
    "    \n",
    "def PGD(x0, f, grad_f, R, prox_R, step, lmbd=0., n_iter=50,\n",
    "         x_true=coefs, verbose=True):\n",
    "    \"\"\"Proximal gradient descent algorithm\n",
    "    \"\"\"\n",
    "    x = x0.copy()\n",
    "    x_new = x0.copy()\n",
    "    n_samples, n_features = A.shape\n",
    "\n",
    "    # estimation error history\n",
    "    errors = []\n",
    "    # objective history\n",
    "    objectives = []\n",
    "    # Current estimation error\n",
    "    err = norm(x - x_true) / norm(x_true)\n",
    "    errors.append(err)\n",
    "    # Current objective\n",
    "    obj = f(x) + R(x, lmbda)\n",
    "    objectives.append(obj)\n",
    "    if verbose:\n",
    "        print(\"Lauching PGDsolver...\")\n",
    "        print(' | '.join([name.center(8) for name in [\"it\", \"obj\", \"err\"]]))\n",
    "    for k in range(n_iter + 1):\n",
    "\n",
    "        #### TODO ####\n",
    "        \n",
    "        obj = f(x) + R(x, lmbda)\n",
    "        err = norm(x - x_true) / norm(x_true)\n",
    "        errors.append(err)\n",
    "        objectives.append(obj)\n",
    "        if k % 10 == 0 and verbose:\n",
    "            print(' | '.join([(\"%d\" % k).rjust(8), \n",
    "                              (\"%.2e\" % obj).rjust(8), \n",
    "                              (\"%.2e\" % err).rjust(8)]))\n",
    "    return x, objectives, errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TO BE COMPLETED\n",
    "\n",
    "def APGD(x0, f, grad_f, R, prox_R, step, lmbd=0., n_iter=50,\n",
    "         x_true=coefs, verbose=True):\n",
    "    \"\"\"Accelerated Proximal gradient descent algorithm\n",
    "    \"\"\"\n",
    "    x = x0.copy()\n",
    "    x_new = x0.copy()\n",
    "    # An extra variable is required for APGD\n",
    "    z = x0.copy()\n",
    "    n_samples, n_features = A.shape\n",
    "    # estimation error history\n",
    "    errors = []\n",
    "    # objective history\n",
    "    objectives = []\n",
    "    # Current estimation error\n",
    "    err = norm(x - x_true) / norm(x_true)\n",
    "    errors.append(err)\n",
    "    # Current objective\n",
    "    obj = f(x) + R(x, lmbda)\n",
    "    objectives.append(obj)\n",
    "    t = 1.\n",
    "    t_new = 1.    \n",
    "    if verbose:\n",
    "        print(\"Lauching APGD solver...\")\n",
    "        print(' | '.join([name.center(8) for name in [\"it\", \"obj\", \"err\"]]))\n",
    "    for k in range(n_iter + 1):\n",
    "\n",
    "        #### TODO ####\n",
    "\n",
    "        obj = f(x) + R(x, lmbda)\n",
    "        err = norm(x - x_true) / norm(x_true)\n",
    "        errors.append(err)\n",
    "        objectives.append(obj)\n",
    "        if k % 10 == 0 and verbose:\n",
    "            print(' | '.join([(\"%d\" % k).rjust(8), \n",
    "                              (\"%.2e\" % obj).rjust(8), \n",
    "                              (\"%.2e\" % err).rjust(8)]))\n",
    "    return x, np.array(objectives), np.array(errors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Algorithms comparison and numerical experiments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Some definitions before launching the algorithms\n",
    "x0 = np.zeros(n_features)\n",
    "n_iter = 200\n",
    "lmbd = 1e-2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>Questions</b>:\n",
    "     <ul>\n",
    "      <li>Compute a precise minimum and a precise minimizer of the linear regression with L2 \n",
    "  penalization problem using the parameters give above. This can be done by using APGD with \n",
    "  1000 iterations.</li>\n",
    "    <li>Compare the convergences of PGD and APGD, in terms of distance to the minimum and \n",
    "  distance to the minimizer. Do your plots using a logarithmic scale of the y-axis.</li>\n",
    "    </ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <b>EXERCISE</b>:\n",
    "     <ul>\n",
    "      <li>Compare the solution you obtain with PGD and APGD with the true parameter `coefs` of\n",
    "  the model. This can be done with 3 `plt.stem` plots. Do this for one problem (linear or logistic), but for both regularizers. Comment.</li>\n",
    "    <li> Pick one problem (logistic or linear) and regularizer (L1 or L2), and study the influence of the correlation \n",
    "  of the features on the performance of the optimization algorithms. Explain.</li>\n",
    "    <li>Pick a different problem (logistic or linear) and a different regularization (L1 or L2), and study the influence of the level of  \n",
    "  penalization ($\\lambda$) on the performance of the optimization algorithms. Explain. Analyse and explain the behavior of the solution vector when $\\lambda$ increases.</li>\n",
    "    </ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = m"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "interpreter": {
   "hash": "ddeb6e512e1ef22841d21b181b577eb7b6fd883489eaad141b98586281502ceb"
  },
  "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.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}