Part5-Proximal-Lasso.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 5: Convexity and Constrained Optimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import importlib\n",
    "\n",
    "import numpy as np # linear algebra\n",
    "import pandas as pd #dataframe library\n",
    "import sklearn # machine learning library\n",
    "assert sklearn.__version__ >= '0.2'\n",
    "from sklearn import datasets\n",
    "import random\n",
    "np.random.seed(seed=42) # to have same results each time\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",
    "import useful_functions\n",
    "importlib.reload(useful_functions)\n",
    "from useful_functions import generate_dataset, unison_shuffled_copies, get_train_test_datasets, gradient_descent\n",
    "from useful_functions import plot_convergence_rate, plot_evolution_on_test_error, plot_evolution_on_test_error_with_bias\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 1: Add an l2 regularization term to your objective function from Part 1 or Part 3. Compare the solution of the unregularized problem to those obtained while solving the problem with a) a small value for the regularization parameter and b) a large value for the regularization parameter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The objective function of part 1 was trained using a small ridge parameter. We direct the reader towards part 1 for more information. \n",
    "Here we just note that adding a small ridge penalty (2nd figure) leads to better results than no ridge (1st figure). In any case, the linear model still underfits the task at hand."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#% of rows to include in training dataset\n",
    "TRAINING_RATIO = 0.75\n",
    "\n",
    "dataset = sklearn.datasets.fetch_california_housing(data_home=None, download_if_missing=True, return_X_y=False, as_frame=False)\n",
    "A = dataset['data']\n",
    "Y = dataset['target']\n",
    "A_shuffled, Y_shuffled = unison_shuffled_copies(A, Y)\n",
    "A_train, Y_train, A_test, Y_test = get_train_test_datasets(A_shuffled, Y_shuffled,TRAINING_RATIO)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAERCAYAAAAudzN9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAl1UlEQVR4nO3de5xdVX338c93ZpJMyAXIRYQohFBSMBQsTa1WsVjRikIFQlseQAEtAdFXHwWr9hFNBEVBabUi1wopClS5thTvGoXipYxcxEgIJhAEQUIMITdy/T1/rH1m9jnZM8yeOZOzM/N9v177NXP276x91hrC/GbttfZaigjMzMyqpq3VFTAzMyviBGVmZpXkBGVmZpXkBGVmZpXkBGVmZpXU0eoKDBdTpkyJ6dOnt7oaZmY7lZ///OfPRsTUopgTVJNMnz6drq6uVlfDzGynIml5bzHf4jMzs0pygjIzs0pygjIzs0pygjIzs0pygjIzs0pygjIzs0pygjIzs0pygmqx7z/0Oy774dJWV8PMrHKcoFps4cPPcNVdy1pdDTOzynGCajGhVlfBzKySnKAqwLsam5ltzwmqxSRwejIz254TVIsJcAfKzGx7TlAtJnkMysysSEsTlKT3SeqStFHSgtz5V0v6rqTfS1oh6UZJe+biknShpJXZcZFyv+klTZe0UNJ6SYslHdHwuSdKWi5pnaTbJE3KxcZIulrS85KelnT2EP8YPAZlZlag1T2o3wKfBK5uOL87cCUwHdgHWANck4vPBY4BDgEOBo4CzsjFbwDuAyYDHwVukjQVQNIs4ArgHcAewHrg0lzZ+cD+2ee+AfiQpLcMppEvxunJzGx7LU1QEXFLRNwGrGw4/82IuDEino+I9cAlwGtzbzkFuDginoiIJ4GLgVMBJM0EDgXmRcSGiLgZeBCYk5U9Cbg9Iu6MiLXAx4DjJE3I4u8Ezo+IVRHxEHBV7dpDQcIZysysQKt7UP31emBR7vUs4IHc6weyc7XYsohY00e8u2xELAU2ATMl7Q7s1ce1m07I+cnMrEDlt3yXdDDwceDtudPjgdW516uB8dk4VGOsFp/WS9lafEIWg+2vPYECkuaSbjey995796M1RdcYUDEzs2Gv0j0oSX8AfBP4vxFxVy60FpiYez0RWBtptkFjrBZf00vZfHxt7nVR2ToRcWVEzI6I2VOnTu1fo4qvM+CyZmbDVWUTlKR9gO+RxoO+0hBeRJogUXMIPbcAFwEzcmNKRfHuspJmAGOAJRGxCniqj2s3nYegzMyKtXqaeYekTqAdaJfUmZ2bBvwA+FJEXF5Q9FrgbEnTJO0FnAMsAIiIJcD9wLzseseSZvrdnJW9Djha0mGSxgHnAbfkxqyuBc6VtLukA4DTa9ceCpIf1DUzK9LqMahzgXm51ycDnyB1KmaQkkx3PCJqY0RXZPEHs9f/lp2rOYGUVFYBjwPHR8SK7BqLJJ1JSlSTSb2003Jl5wGXAcuBDcCFEfGtwTa0N35Q18ysmDz+0RyzZ8+Orq6u0uUu+MZDXPuTx1h8/pFDUCszs2qT9POImF0Uq+wY1EjhtfjMzIo5QbWaVzM3MyvkBNVi3rDQzKyYE1QVuAtlZrYdJ6gWSxsWOkOZmTVygmoxT5IwMyvmBNVifgzKzKyYE1QFuANlZrY9J6gWE/JisWZmBZygWkx+DsrMrJATVIt5koSZWTEnqFbzLAkzs0JOUGZmVklOUC1W6z95ooSZWT0nqBar3eFzfjIzq+cE1WJeLNbMrJgTVEW4A2VmVs8JqsV6bvE5RZmZ5TlBtVj3JImW1sLMrHqcoFrMj0GZmRVzgqoI3+EzM6vnBNViyrpQ3rTQzKyeE1RFuAdlZlbPCarFPAZlZlbMCcrMzCrJCarFaitJ+BafmVk9J6gW635Q15MkzMzqOEG1WM9q5i2thplZ5ThBtZgnSZiZFeso82alh3aOAPYHJsN2S3FHRJzfpLqNKO5AmZnV63eCkrQ/cBtwANsnppoAnKBK6Jkk4RRlZpZX5hbfF4H9gA8Ds4F9C44ZZT5c0vskdUnaKGlB7vxoSTdJekxSSDq8oZwkXShpZXZcJPXcLJM0XdJCSeslLZZ0REP5EyUtl7RO0m2SJuViYyRdLel5SU9LOrtMm8rqmSRhZmZ5ZW7xvQ74fER8romf/1vgk8BfAWMbYv8DfB64saDcXOAY4BDS7/bvAsuAy7P4DcBPgLdmx02S9o+IFZJmAVcAbwPuBa4ELgVOyMrOJ93C3Ad4KbBQ0q8i4luDa6qZmZVRpge1CXi0mR8eEbdExG3AyobzmyLi8xHxP8DWgqKnABdHxBMR8SRwMXAqgKSZwKHAvIjYEBE3Aw8Cc7KyJwG3R8SdEbEW+BhwnKQJWfydwPkRsSoiHgKuql17KPkOn5lZvTIJ6tvAa4eqIiXNAh7IvX4gO1eLLYuINX3Eu8tGxFJS8p0paXdgrz6u3XTyPT4zs0JlEtTZwGsknSNp9FBVqJ/GA6tzr1cD47NxqMZYLT6hl7L5+Pjc66KydSTNzcbQulasWFG6EZDfsNAZyswsr0yCuhuYCFwErMsmGSxrOJYOTTW3szarS81EYG2kqXCNsVp8TS9l8/G1uddFZetExJURMTsiZk+dOrV0I8DPQZmZ9abMJInHqc6NqEWkCRL/m70+JDtXi82QNCF3m+8Q4PqGsgBImgGMAZZExBpJT2Xx7xZce8h4DMrMrF6/E1REHN7sD5fUkdWhHWiX1AlsiYgtksbQcwdsdBbbmPWSrgXOlvQNUtI8hzQNnohYIul+YJ6kc4EjgYPpmSRxHfATSYeRZvGdB9ySS2bXAudK6gL2AE4HTmt227t/BtlX5yczs3qtXuroXGAD8BHg5Oz7c7PYw9nraaQJGhtIU78hTRO/nTQ775fAHdm5mhNIz2qtAj4DHB8RKwAiYhFwJilRPUMaXzorV3YesBRYDvwI+OxQTjHv3lHXXSgzszqlljoCkLQf8HZ6HspdBvxnNhuulIiYT3ruqCg2vY9yAXwoO4rijwGH91H+enpu+TXGNgLvyo4h5zEoM7NiZdfiO5/U22lvCF0k6YKI+HjTajbCuP9kZlav37f4JL0L+CjwM+BY0moL+5NWdPgJ8FFJQzZWM1x5uw0zs2JlelDvJSWnwyNiS+780myywl3A+4Brmli/4a82BuU+lJlZnTKTJA4E/qMhOQGQnfuP7D1WQvcQlPOTmVmdsmvxje8jPiF7j5XgSRJmZsXKJKh7gDMk7dEYkPQS0grjP2tWxUYad6DMzOqVGYM6H/g+8JCkLwO/ys7PIj3IOoG0UriV0LNhYYsrYmZWMWVWkrhT0nHAJaSVG/IeB06JiLuaWbmRoGcxc2coM7O8Us9BRcTtku4A/oS0g65Iqy7cGxHbhqB+w56HoMzMipVeSSJLRPdkhzWJb/GZmdVr9Vp8I573KzQzK9ZrD0rSo8A24ICI2CxpWT+uFxGxX9NqNwL0TJJwijIzy+vrFt9y0h/2td+cVdoPavjwIJSZWaFeE1Tj/k9DsR+U9XAHysysXpnFYveWNLaP+FhJezenWiOHO1BmZsXKTJJ4lLSKeW/+OnuPldCzYWGLK2JmVjFlEtSL/bHfhseoSnMPysysWNlp5n0loAOB5wZelZHNK0mYmdXr80FdSacAp+ROnSvp9IK3TgIOAm5tYt1GhO7noJyfzMzqvNhKEruRljSC1HuaCuzS8J4A1gJXk3bctRL8oK6ZWbE+E1REfAH4AoCkbcD7I+L6HVGxkUIehTIzK1RmNXMvizSEvJKEmVm9Ms9B/bGk9/YRf6+kVzalViOIb/GZmRUr0yuaB7ytj/iRwMcHV52Ryx0oM7N6ZRLUnwI/6iP+I+BVg6vOyFN7UNd9KDOzemUS1BTg933En8veYyV4ioSZWbEyCeoZYFYf8YPoO4FZH3yLz8ysXpkE9T3g7yVtl6QkvQJ4d/YeK8GTJMzMipXZ8v2TwHHAPZKuBu4n/V79Y+BdwCbg/GZXcLjr2bCwxRUxM6uYMs9BLZX0RmABcFZDeBFwWkQ80sS6jQjyIJSZWaEyPSgiogs4KHveaX/SGP/DEfHAENRtRPFisWZm9UolqJqIuJ90i88GqXuSufOTmVmdAS1fJGkXSS/PdtmtO0pe532SuiRtlLSgIfZGSYslrZe0UNI+uZgkXShpZXZcpJ4HipA0PSuzPrvGEQ3XPlHScknrJN0maVIuNkbS1ZKel/S0pLNL/4BK8GrmZmbFyix11CbpI5KeBNYAj5F20G08yvgtafLF1Q2fNQW4BfgYaSuPLuBrubfMBY4BDgEOBo4CzsjFbwDuAyaTVli/SdLU7NqzgCuAdwB7AOuBS3Nl55NuX+4DvAH4kKS3lGxXCR6EMjMrUuYW32eAD5ImRNwMrBzsh0fELQCSZgMvy4WOAxZFxI1ZfD7wrKQDImIxaY+qiyPiiSx+MXA6cLmkmcChwJsjYgNws6T3A3OAy4GTgNsj4s6s7MeAhyRNiIg1wDtJEz5WAaskXQWcCnxrsO3ti8egzMzqlUlQJwPfioi3DlVlcmYB3RMvImKdpKXZ+cWN8ez7Wbmyy7Jk01v8x7lrL5W0CZgpaRmwV8G1jymqpKS5pN4ce+9d6u5m7hq1egyouJnZsFVmDGp34D+HqiINxgOrG86tBib0El8NjM/GocqWzcfH514Xla0TEVdGxOyImD116tQ+G9Qb3+AzMytWJkE9COw5VBVpsBaY2HBuImnsqyg+EVgbaVOlsmXz8bW510Vlm05+EMrMrFCZBPUJ4ExJLx+qyuQsIk2AAEDSOGC/7Px28ez7fGyGpAl9xPPXngGMAZZk405P9XHtIeNbfGZm9cqMQf0JsBz4laRbSTP2tja8JyKi38sdSerI6tAOtEvqBLYAtwKflTQHuIO0z9QvsgkSANcCZ0v6Bmm5pXOAL2YVWCLpfmCepHNJ+1QdTJokAXAd8BNJhwH3AucBt+TGrK4FzpXURZrldzpwWn/bVFbPZhvOUGZmeWUS1Pzc9yf38p6g3Hp855I2Qsxf9xMRMT9LTpcAXwV+BpyQe98VwAzSbUeAf8vO1ZxAWpJpFfA4cHxErACIiEWSziQlqsmkBW7zCWgecBkpGW8ALoyIIZvB50kSZmbFyiSofZv94RExn/rEl499Dzigl1gAH8qOovhjwOF9fO71wPW9xDaSFr99V2/lm8mrmZuZFSuzWOzyoazISCXP4zMzKzSgpY6s+cL3+MzM6vS7B5XtAfViIiLePYj6jDy+xWdmVqjMGNSp/XhPkHbWtX7yauZmZsX6fYsvItoaD2AU8IfAVcBPSatNWAlt6plobmZmPQY1BhURWyPikYg4g7R47IXNqdbIUUtQ25yfzMzqNHOSxDfpeRjW+qkt60Btc4YyM6vTzAQ1mZ7FVq2f5B6UmVmhAW35nidpN+AI4APAzwd7vZGmZyUJZygzs7wy08y30ftIvoDfA0O6PfpwVBuDcnoyM6tXpgd1Ldv/Hg1SYloC3NCwSaD1Q/cYlHtQZmZ1ek1Qkl4PPJRbZPXUHVWpkcRjUGZmxfqaJLEQeFPthaRlkv566Ks0srgHZWZWrK8EtZG0kV/NdDxLr+m6x6CcoMzM6vQ1BrUEOEXSvaR9lQAmS9q7rwtGxOPNqtxI0P2g7rYWV8TMrGL6SlCfJO2ZdG/2OoDPZ0df2gddqxFEvsVnZlao1wQVETdJeoC08d+epJ1mbwN+sUNqNkJ4qSMzs2J9TjOPiEeARwAkzQduznajtSZpy0YBPQZlZlavzI663txwCLgHZWZWzEmnxTzN3MysmBNUi/U8qOsEZWaW5wTVYj3PQbW4ImZmFeME1WK+xWdmVswJqsU8ScLMrJgTVIv5QV0zs2L9TlCSXiXp9IZzb5f0oKQnJV3Q/OoNf16Lz8ysWJke1DygezXzbE2+G4CXAquBD0s6rbnVG/56elCtrYeZWdWUSVCHAHfnXp9A2kn3lRHxCuA7wNwm1m1E8Cw+M7NiZRLUZODp3Ou/Au6MiCez1/8F7N+sio0UHoMyMytWJkE9B+wBIGkM8Grgzlw8gLFNq9kI4TEoM7Ni/V6LD7gf+HtJ3wOOBTqBb+fi+wK/a17VRgZPMzczK1amB3U+aduN/wX+H/C9iOjKxY8CftbEuiHpQEk/kLRa0q8lHZuLvVHSYknrJS2UtE8uJkkXSlqZHReptqZQik/PyqzPrnFEw+eeKGm5pHWSbpM0qZntyvODumZmxfqdoCLix8ChwPuBU4GjazFJk0mTJC5rVsUkdQD/Cfw3MIk0AeOrkmZKmgLcAnwsi3UBX8sVnwscQ5rYcTApeZ6Ri98A3EcaV/socJOkqdnnzgKuAN5BuqW5Hri0We1qJPegzMwKlbnFR0QsIW0F33h+JfCBZlUqcwCwF/AvkQZofiDpblLi+A2wKCJuhO69qp6VdEBELAZOAS6OiCey+MXA6cDlkmaSEu2bI2IDcLOk9wNzgMuBk4DbI+LOrOzHgIckTYiINU1uY3cPymNQZmb1yjyo2y5pl4Zzu0k6R9KnJB3U5Lqpl3MHAbOAB2onI2IdsDQ7T2M8+z4fW9aQbBrj+WsvBTYBM7erjDRXUpekrhUrVvS/ZTltXs3czKxQmTGoK0jjTwBIGgX8D/BZ4J+AeyS9sol1Www8A/yjpFGS3gz8BbALMJ70cHDeamBC9n1jfDUwPhuHKlu2Md4tIq6MiNkRMXvq1Kll2tbNkyTMzIqVSVCvIz3rVHM88ArgvcCfk2bwfaRZFYuIzaRxpLeRnr86B/g68ASwFpjYUGQiUOsVNcYnAmuzW4VlyzbGm8rPQZmZFSuToPYEHs29fhtpHOiyiPgpcCXwmmZWLiJ+ERF/ERGTI+KvgBmkXtwi0gQIACSNA/bLztMYz77Px2ZImtBHPH/tGcAYCsbemsErSZiZFSuToAS0514fDizMvX4KeEkT6tTzgdLBkjol7SLpg6QkuQC4FThI0hxJncDHgV9kEyQArgXOljRN0l6k3tcC6J7ocT8wL7v2saSZfjdnZa8DjpZ0WJb4zgNuGYoJEpCbZu57fGZmdcokqEdJyxsh6bWkZJFPUHux/djNYL2DlPieAd4IvCkiNkbECtKsu08Bq4A/I60NWHMFcDvwIPBL4I7sXM0JwOys7GeA47NrEhGLgDNJieoZ0tjTWU1uVzePQZmZFSszzfwa4J8l/RKYRvrlnV9J4s9IExuaJiL+EfjHXmLfI01FL4oF8KHsKIo/RuoB9va51wPXl6vtwHgMysysWJkHdT9P2nJjI+kh12MjYj10P6j7auAbQ1DHYU0Skp+DMjNrVPZB3fNJSx41nl9Jk8efRpI2ybf4zMwaDHjLd0lTsiWHbJDaBFvdgzIzq1MqQUnaS9K/S3qO9NzT7yStkrRA0rQhqeEIkHpQTlBmZnn9vsWXbfH+U9IW7/fT89zQK4B3Am+S9OqI+E2zKzncdbSJrVudoMzM8sqMQZ0P7A4cFRF1kyEkHUlaXfx80krnVkJ7m9jiQSgzszplbvG9Gbi0MTkBRMQ3SVttvKVZFRtJRrW3sWXbtlZXw8ysUsokqN2BR/qIPwLsNqjajFDtbWKre1BmZnXKJKgn6OPhVuD12XuspI42sdljUGZmdcokqBuBv5H0aUm71k5KmijpAuBvqd/V1vqpo73NPSgzswZlJ0kcBnwY+KCk32bn9yItIns38MnmVm9k6PAkCTOz7ZRZ6mg9acPAM4DvAOuA9aT1+OYCb8i2ULeS2tvElq2eJGFmlld2qaOtwFXZYU3S0d7mHpSZWYMBL3VkzdPhWXxmZtvptQcl6Z0DuWBEXDvw6oxM7W1is2/xmZnV6esW3wIgSDvp9leQdrO1Eka1uwdlZtaorwT1hh1WixHOSx2ZmW2v1wQVET/akRUZyTra2li/aUurq2FmVimeJFEBHb7FZ2a2HSeoCvCDumZm23OCqoD0oK4TlJlZnhNUBYzuaGeTp5mbmdVxgqqAzo42Xti8tdXVMDOrFCeoCugc1e4EZWbWwAmqAjpHtfHCZt/iMzPLc4KqgM5R7bywZSsRnihhZlbjBFUBYzraiMATJczMcpygKqBzVDuAb/OZmeU4QVXAmCxBbfRECTOzbk5QFTA2S1AbnKDMzLo5QVXAhM60Zu+aF7xgrJlZTaUTlKTpkr4haZWkpyVdIqkji71R0mJJ6yUtlLRPrpwkXShpZXZcJEkN112YlV0s6YiGzz1R0nJJ6yTdJmnSULZz17GjAFi9YfNQfoyZ2U6l0gkKuBR4BtgTeCXwF8BZkqYAtwAfAyYBXcDXcuXmAscAhwAHA0cBZ+TiNwD3AZOBjwI3SZoKIGkWcAXwDmAPYH1WjyHjBGVmtr2qJ6h9ga9HxAsR8TTwLWAWcBywKCJujIgXgPnAIZIOyMqdAlwcEU9ExJPAxcCpAJJmAocC8yJiQ0TcDDwIzMnKngTcHhF3RsRaUhI8TtKEoWqkE5SZ2faqnqC+AJwgaRdJ04Aj6UlSD9TeFBHrgKXZeRrj2ff52LKIWNNHPH/tpcAmYGZj5STNldQlqWvFihUDbuSkcaMBeHbNxgFfw8xsuKl6gvoRKWE8DzxBupV3GzAeWN3w3tVArZfTGF8NjM/GocqWbYx3i4grI2J2RMyeOnVq/1vVoHNUO3tMHMPy368f8DXMzIabyiYoSW3At0ljTeOAKcDuwIXAWmBiQ5GJQK1X1BifCKyNtJZQ2bKN8SGxz+RxLF+5big/wsxsp1LZBEWa/PBy4JKI2BgRK4FrgLcCi0gTIACQNA7YLztPYzz7Ph+b0TCm1BjPX3sGMAZY0pxmFdtn0i4sX+kelJlZTWUTVEQ8CzwKvEdSh6TdSJMfHgBuBQ6SNEdSJ/Bx4BcRsTgrfi1wtqRpkvYCzgEWZNddAtwPzJPUKelY0ky/m7Oy1wFHSzosS3znAbc0jFk13fQp43hmzUZPlDAzy1Q2QWWOA94CrAB+DWwBPhARK0iz7j4FrAL+DDghV+4K4HbS7LxfAndk52pOAGZnZT8DHJ9dk4hYBJxJSlTPkMaezhqa5vU4dO/dAfjfR38/1B9lZrZTkLd4aI7Zs2dHV1fXgMtv3LKVQz7xHf5u9sv5xNsPamLNzMyqS9LPI2J2UazqPagRY0xHO6/ffyp3PPg0W7zthpmZE1SVzPmTl/Hs2o3c+cjAn6kyMxsunKAq5A1/+BImjxvNV36yvNVVMTNrOSeoChnd0capfz6dhQ+v4Fe/fb7V1TEzayknqIp552umM250O19a+OtWV8XMrKWcoCpm111G8a7X7csdDz5F12Oecm5mI5cTVAW95/D9eOnETub91yLP6DOzEcsJqoJ2Gd3BvKNfwaLfPs+/fv+RVlfHzKwlnKAq6sg/2pM5h76MSxb+moUPP9Pq6piZ7XBOUBV23ttnceCeEznrq/dy/2+ea3V1zMx2KCeoChs3poNrTvtTpkwYzcn/9jPu8gO8ZjaCOEFV3EsmdPL1M17Dy3Yfy2nX3MMlP3jEEyfMbERwgtoJ7LnrWG488zW85aCX8rnvLOHoS+5m4eJn8EK/ZjacOUHtJCZ0juKSEw/l0pMOZe3GzZy24B6O/MJdXHP3o6xcu7HV1TMzazpvt9Ekg91uo4xNW7Zx871PcP3PHufBJ1cDcNC0ibz2D6bwR9N25cA9JzJ98jja27RD6mNmNlB9bbfhBNUkOzJB5T301PN8/6Hfcecjz3Lf46vYvDX99xzd3saeu3Wy165j2Wu3sUwZP5qJY0exa+4YO7qdzo52xoxqY0xHG2M62tPXUW20t4l2ifY2ITnRmdnQcILaAVqVoPI2btnKI79by0NPPc+vV6zlqede4LfPbeDJ5zbw+3Wb2LhlYJMrJGiXaMslrTaRklibaFM68u8HUPdrbXe9+vep7nVRWTUEeru2me14Fxz7R7xq30kDKttXguoYVK2sUsZ0tHPQtF05aNquhfEXNm9l9YbN3ccLm7eycfM2Nm7Zlr7fso2NW9LXrdui+9gW2fcRbNsWbN3Gdudqf+cE6Zue1zS8rg/0xHv+UNq+TP17ovGNZtZS48a0D8l1naBGkM5R7XSOamePiZ2troqZ2YvyLD4zM6skJygzM6skJygzM6skJygzM6skJygzM6skJygzM6skJygzM6skJygzM6skL3XUJJJWAMsHWHwK8GwTq7MzcJtHBrd5ZBhMm/eJiKlFASeoCpDU1dtaVMOV2zwyuM0jw1C12bf4zMyskpygzMyskpygquHKVlegBdzmkcFtHhmGpM0egzIzs0pyD8rMzCrJCcrMzCrJCcrMzCrJCaqFJE2SdKukdZKWSzqx1XUaLEljJH05a88aSfdJOjIXf6OkxZLWS1ooaZ9cTJIulLQyOy6SpNa0pDxJ+0t6QdJXc+eGc3tPkPRQ9u93qaTDsvPDss2Spkv6hqRVkp6WdImkjiw2LNos6X2SuiRtlLSgITbgNmY/u4VZ2cWSjuhXhSLCR4sO4Abga8B44HXAamBWq+s1yDaNA+YD00l/AB0FrMleT8na+DdAJ/BZ4Ke5smcADwMvA6YBvwLObHWbSrT9O8BdwFez18O2vcCbSCunvDr77zwtO4Zzm78BLMja9VLgQeAfhlObgeOAY4DLgAW584NqI/AT4J+BscAc4Dlg6ovWp9U/kJF6ZL/INwEzc+e+Anym1XUbgrb+IvtHORf4ccPPYANwQPb6x8DcXPzd+f8JqnwAJwBfJyXnWoIazu39MfDugvPDuc0PAW/Nvf4scMVwbDPwyYYENeA2AjOBjcCEXPyu/iRp3+JrnZnA1ohYkjv3ADCrRfUZEpL2ILV1EaltD9RiEbEOWEpPm+vi7CQ/D0kTgfOAcxpCw7W97cBsYKqkX0t6IrvdNZZh2ubMF4ATJO0iaRpwJPAthnebawbTxlnAsohY00u8V05QrTOe1GXOWw1MaEFdhoSkUcB1wL9HxGJevM2N8dXA+Krer885H/hyRPym4fxwbe8ewCjgeOAw4JXAHwPnMnzbDPAj0i/V54EngC7gNoZ3m2sG08YB/65zgmqdtcDEhnMTSeM1Oz1JbaRblpuA92WnX6zNjfGJwNrI7glUkaRXAkcA/1IQHnbtzWzIvn4xIp6KiGdJ4wtvZZi2Ofv3/G3gFtLtrSnA7sCFDNM2NxhMGwf8u84JqnWWAB2S9s+dO4R0K2ynlv3V9GXSX9pzImJzFlpEamPtfeOA/ehpc12cnePncThpAsjjkp4GPgjMkXQvw7O9RMQqUg+i6BfssGwzMAl4OXBJRGyMiJXANaSkPFzbnDeYNi4CZkia0Eu8d60ejBvJB/AfpJl844DXMgxm8WXtuhz4KTC+4fzUrI1zSDOBLqR+JtCZpIHoacBe2T/gSs52ytV5F9KMrtrxOeCmrK3Drr25up8H3AO8hNSTuIt0q3M4t3kZ8BGgA9gNuJV0C3vYtDlrWyfwadIdkM7s3KDamP0++FxW9lg8i6/6B+mvstuAdcDjwImtrlMT2rQP6S/rF0hd+9pxUhY/AlhMuk30Q2B6rqyAi4DfZ8dFZOtF7iwHuVl8w7m9pDGoS7NfNE8D/wp0DvM2vzJrzyrS5nw3Ai8ZTm3O/v1GwzF/sG0k3WX4YVb2YeCI/tTHi8WamVkleQzKzMwqyQnKzMwqyQnKzMwqyQnKzMwqyQnKzMwqyQnKzMwqyQnKrMkkHS4pJJ3a6roMhKRTs/of3uq62MjmBGU2xLLN2uZn6/ZVQpZE50vardV1MeuNH9Q1a7JsYdHRwOaI2Jr1RBYCp0XEghZWrZuk+cA8YN+IeKwh1k5aKWJTRGzb8bUzSzpaXQGz4Sb7pf7Cjvo8SROifq+dQYmIrcDWZl3PbKB8i8+syfJjUNk41MIsdE12PiT9MPd+SXqPpJ9LWi9pjaSFkt7QcN3pWdn5kv4ue/8G4ItZ/ABJl0palF1jffae0xuus4DUewJ4NFen+Vm8cAxK0hRJX5L0G0mbsq9fkjS54X218n8p6YOSlkraKGmJpFMG9cO1EcU9KLOhdSdwAfD/gCtJq34D/C73nq8A/4e0Cvo1wBjgJOC7ko6LiP9quOYxwD8Al5FWjn8+O3848Hrgv4FHSavk/w1wpaQpEfHp7H1XkPbjORb4AGnhU4Bf9NYISbuStvX+A+Bq4F7SJoXvAf5S0qsKenEXAGOzz9uYvXeBpF9HxN29fZZZt1avnuvDx3A7SIkigFOLXje899gsNrfhfAdpx9ZH6Rkrnp69dzNwYMG1xhWcayOtIr0aGJU7Pz+71vSCMqdmscNz5z6VnTur4b3vzc6fX1D+PmB07vw0UqK6odX/jXzsHIdv8Zm11smknUVvy26hTZE0hbTf0O2kpLR/Q5k7IuKhxgtFxLra95I6s1tvk4DvkHpMBwyinscCK0i9wLwrSD2wYwvKXBoRm3L1e5K0UWdje8wK+RafWWsdCEyg/pZfoz1Iv9hrlhS9SdJ4Us/ob0m7vzbafWBVBGBfoCsituRPRsQWSQ8DhxaUWVZwbiVpzzCzF+UEZdZaIvVMTuzjPb9seL2+l/ddDxxF6uXcSdo4bgtpW/IPsOMnRfU2E1A7tBa203KCMht6fT1s+Agwk7R99tqBfkD2wO1RwFci4syG2BEl61RkGfCHkjryvShJHaT6F/WWzAbFY1BmQ6+WeCYVxK4l/X/46YIYkvbo52fUeit1vRNJewJ/X7JORW4DphZc6/Ts/K39vI5Zv7kHZTb0fkWaCHGWpPXAc8AzEfGDiLhJ0jXA+yQdSpoi/izwMuA1pGndM17sAyJijaTvACdnz0bdQxrrOYM0E3ByQ5GfZl8vlHQd6cHiX0ZE4+3EmotIU9a/lNXzPtI083cDD2dxs6ZyD8psiEXEBuAE0vNKnwduAD6ei78LeCewDfgn0oO3p5B6Of9U4qNOJj2jdDRwCel5qY8CXyqo093Ah4H9gKuyOh3fRxtWA68lzdp7K/Cv2dfLgddFE1eyMKvxWnxmZlZJ7kGZmVklOUGZmVklOUGZmVklOUGZmVklOUGZmVklOUGZmVklOUGZmVklOUGZmVklOUGZmVkl/X/HKyeqImgFvAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "nb_of_iterations = 1000\n",
    "step_size = 1e-5\n",
    "ridge_parameter = 0\n",
    "\n",
    "iteration, loss_function_array = gradient_descent(A_train, Y_train, nb_of_iterations, step_size, ridge_parameter, add_bias = False)\n",
    "plt.plot(iteration, loss_function_array)   \n",
    "plt.xlabel(\"iteration\", rotation=0, fontsize=18)\n",
    "plt.ylabel(\"loss function\", rotation=90, fontsize=18)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAERCAYAAAA9oHOJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAq6ElEQVR4nO3de7xXVZ3/8df7AILcROBoShFikoY/MYfp7kST028yKy80w09L7IZmTb/UyunXBdPG0nKmJrK0i5dKG69MTpblxIjdZiSVjERMEFMjjoRHDnfO+fz+WPt72Gef7zmcL5zv2Ycv7+fjsR/w3Z+99nctxPNh7bX2WooIzMzMBlpT2RUwM7N9kxOQmZmVwgnIzMxK4QRkZmalcAIyM7NSDC27AnuLiRMnxpQpU8quhpnZXuXXv/71MxHRXC3mBNRHU6ZMYcmSJWVXw8xsryJpdU8xP4IzM7NSOAGZmVkpnIDMzKwUTkBmZlYKJyAzMyuFE5CZmZXCCcjMzErhBFRnf2zdzD//+BFWtrSVXRUzs0HFCajO/vTcVv71p7/n8XUby66Kmdmg4gRUZ8p+9b5/ZmZdOQHVmbIM5ARkZtaVE1CdqbMPZGZmeU5AA8QdIDOzrpyA6mznIzinIDOzvFITkKQpku6UtF7SGkkLJA3NYiMlXSnpGUmtkhbnyknSZZLWZcflklS47yJJmyQtl3RC4XtPl7Ra0kZJCyWNr3dbnX7MzLoquwd0JbAWOAQ4FngtcG4WuxoYDxyV/Xpertw84GRgBnAMcBJwdi5+I/AAMAH4OHCLpGYASdOBq4B3AAcDm7J61IUnIZiZVVd2AjoMuCkitkTEGuBHwHRJLwbeAsyLiJaIaI+IX+fKzQWuiIgnI+Ip4ArgLABJ04DjgPkRsTkibgUeAk7Lyp4B3BERiyOiDfgkcKqkMfVo4M5JCM5AZmZ5ZSegLwFzssdtk4A3kpLQy4HVwKezR3APSTotV246sDT3eWl2rhJbGREbeol3lo2Ix4BtwLRi5STNk7RE0pKWlpbdaqB7QGZm1ZWdgO4hJYTngCeBJcBC4PnA0UArcCjwAeA6SUdl5UZnsYpWYHQ2DlSMVeJjeihbjHeKiKsjYmZEzGxurrql+S51JqDdKm1m1rhKS0CSmoC7gNuAUcBE4EDgMmAzsB34TERsi4h7gEXAG7LibcDY3O3GAm2RppoVY5X4hh7KFuP9qvIIzj0gM7OuyuwBjQdeACyIiK0RsQ64BjgR+M0uyi4jTUComJGdq8SmFsZ0ivHOspKmAsOBFbvZjl7J76GamVVVWgKKiGeAVcD7JA2VNI40uWApsBh4AvhYFns1MIvUYwK4Hjhf0iRJhwIXANdm910BPAjMlzRC0imkmXK3ZmW/C7xZ0vGSRgEXA7cVxoz6v71+CGdm1kXZY0CnAn8LtAC/B3YA50XEduCtpN5QK/B14MyIWJ6Vuwq4gzS77bfAD7JzFXOAmcB64HPA7IhoAYiIZcA5pES0ljT2cy514sVIzcyqG1rml0fEg6SeTbXYMuCVPcQC+Gh2VIs/3tN9s/gNwA211HV3eRKCmVl1ZfeA9gGVSQhOQWZmeU5AdeZJCGZm1TkB1ZnHgMzMqnMCqrPKGqmeBWdm1pUTUJ25B2RmVp0TUJ15DMjMrDonoAHiHpCZWVdOQHXWuRZcyfUwMxtsnIDqzFtym5lV5wQ0QJx+zMy6cgKqs85JCM5AZmZdOAHVmd8DMjOrzgmozvwekJlZdU5AdebVsM3MqnMCqjPhN1HNzKpxAhogfgRnZtaVE1Cd7XwE5wxkZpbnBFRnnoRgZladE1C9eRKCmVlVTkB11jkJwV0gM7MunIDqzNOwzcyqcwKqM48BmZlVV2oCkjRF0p2S1ktaI2mBpKGFa+ZLCkkn5M5J0mWS1mXH5dLOrd+y+y6StEnS8nzZLH66pNWSNkpaKGl8HdsIeDVsM7OisntAVwJrgUOAY4HXAudWgpIOB2YDfyyUmwecDMwAjgFOAs7OxW8EHgAmAB8HbpHUnN1zOnAV8A7gYGBTVo+68GuoZmbVlZ2ADgNuiogtEbEG+BEwPRdfAFwIbCuUmwtcERFPRsRTwBXAWQCSpgHHAfMjYnNE3Ao8BJyWlT0DuCMiFkdEG/BJ4FRJY+rSwoz7P2ZmXZWdgL4EzJE0UtIk4I2kJISktwHbIuLOKuWmA0tzn5eyM3FNB1ZGxIZe4p1lI+IxUoKbVvwSSfMkLZG0pKWlZXfal9uQbreKm5k1rLIT0D2khPAc8CSwBFgoaTRwKfChHsqNBlpzn1uB0dk4UDFWiY/poWwx3ikiro6ImRExs7m5ua9t6sJbcpuZVVdaApLUBNwF3AaMAiYCBwKXAZ8Gvh0Rq3oo3gaMzX0eC7RFGukvxirxDT2ULcb7l7fkNjOrqswe0HjgBcCCiNgaEeuAa4ATgdcDH8xmxq3JrrtJ0oVZ2WWkCQgVM7JzldjUwphOMd5ZVtJUYDiwoj8bt/P+9birmdneb+iuL9kpe8R1AnAEaYZZ8cdrRMQlfblXRDwjaRXwPklfID0am0san/kgMCx3+X3A+cAPs8/XA+dLupP0dOsC4MvZfVdIehCYL+kTpHGlY9g5CeG7wC8lHQ/cD1wM3FYYM+o3fg/IzKy6PicgSUcAC4Ej6Xl2cQB9SkCZU4Evkma6tQOLgPOy3lD+u9uB9dmsNUjTqKeSZrcBfCM7VzEHuBZYDzwBzI6IFoCIWCbpHFIimgDcDbyzhjrXxFtym5lVV0sP6MvA4aRk8VNgXe+X71pEPAjM6sN1UwqfA/hodlS7/vHe7hsRNwA39LWee8I9IDOz6mpJQK8BvhgRX6hXZRqRx4DMzKqrZRLCNqCnWWm2C+4AmZl1VUsCugt4db0q0qg63wNyBjIz66KWBHQ+8EpJF0jar14VajTektvMrLpaxoB+Tnph9HLgc5KeJs1cy4uIOLy/KtdI3AMyM+uqlgT0BB7KqJknIZiZVdfnBBQRs+pYj4a1cwzIudvMLK/sxUgbnlfDNjOrrqaleKBzk7i3klYiAFgJ/Hu2rYEVdL6IWmotzMwGn1rXgrsE+EdgSCF0uaRLI+JT/VazBiEPApmZVdXnR3CS3kXa3vq/gVNIC5IeQdoa+5fAxyXVbU21vZ0fwZmZdVVLD+j9pOQzKyJ25M4/lq1KfS/wAdKWCpbZ+QjOGcjMLK+WSQhHAd8rJB8AsnPfy66xHE9CMDOrrta14Eb3Eh+TXWM5O7djMDOzvFoS0H3A2ZIOLgYkHQTMIz2is2rcBTIz66KWMaBLgP8EHpb0TeB32fnppA3dxgBn9G/1GoPkHpCZWVEtKyEslnQqsIC0BXbeE8DciLi3PyvXKIQ7QGZmRTW9BxQRd0j6AfAXwGGkn62PAfdHREcd6tcQJHkWnJlZQc0rIWSJ5r7ssD7wq6hmZt15LbgB4kdwZmZd9dgDkrQK6ACOjIjtklb24X7eD6gKT0IwM+uut0dwq0k/Nys/O70f0G4Scg/IzKygx0dwETErIl5XWfkg97nXo5YvlzRF0p2S1ktaI2mBpKGSXiHpJ5L+LKlF0s2SDsmVk6TLJK3LjsuVW/Uzu+8iSZskLZd0QuF7T5e0WtJGSQslja+l3jWTl+IxMyuqZTHSyZL27yW+v6TJNX7/lcBa4BDgWOC1wLnAgcDVwBTghcAGuq4xN4+0COoM4BjgJODsXPxG4AFgAmkB1VskNWf1nA5cBbwDOBjYlNWjbgTuO5qZFdQyCWEVaRXsnrwlu6YWhwE3RcSWiFgD/AiYHhE/jIibI+K5iNhEevfo1blyc4ErIuLJiHgKuAI4C0DSNOA4YH5EbI6IW4GHgNOysmcAd0TE4ohoAz4JnCppTI117zOPAZmZdVdLAtrVbOImav85+yVgjqSRkiYBbyQloaK/ApblPk8HluY+L83OVWIrI2JDL/HOstlGetuAacUvlTRP0hJJS1paWmpqWJf7IG/JbWZWUOs07N5+ih4FPFvj/e4hJYTngCeBJcDC/AWSjgE+BXwkd3o00Jr73AqMzsaBirFKfEwPZYvxThFxdUTMjIiZzc3NfW9VgfekMzPrrtcXUSXNJT3uqviEpPdWuXQ8cDRwe1+/WFITcBdpPOZVpMTwLeAy4KPZNS8Cfgj838IyP23A2NznsUBbRISkYqwS39BD2WK8LtwBMjPralc9oHGkcZrDSL2f5tznyjElu8+3SBMI+mo88AJgQURsjYh1pIkGJwJIeiFwN3BJRHy7UHYZaQJCxQx2PqJbBkwtjOkU451lJU0FhgMraqh7TYTHgMzMinpNQBHxpYg4LCIq6759qPI5d0yNiGMiYl5E9HmgJCKeIU1aeF829Xocqbe1NBsP+inwlYj4WpXi1wPnS5ok6VDS4qjXZvddATwIzJc0QtIppJlyt2Zlvwu8WdLxkkYBFwO3FcaM+pXk94DMzIr6PAYUEU0RcUM/f/+pwN8CLcDvgR3AecB7gKmkJNJWOXLlrgLuIM1u+y3wg+xcxRxgJrAe+Bwwu5IcI2IZcA4pEa0ljf3U0nOrWeoBOQOZmeX1eTFSSS8FXhURX+kh/n7g5xHxYF/vmV07q0ro09nRU7kgjRN9tIf44z3ctxK/AejvZNozeQzIzKyolllw84E39RJ/I2m2mhV4EpyZWXe1JKC/JE2b7sk9wMv2rDqNKY0BuQtkZpZXSwKaCPy5l/iz2TVW4JUQzMy6qyUBrWXnagLVHE3vCWqf5UdwZmbd1ZKA7gbeky3m2YWklwDvzq6xKvwEzsysq1q25P4Madr0fZK+RXrXJoCXAu8irad2SX9XsBFI8jRsM7OCPiegiHhM0utJL3wW35tZBrwzIh7tx7o1DOEekJlZUS09ICJiCXC0pGOBI0g/Wx+JiKW9FtzHeRKCmVl3NSWgiuwF0gf7tSYNzUvxmJkV7VYCkjSStNtotwleEfHEnlaq0TQJvwdkZlZQy1I8TaSlb/4BeF4vlw7Z00o1miFNosMJyMysi1p6QJ8DPkyacHArsK4uNWpATRIdzj9mZl3UkoDeDvwoIk6sV2UalYR7QGZmBbW8iHog8O/1qkgja/J+QGZm3dSSgB4CDqlXRRpZk3tAZmbd1JKAPg2cI+kF9apMo/IYkJlZd7WMAf0FsBr4naTbSdtptxeuiYjwcjwFHgMyM+uulgR0Ue73b+/hmsDrwXXTJNHhLpCZWRe1JKDD6laLBuf3gMzMuqtlMdLV9axII5PHgMzMuqllEoLtJi/FY2bWXS1L8XyrD5dFRLx7D+rTkDwLzsysu1p6QGf18egzSVMk3SlpvaQ1khZIGprFXi9puaRNkhZJemGunCRdJmlddlwuSYX7LsrKLpd0QuF7T5e0WtJGSQslja+l3rXye0BmZt31OQFFRFPxAIYBLwa+DvyKtFpCLa4E1pJecD0WeC1wrqSJwG3AJ4HxwBLg33Ll5gEnAzOAY4CTgLNz8RuBB0grdn8cuEVSM0C2pfhVwDuAg4FNWT3qxmNAZmbd7dEYUES0R8SjEXE2aXHSy2q8xWHATRGxJSLWAD8CppO2/l4WETdHxBbSFPAZko7Mys0FroiIJyPiKeAKst6XpGnAccD8iNgcEbeSVnE4LSt7BnBHRCyOiDZSkjtV0pjd+TPoC48BmZl115+TEH7Izh/yffUlYI6kkZImAW9kZxLq3GU1IjYCj2XnKcaz3+djKyNiQy/x/L0fA7YB04qVkzRP0hJJS1paWmps2k5NEu3uApmZddGfCWgCMLrGMveQEsJzwJOkR20Ls/u0Fq5tBSq9lGK8FRidjQPVWrYY7xQRV0fEzIiY2dzc3PdWFTT5PSAzs272OAFJGidpNnAe8OsayjUBd5HGekYBE0ljSJcBbcDYQpGxQKVXU4yPBdoiPeeqtWwx3u/SJIR63d3MbO/U5wQkqUNSe/Egjf3cBHQA59fw3eOBFwALImJrRKwDrgFOJG16NyP33aOAw7PzFOPZ7/OxqYUxnWI8f++pwHBgRQ11r0najsEZyMwsr5aleK4nrfWWF8CfST+8byyMu/QqIp6RtAp4n6QvkB6NzSWNz9wOfF7SacAPgE8Bv4mI5bm6nC/pzqwOFwBfzu67QtKDwHxJnyCNKx3DzvGp7wK/lHQ8cD9wMXBbLXWvld8DMjPrrscEJOmvgIcjogUgIs6qw/efCnwRuJC0svYi4LyIaMmSzwLgO8B/A3Ny5a4CppJmtwF8IztXMQe4FlgPPAHMzrVjmaRzSIloAnA38M46tK2TV8M2M+uutx7QItK7MjcASFoJfCgivt9fXx4RDwKzeojdDRzZQyyAj2ZHtfjjPd03i99A1q6B4B6QmVl3vY0BbSWNjVRMofZZbobfAzIzq6a3HtAKYK6k+0mPsgAmSJrc2w0j4on+qlyj8HtAZmbd9ZaAPkN6THV/9jlI4zVf3MU9h+xxrRpMeg+o7FqYmQ0uPSagiLhF0lLSWMohwHzSS6K/GZCaNRA/gjMz667XadgR8SjwKICki4BbswF8q0GahOAEZGaWV8uOqN68bjd5FpyZWXdOKgPA7wGZmXXnBDQA0lI8ZdfCzGxwcQIaAN4R1cysOyegAeD3gMzMunMCGgBNTX4EZ2ZW5AQ0APwIzsysu1r2A3qZpPcWzr1V0kOSnpJ0af9XrzH4PSAzs+5q6QHNB95S+ZCtCXcj8DzSltYXSqrrtgZ7K0m0d5RdCzOzwaWWBDQD+Hnu8xxAwLER8RLgx8C8fqxbwxg2RLR3OAOZmeXVkoAmAGtyn/83sDginso+fx84or8q1kiGNIkd7X4EZ2aWV0sCehY4GEDScOAVwOJcPID9+61mDWTYkCZ2eBq2mVkXfV4LDngQeI+ku4FTgBHAXbn4YcCf+q9qjWNIk9jhR3BmZl3UkoAuIY3z/A9p7OcnEbEkFz8J+O9+rFvDGNYk94DMzApqWQ37F5KOI439tALfq8QkTSAlp9v7vYYNYEhTExHQ3hEMaVLZ1TEzGxRq6QEREStIW3UXz68DzuuvSjWaoUNS0tne3sGQJm8Ya2YGtb2IOkTSyMK5cZIukPRPko6u5YsltRWOdklfzsX/TtLDkjZI+p2kk3MxSbpM0rrsuFyScvEpkhZJ2iRpuaQTCt99uqTVkjZKWihpfC11r9XQrNfj9eDMzHaqZRbcVaTxHwAkDQN+Bnwe+Bhwn6Rj+3qziBhdOUiz6zYDN2f3ngR8BzgfGAt8BLhB0kFZ8XnAyaR3k44hjT+dnbv9jcADpKnjHwdukdSc3Xt61pZ3ZN+7Cbiyr/XeHUOHpD9mT8U2M9uplgT0GtK7PhWzgZcA7wdeRZoB94+7WY/ZwFrg3uzz84FnI+KHkfwA2AgcnsXnAldExJPZe0hXAGcBSJoGHAfMj4jNEXEr8BBwWlb2DOCOiFgcEW3AJ4FTJY3ZzbrvUqUH5JlwZmY71ZKADgFW5T6/CVgWEV+NiF8BVwOv3M16zAWuj+hcMG0J8LCkt2SP/k4GtgK/yeLTgaW58kuzc5XYyojY0Eu8s2xEPAZsA6btZt13qTIG5JlwZmY71TIJQUB+BH0WcFvu8x+Bg6hRtqbca4F3V85FRLuk64EbSO8bbQPeFhEbs0tGk2biVbQCo7NxoGKsEp/UQ9lKvFsPSNI8suWFJk+eXGvTOg1ryh7BOQGZmXWqpQe0ijQFG0mvJvWIFuXih9L9B3tfnAn8LCI6e1fZpIHLSUluP1KC+kZujKmNNDZUMRZoy3pQxVglvqGHssV4p4i4OiJmRsTM5ubm2luWqUy93uEVSc3MOtWSgK4B3irpt8B/kMZs8ishvBxYvht1OBO4rnDuWNI6c0sioiMi7iO95FqZzbaMNAGhYkZ2rhKbWhjTKcY7y0qaCgynyvTy/uJHcGZm3fU5AUXEF0lbMmwlzTA7JSI2QeeLqK8A7qzlyyW9ivRo7OZC6D7g+EqPR9JLgePZOQZ0PXC+pEmSDgUuAK7N6rmCtGzQfEkjJJ1Cmil3a1b2u8CbJR0vaRRwMXBbYcyoXw1t8iw4M7OiWl9EvYS0JE/x/Dp2Y/yHNPmg2w//iLhH0kWk6dMHAy3ApRHx4+ySq4CppNltAN/IzlXMISWk9cATwOyIaMnuvUzSOaRENAG4G6jrPkY7e0B+BGdmVlFTAsqTNBEgIp7Z3XtExNm9xBYAC3qIBfDR7KgWf5w0ftTTvW8gTXAYEMM6V0JwD8jMrKKWMSAkHSrpOknPkt77+ZOk9ZKuzV4etSqGD02TB7duby+5JmZmg0efe0DZdOlfkbbgfpCdg/ovIU0k+BtJr4iIP/R3Jfd2w4emPL9lhx/BmZlV1Lodw4HASRHRZbKBpDeS3gm6hGxFAttpxDD3gMzMimp5BPcG4Mpi8gGIiB8CXwX+tr8q1khGDHMPyMysqJYEdCDwaC/xR4Fxe1SbBlUZA9riHpCZWadaEtCT9DKzDPir7BorGJ71gPwIzsxsp1oS0M3A2yR9VtIBlZOSxkq6FPg74N/6u4KNoHMMyI/gzMw61ToJ4XjgQuDDkp7Ozh9KWqT058Bn+rd6jWGEH8GZmXVTy1I8m0iLgp4N/Ji0P88m0npw84DXRcTmelRybzdsiBjSJDZtcwIyM6uodSmeduDr2WF9JIkD9h9G6+btZVfFzGzQqGklBNt940YO41knIDOzTj32gCSduTs3jIjrd786jWvc/sN4dtO2sqthZjZo9PYI7logSDuh9lWQtkqwgnEj9+NPz20puxpmZoNGbwnodQNWi33AuJHDeGRN3bYcMjPb6/SYgCLinoGsSKMbt/9+noRgZpbjSQgDZNzIYbRt3cE2v4xqZgY4AQ2YA0cOA+DZzZ6IYGYGTkAD5qCxIwBY0+qJCGZm4AQ0YCaPHwnA6nWbSq6Jmdng4AQ0QCoJ6Ik/OwGZmYET0IAZNXwoE0fvxxPuAZmZAU5AA+qFE0ax6pmNZVfDzGxQKC0BSWorHO2SvpyLj5R0paRnJLVKWpyLSdJlktZlx+WSlItPkbRI0iZJyyWdUPju0yWtlrRR0kJJ4weizUc+bwwP//E5ImIgvs7MbFArLQFFxOjKARwMbCZteldxNTAeOCr79bxcbB5wMjADOAY4ibRNRMWNwAPABODjwC2SmgEkTQeuAt6Rfe8m4Mp+bl5V0w89gA1bd/Dkeu9aYWY2WB7BzQbWAvcCSHox8BZgXkS0RER7RPw6d/1c4IqIeDIingKuAM7Kyk4DjgPmR8TmiLgVeAg4LSt7BnBHRCyOiDbgk8CpksbUu5EvOXQsAMuefq7eX2VmNugNlgQ0F7g+dj6bejmwGvh09gjuIUmn5a6fDizNfV6anavEVkbEhl7inWUj4jFgGzCtvxrTkxcfPIYmwe+ebq33V5mZDXqlJyBJk0k7rV6XO/184GiglbTl9weA6yQdlcVHZ7GKVmB0Ng5UjFXiY3ooW4zn6zZP0hJJS1paWmptWjf77zeEFx00mqVPOgGZmZWegIAzgZ9FxKrcuc3AduAzEbEtWxh1EfCGLN4GjM1dPxZoy3pQxVglvqGHssV4p4i4OiJmRsTM5ubm2ltWxcwp47l/9XraOzwRwcz2bYMlAV1XOPebXZRZRpqAUDEjO1eJTS2M6RTjnWUlTQWGAytqq/bu+cspB7Jh6w6Wr/E4kJnt20pNQJJeBUyi6+w3gMXAE8DHJA2V9GpgFnBXFr8eOF/SJEmHAheQNtAjIlYADwLzJY2QdAppptytWdnvAm+WdLykUcDFwG2FMaO6+cspacb3fav+PBBfZ2Y2aJXdA5pLlR/+EbEdeCtwIml85uvAmRGxPLvkKuAO0uy23wI/yM5VzAFmAuuBzwGzI6Ilu/cy4BxSIlpLGvs5tx6Nq+b5B47k0ANGcN/j6wfqK83MBqXedkStu4g4u5fYMuCVPcQC+Gh2VIs/Tuox9XTvG4Abaqhqv3rZYeO599Fn6OgImppq2fHczKxxlN0D2ie99sXNrNu4jd885dlwZrbvcgIqwWunHYQEP12+tuyqmJmVxgmoBONH7cdxkw9kkROQme3DnIBK8tdHHsRDT7Xy9LNeF87M9k1OQCV50/86BIDvL3265JqYmZXDCagkUyaO4rjJ41j4wFNlV8XMrBROQCU65aWTWL5mA8u8OKmZ7YOcgEr05hmHMmJYE9f/YnXZVTEzG3BOQCUaN3I/Tjvu+dz+4FM807a17OqYmQ0oJ6CSves1h7FtRwfX/eLxsqtiZjagnIBKdnjzaN50zCF882erWLthS9nVMTMbME5Ag8BH3vBitu3o4F9+MiA7QpiZDQpOQIPAlImjmPuqKdz4P3/gVyvXlV0dM7MB4QQ0SFzwhmlMHj+Sj97yG1o3bS+7OmZmdecENEiM3G8o//L3M/hj62Y++L0H2NHeUXaVzMzqygloEPmLF47n0285mntWtHD+TUudhMysoZW6IZ11d/rLJ/Ps5m1c/qNH2LK9nX/++2MZPdz/mcys8bgHNAidO+tFzH/zS/jP5Wt564Kf8cAT3r7bzBqPE9Ag9c5XH8a33/0yNm5t59Sv/oKP3LyUx5/ZWHa1zMz6jSKi7DrsFWbOnBlLliwZ8O9t27qDL/5kBdf/ajU72jt43YsP4q0vncTrjzyIUX40Z2aDnKRfR8TMqjEnoL4pKwFVrN2whWt//ji33f8Ua57bwtAmcewLxvHyqeM56pCxHPm8MUyZMIqhQ9ypNbPBwwmoH5SdgCraO4L/WfVn7lnRwi9XruO3T7XS3pH+Gw5tEgePHcEhB4zgeQeM4OCxIxg7Yhhj9x/KmBHDGDsi/TpiWBP7DW1i+NAhDB/axPChOz8PGyKGNAlJJbfUzBpBbwmotGc4ktoKp/YHroyIfyhcNx+4CPibiLg7Oyfgc8B7ssu+CVwYWTaVNAW4Bng58ATwgUrZLH468FlgIvAT4F0R8ef+bF+9DGkSrzx8Aq88fAIAW7a38/u1bTyyZgO/b2ljTesW/ti6mWVPP8dPl69l07b2PfquJkGTRJMqialyvnKkz/l01VPyqpzOh/Mlu56vfr8ud1bv15pZ/5g1rZlPnPSSfr9vaQkoIkZXfi9pFPAn4Ob8NZIOB2YDfywUnwecDMwAgpREVgJfy+I3Ar8ETsyOWyQdEREtkqYDVwFvAu4HrgauBOb0Y/MGzIhhQzh60gEcPemAqvHt7R20bdnBc1u2s2HLDp7bvJ2tOzqyo51t2e8rv25v76Ajgo6OoCOgPaLr544gIrLz0NERnT0wSP8xOn+f+xCVSJdz+Wv7co/u13fpv7szb1YXh4zbvy73HSyj2LOBtcC9hfMLgAtJCSJvLnBFRDwJIOkK4L3A1yRNA44D3hARm4FbJX0IOI2UoM4A7oiIxVnZTwIPSxoTERvq0bgyDRvSxIGj9uPAUfuVXRUzsy4Gy4j1XOD6yP0zWNLbgG0RcWeV66cDS3Ofl2bnKrGVhWRSjHeWjYjHgG3AtD1thJmZ9V3pPSBJk4HXAu/OnRsNXAq8oYdio4HW3OdWYHQ2NlSMVeKTeihbiY+pUrd5pMd9TJ48uQ+tMTOzvhoMPaAzgZ9FxKrcuU8D3y6cy2sDxuY+jwXash5UMVaJb+ihbDHeKSKujoiZETGzubm5T40xM7O+GSwJ6LrCudcDH5S0RtIa4AXATZIuzOLLSBMQKmZk5yqxqZLG9BLvLCtpKjAc8G5wZmYDqNQEJOlVpEdjNxdCrweOBo7NjqeBs4GvZPHrgfMlTZJ0KHABcC1ARKwAHgTmSxoh6RTgGODWrOx3gTdLOj6bfXcxcFsjTkAwMxvMyh4DmkuVH/4R0WVbUEntwPqIqLw7dBUwFXgo+/yN7FzFHFJCWk96D2h2RLRk914m6RxSIpoA3A28sx/bZGZmfeCVEPposKyEYGa2N+ltJYTBMAZkZmb7IPeA+khSC7B6D24xEXimn6qzN9jX2gtu877Cba7NCyOi6jRiJ6ABImlJT93QRrSvtRfc5n2F29x//AjOzMxK4QRkZmalcAIaOFeXXYEBtq+1F9zmfYXb3E88BmRmZqVwD8jMzErhBGRmZqVwAjIzs1I4AdWRpPGSbpe0UdJqSaeXXac9JWm4pG9m7dkg6QFJb8zFXy9puaRNkhZJemEuJkmXSVqXHZdnezjtNSQdIWmLpO/kzjVsmyXNkfRw9nf4MUnHZ+cbss2Spki6U9L6bDX+BZKGZrG9vs2SPiBpiaStkq4txHa7fdmf26Ks7HJJJ/SpQhHho04HcCPwb6RN8F5D2vhuetn12sM2jQIuAqaQ/gFzEmkvpSmkt6VbgbcBI4DPA7/KlT0beAR4PmkV9N8B55Tdphrb/2PS1vHfyT43bJuBvyGt/vGK7L/1pOxo5DbfSVrIeATwPNKCxx9slDYDpwInA18Frs2d36P2Ab8E/hnYHzgNeBZo3mV9yv4DadQj+0G9DZiWO/dt4HNl160Obf1N9pduHvCLwp/BZuDI7PMvgHm5+Lvzf8kH+0FaZf0mUgKuJKCGbXNW93dXOd/IbX4YODH3+fOklfYbqs3AZwoJaLfbB0wDtgJjcvF7+5KA/QiufqYB7ZH2J6pYCkwvqT51IelgUluXkdq2tBKLiI3AY+xsc5c4e9Gfh6SxpL2jLiiEGrLNkoYAM4FmSb+X9GT2OGp/GrTNmS8BcySNlDQJeCPwIxq7zbBn7ZsOrIyu2+r0qf1OQPUzmtSlzWsFxlS5dq8kaRhpX6XrImI5u25zMd4KjB6Mz8qruAT4ZkT8oXC+Udt8MDAMmA0cT9oY8qXAJ2jcNgPcQ/rB+RzwJLAEWEhjtxn2rH27/bPOCah+2oCxhXNjSeMlez1JTaRHituAD2Snd9XmYnws0BZZn32wknQscALwL1XCDdlm0uMXgC9HxB8j4hnSM/4TadA2Z3+n7wJuIz2CmggcCFxGg7Y5Z0/at9s/65yA6mcFMFTSEblzM0iPqvZq2b96vkn6V/JpEbE9Cy0jtbFy3SjgcHa2uUucvefPYxZpksUTktYAHwZOk3Q/DdrmiFhP6gFU+wHakG0GxgMvABZExNZIOzNfQ0q6jdrmij1p3zJgqqQxPcR7VvZgWCMfwPdIM+FGAa+mAWbBZe36GvArYHThfHPWxtNIM2kuo+tMmnNIg7yTgEOzv6CDbqZQlfaOJM2IqhxfAG7J2tuQbc7qfjFwH3AQqSdwL+lRZCO3eSXwj8BQYBxwO+kxc0O0OWvXCOCzpCcYI7Jze9S+7OfBF7Kyp+BZcOUfpH9RLQQ2Ak8Ap5ddp35o0wtJ/yreQup6V44zsvgJwHLSI5z/Aqbkygq4HPhzdlxOth7h3nSQmwXXyG0mjQFdmf0wWQP8KzCiwdt8bNae9aQN2G4GDmqUNmd/d6NwXLSn7SM9IfivrOwjwAl9qY8XIzUzs1J4DMjMzErhBGRmZqVwAjIzs1I4AZmZWSmcgMzMrBROQGZmVgonILMaSZolKSSdVXZddoeks7L6zyq7LrZvcwIy20PZZlwXZWvGDQpZkrxI0riy62LWE7+IalajbNHK/YDtEdGe9SQWAe+MiGtLrFonSRcB84HDIuLxQmwIaZWDbRHRMfC1M0uGll0Bs71N9kN7y0B9n6Qx0XWvlT0SEe1Ae3/dz2x3+RGcWY3yY0DZONCiLHRNdj4k/Vfuekl6n6RfS9okaYOkRZJeV7jvlKzsRZL+Prt+M/DlLH6kpCslLcvusSm75r2F+1xL6v0ArMrV6aIsXnUMSNJESV+R9AdJ27JfvyJpQuG6Svm/lvRhSY9J2ipphaS5e/SHa/sU94DM9sxi4FLg/wFXk1aMBvhT7ppvA/+HtIL2NcBw4AzgJ5JOjYjvF+55MvBB4Kuklcefy87PAv4K+A9gFWmV9bcBV0uaGBGfza67irQfyynAeaRFNSFtnV6VpANI2y6/CPgWcD9pA7r3AX8t6WVVemGXAvtn37c1u/ZaSb+PiJ/39F1mncpendWHj73tICWCAM6q9rlw7SlZbF7h/FDSbpur2DkWOyW7djtwVJV7japyrom0CnErMCx3/qLsXlOqlDkri83Knfun7Ny5hWvfn52/pEr5B4D9cucnkRLRjWX/N/Kxdxx+BGdWX28n7Qy5MHvENVHSRNJeM3eQks4RhTI/iIiHizeKiI2V30sakT0aGw/8mNTjOXIP6nkK0ELqxeVdRepBnVKlzJURsS1Xv6dIGzEW22NWlR/BmdXXUcAYuj6SKzqY9IO7YkW1iySNJvVs/o60c2fRgbtXRQAOA5ZExI78yYjYIekR4LgqZVZWObeOtGeU2S45AZnVl0g9i9N7uea3hc+berjuBuAkUi9lMWljsB2kLaPPY+AnFfU0k04DWgvbazkBme253l6mexSYRtreuG13vyB7ofQk4NsRcU4hdkKNdapmJfBiSUPzvSBJQ0n1r9bbMdsjHgMy23OVxDK+Sux60v9nn60SQ9LBffyOSm+jS+9C0iHAe2qsUzULgeYq93pvdv72Pt7HrM/cAzLbc78jTTQ4V9Im4FlgbUT8NCJukXQN8AFJx5GmUD8DPB94JWna89RdfUFEbJD0Y+Dt2btB95HGWs4mzaSbUCjyq+zXyyR9l/Ti7G8jovi4r+Jy0pTur2T1fIA0DfvdwCNZ3KxfuQdktociYjMwh/S+zheBG4FP5eLvAs4EOoCPkV4snUvqpXyshq96O+kdnTcDC0jvC30c+EqVOv0cuBA4HPh6VqfZvbShFXg1adbbicC/Zr9+DXhN9ONKDGYVXgvOzMxK4R6QmZmVwgnIzMxK4QRkZmalcAIyM7NSOAGZmVkpnIDMzKwUTkBmZlYKJyAzMyuFE5CZmZXi/wNfghr8VyMvnQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "nb_of_iterations = 1000\n",
    "step_size = 1e-5\n",
    "ridge_parameter = 0.2\n",
    "\n",
    "iteration, loss_function_array = gradient_descent(A_train, Y_train, nb_of_iterations, step_size, ridge_parameter, add_bias = False)\n",
    "plt.plot(iteration, loss_function_array)   \n",
    "plt.xlabel(\"iteration\", rotation=0, fontsize=18)\n",
    "plt.ylabel(\"loss function\", rotation=90, fontsize=18)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 2 : Add a l-1 regularization term to your objective function from Part 1 or Part 3 and solve the resulting problem. Can you find a value of the regularization parameter that yields a sparse solution? Does it provide a good value for the data-fitting term?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Regularization is obtained by introducing a penalty. It is often called LASSO regression, and is defined as $$ \\min_{x} f(x) = \\lVert Ax-y \\rVert_{2}  ^2+\\lambda \\lVert x \\rVert_{1}$$ with $\\lambda$ > 0 the regularization parameter, $x$ the paramer to optimize, $y$ the targe value and $A$ the matrix containing the data.\n",
    "\n",
    "LASSO optimsiation usually yields sparse solution for $\\lambda$ high enough. In our case for a $\\lambda = 3000 $ (high value as we have lot of samples and we don't apply any division by the number of rows), we reach a sparse optimal parameter. As for the data fitting term, the results are quite similar to the ones without regularization (consequence of underfitting and capping of y values in the dataset, see part 1 question 4 for more explanations). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAacAAAESCAYAAABZ6BpeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAtn0lEQVR4nO3deZxU1Zn/8c+3F2iwaWRpUDSCCESFiEsnOjFGHbck6sRtMrjFJYrGOL/J6GSSGTUaNSYmcSaL45qocaLGjKCJo6PRiJoYt8YNEQRBQVygRQSaXXh+f9xbUBTdDdVd3VV0fd+vV7266px7bj2nTXj63HvuOYoIzMzMSklFsQMwMzPL5eRkZmYlx8nJzMxKjpOTmZmVHCcnMzMrOU5OZmZWcqqKHUB3MXDgwBg2bFixwzAz26pMnjz5g4iozy0vanKSdD5wOvAp4K6IOD0t7wHcCTQAQ4GDI+LxrHYCfgiclRb9Cvh2pA9tSRoG3ArsC8wFzo+IR7PanwT8ABgIPAKcGREfpnU9geuBE4DlwI8i4j8215dhw4bR2NjYjt+CmVn5kjSnpfJiX9Z7F7gSuKWFur8ApwDvt1A3HjgGGAvsARwFnJNVfxfwIjAAuAi4R1I9gKTRwI3AqcBgkgR0XVbby4CRpEkR+FdJX2hP58zMrH2KmpwiYmJE3AcszClfHRE/jYi/AGtbaHoacE1EzIuId4BrSEZgSBoF7A1cGhErImICMAU4Pm17MnB/RDwZEc3AJcBxkvqk9V8FroiIRRExDbg5c24zM+saxR45tddo4OWszy+nZZm62RGxtI369W0jYhawGhglqR8wpI1zb0TSeEmNkhqbmpo60B0zM8u2tSanWmBx1ufFQG16Lyq3LlPfp5W22fW1WZ9baruRiLgpIhoioqG+fpP7eWZm1k5ba3JqBuqyPtcBzemEiNy6TP3SVtpm1zdnfW6prZmZdYGtNTlNJZkMkTE2LcvUDc+6h9RS/fq2koYDPYEZEbEIeK+Nc5uZWRcoanKSVCWpBqgEKiXVSKpK63qmdQA90jqln28HLpC0g6QhwIXAbQARMQN4Cbg0bXMsyYy+CWnbO4CjJR0gaRvgcmBi1j2q24GLJfWTtCtwdubcZmbWNYo9croYWAF8h2Ta+Iq0DOD19PMOwMPp+6Fp3Y3A/SSz8F4FHkjLMsaRPCO1iOR5qBMiogkgIqYC55IkqQUk95POy2p7KTALmAM8Afw4Ih4qVIdz/WnafK57/I3OOr2Z2VZJ3mywMBoaGqI9D+Fe9oepTHhhHlMuO6ITojIzK22SJkdEQ255sUdOZa+uVzXNqz5m3Tr/kWBmluHkVGR9e1UTAUtXflzsUMzMSoaTU5H17VUNwOIVa4ociZlZ6XByKrK6mmTtXScnM7MNnJyKzCMnM7NNOTkVWd/eSXJastLJycwsw8mpyDxyMjPblJNTkdXVODmZmeVyciqy3j0qqaqQk5OZWRYnpyKTRN9e1SxxcjIzW8/JqQT07VXtkZOZWRYnpxJQ5+RkZrYRJ6cSUOfLemZmG3FyKgG+rGdmtjEnpxLQt1cVS7zwq5nZek5OJSAzcvLeWmZmCSenEtC3VzVr1wXLVq8tdihmZiXByakEeJUIM7ONOTmVgMz6ep6xZ2aWcHIqAV781cxsY05OJaDOycnMbCNOTiXAIyczs405OZWAOt9zMjPbSFGTk6TzJTVKWiXptpy6QyRNl7Rc0iRJQ7PqLpO0RlJz1mt4Vv2wtM3y9ByH5pz7JElzJC2TdJ+k/ll1PSXdImmJpPclXdCJvwIA+vSsQnJyMjPLKPbI6V3gSuCW7EJJA4GJwCVAf6ARuDun7d0RUZv1mp1VdxfwIjAAuAi4R1J9eu7RwI3AqcBgYDlwXVbby4CRwFDgYOBfJX2h411tXUWFqKvxEkZmZhlFTU4RMTEi7gMW5lQdB0yNiP+JiJUkCWOspF03d05Jo4C9gUsjYkVETACmAMenh5wM3B8RT0ZEM0kCPE5Sn7T+q8AVEbEoIqYBNwOnd6SfW8Lr65mZbVDskVNrRgMvZz5ExDJgVlqecbSkDyVNlfT1nLazI2JpVtnLWW1zzz0LWA2MktQPGJJdn9N2I5LGp5clG5uamvLt40bqelU5OZmZpUo1OdUCi3PKFgOZ0c3vgN2AeuBs4LuSTtzCtm3V12Z9bqntRiLipohoiIiG+vr6zfWpTR45mZltUKrJqRmoyymrA5YCRMRrEfFuRKyNiL8CPwNO2JK2m6lvzvrcUttO07dXtVcmNzNLlWpymgqMzXyQtA2wS1rekgCU1XZ41j0k0nNNzarPPvdwoCcwIyIWAe9l1+e07TQeOZmZbVDsqeRVkmqASqBSUo2kKuBeYIyk49P67wKvRMT0tN2XJfVT4jPA/wN+DxARM4CXgEvT8x0L7AFMSL/2DpL7VQekSe9yYGLWParbgYvT8+9Kctnwts7+XXirdjOzDYo9croYWAF8BzglfX9xRDSRzK77PrAI2BcYl9VuHPAGyeW224GrI+LXOfUNadsfAiek5yQipgLnkiSpBST3k87LanspyeSLOcATwI8j4qHCdblldTXVrP54HSvXeNsMM7OqYn55RFxGMk28pbpHgRanjkfEiS2VZ9W/BRzURv2dwJ2t1K0CzkxfXSZ7ZfKa6squ/Gozs5JT7JGTpby+npnZBk5OJcLJycxsAyenEuFtM8zMNnByKhEeOZmZbeDkVCK8VbuZ2QZOTiWiriaZOLl4hVeJMDNzcioRVZUV1Pb04q9mZuDkVFLqapyczMzAyamk1PWqZslKJyczMyenEuLFX83MEk5OJaRvr2rP1jMzI8+19SQJOBQYCQxgwzYVGRERVxQotrLjkZOZWWKLk5OkkcB9JIux5ialjACcnNrJ22aYmSXyGTn9gmTDv28DjwELOyWiMta3VzXLV69lzdp1VFf6iquZla98ktPngJ9GxE86K5hyl71KxIDankWOxsysePL583w18GZnBWJeX8/MLCOf5PQwsH9nBWJOTmZmGfkkpwuAv5F0oaQenRVQOavrlVlfz8nJzMpbPvecngK2AX4E/FDSu8DanGMiInYpVHDlZv09p5Ve/NXMyls+yWkuyVRx6yTecNDMLLHFySkiDurEOAyoq/GeTmZm4OWLSkpNdSU9qyo8cjKzspfX8kUAknYBvgwMT4tmA7+PiFmFDKxceX09M7M8R06SrgCmAz8BzktfPwFel3R5vl8u6XxJjZJWSbotp+4QSdMlLZc0SdLQrDpJulrSwvT1o3Tdv0z9sLTN8vQch+ac+yRJcyQtk3SfpP5ZdT0l3SJpiaT3JV2Qb786wuvrmZnlkZwknQlcBDwLHEuy+OtI4BjgaeAiSWfk+f3vAlcCt+R810BgInAJ0B9oBO7OOmR8+r1jgT2Ao4BzsurvAl4kWZz2IuAeSfXpuUcDNwKnAoOB5cB1WW0vS/s1FDgY+FdJX8izX+3m5GRmlt/I6RskiemgiPh9RMxKX38g+Uf8OeD8fL48IiZGxH1suk7fccDUiPifiFhJkjDGSto1rT8NuCYi5kXEO8A1wOkAkkYBewOXRsSKiJgATAGOT9ueDNwfEU9GRDNJAjxOUp+0/qvAFRGxKCKmATdnzt0VnJzMzPJLTrsBv42ITR7CSct+mx5TCKOBl7POvwyYlZZvUp++z66bHRFL26jPPvcskqWZRknqBwxp49wbkTQ+vSzZ2NTUlFcHW+OVyc3M8l9br7aN+j7pMYVQCyzOKVucfkdL9YuB2vS+U75ts+trsz631HYjEXFTRDREREN9fX2bHdpSnhBhZpZfcnoeOEfS4NwKSYNI7gM9W6C4moG6nLI6YGkr9XVAc0REO9pm1zdnfW6pbaer61XN0lUfs26dn3c2s/KVT3K6AtgemCbpx5LOSF8/AaYB25FMbiiEqSSTHQCQtA3JXlJTW6pP32fXDc+6h9RSffa5hwM9gRkRsQh4r41zd7q+vaqJgKVewsjMytgWJ6eIeJJkosJS4ELgV+nrgrTsuIj4cz5fLqlKUg1QCVRKqpFUBdwLjJF0fFr/XeCViJieNr0duEDSDpKGpPHclsY5A3gJuDQ937EkM/ompG3vAI6WdECa9C4HJmbdo7oduFhSv3QCxtmZc3cFr0xuZpbnQ7gRcb+kB4B9gJ1JtmufBbwQEeva8f0XA5dmfT4F+F5EXCbpeOBa4DcklwvHZR13I8lDwFPSz79MyzLGkSSURSRrAp4QEU1pH6ZKOpckSQ0AHgWyp8BfClwPzAFWAFdHxEPt6Fu71NUk/0mWrHRyMrPypeQ2jXVUQ0NDNDY2dvg8z85eyD/c9Ax3nLUv+48YWIDIzMxKl6TJEdGQW+619UpM396+rGdm1uplPUlvAuuAXSNijaTZW3A+7+fUQb7nZGbW9j2nOST7N2Wu+3k/py6Q2TbDycnMylmrySl3/ybv59Q1eveopKpCfhDXzMpaPgu/7iSpVxv1vSTtVJiwypckr69nZmUvnwkRb5KsRt6av0uPsQ5ycjKzcpdPctJm6ivwPamC8OKvZlbu8p1K3lby2Q34qP2hWEadF381szLX5goRkk4j2Tsp42JJZ7dwaH9gDMmyQ9ZBfXtV8/aHy4sdhplZ0Wxu+aJtSZYpgmTUVA/0zjkmsxL4LSS7zloH9e1V5ct6ZlbW2kxOEfEz4GcAktYB34yIO7sisHKWmRARESRbVJmZlZctXvg1IrzUURfp26uateuCZavXUtszr7V5zcy6hXyec9pL0jfaqP+GpD0LElWZy6wS4UkRZlau8hkNXQoc2Ub9F0n2XbIO8vp6Zlbu8klOnwaeaKP+CeAzHQvHwMnJzCyf5DQQ+LCN+o/SY6yD6pyczKzM5ZOcFgCj26gfQ9vJy7aQR05mVu7ySU6PAmdJ2iRBSdod+Fp6jHVQZuTkCRFmVq7ymad8JXAc8LykW4CXSB7A3Qs4E1gNXFHoAMtRn55VSE5OZla+8nnOaZakQ4DbgPNyqqcCZ0TEzALGVrYqKkRdjRd/NbPyldcTnhHRCIxJn2caSbJS+esR8XInxFbWvG2GmZWzdi0/EBEvkVzWs05S5/X1zKyMtSs5SeoNDKCFPZ4iYm5Hg7Jk5LRk5cfFDsPMrCjyWb6oQtJ3JL0DLAXeItn5NvdVMJJ2k/SYpMWS3pB0bFo+TFJIas56XZLVTpKulrQwff1IWSuopu0nSVouabqkQ3O+9yRJcyQtk3SfpP6F7NeW8GU9Mytn+Yycfgj8C8nkhwnAwk6JKCWpCvg9cANwGHAgcL+kvUhmBgJsGxEtDS/GA8cAY0lmFD4CzE7PBXAX8DTwpfR1j6SREdGUTpW/kWSppheAm4DrgHGF7mNbnJzMrJzlk5xOAR6KiC91VjA5dgWGAP8ZEQE8Jukp4FTgV5tpexpwTUTMA5B0DXA2cIOkUcDewOERsQKYIOmbwPEkyetk4P6IeDJtewkwTVKfiFha6E62xlu1m1k5y+ch3H4kI5mu0tJGRiJZiSJjjqR5km6VlL100mggewbhy2xY3WI0MDsn0eTWr28bEbNIRmqjNglGGi+pUVJjU1PTFnZry9TVVLP643WsXLO2oOc1M9sa5JOcpgDbd1YgLZhOsmTStyRVSzqc5NJeb+ADkoVohwL7AH2AO7La1gKLsz4vBmrT+065dZn6Pq20za1fLyJuioiGiGior6/Pv4dt8BJGZlbO8klO3wPOlfSJzgomW0SsIblvdCTwPnAh8DtgXkQ0R0RjRHwcEfOB84HDJdWlzZuBuqzT1QHN6eXB3LpM/dJW2ubWd4nt+9YAMG/Riq78WjOzkpDPPad9gDnAa5LuJZmZl3vNKSKiYEsYRcQrJKMlACT9Ffh1S4dmDkl/TiWZDPFc+nlsWpapG55zD2kscGdO28x3Dgd6AjM61Jk8jRyUDNRmzl/KPkP7deVXm5kVXT7J6bKs96e0ckxQwPX1JO1BkhQqSJZM2h64TdK+JFt0zCS5F/Zz4PGIyFyOux24QNKDaUwXAr8AiIgZkl4CLpV0MckmiXuQTIiA5PLg05IOIJmtdzkwsSsnQwDs2K8XNdUVzFzQ3JVfa2ZWEvJJTjt3WhStOxU4C6gG/gwcFhGr0tHMVcAgYAnJVPETs9rdCAwnuU8G8Mu0LGMcyRqBi4C5wAkR0QQQEVMlnUuSpAaQrLR+Rmd0ri0VFWLEoFpmzO/SnGhmVhKU3IaxjmpoaIjGxsaCnvOCu1/i6dkLefrfDinoec3MSoWkyRHRkFuez4QI62IjBtfy3uKVLFnpGXtmVl62+LJeuofT5kREfK0D8ViWUemkiDcWNLP3Tp4UYWblI597TqdvwTFBsiOuFcDIwbVAMmPPycnMyskWX9aLiIrcF8lEhU8CNwPPkMycswLZsV/vZMbefM/YM7Py0qF7ThGxNiJmRsQ5JAvBXl2YsAygskLsUl/r6eRmVnYKOSHi/9jwrJAVyKjBfZjp6eRmVmYKmZwGkKxLZwU0YlAt7y5eyVLP2DOzMtLh5CRpW0knAP8MTO54SJZt1OANM/bMzMpFPlPJ17FhDbtNqoEPgQsKEZRtMHJQOmNvQTN7ecaemZWJfKaS386mySlIktIM4K6uXn+uHHyif296VlX4vpOZlZVWk5OkzwPTstacO72rgrINPGPPzMpRW/ecJgGHZT5Imi3p7zo/JMs1cnCtn3Uys7LSVnJaRbKPUcYwPBuvKEYN7sM7H62gedXHxQ7FzKxLtHXPaQZwmqQXSLaWABggaae2ThgRcwsVnCVGpJMi3ljQzJ6f2La4wZiZdYG2ktOVJLvDvpB+DuCn6astlR2OyjaSmU4+c/5SJyczKwutJqeIuEfSy8BBJDvQXgrcB7zSJZHZejv1702PKu+Ka2blo82p5BExk2QrdCRdBkyIiDu7IC7Lsn7GnqeTm1mZ2OLnnNJVyK1IRg6qZfKcRZs/0MysG3DC2UqMGlzLOx+tYJln7JlZGXBy2kqMSHfFndXk+05m1v05OW0lRqW74s7ww7hmVgacnLYSO/XvTY/KCmYu8KQIM+v+nJy2ElWVFQyv38bLGJlZWSjp5CRpN0mPSVos6Q1Jx2bVHSJpuqTlkiZJGppVJ0lXS1qYvn4kSVn1w9I2y9NzHJrzvSdJmiNpmaT7JPXvmh63beTgPh45mVlZ2OLkJOkzks7OKfuypCmS3pF0VSEDk1QF/B74X6A/MB74jaRRkgYCE4FL0rpG4O6s5uOBY4CxwB7AUcA5WfV3AS+S7N57EXCPpPr0e0cDNwKnAoOB5cB1hexbe40aVMvbH65g+WrP2DOz7i2fkdOlwPpVydM19u4CtgMWA9+WdEYBY9sVGAL8Z0SsjYjHgKdIksZxwNSI+J+IWAlcBoyVtGva9jTgmoiYFxHvANcAp6dxjwL2Bi6NiBURMQGYAhyftj0ZuD8inoyIZpIEeJykPgXsW7uMTCdFzFqwrMiRmJl1rnyS01iS5JAxjmQH3D0jYnfgjyQjlkJRK2VjgNHAy5nCiFgGzErLya1P32fXzc7ZGDG3Pvvcs4DVwKhNgpHGS2qU1NjU1LTlPWunkekaezO8UoSZdXP5JKcBwPtZn48AnkxHJgB/AEYWKjBgOrAA+JakakmHAwcCvUm27licc/xiIDO6ya1fDNSm953ybZtbv15E3BQRDRHRUF9fn0/f2mXo+hl7nhRhZt1bPsnpI5J7MEjqCewHPJlVH0CvQgUWEWtI7hsdSZIULwR+B8wDmoG6nCZ1QGZIkVtfBzRHRLSjbW590WyYsVf0UMzMOlU+yekl4CxJ+5Dch6kBHs6q3xmYX7jQICJeiYgDI2JARBwBDAeeA6aSXGYEQNI2wC5pObn16fvsuuE595By67PPPZxk08UZhepXR4wY5C3bzaz7yyc5XUGydcZzwL8Dj0ZEY1b9UcCzBYwNSXtIqpHUW9K/pN9/G3AvMEbS8ZJqgO8Cr0TE9LTp7cAFknaQNIRk1HUbQETMIEm0l6bnPpZkRt+EtO0dwNGSDkiT3uXAxJx7VEUzclAf3l60nBWr1xY7FDOzTrPFySki/koyy+2bJDPfjs7USRpAMiHi+sKGx6nAeyT3ng4BDouIVRHRRDK77vsku/TuSzJBI+NG4H6SWXivAg+kZRnjgIa07Q+BE9JzEhFTgXNJktQCkntN5xW4X+02anAtEV5jz8y6NyW3YayjGhoaorGxcfMHdtDspmb+9ponuPKYMZyy39DNNzAzK2GSJkdEQ255Pg/hVkrqnVO2raQLJX1f0phCBGpt23ngNgwb0Js/vlbQ23tmZiUln3tON5LcbwJAUjXwF+DHwL8Bz0vas6DR2SYkccTo7Xh61gcsXrGm2OGYmXWKfJLT50ieZco4Adgd+AbwWZKZet8pXGjWmiPGbMeatcGk6QuKHYqZWafIJzltD7yZ9flIkiWEro+IZ4CbgL8pZHDWsj133JZBfXry8NT3N3+wmdlWKJ/kJKAy6/NBwKSsz+8BgwoQk21GRYU4fPRgHn+9iZVrPKXczLqffJLTmyRLFiFpf5KRVHZyGsKmy/5YJ/nC6O1ZsWYtT87o/DX9zMy6Wj7J6Vbgy5JeJdnGYgEbrxCxL8l6eNYF9h3en769qnnIl/bMrBvK5yHcn5Jsm7GKZC+kYyNiOax/CHc/4MFOiNFaUF1ZwSG7DeJP0xawZu26YodjZlZQee2EGxFXRMQ+EfG36SSITPnCiBgUET8qfIjWmiNGb8fiFWt47s0Pix2KmVlBtXubdkkD0x1prUg+P7KeXtWVPPSqL+2ZWfeSV3KSNETSryV9RPJc03xJiyTdJmmHTonQWtWrRyUHjqrnj6+9z7p1XobKzLqPfJYv2gloJFmMdTZwZ/qaDXwVeE7SJzojSGvdF8Zsx/wlq3hp3kfFDsXMrGDy3TKjH3BUROwdEaemr31IHsjtnx5jXejgXQdRVSE/kGtm3Uo+yelw4LqI2GRGXkT8H8l2GV8oVGC2Zfr2quazIwby8Kvv4xXmzay7yCc59QNmtlE/E9i2Q9FYuxwxejBvLVzOjPne48nMuod8ktM8kiWLWvP59BjrYoftPhgJX9ozs24jn+T0P8DfS/qBpL6ZQkl1kq4CvgLcXegAbfMG9alhn536eUq5mXUb+U6IeBr4NvCBpDmS5gALSbbK+CtwZeFDtC1xxOjteO29Jbz94fJih2Jm1mH5LF+0HDgQOAf4I7AMWE6yvt544OCIWNEZQdrmHTF6OwAenPJekSMxM+u4qnwOjoi1wM3py0rITgN68+lh/bj96Tmcsf/O9Khq9+IfZmZF53/BupFvHDyCdz5awX0vvlPsUMzMOqTVkZOkr7bnhBFxe/vDsY44cFQ9n9qhL9c9/gbH7b0DVZX+28PMtk5tXda7DQiSHXC3VAAFS06ShgHXkWz/vgq4B/gmsCPJ5ofLsg6/OiKuSNsJ+CFwVlr3K+DbkT6lmp73VpI9qOYC50fEo1nfexLwA2Ag8AhwZkSU/NLfkvjGwSM49zeTeWDKe3x5Ty93aGZbp7aS08FdFkXrriPZ1HB7kgd8HwHOA/6Q1m8bER+30G48cAwwliRhPkKyBuANaf1dJDMPv5S+7pE0MiKaJI0GbiRZkukF4KY0jnEF7lunOHz3wYwaXMt/TXqDo/cYQkVFPn9bmJmVhlaTU0Q80ZWBtGJn4NqIWAm8L+khYDQbklNrTgOuiYh5AJKuAc4GbpA0CtgbODydXThB0jeB40mS18nA/RHxZNr2EmCapD4RsbTgPSywiopk9PRPv32JP742ny+M2a7YIZmZ5a3Ub0r8DBgnqXe6JccXgYey6udImifp1py9pUYDL2d9fjkty9TNzkk0ufXr20bELGA1MCo3OEnjJTVKamxqampfDzvBUXsMYdiA3lw7aabX2zOzrVKpJ6cnSJLFEpKlkRqB+4APgE8DQ4F9gD7AHVntaoHFWZ8XA7Xpvajcukx9n1ba5tavFxE3RURDRDTU19fn27dOU1khzjtoBK++s4THZ5RO0jQz21Ilm5wkVZA84DsR2IZkckI/kokPzRHRGBEfR8R84HzgcEl1afNmoC7rdHVAczohIrcuU7+0lba59VuFY/bagR227cW1j73h0ZOZbXVKNjmR7A/1CZJ7TqsiYiHJDLsvtXBs5l/fzN3/qSSTITLGpmWZuuGS+rRRv76tpOFAT2BG+7vS9XpUVXDOgcOZPGcRz8wu+YmGZmYbKdnkFBEfkEwX/7qkKknbkkx0eFnSvpI+KalC0gDg58DjEZG5HHc7cIGkHSQNAS4kmRpPRMwAXgIulVQj6VhgD2BC2vYO4GhJB0jaBrgcmLg1TIbI9ZWGT1DfpyfXTmprpxMzs9JTsskpdRzJBoZNwBvAx8A/A8NJJkYsBV4leQbqxKx2NwL3A1PS+gfSsoxxQAOwiOR5qBMiogkgIqYC55IkqQUk95rO65TedbKa6krGHzCcp95YyAtzFxU7HDOzLSbfjyiMhoaGaGxsLHYYm1i26mM+d/VjjNmhL7ef+RmSOSFmZqVB0uSIaMgtL/WRk3XQNj2r+Me/HcmfZ37A/77iFcvNbOvg5FQGTvvsMPbYsS/fu38qi5evKXY4Zmab5eRUBiorxFXHfopFy9fwg/+bVuxwzMw2y8mpTIzZoS9nfW5nfvv82zw7e2GxwzEza5OTUxn5p0NHsmO/XvzbvVNY9fHaYodjZtYqJ6cy0rtHFVceM4bZTcu4btKsYodjZtYqJ6cyc9AnB/HlPYdw/eOzeGPBVvdcsZmVCSenMnTJUbvTq0cl/z7xVdat83NuZlZ6nJzK0MDanlz0pd147q0P+V3j28UOx8xsE05OZervG3Zkv+H9+f6D05i7cHmxwzEz24iTU5mSxI9PGIuA8+6czMo1nr1nZqXDyamMfaJ/b/7jK3vy6jtL+N79rxU7HDOz9Zycytyhuw/m6wftwl3PzWXC5HnFDsfMDHByMuDCw0ax3/D+XHTfFKa/v6TY4ZiZOTkZVFVW8PMT96JPTTVf/80LLF3pxWHNrLicnAyAQX1quPbEvZj74XK+M2EK3ufLzIrJycnW23f4AP71iE/ywJT3uPWpt4odjpmVMScn28j4zw/nsN0Hc9WD0/jfV94tdjhmVqacnGwjkrjmK2PZa6dtOf/OF7npyVm+xGdmXc7JyTZRV1PNf39tX4781PZc9eB0LvvDVNZ6DT4z60JVxQ7ASlNNdSW/OHEvhmxbw81/fpN3F6/k5+P2olePymKHZmZlwCMna1VFhbjoyN257OjdeXTafMbd/AwfNK8qdlhmVgacnGyzTt9/Z244ZR+mv7eE4677K8+/9WGxQzKzbq6kk5OkYZIelLRI0vuSrpVUldYdImm6pOWSJkkamtVOkq6WtDB9/UiScs47KW07XdKhOd97kqQ5kpZJuk9S/67rdWk6YvR23DV+P1asWcvf3/A0J938DM/MXljssMysmyrp5ARcBywAtgf2BA4EzpM0EJgIXAL0BxqBu7PajQeOAcYCewBHAedk1d8FvAgMAC4C7pFUDyBpNHAjcCowGFiexlH29t6pH09+62AuPnI3Zi5oZtxNz/CVG5/mqTc+8Iw+MysolfI/KpKmARdGxIPp5x8DdcBk4PSI+Gxavg3wAbBXREyX9Ffgtoi4Ka3/GnB2ROwnaRQwBRgYEUvT+j8Dd0TEDZKuAoZFxElp3S7ANGBA5viWNDQ0RGNjY2f8GkrSyjVr+e1zc7n+iVnMX7KKfYb248z9d+aQ3QZRU+1JE2a2ZSRNjoiG3PJSHzn9DBgnqbekHYAvAg8Bo4GXMwdFxDJgVlpObn36Prtudk6iya3PPvcsYDUwKjc4SeMlNUpqbGpqancnt0Y11ZWcvv/OPPGtg7nimDG899EKvnHnCzRc+SgX3P0Sk15fwJq164odppltpUp9KvkTwNnAEqAS+DVwH8llutxssBjok76vTT9n19Wm951y6zL1O7TSNvfc66Ujs5sgGTltYZ+6lZrqSk7dbygnfWYnnpm9kD+89C7/9+p7THzxHfpv04MvjtmOo/YYwqeH9aOqstT/FjKzUlGyyUlSBfAwyf2fz5IkjVuAq4Fmkst72eqAzGgot74OaI6IkJRv29x6a0Flhdh/xED2HzGQy48ZzROvN/GHl99lwgvzuOPZufTrXc0huw3m8N0Hc8DIej8vZWZtKtnkRDLR4RPAtRGxClgl6VbgSuDnwGmZA9N7TrsAU9OiqSSTIZ5LP4/NqRsuqU/Wpb2xwJ05bTPnHg70BGYUtHfdWM+qSg4fvR2Hj96OZas+5okZTfxx6vs8PPV97pk8j5rqCj4/sp7Ddh/MgZ+sZ1CfmmKHbGYlptQnRMwmuWz2E5KR060ks+e+CbwBnAk8AHwPODAi9kvbnQv8E3AoEMAjwC8i4oa0/hngL8DFJPexbgVGRkRTOlvvaeBI4AWSkVtVRIxrK9ZymxDRHmvWruPZ2R/yx9fe55HX5vPe4pUAfGqHvhz8yXoO2nUQY3fclsoKbeZMZtZdtDYhotST057AT0lGMmuBScA3ImJB+mzStcBQ4FmS2Xtvpe1EcvnvrPRUvwS+HWlnJQ0DbgP2Beam53w063tPAn5IMtX8UeCMiGjzyVMnp/xEBFPfXcITM5qYNH0BL8xdxLqAfr2r2X/EQIYN2IbBdT0ZVFfDdnU1DK6rYWBtD9+3MutmtsrktDVxcuqYj5av5s8zP2DS6wt4dvaHvL9k5SaLzUpQXVlBpURlhahQsotvhZL3FRISiGR1dYn0c/IeNtRl3m84eYtvc76/80d0HjPa1uiq4z7Fp4e1b62C1pJTKd9zsjKybe8eHD12CEePHQLA2nXBwmWrmL94FfOXrGT+0pUsWLKKVR+vY10Ea9dteH28LogIImBdBAFEJKOz5H2S5DLlmfcZ2X+gtfqnWhf8DRdd8SVmnaBXJzzb6ORkJamyQgzqU8OgPjV8ir7FDsfMupgv4JuZWclxcjIzs5Lj5GRmZiXHycnMzEqOk5OZmZUcJyczMys5Tk5mZlZynJzMzKzkePmiApHUBMxpZ/OBJDv5lhv3u/yUa9/d79YNjYj63EInpxIgqbGltaW6O/e7/JRr393v/PmynpmZlRwnJzMzKzlOTqXhpmIHUCTud/kp176733nyPSczMys5HjmZmVnJcXIyM7OS4+RkZmYlx8mpiCT1l3SvpGWS5kg6qdgxdQZJ50tqlLRK0m05dYdImi5puaRJkoYWKcyCk9RT0q/S/7ZLJb0o6YtZ9d2577+R9J6kJZJmSDorq67b9jtD0khJKyX9Jqus2/Zb0uNpf5vT1+tZde3qt5NTcf0XsBoYDJwMXC9pdHFD6hTvAlcCt2QXShoITAQuAfoDjcDdXR5d56kC3gYOBPqS9PN3koaVQd9/AAyLiDrg74ArJe1TBv3O+C/g+cyHMun3+RFRm74+CR3rt2frFYmkbYBFwJiImJGW/TfwTkR8p6jBdRJJVwI7RsTp6efxwOkR8dn08zYkS53sFRHTixZoJ5L0CvA9YABl0ndJnwQeB/4J2JZu3m9J44DjgNeAERFxSnf/37qkx4HfRMQvc8rb3W+PnIpnFLA2k5hSLwPdceTUmtEkfQYgIpYBs+imvwNJg0n+u0+lDPou6TpJy4HpwHvAg3TzfkuqAy4HLsyp6tb9Tv1A0geSnpJ0UFrW7n47ORVPLbA4p2wx0KcIsRRL2fwOJFUDdwC/Tv9i7PZ9j4jzSPpzAMmlnVV0/35fAfwqIt7OKe/u/f42MBzYgeTB2/sl7UIH+u3kVDzNQF1OWR2wtAixFEtZ/A4kVQD/TXJ/8fy0uCz6HhFrI+IvwI7A1+nG/Za0J3Ao8J8tVHfbfgNExLMRsTQiVkXEr4GngC/RgX47ORXPDKBK0sissrEkl3zKxVSSPgPrr0fvQjf6HUgS8CuSSS/HR8SatKrb9z1HFRv61137fRAwDJgr6X3gX4DjJb1A9+53SwIQHei3k1ORpNdeJwKXS9pG0v7Al0n+wu5WJFVJqgEqgUpJNZKqgHuBMZKOT+u/C7zSHW4QZ7ke2A04OiJWZJV3275LGiRpnKRaSZWSjgBOBB6jG/eb5HLWLsCe6esG4AHgCLpxvyVtK+mIzP+vJZ0MfB54mI70OyL8KtKLZGrlfcAyYC5wUrFj6qR+Xkbyl1T267K07lCSG+YrSGZ0DSt2vAXs99C0rytJLm9kXid3574D9cATwEfAEmAKcHZWfbfsdwu/h8tIZrB1636n/72fJ7lU9xHwDHBYR/vtqeRmZlZyfFnPzMxKjpOTmZmVHCcnMzMrOU5OZmZWcpyczMys5Dg5mZlZyXFyMrNNSHorXWnarCicnMzMrOQ4OZmZWclxcjIzs5Lj5GTWRST1lPTvkqZKWinpI0n3S9or57iDJIWk0yX9o6QZ6fEzJP1jK+f+vKRHJC2WtELSC5K+1sqxIyTdKmmepNWS3pX0e0n7tHDsrpIekLQ0Pfc9krbLOaa/pP+UNCuNc6GkyZK+1ZHfl5U3r61n1gXSzQb/CHyWZOX5yUBf4Gxge+DzEdGYHnsQMAl4AdgOuJFkUc0TgU+TLJr7vaxzH02y+vP7JCtjLwXGAZ8BroqIi7KObQD+BFSTbOXxKskCxAcCD0TEL9Lj3gLWkGwKdy/JbqZjgXOARyPi8Kxz/olkFeob0+N6A7sCO0XEkR37zVm5cnIy6wKS/hn4D+ALEfFwVnkdSYKYHREHpWUHkSSnZmC3iJiXlvcA/gLsBewcEfMkVQKzSRLd7hHxbtaxk4D9gF0jYma6t9QUYATwmYh4JSfGiohYl75/i2RV9X+IiN9lHfNfwHlpXNMl9SVZifr6SHa+NSsIX9Yz6xqnkGwbMFnSwMwL6AE8AnxOUq+cNndkEhNARKwm2WW1Cjg6Ld4H2Am4JZOYso79Mcn/x7+cFu8JjAZuzU1MaZt1OUXvZiem1GPpzxHpzxUk26/vK2lY6903y09VsQMwKxO7Ab2ApjaOGQi8nfV5WgvHvJb+HJ7+3Dn92dLOoq/mHJvZdfnFNiPdYHYLZQvTnwMgSYKSvgn8DHhT0mskCey+iPjTFn6P2SacnMy6RuaS2gVtHJObuFq65q7NfN5cDK2dtyVrt+BcRMQNkn4PHEly7+oE4HxJd0fEuDziM1vPycmsa8wk2TH0sRYun7Vm9xbKdkt/ZkY1s9Kfo9tonzn29fTnXi0c2yER8R7wS+CX6X2w/wZOlHRNRDxf6O+z7s/3nMy6xu0kM+9aHDlJGtxC8cmSdsw6pgfwzyQjmv9Ni18A5gJnZE/xTmcHfotklPT7tPhlkst/Z0raJJmlEybyIqm3pN7ZZRGxFsjc0+qf7znNwCMns67yM+Aw4MeS/pbkvswSkskMhwArgYNz2swAnpV0A8n08JNIppJfERFvQ5IIJJ1PMt37eUmZqeT/QDJT76qImJkeG5LOIJlK/pykzFTybUkuxz0E/CLPfo0CnpB0b3quRSSju68DbwJ/zvN8ZoCTk1mXiIg1ko4kmYZ9KpB5Tuld4Dng1y00+wVQB/wjSRKbC3wzIn6Wc+77JR0CXEwyWupBMpni7Ij4Zc6xz0v6NHAJ8BXgXOCDNIan2tG1t4FbSBLrMUBP4B3gZuDqiFjejnOa+Tkns1KT9ZzTGRFxW1GDMSsS33MyM7OS4+RkZmYlx8nJzMxKju85mZlZyfHIyczMSo6Tk5mZlRwnJzMzKzlOTmZmVnKcnMzMrOT8f+B2nBmgGlF8AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal parameter is:\n",
      "[ 0.70799865  0.11092323  0.          0.          0.          0.\n",
      " -0.01320565  0.        ]\n"
     ]
    }
   ],
   "source": [
    "#GD implementation with lasso \n",
    "n_epochs = 50\n",
    "step_size = 1e-5\n",
    "lasso_parameter = 3000\n",
    "D = A_train.shape[1] # dimension of pbm \n",
    "\n",
    "x = np.random.rand(D) # random initialization\n",
    "\n",
    "loss_function_array_lasso = []\n",
    "epochs = []\n",
    "\n",
    "for epoch in range(n_epochs):\n",
    "    #here only one iteration per epoch\n",
    "    epochs.append(epoch)\n",
    "    \n",
    "    loss_function_train = np.linalg.norm(np.dot(A_train,x)-Y_train)**2 + lasso_parameter*np.linalg.norm(x, ord=1)\n",
    "    loss_function_array_lasso.append(loss_function_train)\n",
    "    gradients = 2*A_train.T.dot(np.dot(A_train,x)-Y_train) #gradient of the smooth part \n",
    "\n",
    "    #ista \n",
    "    threshold = step_size*lasso_parameter\n",
    "    for i in range(D): #component wise\n",
    "            vali = x[i]-step_size*gradients[i]\n",
    "            if vali < -threshold:\n",
    "                x[i] = vali+threshold\n",
    "            elif vali > threshold:\n",
    "                x[i] = vali-threshold\n",
    "            else:\n",
    "                x[i] = 0\n",
    "\n",
    "plt.plot(epochs, loss_function_array_lasso)   \n",
    "plt.xlabel(\"epochs\", rotation=0, fontsize=18)\n",
    "plt.ylabel(\"loss function\", rotation=90, fontsize=18)\n",
    "plt.show()\n",
    "print(\"Optimal parameter is:\")\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we provide a plot where we display error on test set for various values of $\\lambda$. As we can see the error is severy reduced when we apply some regularization but then the error achieved plateau and remains high (on average we are 36% off from the target)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lasso: 36.54963976828266%\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$E$')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEVCAYAAAAPRfkLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAfGklEQVR4nO3df5RcZZ3n8fe3qqu6u7qT7g7dCYMBIvJLwQOszUHFRUYYFRR1jUdUdMWdc+LKsivjj6POgpN1UGdUVtZdZciIi4PK4A90BMQfqFFBdmaCoEMcfhNCZEg6kF/d6fSPqu/+caurK01VuqrSXU+lns/rnHtu9b23bn2fVPp++3me+zzX3B0RERGAVOgARESkdSgpiIhIiZKCiIiUKCmIiEiJkoKIiJQoKYiISElH6AAOxuDgoK9atSp0GCIih5R77rlnu7sPVdp3SCeFVatWsWHDhtBhiIgcUszsiWr71HwkIiIlSgoiIlKipCAiIiVKCiIiUqKkICIiJUoKIiJSEnVSKBSc6XwhdBgiIi0j6qTwN798lD/5/C9DhyEi0jKiTgrrHxzh8e1j7JvKhw6lddxxR7KISJQO6RHNB2M6X+BftuwCYPvoBCsHcoEjahFXXpmszz03bBwiEkS0NYWHt40yXqwhjOyZCByNiEhriDYp3PfkztJrJQURkUS8SWHzTjpSBsD20cnA0YiItIZok8Jvt+zkjGOWAaopiIjMiLKjeWximoe27uHSk47j90/tZmR0X+iQWse114aOQEQCirKmMDFd4D+d+XxeefwQQ0s6VVMod8IJySIiUYqyprCsJ8vlr38RAIO9nepTKHfLLcn6ggvCxiEiQUSZFMoNLenk3s07Q4fROq66KlkrKYhEKcrmo3JDvUnzkbuHDkVEJDglhSWdjE/lGZvUVBciItEnhcHeTgC2q7NZRERJYWhJkhRGRpUURETU0TyTFFRTSNxwQ+gIRCSg6JNCqflINYXEkUeGjkBEAoq++WhZT5aUqaZQctNNySIiUYq+ppBOGYf1alRzyTXXJOsLLwwbh4gEEX1NAWbHKoiIxK5pScHMOs3sOjN7wsz2mNm9ZnZecd8qM3MzGy1brmhWbINLOtWnICJCc5uPOoAngVcCm4HzgW+a2YvLjul39+kmxgQkNYVHtu5p9seKiLScptUU3H3M3de6+yZ3L7j7rcDjwEuaFUM1Q0uSSfE01YWIxC5YR7OZrQCOBzaWbX7CzBz4CfBhd9/ejFgGe7NM5gvsHp+mL5dpxke2rm9/O3QEIhJQkI5mM8sAXwe+6u4PANuB04GjSWoOS4r7K713jZltMLMNIyMjCxLP7KhmPWyHwcFkEZEoNT0pmFkKuAGYBC4FcPdRd9/g7tPuvrW4/dVmtnTu+919nbsPu/vw0NDQgsQ0kxS26Q4kuP76ZBGRKDW1+cjMDLgOWAGc7+5TVQ6dady3ZsQ11KupLkpmEsLFF4eMQkQCaXafwjXAC4Fz3X18ZqOZnQHsBB4GBoAvAOvdfVczgpqpKegJbCISu2aOUzgaeC9wKvB02XiEi4BjgB8Ce4D7gQng7c2Kra87QyZtqimISPSaVlNw9yc4cHPQjc2KZS4z06hmERE0zUXJ4JJOPVNBRKIX/YR4M4Z6O/m3XbollR/8IHQEIhKQagpFQ6opJHK5ZBGRKCkpFA32dvLM6AT5QuRTXXzpS8kiIlFSUiha1pOl4LB7vNrQiUh885vJIiJRUlIoWtaTBeDZvRqrICLxUlIoGigmhR1jSgoiEi8lhaJluWJNQUlBRCKmpFC0rLdYU1DzkYhETOMUimZrCpF3NK9fHzoCEQlINYWi7myarkxKNQURiZqSQplluaz6FD73uWQRkSgpKZQZ6Mnq7qNbb00WEYmSkkKZZT1ZjVMQkagpKZQZyKmmICJxU1Ios6xHfQoiEjfdklpmIJdl975ppvIFMulI82V3d+gIRCQgJYUyy3oyQDKAbfmSrsDRBHL77aEjEJGAIv1zuLLZ+Y8iH8AmItFSUiij+Y+Av/zLZBGRKCkplCnVFGK+LfWnP00WEYmSkkKZ0jMVYq4piEjUlBTK9OeKHc1KCiISKSWFMp0daXo7O9ixVx3NIhIn3ZI6R193hl0xP6f5sMNCRyAiASkpzLG0O8Ou8Yibj77zndARiEhATWs+MrNOM7vOzJ4wsz1mdq+ZnVe2/xwze8DM9prZz83s6GbFVq4/9pqCiEStmX0KHcCTwCuBPuAK4JtmtsrMBoGbi9uWARuAm5oYW0lfd4adMfcpfOxjySIiUWpa85G7jwFryzbdamaPAy8BDgM2uvu3AMxsLbDdzE509weaFSOoT4G77w4dgYgEFOzuIzNbARwPbAROAn47s6+YQB4tbp/7vjVmtsHMNoyMjCx4XP25yJOCiEQtSFIwswzwdeCrxZpAL7BrzmG7gCVz3+vu69x92N2Hh4aGFjy2pd0ZJqYL7JvKL/i5RURaXdOTgpmlgBuASeDS4uZRYOmcQ5cCe5oYGpA0HwGqLYhIlJqaFMzMgOuAFcBqd5+58m4ETik7rgd4QXF7U82Mao42KaxcmSwiEqVmj1O4BnghcK67j5dt/y7wWTNbDdwGfBz4XbM7mWG2phDtHUhf+1roCEQkoGaOUzgaeC9wKvC0mY0Wl4vcfQRYDXwS2AGcAbytWbGVU/ORiMSsmbekPgHYAfbfAZzYrHiq6e9OZkqNNilcdlmyvvrqkFGISCCa5mKO2eajSKe6uO++0BGISECaJXWOJV0dmMHuWGsKIhI1JYU5UiljaZcGsIlInJQUKujrzrBTSUFEIqQ+hQqinuri+ONDRyAiASkpVBD1pHjr1oWOQEQCUvNRBUu7M+yKdfCaiERNSaGCqB+0s2ZNsohIlNR8VMFM85G7k0zXFJGHHgodgYgEpJpCBX3dGaYLztikps8WkbgoKVQwkEumuoh2VLOIREtJoYK+XOQzpYpItNSnUMFMTWFHjDWFU08NHYGIBKSkUMFAsaawI8aagmZHFYmamo8q6FefgohESkmhgv6Y+xTe+c5kEZEoqfmogkw6RW9nR5x9Clu2hI5ARAJSTaGK/lwmzpqCiERNSaGKgVw2zpqCiERNSaEK1RREJEbqU6hiIJflyWf3hg6j+V72stARiEhASgpV9OcycY5T+PSnQ0cgIgGp+aiK/lyW3fumyBc8dCgiIk2jpFDFQC6DO/E9V2H16mQRkSip+aiK8plSl/VkA0fTRM88EzoCEQmoqTUFM7vUzDaY2YSZXV+2fZWZuZmNli1XNDO2ufpinv9IRKLV7JrCU8CVwGuA7gr7+919urkhVaZnKohIjJqaFNz9ZgAzGwZWNvOz6xX1TKkiEq1W61N4wswc+AnwYXffHiqQaGdKPeec0BGISECtkhS2A6cD9wGHAV8Evk7SzLQfM1sDrAE46qijFi2gpV0dpFMW36jmK4J25YhIYAfsaDazW8ysd7GDcPdRd9/g7tPuvhW4FHi1mS2tcOw6dx929+GhoaFFi8nM6OvOaP4jEYnKfHcfnQ/kZn4ws5vM7LCyn1OVLtwLYGbEmC3CuWsW5fxH552XLCISpfmaj+ZelM8H+oCZm9mHgD/UcJ7kZGYdxWPTQNrMuoBp4CXATuBhYAD4ArDe3XfVct7FEuVMqePjoSMQkYAWYpxCPee4HBgHPgq8s/j6cuAY4IfAHuB+YAJ4+wLEdlAGYp3/SESitRAdzTVPDuTua4G1VXbfuACxLKj+XJbfP7U7dBgiIk1Ty1/57zGzlxabeqCOJHCo6+9WTUFE4jJfTWE98BHg0yRt/x3AX5vZXcBvgG2LGl1gAz1Zxqfy7JvK05VJhw6nOV7/+tARiEhAB0wK7v4qADM7hqQzeGa5Alg2c9hiBhhSf3FU8869UxzeF0lS+NCHQkcgIgHV1Kfg7o8BjwHfmtlmZquAYeDfLUpkLaA0/9H4JIf3dc1ztIjIoa/hjmZ33wRsAr69UMG0mpmawo6xiPoVzj47Wa9fHzIKEQlED9k5gP7uSOc/EpFoKSkcwECPZkoVkbgoKRzATJ9CdKOaRSRaSgoH0JVJ05VJxfecZhGJVqtMnd2yBnJZdoxFVFN461tDRyAiASkpzKMvtlHNl1wSOgIRCUjNR/MYyGXjuvto795kEZEoKSnMY6AnsgftnH9+sohIlJQU5tGfy6qjWUSioaQwj/7u5Olr7m07xZOISImSwjwGclmmC86eienQoYiILDolhXmUZkqNaf4jEYmWbkmdR/lMqUeRCxxNE1x8cegIRCQgJYV5RDf/kZKCSNTUfDSPvthmSt2+PVlEJEqqKcxjoPRMhUiSwlvekqz1PAWRKKmmMI++7siaj0QkakoK8+hIp1ja1aEBbCISBSWFGgz0ZOOa6kJEoqWkUIP+2GZKFZFoqaO5Bv25iGoK73tf6AhEJKCm1hTM7FIz22BmE2Z2/Zx955jZA2a218x+bmZHNzO2AxnIRTRT6oUXJouIRKnZzUdPAVcCXynfaGaDwM3AFcAyYANwU5Njq6o/l2VnLM1HTz6ZLCISpaY2H7n7zQBmNgysLNv1ZmCju3+ruH8tsN3MTnT3B5oZYyUDuSx79k0znS/QkW7zbph3vStZa5yCSJRa5Qp3EvDbmR/cfQx4tLh9P2a2ptgEtWFkZKQpwZUmxdNtqSLS5lolKfQCu+Zs2wUsmXugu69z92F3Hx4aGmpKcKWkEEu/gohEq1WSwiiwdM62pcCeALE8x8xMqbotVUTaXaskhY3AKTM/mFkP8ILi9uBK02crKYhIm2tqR7OZdRQ/Mw2kzawLmAa+C3zWzFYDtwEfB37XCp3MMNt8FMVtqR/8YOgIRCSgZtcULgfGgY8C7yy+vtzdR4DVwCeBHcAZwNuaHFtVUfUpXHBBsohIlJp9S+paYG2VfXcAJzYznlr1dnbQkbI4+hQefDBZn3BC2DhEJAhNc1EDM4tnANt735usNU5BJEqt0tHc8gZymTiaj0QkakoKNeqPaf4jEYmWkkKNomk+EpGoKSnUKKqZUkUkWuportFALDWFyy8PHYGIBKSkUKP+XJaJ6QLjk3m6s+nQ4Syec88NHYGIBKTmoxpFM6r5vvuSRUSipJpCjQbKksIR/d2Bo1lEl12WrDVOQSRKqinUqF+T4olIBJQUaqSZUkUkBkoKNRqIpU9BRKKmpFCjvphmShWRaKmjuUadHWly2XT7z5T6qU+FjkBEAlJSqMNALtv+zUcvf3noCEQkIDUf1aE/l2FXu9cUfv3rZBGRKKmmUIcoagp//ufJWuMURKKkmkId+nIZ3ZIqIm1NSaEOmilVRNqdkkIdBnJZdo1PUSh46FBERBaFkkId+nNZCg579k2HDkVEZFGoo7kO/d2zo5pnBrO1nauvDh2BiASkpFCHgZ7ZpLCKnsDRLJJTTw0dgYgEpOajOkQxU+oddySLiERJNYU6lGZKHW/jO5CuvDJZ6wlsIlFqqZqCma03s31mNlpcHgwdU7nSTKljbVxTEJGotVRSKLrU3XuLywmhgym3pCuDmWZKFZH21YpJoWWlU0Zfd6b9Z0oVkWi1YlL4tJltN7O7zOzs0MHMFcX8RyISrVbraP4I8HtgEngbcIuZneruj84cYGZrgDUARx11VNMD7M9l2DXexjWFa68NHYGIBNRSNQV3/0d33+PuE+7+VeAu4Pw5x6xz92F3Hx4aGmp6jG1fUzjhhGQRkSi1VFKowAELHUS5/u5Me999dMstySIiUWqZ5iMz6wfOAH4BTAMXAmcBl4WL6rn6c9n2vvvoqquS9QUXhI1DRIJomaQAZIArgROBPPAA8CZ3b7mxCmOTeSanC2Q7Wr2iJSJSn5ZJCu4+ApweOo759PfMjmpevqQrcDQiIgtLf+rWaWZUc1vPfyQi0VJSqFN/d1JT2DHWxv0KIhKtlmk+OlT0z8x/1K41hRtuCB2BiASkpFCngWKfwq52nSn1yCNDRyAiAan5qE4D7V5TuOmmZBGRKKmmUKfuTJpsOtW+o5qvuSZZX3hh2DhEJAjVFOpkZvTnMuxs51HNIhItJYUGtP38RyISLSWFBjxvoJvNz+4NHYaIyIJTUmjAcct7eWz7GNP5QuhQREQWlDqaG3Ds8l4mpws8uWOc5w/2hA5nYX3726EjEJGAVFNowHErlgDw8NY9gSNZBIODySIiUVJSaMCxy3sBeHjbaOBIFsH11yeLiERJSaEBvZ0dHNHXxSNKCiLSZpQUGnTsiiU8vK0Nm49EJGpKCg06bnkvj2wbpVDw0KGIiCwY3X3UoOOW97JvqsA1v3iUFwz1sHxpF8uXdDK0pJPOjnTo8EREGqKk0KDhVQPksmk++6PnPi10IJdhxdIuhpZ0smJpFyuWdnLc8iW89uTD6cooYYhI6zL3Q7f5Y3h42Dds2BDs8/MF59mxSbbu3sfIngm27t7HtrL1tpn1ngnyBWewt5N3v+xo3vnSo0tTcLecvcWR2rlc2DhEZNGY2T3uPlxpn2oKByGdMoaKTUYHUig4v370Gf72V49x1U8e4ovrH+E/nLaS95y5iuOLYx5ahpKBSNSUFJoglTJecdwgrzhukAef3sNX7nyc7/xmCzf+02Zecewg7zlzFX98wnJSKQsdKnzpS8n6kkvCxiEiQaj5KJBnxya58Z8283d3b2Lr7gmeP9jDu192NG8ZPpLezoC5+uyzk/X69eFikLpM5wuMTeYZn8wzNjnN3ok8eyen2TuZZ29p2zR7p/LsnUh+To7NJ9snZ4+fqjCfV7UrRKVLh1c5utplpp7LT7Vr1aLGV8ex1SKpHEeVM9RRxhc/r48b/vSMasEckJqPWtCyniz/5Y+PZc1Zx3D7/U/zlTsfZ+0tv+eqHz/E6pes5K3DR/KiI5aGDlMW0EJevEvHT+aZnK59YsaUQU+2g1xnmly2g1w2TS6bpj+X5Yj+NNmOynepV6vDmlXeU3FrlZNYlR2VTl09jtrPXfXYqhX1Os5R7QwVy1J7uSude+XA4jT1KikElkmneMMpR/CGU47g3s07uP7Xm/jGP27m+l9v4sXP6+Otpx/JG045gr7uTOhQ216h4ExMFxifyrNvKs/4VHIBL39d2jeZb/rFuyfbQX8uy/MGyi/osxf2ns7ZbT3ZNN3Fbd2Z2X2dHamqF3IRUPNRS9q5d5Lv3fsHbtqwhX/9t910dqQ489hBXnXicl514nKO6O9evA9voeaj6XyBiekCk9Oz68l8nn1TBSbz+28fn8qzr3jRLr+oz24r7HeBf+7+5Lz1Sqds9qKc7UguxKWLefkFev8Ldfm+XOfshX3mIq+LtyymQ6b5yMyWAdcBrwa2Ax9z92+Ejar5+nNZLj7z+bz75au4/w+7ufneLfz0X7fxswe2AXD40i6OGeqhP5chm06RTqVIpyitO1IpUmZ0pC1Zp4xUKlmnU8m2at6wcxyA769/tOJ+xykUnOmCky9f5518ofDc7YXi9vz+25MLej65qOcLTOx3oU/+qj7YweKdHSm6s2m6M8nSlUnTlUm2DeQydBW3zxzTVfa6O5OmM5Oqur8rk9bFW9pSS9UUzOxGkqk3/hQ4FbgNeLm7b6x0fLvWFCpxdx4dGWX9gyNsfGo3m5/dy869k0wVL7b5gpP3steF525faClLElB6JuGkZxPPftvLlpmfsx0pOjuSNuzkdbJk0yk6M2my6dntlY7NdqToTKfozKTIptN0Z1P7XeS7OtKtcTeXSAs6JGoKZtYDrAZOdvdR4E4z+z7wLuCjQYNrAWbGscuXcOzyxsY1uDsFh+lCgULhQJ1q80unjLSZLroibahlkgJwPJB394fKtv0WeGWgeNqKmZE2SKc0zYaIVNdKs6T2ArvmbNsF7PensZmtMbMNZrZhZGSkacGJiMSglZLCKDD3xvylwH4PLXD3de4+7O7DQ0NDTQtORCQGrZQUHgI6zOy4sm2nABU7mUVEZOG1TFJw9zHgZuATZtZjZmcCbwRuCBuZiEg8WiYpFF0CdAPbgBuB91W7HVVERBZeK919hLs/C7wpdBwiIrFqtZqCiIgE1FIjmutlZiPAEw2+fZBkKo12oLK0pnYpS7uUA1SWGUe7e8XbNw/ppHAwzGxDtWHehxqVpTW1S1napRygstRCzUciIlKipCAiIiUxJ4V1oQNYQCpLa2qXsrRLOUBlmVe0fQoiIvJcMdcURERkDiUFEREpUVIQEZGStk0KZrbMzL5rZmNm9oSZveMAx/6ZmT1tZrvM7Ctm1tnMWOdTa1nM7GQz+5GZbTezluwsqqMs7zaze8xst5ltMbPPmFlLTctSR1neZmYPFv9/bTOzr5rZ3Gnig6rn96XsPT8zM2+l76WO7+RiM8ub2WjZcnZzoz2wOq9hx5jZrWa2p/j7/5lGP7dtkwLwRWASWAFcBFxjZifNPcjMXkPyuM9zgFXAMcD/aF6YNampLMAU8E2SZ1y3qlrLkgMuIxm1eQbJ9/OhJsVYq1rLchdwprv3kfz/6gCubFqUtam1LACY2UW02NxpRfWU42537y1b1jcryBrVeg3LAj8BfgYcDqwEvtbwp7p72y1AT/Ef8/iybTcAf1Xh2G8Anyr7+Rzg6dBlaKQsZfuPTb7a8PEfbFnKjvsAcEvoMhxsWUieMPh3wA9Cl6HRsgB9JM8/eSngQEfoMtRbDuBi4M7QMS9QWdYAv1qoz27XmkK15z1X+ovhpOK+8uNWmNlhixhfPeopS6s7mLKcRWs9cKmuspjZK8xsF8mTBFcDVy96hLWr93v5FHAN8PRiB1anestxWrGp5SEzu6KVmsGorywvBTaZ2e3F8qw3sxc3+sHtmhRqet5zlWNnXlc6NoR6ytLqGiqLmb0HGAY+t0hxNaKusrj7nZ40H60EPgtsWtTo6lNzWcxsGDgT+N9NiKte9XwnvwROBpaTJOm3Ax9e1OjqU09ZVgJvA74AHAHcBvxDsVmpbu2aFGp63nOVY2deVzo2hHrK0urqLouZvQn4K+A8d2+l2S0b+l7c/Q/AD4G/X6S4GlFTWcwsBXwJeL+7TzcptnrU/J24+2Pu/ri7F9z9X4BPAG9pQoy1quf/1zhJU9jt7j5J8sfTYcALG/ngdk0K9TzveWNxX/lxW939mUWMrx7t9OzquspiZq8F/ha4oPiL20oO5nvpAF6wKFE1ptayLCWpsd1kZk8D/1zcvsXM/v3ihzmvg/lOHLBFiaox9ZTldyTxL4zQHSqL2FHz9ySP9Owhqe7uAk6qcNxrSdpGXwQMkPTgz9vx2aJlMaCrWBYvvu4MHX+DZXkV8AxwVuiYF6AsFwFHFb+fo4FfADeHjr/eshTjP7xsOb34/+x5QDZ0Ger8Ts4DVhRfnwjcD/xF6PgbLMsJwF7gXCAN/BnwaKPfSfCCL+I/6DLge8AYsBl4R3H7USRVs6PKjv0AsBXYDfzfFryQ1lQWkltqfc6yKXT8DZbl58B0cdvMcnvo+BssyyeBLcXjtpBMZHZY6PgbKcuc98z8f2uJu4/q/E4+V/ydHwMeI2k+yoSOv9HvBHgz8EjxGraeCsmj1kUT4omISEm79imIiEgDlBRERKRESUFEREqUFEREpERJQURESpQURESkRElBRERKlBREFpGZWXG9tvxnkValwWsii8jMLgUmSKZCngZ+7O4/DxuVSHWqKUhUzGytmU2a2fdq+avdzM42s02Nfp67/x+SObX+G3BbIwnhYGMQqYeSgsTmfwKXAG9k/9lxS8zsb8zs87WczMxeaGb7zOzeKvsvAXaSzHX/ulZ7DrDIXEoKEhV33w1cT9Kk85ynUxVrDxcA/1DjKf8XSaI5qcpDTa5x93XAmLt/jGSGVJGWpaQgMcqQzDJ5coV9p5NMOX7nfCcxszeTTFX8meI5XzT3GC922rn72vKfy86xxcw+MGfbi4u1j+ecT2SxKSlIjD5O8mSqSs+xfRNJ2/8BnyxmZt0k0y9/yN13Ak8BpzUQy90kiajc1cCX3f33DZxP5KAoKUhUzOxkkoeQ3ELlmsIbSeawn89HSR6BONOXsBE4tYGQ9ksKxcePngb8RQPnEjloSgoSjWJ/wTqSR3x+GTjSzPrK9h8LHAP8aJ7zrAL+K/DfyzbfT2NJ4f8BLzCzZWbWSVL7+IS3zuNgJTIdoQMQaaL/DBwJvIak+QiS2sJdxddvAn7q7mPznOfzJLeZPl52V6sBo2Zmc/sN5nEPMEny7OPTSMYyfLGO94ssKCUFiYKZ/RHwaeA/uvseYI+Z7SbpV5hJCm8EvjrPeV4NvBJ4CcnFfMaJwLeA55M83rEm7j5RvJ31AuDdJI9cnKr1/SILTUlBYvEF4Cfu/v2ybfdT7FcwsyHgpcBbqp3AzDIkt6Be5e6/mbNvc/HladSRFIruBt5fjO/WWt5Q7Hu4hOSZ4juADwI3ufuX6/xskf2oT0Hanpm9DvgTklHF5X7H7B1IFwD/7O5bD3Cq9wNDJAlmP8XxD0/TWL/CfUAB+MA8x5U7C3gt8CrgHcDrgGPMrKuBzxcpUU1B2p673wb0V9j+vrIf573ryN0/R9IRXG3/HzUWIRcB17r7xjreM+nuBTN7DHghMAXsRb/TcpBUUxBJ3AXc2KwPM7OUma0ws4+Q1FYur/MUj5jZT0l+h38M/ApIu/voAocqkdFfFSKAu3+myR95FvAz4EFgtbvvqOfNxb6D8v6Dry1gbBIxJQWRA9tEMsJ4Qbn7emqvqS9KDCKV6HkKIiJSoj4FEREpUVIQEZESJQURESlRUhARkRIlBRERKVFSEBGREiUFEREpUVIQEZESJQURESn5/77Afkc/GQNqAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Code from Gabriel Peyré\n",
    "\n",
    "q = 200\n",
    "D = A_train.shape[1] # dimension of pbm \n",
    "lmax = np.linalg.norm(A_train,2)**2\n",
    "lambda_list = lmax*np.linspace(.6,1e-3,q)\n",
    "niter = 500\n",
    "step_size = 1e-5\n",
    "\n",
    "#Soft threholding operator\n",
    "def Soft(x,s): return np.maximum( abs(x)-s, np.zeros(len(x))  ) * np.sign(x)\n",
    "\n",
    "#ISTA \n",
    "C = A_train.transpose().dot(A_train)\n",
    "u = A_train.transpose().dot(Y_train)\n",
    "def ISTA(x,lasso_parameter,step_size): return Soft( x-step_size*( C.dot(x)-u ), lasso_parameter*step_size )\n",
    "\n",
    "\n",
    "X = np.zeros( (D,q) )\n",
    "E = np.zeros( (q,1) )\n",
    "x = np.zeros(D)\n",
    "\n",
    "for iq in np.arange(0,q):\n",
    "    Lambda = lambda_list[iq]\n",
    "    # ISTA #\n",
    "    for i in np.arange(0,niter):\n",
    "        x = ISTA(x,Lambda,step_size)\n",
    "    X[:,iq] = x.flatten() # bookkeeping\n",
    "    E[iq] = np.linalg.norm(A_test.dot(x)-Y_test) / np.linalg.norm(Y_test)\n",
    "\n",
    "# find optimal Lambda\n",
    "i = E.argmin()\n",
    "lambda0 = lambda_list[i]\n",
    "xSparse = X[:,i]\n",
    "print( 'Lasso: ' + str(E.min()*100) + '%' )\n",
    "# Display error evolution.\n",
    "plt.clf\n",
    "plt.plot(lambda_list/lmax, E)\n",
    "plt.plot( [lambda0/lmax,lambda0/lmax], [E.min(),E.max()], 'r--')\n",
    "plt.axis('tight')\n",
    "plt.xlabel('$\\lambda/|A^* y|_\\infty$')\n",
    "plt.ylabel('$E$')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.12 ('base')",
   "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.12"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "965fe59137d5334f034bf4d21dcfaab0895405927d7ec8dd4d8796ddf51dcb72"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}