gradient_descent.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from numpy.random import randn, randint\n",
    "from scipy.linalg import norm\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate random data matrix\n",
    "random.seed(0)\n",
    "n,d = 6,4\n",
    "X = randn(n,d)\n",
    "\n",
    "# optional: give it linearly dependent columns\n",
    "# X[:,3] = X[:,2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# form data from noisy linear model\n",
    "wtrue = randn(d)\n",
    "y = X.dot(wtrue) + .1*randn(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# look at least squares objective as a function of w\n",
    "def f(w):\n",
    "    return norm(y - X.dot(w))**2\n",
    "\n",
    "def gradf(w):\n",
    "    return 2*X.T @ (X @ w) - 2*X.T @ y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fa3f7cf6a90>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEGCAYAAACevtWaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3hUZfbA8e+BEIqAtCDdBClCkBpExEYVkaWKoghBKaIugguuIuruT+yoa8GGCARFECtYVkUEWVkWCZ1QJFQDQTEUQTq8vz/OIAGTkDKTOzM5n+eZJzNzp5y5hJM7733fc8Q5hzHGmPBSyOsAjDHG+J8ld2OMCUOW3I0xJgxZcjfGmDBkyd0YY8JQhNcBAFSoUMFFR0d7HYYxxoSUJUuW/Oqci8poW1Ak9+joaBITE70OwxhjQoqIbM1smw3LGGNMGLLkbowxYciSuzHGhKGgGHM3xgTWsWPHSElJ4fDhw16HYnKhWLFiVKtWjSJFimT7OZbcjSkAUlJSKFWqFNHR0YiI1+GYHHDOkZaWRkpKCjExMdl+ng3LGFMAHD58mPLly1tiD0EiQvny5XP8rcuSuzEFhCX20JWbf7uQTu579sDw4bB3r9eRGGNMcAnpMffkZBg3Dvbtg0mTvI7GGGOCR0gfuTdvDqNGweTJMGuW19EYY0zwCOnkDvDww9CoEQweDL/+6nU0xphzGTFiBPXr12fQoEFcffXVnDhxwuuQADh69ChXXXUVx48f9zoUvwj55B4ZCVOmwO7dcPfdXkdjjMnKpk2bWLBgAWvWrKFx48b06NGDwoUL5+q15s2bR//+/f0WW2RkJG3btuW9997z22t6KaTH3E9p2BD+7//gwQehe3fo3dvriIwJXsOHw/Ll/n3Nxo3hhReyfsz69etp164dx48fp0mTJgB8/PHHAGzevJn+/fvz3XffsXTpUpo1a8avv/5KmTJlaNSoET/88AMlSpTIUUxJSUkMGzaMbdu20bdvX3755Rf69etHoUKFGDlyJHPnzgVg9erVDBo0iIULF9KtWzdGjRpFnz59cr4TgkxYJHeA++6DmTPhrrvgyiuhalWvIzLGpFe3bl3i4+OJjo6mX79+1KhRg1OlvsuUKcP+/fsBePnll7nsssvYs2cPCxYsoH379jlO7IcPH6ZXr168//771KxZk4svvphmzZrRvHlzDh48yI8//vjHYx955BHGjBkDQIMGDVi8eLF/PrDHzpncRWQi0Bn4xTnX4KxtI4GxQJRz7lfffaOAAcAJ4B7n3Fd+jzoDERE6PNO4MQwYAP/+N9i0XmP+7FxH2IG0atUqunbt+sdR+Snnn38+Bw8eJC0tjdTUVFq1asWePXsYP348zz///Bmv0aJFC44cOcKBAwfYvXs3jRs3BuDpp5/m2muvBeCbb76hSZMmxMbGAjqePmLECABKlChBsWLF2Lt3L5s2bWLPnj20a9cOgMKFCxMZGcn+/fspVapUwPdHIGVnzH0y0PHsO0WkOtAe2JbuvvpAbyDW95xXRSR3A2q5UKcOPPssfPUVvP56fr2rMSa7kpKSiI2NpXjx4mesuCxUSFPRm2++yYABAyhVqhQrV67kxIkT1KlT54zXWLRoEcuXL2fChAl06dKF5cuXs3z58j8SO8CyZcto2rQpADt27KBkyZK0atXqj+3169dn3bp1PPzwwzz22GNnvP6RI0coVqyY3z97fjtncnfOzQd2Z7DpX8DfAZfuvq7AdOfcEefcZiAZuNQfgWbXnXdC+/YwciRs2JCf72yMycr+/fspUqQIJUqUoGzZspw4ceJPCX7WrFl0796d0qVL8+yzzzJkyJBcvVfRokVJSUkBYNSoURw9evSM7bGxsUyaNAnn3BlJPy0tjaioqBwV6ApWuZotIyJdgO3OuRVnbaoK/JTudorvvnwjoguaIiOhb18Ik1lNxoS81atX06DB6ZHdDh068P333/9xOzIykuuuu46IiAhKly7N77//TufOnXP1Xrfccgvz58+nbt26NGrUiJYtWzJ8+PA/tsfGxvLmm2/+MdZ+yty5c+nUqVOu3jPoOOfOeQGigdW+6yWARcD5vttbgAq+668At6Z73ltAz0xeczCQCCTWqFHD+dv06c6Bc//8p99f2piQs2bNGq9D+JOlS5e6W2+91eswztC9e3e3bt06r8PIUEb/hkCiyyRv5+bI/SIgBlghIluAasBSEamEHqlXT/fYasCOTP6ojHfOxTnn4qKiMuzvmic33QR9+sCYMbBokd9f3hiTR02aNKF169ZBtYipW7du1K1b1+tQ/CLHyd05t8o5V9E5F+2ci0YTelPn3E5gFtBbRIqKSAxQG/jBrxHnwLhxUKUK3HorHDjgVRTGmMzcfvvtuV7E5G+RkZH069fP6zD85pzJXUSmAQuBuiKSIiIDMnuscy4JmAGsAb4E7nbOefZnuUwZnR65cSP87W9eRWGMMfnvnPPcnXM3n2N79Fm3Hwcez1tY/nPNNfD3v8PTT0OnTtCtm9cRGWNM4IV8bZnsePRRaNoUBg6EHRmeATDGmPBSIJJ7ZCS8+y4cPAj9+8PJk15HZIwxgVUgkjtA3bq67Hr2bG+XXxtjTH4oMMkdYNAgHXN/4AFYutTraIwxJnAKVHIXgQkToGJFuPlmmx5pjAlfBSq5A5QvD++8o3Vnhg3zOhpjTEExadIkhgwZQkxMDEOGDOGNN94I6PsVuOQOOj3ywQdh4kQIk6YrxoSM3LTZ87oFnj9ivu2223jiiSdo27Ytr7/+OnfccUdAYy6QyR3gH/+Ali219+rmzV5HY0zBkNs2e9lpgZfbtnuxsbFUrlyZWrVq/XEpVaoUo0eP9nvMS5Ys+aMUcaAV2ORepIhOjxTRtnxnVQQ1xvjZ+vXrufrqq9m6dStNmjRhwoQJdO3aFdA2e1dffTUAS5cuRURIS0vjxIkTNGjQgIMHD9KtWzemTp2a6/dPSkqiXbt21KlThzFjxjB06FAWL17MwIED6du3L8nJySQnJ7NhwwYuuOACBgwY4PeYlyxZQrNmzc643bp16z9ur169mpYtW+b6M6YXNm32ciM6Wk+w9uoFDz0EzzzjdUTG5AOPmqjmtc1eXlrgZdV2r2bNmjRt2pQnnniCiIgI5s2bR0xMDDVr1gTwa8wrVqxgWLqTffXq1cu05V9eFdgj91NuuAGGDIGxY7U1nzEmcFatWkWjRo1y1GbvzjvvBM5sgZdeixYtaNy4MQMHDmTWrFk0btyYxo0b89VXpzt8pm+7V7x48TPa7pUvX56WLVvy2WefAfDWW28xcODAgMQ8bdo0ihcv/sdrpG/5t3Tp0jNa/uVVgT5yP+X552HBAujXTw9orLm2CWseruI71WbvyJEjWbbZW7NmTYZt9jJqgbfIV9N73rx5TJ48mcmTJ//pfc/Vdm/QoEG88MILtG7dmvnz5zNx4sSAxpzeqZZ/Y8aM+VPLv7wo8EfuAMWLw4wZcOgQ3HKLdW8yJhDy2mYvLy3wztV2r02bNvz4448899xz9OrVi8jIyHyLObOWf3llyd3n4ovhtddg/nz4v//zOhpjwk9e2+zlpQXeudruiQj9+/fn8ccfP2NIJj9izqzlX55l1qIpPy/NmjXLW/8pP7rtNudEnPv6a68jMcZ/wqHNXqBb4B04cMCtWrUqy8d4GXN+tNkLa+PGQf362qLPygMbEzg5abOXHy3wzjvvvDOO0jMSbDFnRTT5eysuLs4lJiZ6HcYf1q6FuDho1gy+/RYi7LSzCXFr166lXr16Xodh8iCjf0MRWeKci8vo8dlpszdRRH4RkdXp7hsrIutEZKWIfCwiZdJtGyUiySKyXkSuzcNn8Uy9ejB+PPznP+BbpGaMMSElO8Myk4GOZ903G2jgnGsI/AiMAhCR+kBvINb3nFdFJDi63+ZQnz5wxx26sOnTT72Oxhhjcuacyd05Nx/YfdZ9XzvnTk0Y/B9QzXe9KzDdOXfEObcZSAYu9WO8+eqFF6BJE53/bvVnTKgLhiFYkzu5+bfzxwnV24FTazurAj+l25biuy8kFSsGH3yg13v2hHTTW40JKcWKFSMtLc0SfAhyzpGWlpblQqiM5OlUoYiMBo4DpyrjSEaxZfLcwcBggBo1auQljICqWROmTIEuXWDoUHjzTa8jMibnqlWrRkpKCrt27fI6FJMLxYoVo1q1aud+YDq5Tu4iEg90Btq604cDKUD1dA+rBmQ4odA5Nx4YDzpbJrdx5Ie//EXrvz/xBFx+Odx2m9cRGZMzRYoUISYmxuswTD7K1bCMiHQE7ge6OOcOpts0C+gtIkVFJAaoDfyQ9zC99+ij0LYt3Hmn9V81xgS/7EyFnAYsBOqKSIqIDADGAaWA2SKyXEReB3DOJQEzgDXAl8Ddzrlzz/YPAYULw7Rp2n+1Rw9IS/M6ImOMyZwtYsqhxYvhiiu0Vd8XX2jSN8YYL+RpEZM5U/Pm8Mor8PXX8MgjXkdjjDEZs+SeCwMH6uWJJ+Cjj7yOxhhj/sySey6NGwctWkB8PKxZ43U0xhhzJkvuuVS0KHz4IZQsCd26wd69XkdkjDGnWXLPg6pVdQXrli1aiyYbVUCNMSZfWHLPo1at4KWXdObMQw95HY0xxiirVO4HQ4ZoY+2nnoKGDeHmm72OyBhT0NmRu5+89BJceSXcfjssWeJ1NMaYgs6Su59ERur4e8WKeoJ1506vIzLGFGSW3P2oYkWYORN279YEbyWCjTFeseTuZ40bw9tvw6JFMGgQBEF1B2NMAWTJPQB69IAxY+Cdd+Dpp72OxhhTENlsmQAZPVpXrj74IFx8sQ7TGGNMfrEj9wARgbfegksv1QVOVgPeGJOfLLkHUPHieoK1QgXt5rR9u9cRGWMKCkvuAXbBBfDZZ7B/vyb4Awe8jsgYUxBYcs8Hl1wC06fDihVwyy1Wg8YYE3jZabM3UUR+EZHV6e4rJyKzRWSD72fZdNtGiUiyiKwXkWsDFXio6dRJV7F++in87W9eR2OMCXfZOXKfDHQ8674HgDnOudrAHN9tRKQ+0BuI9T3nVRGxRnQ+d98N996rSf6ll7yOxhgTzs6Z3J1z84HdZ93dFUjwXU8AuqW7f7pz7ohzbjOQDFzqp1jDwtix0L27JvlZs7yOxhgTrnI75n6Bcy4VwPezou/+qsBP6R6X4rvP+BQurIub4uKgd29dyWqMMf7m7xOqksF9GS7AF5HBIpIoIom7du3ycxjBrUQJHXuvXFln0CQnex2RMSbc5Da5/ywilQF8P3/x3Z8CVE/3uGrAjoxewDk33jkX55yLi4qKymUYoatiRfjySzh5Eq67DgrY3zdjTIDlNrnPAuJ91+OBmenu7y0iRUUkBqgN/JC3EMNX7dp6BJ+SAp07w++/ex2RMSZcZGcq5DRgIVBXRFJEZADwFNBeRDYA7X23cc4lATOANcCXwN3OOZvVnYWWLXUOfGIi3HgjHDvmdUTGmHAgLghq0sbFxbnExESvw/DU+PFwxx3Qvz9MnKi1aYwxJisissQ5F5fRNqsKGSQGD4bUVPjnP6FSJXjySa8jMsaEMkvuQeSRRzTBP/WUnnC9916vIzLGhCpL7kFEBF55BdLStERB+fLQr5/XURljQpEl9yBzapHT7t1w++1QrpzOpDHGmJywqpBBqGhR+OQT7cfaqxd8953XERljQo0l9yBVqpQucoqJ0VWsBXwykTEmhyy5B7EKFWD2bB1779hRe7IaY0x2WHIPclWraoIvUgTat4eNG72OyBgTCiy5h4BatTTBHz4MbdvCTz+d+znGmILNknuIaNAAvv4a9uzRBL9zp9cRGWOCmSX3ENKsGXzxBWzfDu3aWSVJY0zmLLmHmFattJLkxo06Br/77B5ZxhiDJfeQ1KYNzJwJa9dChw6wd6/XERljgo0l9xDVoQN8/DGsXAnXXgv79nkdkTEmmFhyD2GdOsEHH8CyZZbgjTFnsuQe4rp0gfffhyVLLMEbY04L/eQeBM1GvNa1qx7Bn0rwNgZvjMlTcheRe0UkSURWi8g0ESkmIuVEZLaIbPD9LOuvYP8kJUVX+Dz8MCQnB+xtQsGpBL90qc2iMcbkIbmLSFXgHiDOOdcAKAz0Bh4A5jjnagNzfLcDY98+Te5PPKHdplu10n51BfTQtWvX0ydZ27aFX3/1OiJjjFfyOiwTARQXkQigBLAD6Aok+LYnAN3y+B6Zi42Fr76Cbdu0fdGePdqItFIl6N0b/v1vOH48YG8fjK6/HmbNgnXroHVr+PlnryMyxngh18ndObcdeBbYBqQC+5xzXwMXOOdSfY9JBSr6I9AsVa0K998PSUnwww8wcKAWY+nUCapXh/vug9WrAx5GsLj2Wvj8c9i0Ca66SkevjDEFS16GZcqiR+kxQBXgPBG5NQfPHywiiSKSuMtf6+hFoHlzGDcOduyADz+ESy+FF16ASy7R9fsvvVQg1u23aaO1aHbuhCuv1ERvjCk48jIs0w7Y7Jzb5Zw7BnwEXA78LCKVAXw/f8noyc658c65OOdcXFRUVB7CyETRotCjhy7l3LFDE7xzMGwYVKkC3brpAPXRo/5/7yDRqhXMmQO//aYJ3urBG1Nw5CW5bwMuE5ESIiJAW2AtMAuI9z0mHpiZtxD9ICpKk/rSpXq2cdgw+N//NPlXqQJDh2qrozCcVhkXp236Tp7UIRrr6GRMwZCXMfdFwAfAUmCV77XGA08B7UVkA9Dedzt4XHIJPPusDkR//rlOK3nzTR3OadAAnnlGj/TDSIMG8P33ULq0DtfMm+d1RMaYQBMXBEercXFxLtHLQ8o9e2DGDEhIgIULoVAhnSzev7/OLyxe3LvY/Gj7dq1Js3EjTJ+uI1PGmNAlIkucc3EZbQv9Far+ULasTqH8739h/XoYNUoHqG++GSpXhsGDYcGCkB+2qVoV5s+Hxo2hZ0+YMMHriIwxgWLJ/Wx16sBjj8GWLXo2smtXePdduOIK3fboo7otRJUvrx+rQwcYNEjXf4X43yxjTAYsuWemUCEdoE5I0PmEkyfrnPl//ANiYuCaa2DSJNi/3+tIc+y883ShU58+MHq0nk8+ccLrqIwx/mTJPTtKloT4ePj2Wz1qHzNGB7Bvv11Xw/brp4fDJ096HWm2FSkCU6bAyJHwyivQqxccOuR1VMYYf7HknlMXXggPPQQ//qjj8LfeqofB7dpBdLQeCq9f73WU2VKoEIwdq0sAPvlEP0JamtdRGWP8wZJ7bonA5ZfDG29AaqpOP2nQQGvcXHwxtGwJr72mM3GC3LBhOlloyRINe+NGryMyxuSVJXd/KF4cbroJvvhC58+PHQsHDsBdd+mwzY036pz6IC5idsMNOrK0ezdcdpnOCDXGhC5L7v5WubIOZK9cqYfCQ4bA3LnQuTNUqwYjRui2INSqlSb1MmW0ouSMGV5HZIzJLUvugSICTZvCiy/qyddPPtExj5dfhkaNoEkTHez+JcPSO56pXVsTfPPm+mVkzBibKmlMKLLknh8iI0930tixQytTFi4M996rK4u6dNEKlkeOeB0pABUqwDffQN++8Mgj+vPwYa+jMsbkhCX3/FahwulCZatXa4JPTNRB7ypV4O67tSa9x4fLRYvqFP/HH4epU3WYJjXV05CMMTlgyd1LsbFaqGzbNu0a1aEDTJwILVpA/fo688bDThsi8OCD+qVi5Uodqlm82LNwjDE5YMk9GEREQMeOMG2aroYdP17rBIwaBTVqaNKfOhUOHvQkvB49tOxORISWDZ461ZMwjDE5YMk92Jx/vhZ9+f572LDh9IKpW2/VaZUDBmj1r3xeDduokY4WXXqphvK3vwX1zE5jCjxL7sGsVi0tVLZpk06n7NlT5ydefbVu++c/87V/XsWKeqJ16FD417/0C0UB6FhoTEiy5B4KChU6Xahs504tClOzpib+iy7SZD9xovbTC7AiRXSyz+TJOlTTrJke0Rtjgosl91Bz3nk6N/Gbb7SI2eOPa8IfMECHbW69FWbPDniZx/h4Te6FC2s15Nde83yCjzEmnTwldxEpIyIfiMg6EVkrIi1FpJyIzBaRDb6fZf0VrDlLjRo6nWXdOl15FB+vZQ46dNACZ6NG6bYAadpUF+G2b6+VFvr21aoLxhjv5fXI/UXgS+fcxUAjtEH2A8Ac51xtYI7vtgkkES0I89prOhl9xgxttzR2LNSrp1MrX31VC8f4Wbly8OmnOkI0bZpOl1y92u9vY4zJoVwndxEpDVwFvAXgnDvqnNsLdAUSfA9LAKxTZ34qVkyLs3/2mc6Rf+45XV56991a9+aGG7RE8bFjfnvLQoXg4Yd1NGjPHp1RM3GiDdMYkxXn9Hhs1qzAvH5ejtxrAruASSKyTEQmiMh5wAXOuVQA38+KfojT5EalSjpnccUKWLZMx07mz9dSCFWr6urY5cv99nZt2ujLXXaZngLo0ydfzvEaE3J279bJb3fdFbgCfXlJ7hFAU+A151wT4HdyMAQjIoNFJFFEEnfZfLrAa9xY5y9u366HClddpUM1TZroJPbnn4eff87z21SqpEfwY8bAe+/py9tsGmNO+/57/e/42Wf6xXrKlMC8T16SewqQ4pxb5Lv9AZrsfxaRygC+nxmWPXTOjXfOxTnn4qKiovIQhsmRIkXgL3+BDz7Q8flXXtGhnBEj9Gi+c2d4//08VQorXFjXXn33nY7+tGqljbitT6spyI4d00J8V1+ttQT/+1/9Yl0oQHMWc/2yzrmdwE8iUtd3V1tgDTALiPfdFw/MzFOEJnDKldPvhYsWwZo1Wod+2TJtLlKlCtx5J/zvf7kePL/iCh0R6tlTuw+2bg1bt/r5MxgTApKT9f/DmDHacnnZMoiLC+x7isvDWS8RaQxMACKBTcBt6B+MGUANYBvQyzmX5TSNuLg4l5iYmOs4jB+dOKEtmRIStETxoUNQt67+RvbtC9Wr5/glnYN33tFzuiK6CKpfP71uTDhzTjtxjhypR+tvvKHzHfxFRJY45zL8M5Gn5O4vltyD1G+/6RBNQgL85z+ajdu00fn0PXrogqoc2LxZn/qf/0C3bvqLXtFOt5swtX07DBwIX36pzecnTdJmbP6UVXK3Faomc6VLny5UtnGjDhhu2qSH3ZUqwW23wbx52S5iFhOjJXKefVYrHMfG6t8OY8KJc/D223DJJXreadw4+Oor/yf2c7HkbrKnZk0tVJacrL+xN96ohd5bt9b6No88on8AzqFwYT13u2QJREfry/TqFXTdBo3Jle3bdb5Cv366fnD5ch2ODNRJ06xYcjc5U6iQTqN86y2taTN1KtSpA489ppUqr7gC3nwT9u3L8mViY7ViwpNP6szM+vV1XD4IRgmNyTHnYMIE/b3+9luddTx/vv7X8Iold5N7JUrALbfod85t27Rz1O7dMHiwDtuc2pbJHMiICHjgAZ05UKeOnq/t1Mlm1JjQsmGDnooaNEjnr69YAcOH67dUL1lyN/5RrRrcfz8kJemqpQED9ExSx45a4Oz++3W6ZQbq19eTrC+9pD/r19dxeT9WSDDG744c0aKsDRvqAcqbb+pRe+3aXkemLLkb/xLR6mHjxukiqQ8+0KLvzz2n31lPbUtLO+NphQtrE5A1a6BtW7jvPp0HvHChR5/DmCzMm6dH6Q89pOv+1qzRmTFejK1nJohCMWGnaFFdwTRrFuzYoQORx49rFq9cWadTzpx5xiF6jRp618cf6wjP5ZfrlwCrUGGCQWqqtkxo3VqP3D//XGd8VanidWR/Zsnd5I+KFXUgctkyHZQcOlTXX3frpv8zhg2DpUvBOUT07rVr9Qh+yhQdkx83zvq2Gm8cO6ZfPuvU0WQ+erSWtu7UyevIMmfJ3eS/hg31f0pKilZPat0aXn9dh28aNtQB99RUSpaEZ56BlSt109ChWojsm2+8/gCmIPniC52zPnKkThRLStLJYSVKeB1Z1iy5G+9ERMD112vN0507tbh1yZJ6uF6tmh4Wvfce9WIOM3u2Tqv//Xft/NS1K/z4o9cfwISzpCS47jr9FXVOm9J89pnO+A0FltxNcChbFoYM0TOo69bpHMnVq6F3b6hUCRlyBz0qL2RNkuOpp3RWQv36ukDEFkAZf0pN1WmNDRvqr+Pzz8OqVXriNJTqIVlyN8Gnbl2dY7Zli47BdOmiK5wuv5xijepy//HH2TRvG4MHa32aWrW02p71bzV5sW+fdhSrVUvLKd1zjy66vvdeLfoVaqxwmAkN+/frtMqEBC1/IALXXMOODvGMWNCT6Z+VJCpKp6bdcYdO1DEmOw4d0rYGTz6pM7RuukmPLS66yOvIzs0Kh5nQV6rU6UJlmzZpnZutW6kyqj/T5lbil07x9K36LcOHnaR2bT0/e/So10GbYHb4MLz8sibx++7T3r9LlsD06aGR2M/FkrsJPTExWqgsOVmXtN58M1Hff8Jzy9vye8UYHj35EM/duYHatXXY5sgRrwM2weTQIZ1WW6uWDr3UqaPHDP/+NzRt6nV0/mPJ3YQukdOFynbuhHffpXiTesSnPskG6vDp7stZOuQNmsTs5cUX4eBBrwM2XjpwQGfZxsTotNqYGD0xP2+etr4LNzbmbsLPjh3wzju4hARkzRqOSlE+dl35pHQ89Yd34K57Iihf3usgTX75+Wcdfnn1VdizR6fSPvSQzlkPddaJyRRMzukgakICx96eRpF9aaRSifci+vB7z3h6PXqJpyVZTWCtWQMvvqjn4I8ehe7d4e9/hxYtvI7MfwKa3EWkMJAIbHfOdRaRcsB7QDSwBbjRObcnq9ew5G4C7uhR+PxzfhuXQIm5nxPhjrOEpiyJjafu/93CVT0qhNQcZpOxkye1yvSLL+rPYsW0ccaIEd7WVg+UQM+WGQasTXf7AWCOc642MMd32xhvRUZC9+6UnvMJET/v4LcxL1K5MgxOGsblN1RmTulufDbgY/b+YlNsQlFamo6n166tC5tXrtQSAT/9pCfVwzGxn0uekruIVAOuByaku7srkOC7ngB0y8t7GON3UVGUfugequxYwpHFK9lw3TCaHPkfnSf24PgFVfim3lCWT0jEnfR+yNJkzjk9GdqnD1StqtMZq1aFadN0/dvo0VChgtdReidPw/0xEJ4AABFjSURBVDIi8gHwJFAKGOkbltnrnCuT7jF7nHNlM3juYGAwQI0aNZpttfY7xkvHj5P86tfsfTGBBptmUowjbIisz09t4rl4zK1UiQvCmq4F1Nat2oA6IUFnw5Ypo2V4Bw/WAl8FSUCGZUSkM/CLc25Jbp7vnBvvnItzzsVFRUXlNgxj/CMiglr3dCJu43sc35bKgr6vc7jo+bT58n4uaF6dxRU6Mu+OaexNPeR1pAXS7t0647V1a22s/vDDWltuyhSdHPXyywUvsZ9Lro/cReRJoC9wHCgGlAY+ApoD1zjnUkWkMjDPOVc3q9eyE6omWG2d/SNbHp1CrYVvU/XENvZRmh8uvBHi42nxt1aUPt/OwgbKnj1aifH99/Xk6LFjp3vt3nqrJvmCLuBTIUXkGk4Py4wF0pxzT4nIA0A559zfs3q+JXcT7NyJk6x7fR77X5lC7LoPOM/9zkYuYmGdfkTc1o9r+kdTqZLXUYa+7ds1oc+cqTXjjh/XI/SbbtJ+602ahFZlxkDL7+ReHpgB1AC2Ab2cc7uzer4ldxNKTv52gI1jP4KEBGr/9C0Ac7mG/8TEE3FTT9p0LUXz5toX1mTt+HFYtEh7qX/xhTbjAqhZUzs09uypbXeDqTdpMLFFTMYEiNuylZ+fe5vIaQmUS0vmd0rwIT35uFQ8Ee1b07ptIdq21eEEO+LUGS5r1uiy/2+/hblztdRuoULQsiX85S96qVfP9ld2WHI3JtCcg4ULOfxGAoXen07kod/YXrg6k0/0ZQr92FuxLldcAVdeCZddBo0b6wKbcHfokLbNXbhQa7x9/73OSQet7dKmDXTsCG3bar8WkzOW3I3JT4cOwaxZuIQE+Oor5ORJkitcxmQXzytpN7GXshQpogm+WTMdR27SBBo0gOLFvQ4+937/XZtnLVuml8REXUx0qqn5RRfpH7crr9RZLzEx3sYbDiy5G+OV1FSYOlUnZa9ejStalO1NuzC7Sjzv7LqWxOUR/PabPlREx5pjY3VYok4dvVx0EVxwQXCMO588qR9p40btYbthgw6zJCXpwqFT6aRsWf2D1aKFXi69FCpX9jT0sGTJ3RivOaeHswkJ8O678OuvcMEFuFv6sL1dPIsONWTVKk2SSUm6OOfYsdNPL1oUatTQS5UqeqlUSVdgRkVBuXJw/vlQurT2NSlePHt/DE6e1C8a+/fr2Pdvv+mwya5dGmJqqs4j37EDtm3T5fzpm6BERuqS/9hY/eZxySWa1GvUsDHz/GDJ3ZhgcvSodoZISIDPPtMs3rgxxMfrfL+KFTl+XFdibtgAmzfrUfHmzZCScjrZpk/+GYmM1HH9woVPX06cOH05fPjc3aoiI/UPSeXKmrCjo+HCC/XbRO3aep/NCvKOJXdjgtWvv2pft4QEHaSOiIDrrtNE37lzps1gT57URT5pafoSaWl61L1vnx6FHz58+nL8uCbzkyf1aP5Uoi9WTI/wixfXI/5Tl/Ll9RtB+fI6vGJH4MHLkrsxoWDNGk3y77yjh+blykHv3promze3LGv+xBpkGxMK6teHp5/Wwe0vv4Rrr4WJE/WMZGwsPPWULuE0JhssuRsTbAoX1sT+7rvaG3b8eD2KHzUKqleHDh10Bo41hTVZsORuTDA7/3wYNEhX/2zYoM0/f/xRK2dVqgQDBsD8+afnIBrjY8ndmFBRqxY8+ihs2qTr9nv2hBkz4OqrdfrKP/+p24zBkrsxoadQIbjmGpg0SYdtpkzR5P7oo/rzqqvgrbf4Y3WUKZAsuRsTys47Twucz56tE+OfeAJ+/hkGDtRhmz594OuvdS6kKVAsuRsTLqpX15Ou69Zppa7+/bWO7rXX6sqjU9tMgWDJ3ZhwI6KlJ199VesHzJihK2DHjtWiNS1a6LbdWbZZMCHOkrsx4axYMejVS8scpKTAs89qMZm779aaAjfcoK2PzlXLwIQcS+7GFBSVKsGIEbBihbY8uusunUbZpQtUrQr33gvLl3sdpfGTXCd3EakuInNFZK2IJInIMN/95URktohs8P20EvzGBBMRLd34r3/pitdZs3SGzauv6v2NGsFzz+lMHBOy8nLkfhwY4ZyrB1wG3C0i9YEHgDnOudrAHN9tY0wwKlJE+9p98IHWs3nlFR3KGTlSO1Nff72O2R8+7HWkJodyndydc6nOuaW+6/uBtUBVoCuQ4HtYAtAtr0EaY/JB+fI6VLNokRYxu+8+HcK56SYdnx8yRGfh2GrYkOCXqpAiEg3MBxoA25xzZdJt2+Oc+9PQjIgMBgYD1KhRo9nWrVvzHIcxxs9OnNBO1gkJ8NFHejK2Th2tVNm3r06/NJ4JaFVIESkJfAgMd85le0mcc268cy7OORcXFRWV1zCMMYFQuDC0b69liHfu1JWvlSrB6NE6d75dO3j7bW2gaoJKnpK7iBRBE/tU59xHvrt/FpHKvu2VgV/yFqIxJiiULg233w7ffadNVB95RGvZ9OunCf+222DePO0KYjyXl9kyArwFrHXOPZ9u0ywg3nc9HpiZ+/CMMUGpZk0tVJacrMn+xhvhww+hdWutb/PII/oHwHgm12PuInIF8B9gFXDqT/WDwCJgBlAD2Ab0cs5luRTOOjEZEwYOHoSPP9bx+W++0ROvrVrp+PyNN2r5YuNX1mbPGJO/UlJ0nD4hQevZFCsG3bppom/f3rpq+4m12TPG5K9q1eCBB3RK5aJFOlb/1Vfa/Lt6dfj73yEpyesow5old2NM4IjApZfq4qjUVF0sFRcHzz8PDRro9XHjIC3N60jDjiV3Y0z+KFpUu0fNmqWrYf/1L51HP3SoLpLq0QNmzrQiZn5iyd0Yk/8qVoThw2HZMl0FO3QoLFig4/JVqsCwYVrcLAjOCYYqS+7GGG81bKiFyrZv19LErVvD669Ds2a67dlndUjH5Igld2NMcIiIOF2obOdOeO01KFlSa9xUqwadOsF771kRs2yy5G6MCT5ly54uVLZunc68Wb0aevfW1bB33AH//a8N22TBkrsxJrjVrQuPPw5btujiqC5ddA59q1a67bHHtDm4OYMld2NMaChUCNq2hSlTdNhm0iTtIPXwwxAdDW3a6KKpAwe8jjQoWHI3xoSeUqWgf3+YOxc2b4ZHH4Vt2/S+SpV0Jey33xboImaW3I0xoS06Wo/eN2yA77+HW26BTz7Ro/yYmNPbChhL7saY8CCi4/Djx+uwzbRpUK8ePPGENhg5tW3vXq8jzReW3I0x4ad4cZ1Z8+WX8NNP8PTTmtTvuEOHbXr3hi++gOPHvY40YCy5G2PCW5UqWqhs9WpYvBgGDdJZN9dfr0XMRo6EVau8jtLvLLkbYwoGES1U9vLLWtvmo4+gRQt48UVdCdu0qV7ftcvrSP3CkrsxpuCJjITu3fXE644dmtRFtN5NlSrQtasm/6NHvY401yy5G2MKtqgouOceWLJEh2eGD4cfftAKlpUrw1//ComJIbcaNmDJXUQ6ish6EUkWkQcC9T7GGOM3DRrA2LF6EvaLL7Rr1IQJ0Ly5bnvmGT3SDwEBSe4iUhh4BbgOqA/cLCL1A/FexhjjdxER2jVq+nSdVvnGG9oD9v779SRsx4461fLQIa8jzVSgjtwvBZKdc5ucc0eB6UDXAL2XMcYETpkyMHiwFipbvx5GjYK1a3WxVKVKOvvm+++DbtgmUMm9KvBTutspvvv+ICKDRSRRRBJ3hcnZaWNMmKtTRwuVbd4Mc+Zoc5Fp0+DKK6F2bS2DsGWL11ECgUvuksF9Z/xZc86Nd87FOefioqKiAhSGMcYEQKFCpwuV7dwJkydDjRrwj39oyYNrrtHCZvv3exdigF43Baie7nY1IDTOQhhjTE6ULHm6UNmWLXpkv2MH3H67Dtv07auLpk6cyNewApXcFwO1RSRGRCKB3sCsAL2XMcYEhwsvhNGjdWx+wQJN7J9+qrNuoqPhwQd1Wz4ISHJ3zh0H/gp8BawFZjjnkgLxXsYYE3RE4PLLtRfszp3aHrBhQ61xc/HFcNll2kZw9+7AheCC4AxvXFycS0xM9DoMY4wJrJ07YepUHaNfvVpXyv71r9ogPBdEZIlzLi6jbbZC1Rhj8kulSjBiBKxcqStihwzR4ZoAiAjIqxpjjMmciBYqa9o0YG9hR+7GGBOGLLkbY0wYsuRujDFhyJK7McaEIUvuxhgThiy5G2NMGLLkbowxYciSuzHGhKGgKD8gIruArXl4iQrAr34Kx58srpyxuHLG4sqZcIzrQudchjXTgyK555WIJGZWX8FLFlfOWFw5Y3HlTEGLy4ZljDEmDFlyN8aYMBQuyX281wFkwuLKGYsrZyyunClQcYXFmLsxxpgzhcuRuzHGmHQsuRtjTBgKueQuImNFZJ2IrBSRj0WkTCaP6ygi60UkWUQeyKfYeolIkoicFJFMpzaJyBYRWSUiy0Uk4P0FcxBXvu4zESknIrNFZIPvZ9lMHhfw/XWuzy7qJd/2lSISuC4LOYvrGhHZ59s3y0XkkXyKa6KI/CIiqzPZ7tX+OldcXu2v6iIyV0TW+v4vDsvgMf7dZ865kLoAHYAI3/WngaczeExhYCNQE4gEVgD18yG2ekBdYB4Ql8XjtgAV8nGfnTMuL/YZ8AzwgO/6Axn9W+bH/srOZwc6Af8GBLgMWJQP/27Ziesa4LP8+l1K975XAU2B1Zlsz/f9lc24vNpflYGmvuulgB8D/TsWckfuzrmvnXPHfTf/B1TL4GGXAsnOuU3OuaPAdKBrPsS21jm3PtDvk1PZjMuLfdYVSPBdTwC6Bfj9MpOdz94VmOLU/4AyIlI5COLyhHNuPrA7i4d4sb+yE5cnnHOpzrmlvuv7gbVA1bMe5td9FnLJ/Sy3o3/pzlYV+Cnd7RT+vCO95ICvRWSJiAz2OhgfL/bZBc65VNBffqBiJo8L9P7Kzmf3Yv9k9z1bisgKEfm3iMQGOKbsCub/g57uLxGJBpoAi87a5Nd9FpQNskXkG6BSBptGO+dm+h4zGjgOTM3oJTK4zy9zPrMTWza0cs7tEJGKwGwRWec74vAyroDss6ziysHL+H1/nSU7nz1gv1NZyM57LkXrixwQkU7AJ0DtAMeVHV7sr+zwdH+JSEngQ2C4c+63szdn8JRc77OgTO7OuXZZbReReKAz0Nb5BqvOkgJUT3e7GrAjP2LL5mvs8P38RUQ+Rr9+5ylZ+SGugOyzrOISkZ9FpLJzLtX39fOXTF7D7/vrLNn57AH7ncpLXOkThHPuCxF5VUQqOOe8LpDlxf46Jy/3l4gUQRP7VOfcRxk8xK/7LOSGZUSkI3A/0MU5dzCThy0GaotIjIhEAr2BWfkVY1ZE5DwRKXXqOnqCOMMz+/nMi302C4j3XY8H/vQNI5/2V3Y++yygn29Gw2XAvlNDSgF0zrhEpJKIiO/6pej/6bQAx5UdXuyvc/Jqf/ne8y1grXPu+Uwe5t99lt9njfN6AZLRcanlvsvrvvurAF+ke1wn9Iz0RnRoIj9i647+9T0C/Ax8dXZs6MyHFb5LUn7Elp24vNhnQHlgDrDB97OcV/sro88ODAGG+K4L8Ipv+yqymA2Vz3H91bdfVqATDC7Pp7imAanAMd/v1oAg2V/nisur/XUFOsSyMl3u6hTIfWblB4wxJgyF3LCMMcaYc7PkbowxYciSuzHGhCFL7sYYE4YsuRtjTBiy5G6MMWHIkrsxxoShoCw/YIzXfAWlXgRqAG+jRc2mOOcWexqYMdlki5iMOYuIFEMLTPUCNgHrgCXOuR6eBmZMDtiRuzF/1g5Y5pxLAvDVdXnO25CMyRkbczfmz5qgR+6ISBXggHNugbchGZMzltyN+bMjnO7w9STa4s6YkGLJ3Zg/exe4SkTWo9UDF4rICx7HZEyO2AlVY4wJQ3bkbowxYciSuzHGhCFL7sYYE4YsuRtjTBiy5G6MMWHIkrsxxoQhS+7GGBOG/h+61xhkKVNi6AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# gradient at w approximates f(w) near w\n",
    "w = randn(d)\n",
    "v = randn(d)\n",
    "alphas = np.linspace(-2,2,100)\n",
    "plt.plot(alphas, [f(w + alpha*v) for alpha in alphas], color='b', \n",
    "     label=r\"$f(w + \\alpha v)$\")\n",
    "plt.plot(alphas, [f(w) + alpha*gradf(w).dot(v) for alpha in alphas], color='r', \n",
    "      label=r\"$f(w) + \\alpha \\nabla f(w)^T v$\")\n",
    "plt.xlabel(r\"$\\alpha$\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fa3f7fd9f70>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEGCAYAAACevtWaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deVhVVffA8e9mdhYBZ3EeEBMVx0yt1LLSHErLyrRy6G3S5sxfWZlvg1Y2l5ZjNthgWtlbWZbmUI6pOKKAoiiCyDxe9u+PjYSKTN7LuVzW53nuI5x7uGcdhMW+++yzltJaI4QQwrW4WR2AEEII+5PkLoQQLkiSuxBCuCBJ7kII4YIkuQshhAvysDoAAH9/f92sWTOrwxBCiApl69atcVrrgMKec4rk3qxZM7Zs2WJ1GEIIUaEopaIu9pxMywghhAuS5C6EEC5IkrsQQrggp5hzL0x2djbR0dFkZGRYHYoQlvHx8aFx48Z4enpaHYqoYJw2uUdHR1OjRg2aNWuGUsrqcIQod1pr4uPjiY6Opnnz5laHIyoYp52WycjIwM/PTxK7qLSUUvj5+cm7V1EmTpvcAUnsotKT3wFRVk6d3IUQwqW99RasWuWQl5bkLoQQVoiKgscfh6++csjLS3IXQggrTJ8OSsHzzzvk5YtN7kqpJkqpNUqpvUqpMKXU5LztzymljimlduQ9ri/wNVOVUuFKqf1KqWsdEnkF9NxzzzF79uxyPcbll19+ya955swZ3nvvvXO22eN1C3r00Udp3749EyZMoF+/fthsNru+flllZWXRt29fcnJyrA5FuJJdu2DxYnjwQWjSxCGHKMnIPQd4VGsdBPQE7ldKtc977g2tdae8xyqAvOduBYKBQcB7Sil3B8RerrTW5ObmWh1GqW3YsOGCbaU9l8KSe2GvW1aHDx9m/fr17Nmzh06dOjFixAjc3cv2I/P7778zbtw4u8Xm5eVF//79+eKLL+z2mkLw9NNQsyZMneqwQxS7zl1rHQPE5H2crJTaCzQq4kuGAp9rrTOBCKVUONAd2FjWIJ//Low9x5PK+uWFat+wJtOHBBe5T2RkJNdddx1XXXUVGzdu5Ntvv2Xy5MkcPXqUjIwMJk+ezMSJE/P3u+KKK9iwYQONGjVixYoVVKlShZkzZ7J48WKaNGlCQEAAoaGhALz++uvMnz8fgPHjxzNlypT8Yw4aNIgrrriCTZs2ERISwl133cX06dOJjY1l6dKldO/e/Zw4L3YMgOrVq5OSklLouaxbt4633nqLrKwsevTowXvvvYe7uzuLFy9m9uzZKKXo2LEj2dnZHDp0iE6dOjFw4EBmzZqV/7oXO5eivicF7d+/nwEDBpCTk0Pnzp0BWL58OQARERGMGzeOP/74g23bthEaGkpcXBy1a9cmJCSEv//+m6pVq5b6/z4sLIzJkydz5MgRxowZQ2xsLHfccQdPPfUUa9asAWD37t1MmDCBjRs3MmzYMKZOncrtt99e6mMJcYF16+D77+Gll6BOHYcdplQ3MSmlmgGdgb+A3sADSqk7gS2Y0X0CJvFvKvBl0RT9x8Cp7d+/nwULFuSPXOfPn0+dOnVIT0+nW7du3HTTTQAcPHiQzz77jHnz5jFq1Ci+/vprgoKC+Pzzz9m+fTs5OTl06dKF0NBQtm7dyoIFC/jrr7/QWtOjRw/69euXn9zCw8P58ssvmTt3Lt26dePTTz/lzz//ZOXKlfz3v//l22+/zY9v69athR6juHPZu3cvX3zxBevXr8fT05P77ruPpUuXEhoaysyZM1m/fj3+/v6cPn2apKQkdu/ezY4dOy54zYudi6+vb6HfkzvuuOOcr2/bti1jx46lWbNm3HnnnQQGBnK2/HPt2rVJTk4G4O2336Znz54kJCSwfv16Bg4cWKbEnpGRwciRI/nyyy9p0aIF7dq1IzQ0lMsuu4wDBw7k7/fss88yY8YMADp06MDmzZtLfSwhLqA1PPkkNGwIDz3k0EOVOLkrpaoDXwNTtNZJSqn3gRmAzvv3NeBuoLCFubqQ15sITAQIDAws8tjFjbAdqWnTpvTs2TP/87feeit/ZHn06FEOHjxI/fr1ad68OZ06dQIgNDSUyMhI4uLiGD58eH4SuvHGGwH4888/GT58ONWqVQNgxIgRrFu3Lj+5N2/enMsuuwyA4OBg+vfvj1KKyy67jMjIyHPiW7duXaHHKO5cfv31V7Zu3Uq3bt0ASE9Pp27duiQmJnLzzTfj7+8PQJ06dUhKuvi7poudy4033ljo96Qwu3btYujQofmj8rNq1apFWloa8fHxxMTE0Lt3bxISEpg7dy6vv/76Oa/Ro0cPMjMzSUlJ4fTp0/nHfeWVV7j22n8v+6xevZrOnTsTHGx+prKysnj00UepWrUqPj4+nDlzhsOHD5OQkMCAAQMAcHd3x8vLi+TkZGrUqHHR74UQxVq+HDZuhLlzoQyDk9IoUXJXSnliEvtSrfU3AFrrkwWenwd8n/dpNFDwCkFj4Pj5r6m1ngvMBejatesFyd9ZnE1aYOZzV69ezcaNG6latSpXXnll/t2D3t7e+fu5u7uTnp4OFH4TitZFn27B13Jzc8v/3M3NrdALeyW90aXguWitGTt2LC+99NI5+7z11lulunGmqHO52PfkfGFhYQQHB5OZmXnO3ZhubuaS0Lx587jnnnvYs2cPO3fuxGaz0aZNm3Ne46+//gLM/9HChQtZuHBhocfavn07Xbp0AeD48eNUr16d3r17A9C+fXv27dvHjBkzePHFF8/5uszMTHx8fC56rkIUKzvbzLEHBcFddzn8cCVZLaOAj4G9WuvXC2xvUGC34cDuvI9XArcqpbyVUs2B1sDf9gvZOomJifj6+lK1alX27dvHpk2bity/b9++LF++nPT0dJKTk/nuu+/yt3/77bekpaWRmprK8uXL6dOnT5liutgxitO/f3+++uorYmNjATh9+jRRUVH079+fZcuWER8fn7+9Ro0a+dMjhR3/Us4lOTkZT09Pqlatiq+vLzab7YIEv3LlSoYPH07NmjWZPXs29957b4lf/3ze3t5ER0cDMHXqVLKysvKfCw4OZsGCBWit8xM+QHx8PAEBAVK8S1yajz6CAwfglVfAw/FlvUpyhN7AGGCXUurspOvTwGilVCfMlEskMAlAax2mlFoG7MGstLlfa+0c69ou0aBBg/jggw/o2LEjbdu2PWe6pjBdunThlltuoVOnTjRt2jQ/6XXp0oVx48blXxgdP358/pRMaV3sGMVp3749L774Itdccw25ubl4enry7rvv0rNnT6ZNm0a/fv1wd3enc+fOLFy4kN69e9OhQweuu+46Zs2adc7xCzuXi03BnG/37t106NAh//NrrrmGP//8M39KxMvLi+uuuw4PDw9q1qxJamoqgwcPLuF350K33XYbQ4cOpW3btkyaNInMzEymTJnCnDlzCA4OZuzYsRfMr69Zs4brr7/+Iq8oRAkkJ8Nzz0HfvnAJP7+lorW2/BEaGqrPt2fPngu2Cde3bds2fccdd1gdxjmGDx+u9+3bZ9nx5XfBBUyfrjVovWmTXV8W2KIvklflDlXhVDp37sxVV13lVDcxDRs2jLZt21odiqioTpyA2bPh5puhR49yO6zT1nMXldfdd99tdQj5vLy8uPPOO60OQ1Rk06dDVpZZ116OZOQuhBCOEhZmLqTedx+0alWuh5bkLoQQjvLEE1CjBjzzTLkfWqZlhBDCEVavNrXaX30V/PzK/fAychdCCHvLzTW12ps2NZUfLSAjdyGEsLfFi2HHDli6FCy6s1lG7kIIYU+pqaakb48eMHq0ZWHIyF0IIexp1iyIiTHt8yxscC4j93JUsEvSpXQyKqx5hqOPWdzx7dmZ6WxXpgcffJD09HSn6sx0VlljPL+zk3R6cjHHjpkLqKNGgZ27lZWWJPdLpMvYoelSOhmVJrkXd8yyxO/IzkwFuzK9/fbbzJ8/v9w7MwUHB9OgQQNatWqV/6hRowbTpk275BjP7+wknZ5czLRpYLPByy9bHYkk9+LMmDGDdu3aMXDgQEaPHs3s2bOJjIwkKCiI++67jy5dunD06FGGDRtGaGgowcHBzJ07N//rZ86cSdu2bRkwYAD79+/P3169evX8jz/55BO6d+9Op06dmDRpEjabLf8YEyZMIDg4mGuuuSa/ZO5TTz2V3xnp8ccfvyDm4o5ZWPyFxQCwePFiOnbsSEhICGPGjLno8c++9uuvv06HDh3o0KEDc+bMyT92Uedz1v79++nXrx9RUVF07tyZ1NRUli5dytChQ4mIiKBfv34AbNu2DaUU8fHx2Gw2OnToQFpaWmn/awFTbnjAgAG0adOGGTNm8OCDDzJu3DjGjBlDeHg44eHhHDx4kHr16nHPPfcUGSNQojiHDRvG0qVL82M4/3NRQW3bZi6kTp4MzZtbHU0FKRw2ebLW/frZ9zF5crFFeTZv3qxDQkJ0WlqaTkpK0q1atdKzZs3SERERWimlN27cmL9vfHy81lrrtLQ0HRwcrOPi4vSWLVt0hw4ddGpqqk5MTNQtW7bUs2bN0lprXa1atfzzHDx4sM7KytJaa/2f//xHL1q0SEdERGh3d3e9fft2rbXWI0eO1EuWLNFaax0REaGDg4MLjbkkxzw//ovFsHv3bt2mTRt96tSpc86xsONXq1Yt/9gpKSk6OTlZt2/fXm/bti3/ay52PgVNmzZNz5s3T2utdWZmpq5Xr57WWuvTp0/rzp07a621HjdunO7Zs6c+ePCgXrFihZ4yZcpF/w/XrFmjx44dW+hz6enpOigoSO/evVunpaXpwMBAPXz4cB0XF6cDAwN1dna21lrr3377TQ8YMKDYGEsaZ05Ojvb398//mvM/P58UDqsAcnO17tNH64AArc+cKbfDUkThMLmgWoQ///yToUOH5vf9HDJkSP5zJenQtGnTpmK7JF2sI1Lfvn1L3MmooJJ2ZiprV6ailKTDVHHnc7YrE3BOZ6bSdGWCknVmulhXJj8/P3r16sX333/PsGHD+Pjjjxk/fnyxMZY0zvM7O0mnJxfw9demN+qHH0KtWlZHA1SU1TIF3t6XJ11El6GSdmgqrquRvkhHpMjIyBJ3MjpfSToplXdXJihZZ6azXZkAqlSpkv99LE1XJihZZ6aiujJNmDCBOXPmcNVVV7F27dr8BuBFxViaOM/v7CSdniqwjAxzw1LHjnDPPVZHk0/m3ItwxRVX8N1335GRkUFKSgo//PBDoftdrENTSbokXawjUlGK64xU2s5MpenKVNTx7dmVCbigM1N5dmW6+uqrOXDgAK+99hojR47Ey8urRDGWJM7zOztJp6cKbs4ciIyEN96AMl74dwRJ7kXo1q0bN954IyEhIYwYMYKuXbtSq5C3XIMGDSInJ4eOHTvyzDPP5E93FOySdNNNNxWa6Ap2ROrYsSMDBw4kJiamyLj8/PzyOyOdf0G1JMcsaQzBwcH5XZlCQkJ45JFHijx+wa5MPXr0KHWHqfO7MsG/nZnAMV2Z1q5dS9u2bQkJCaFXr15MmTIFMO9+xo0bx8yZM8+ZkikuxpLEeX5nJ+n0VIGdOAEzZ8KwYXD11VZHc66LTcaX58OZOzElJydrrbVOTU3VoaGheuvWrRZHVLlY2ZkpJSVF79q1q9j9Shvj+Z2diuv05Cy/C6IQ48Zp7eWl9cGDlhweuaBadhMnTmTPnj1kZGQwduzY/DlaUT4KdmYq61r3sqpWrdoFo/TClCbG8zs7SaenCmzzZli4EJ58stxrtZeE0sVcBCsPXbt21Vu2bDln2969ewkKCrIoIiGch/wuOKHcXOjd28y1HzhgarZbQCm1VWvdtbDnZOQuhBCltXQpbNoECxZYltiLIxdUhRCiNFJSzFRMt27gxP11nXrkrrUu1TprIVyNM0ybivPMnGmqPn79Nbg57/jYaSPz8fEhPj5efrhFpaW1Jj4+Xm5uciYHD8Jrr5kRe69eVkdTJKcduTdu3Jjo6GhOnTpldShCWMbHx4fGjRtbHYY4a8oU01nplVesjqRYTpvcPT09ae4MldWEEALg++9Nw+vZs6F+faujKZbTTssIIYTTyMw0o/Z27SxreF1aTjtyF0IIp/Haa3DoEPz0E+TVGXJ2MnIXQoiiHDkCL74II0bANddYHU2JSXIXQoiiPPyw+feNN6yNo5RkWkYIIS7mp5/gm2/M2vbAQKujKZViR+5KqSZKqTVKqb1KqTCl1OS87XWUUr8opQ7m/etb4GumKqXClVL7lVLXOvIEhBDCITIzzcXT1q3h0UetjqbUSjItkwM8qrUOAnoC9yul2gNPAb9qrVsDv+Z9Tt5ztwLBwCDgPaWU81SwF0KIknjtNXPT0ttvQ4EuYhVFscldax2jtd6W93EysBdoBAwFFuXttggYlvfxUOBzrXWm1joCCAe62ztwIYRwmIgImDHDXES9tmJOPpTqgqpSqhnQGfgLqKe1jgHzBwCom7dbI+BogS+Lztt2/mtNVEptUUptkbtQhRBOZfJk0zLPov7N9lDi5K6Uqg58DUzRWicVtWsh2y4oEKO1nqu17qq17hoQEFDSMIQQwrFWroTvvoPnnoMmTayOpsxKlNyVUp6YxL5Ua/1N3uaTSqkGec83AGLztkcDBb8jjYHj9glXCCEcKC0NHnoIgoPN6L0CK8lqGQV8DOzVWr9e4KmVwNi8j8cCKwpsv1Up5a2Uag60Bv62X8hCCOEgL74IUVHw/vvg6Wl1NJekJOvcewNjgF1KqR15254GXgaWKaXuAY4AIwG01mFKqWXAHsxKm/u11ja7Ry6EEPYUFgazZsHYsdCnj9XRXDKn7aEqhBDlJjcXrrzSJPh9+6CCXAeUHqpCCFGURYtg3Tr46KMKk9iLI7VlhBCVW1wcPP449O4Nd91ldTR2I8ldCFG5Pf44JCbCBx84dU/U0nKdMxFCiNL67TdYuBAeeww6dLA6GruS5C6EqJwyMuDee6FlS3j2WaujsTu5oCqEqJxmzjSFwX75BapUsToau5ORuxCi8tmzB155Be64AwYMsDoah5DkLoSoXHJzYcIEqFEDXn+9+P0rKJmWEUJULh98ABs2mLXtLrKmvTAychdCVB5Hj8JTT8HAgTBmjNXROJQkdyFE5aA13Hcf2Gzw4YegCqtO7jpkWkYIUTksWwbff2/a5zVvbnU0DicjdyGE64uLM82uu3Y19dorARm5CyFc35QpkJAAq1eDR+VIezJyF0K4th9+gKVL4emnoWNHq6MpN5LchRCuKzERJk0ybfOmTbM6mnJVOd6fCCEqpyeegJgY+OYb8PKyOppyJSN3IYRrWr0a5s6Fhx+G7t2tjqbcSXIXQrie5GQYPx7atIEZM6yOxhIyLSOEcD1PPglHjsCff7pkxceSkJG7EMK1/PYbvP++mY65/HKro7GMJHchhOtIToZ77oHWrSvtdMxZMi0jhHAdjz8OUVGwbh1UrWp1NJaSkbsQwjX8/LMpCPboo9C7t9XRWE6SuxCi4jtzxkzHBAVV+umYs2RaRghR8T388L83K/n4WB2NU5CRuxCiYluxAhYuhKlToVs3q6NxGpLchRAVV2ys6YfauTM884zV0TgVSe5CiIpJa5g4EZKSYMmSSlc7pjgy5y6EqJgWLTJTMq+9Zqo+inPIyF0IUfFERpqOSv36mUYc4gLFJnel1HylVKxSaneBbc8ppY4ppXbkPa4v8NxUpVS4Umq/UupaRwUuhKikbDYYM8Y0uF64ENxkjFqYkkzLLATeARaft/0NrfXsghuUUu2BW4FgoCGwWinVRmtts0OsQggBr75qCoItWQLNmlkdjdMq9k+e1notcLqErzcU+Fxrnam1jgDCgcpXSFkI4Rhbt8Kzz8Itt8Dtt1sdjVO7lPczDyilduZN2/jmbWsEHC2wT3TetgsopSYqpbYopbacOnXqEsIQQlQKaWkmodevb6o+KmV1RE6trMn9faAl0AmIAV7L217Yd1sX9gJa67la665a664BAQFlDEMIUWk88ggcOGDm2X19i929sitTctdan9Ra27TWucA8/p16iQaaFNi1MXD80kIUQlR6335rioI99hj07291NBVCmZK7UqpBgU+HA2dX0qwEblVKeSulmgOtgb8vLUQhRKV27JgpCtalC7z4otXRVBjFrpZRSn0GXAn4K6WigenAlUqpTpgpl0hgEoDWOkwptQzYA+QA98tKGSFEmeXmwtixkJEBn34qd6GWQrHJXWs9upDNHxex/0xg5qUEJYQQAMyaBb/+CvPmQdu2VkdTocjqfyGEc/rrL/i//4ORI820jCgVSe5CCOdz5gzceis0agRz58qyxzKQwmFCCOeiNUyaBEePml6otWtbHVGFJMldCOFcPvoIli2D//4XevWyOpoKS6ZlhBDOY+dOU+1x4EB48kmro6nQJLkLIZxDSgqMGmWmYZYskWqPl0imZYQQzuH++015gV9/hXr1rI6mwpM/jUII6y1YAIsXm4qPV11ldTQuQZK7EMJaO3fCfffB1VdLk2s7kuQuhLBOUhLcfLOp8vjpp+DubnVELkPm3IUQ1tAaxo+Hw4fht99knt3OJLkLIazxzjvw5Zfw8svQt6/V0bgcmZYRQpS/DRtM840hQ+Dxx62OxiVJchdClK+TJ00xsKZNzQoZWc/uEDItI4QoPzk5piBYQgKsWiV1YxxIkrsQovxMnQq//w6LFkFIiNXRuDR5PySEKB9ffAGzZ5s17XfeaXU0Lk+SuxDC8Xbtgrvvht694Y03rI6mUpDkLoRwrIQEGD4catUySx+lD2q5kDl3IYTj2Gxw221w5IiZa2/QwOqInMqWyNNU9/GgXf2adn9tGbkLIRzn6afhf/+Dt9+Gyy+3OhqnobVm4foIbp27iZdW7XPIMWTkLoRwjM8+g1dfhXvvNW3zBADpWTaeXr6L5duPMSCoLq+N6uSQ40hyF0LY37Zt5gJqnz7w5ptWR+M0ouJTufeTbew7kcQjA9vwwFWtcHNzTPNvSe5CCPuKiYGhQyEgAL76Si6g5vk57ASPfvkP7m6K+eO6cVXbug49niR3IYT9ZGSYlTGnT8P69VDXsQmsIsix5fLaLwd4//dDXNaoFu/d3oUmdao6/LiS3IUQ9qE1TJgAf/0FX38NnRwzl1yRnEzK4MHPtvN3xGlGd2/C9CHB+HiWT816Se5CCPt45RX45BOYMQNGjLA6GsutD49j8ufbSc208fqoEEZ0aVyux5fkLoS4dF99ZerGjB4N06ZZHY2lbLmat387yJu/HqRlQHU+ndCFNvVqlHscktyFEJdm82YYMwZ69YL580E5ZvVHRRCbnMGUz3ew4VA8wzs34sVhHajmbU2aleQuhCi7I0fgxhuhfn349lvw8bE6IsusO3iKh7/4h5TMbF69uSMjQxujLPxDJ8ldCFE2SUkweDCkpcHq1ZV2ZUy2LZfX81bDtK5bnaXje9C2fvlPw5yv2PIDSqn5SqlYpdTuAtvqKKV+UUodzPvXt8BzU5VS4Uqp/Uqpax0VuBDCQtnZppvS3r1mvj042OqILHH0dBqjPtzI+78fYnT3QFY+cIVTJHYoWW2ZhcCg87Y9BfyqtW4N/Jr3OUqp9sCtQHDe17ynlCqfdT9CiPKhNfznP/Dzz/DhhzBwoNURWWLFjmNc/+Y6wk+m8M5tnXlpxGVU8XKedFfstIzWeq1Sqtl5m4cCV+Z9vAj4HXgyb/vnWutMIEIpFQ50BzbaJ1whhOVeegk+/hj+7/9MiYFKJiUzh2dX7OabbccIberLnFs6lctNSaVV1jn3elrrGACtdYxS6uxkWyNgU4H9ovO2XUApNRGYCBAYGFjGMIQQ5WrRIrPU8Y474IUXrI6m3G2NSuDhL3YQnZDGQ/1b89DVrfBwd87iuva+oFrYpWFd2I5a67nAXICuXbsWuo8Qwon89BOMHw/9+5uReyVa8phjy+WdNeG8/Vs4DWr58MWkXnRrVsfqsIpU1uR+UinVIG/U3gCIzdseDTQpsF9j4PilBCiEcAJbt8JNN5kLp998U6mKgR0+lcIjy/5hx9EzjOjciOeGBlPTx9PqsIpV1vcTK4GxeR+PBVYU2H6rUspbKdUcaA38fWkhCiEsdegQ3HAD+PnBqlVQ0/5dg5yR1polm6K44a0/iYhL5e3RnXn9lk4VIrFDCUbuSqnPMBdP/ZVS0cB04GVgmVLqHuAIMBJAax2mlFoG7AFygPu11jYHxS6EcLQTJ+CaayAnB9asgYYNrY6oXJxIzOCJr3ey9sAp+rT2Z9bNIdSvVbFu0CrJapnRF3mq/0X2nwnMvJSghBBOIDERBg2Ckyfht98gKMjqiBxOa823O44xfUUY2TbNC0ODGdOzqaV3mpaV3KEqhLhQerppuLFnD3z/PXTvbnVEDncqOZP/+3YXP4WdJLSpL7NHhtDcv5rVYZWZJHchxLmys2HUKFi7FpYuNdMyLkxrzcp/jjN9ZRhpWTaeuq4dE/q0wN1B7e/KiyR3IcS/cnNh3DgzWn//fVPC14XFJmXwzIrd/BR2kpAmtXltZEda1XWO8gGXSpK7EMLQGh54AD791NyFeu+9VkfkMFprvtl2jBe+30N6to0nB7VjQp/mTntDUllIchdCmMT+xBNmtP7EE/DUU1ZH5DDRCWlMW76bPw6comtTX165uSMtA6pbHZbdSXIXQsDzz8Ps2XD//fDyy1ZH4xC2XM2iDZHM/nk/AM8Nac+dvZrhVsHn1i9GkrsQld0rr5jkfvfd8NZbLllWYM/xJJ5evosdR89wZdsAZg6/jEa1q1gdlkNJcheiMnvjDTMFM3o0zJ0Lbq4z5wyQnmVjzq8H+GhdBLWrePLmrZ24MaRhhVy3XlqS3IWorN5+Gx55xNSMWbQI3J2nFrk9/LbvJM+uCCM6IZ1bujZh6vXtqF218tTEkeQuRGX03nvw0EMwfDh89hl4Vox6KSURk5jO8yv38L+wE7SqW53PJ/akZws/q8Mqd5Lchahs3n/fXDgdMgQ+/9xlEntWTi4L1kfw5q8HydWaJwa1ZfwVLfDycK2pppKS5C5EZfL222bEPmQIfPmly5Tu3RAex7MrwwiPTWFAUD2mD2VnRzQAABhDSURBVGnvlN2RypMkdyEqizlz4OGHYdgw+OILl0jsx86k898f9vLDrhia1KnCx2O70j+ontVhOQVJ7kJUBrNmmZuTbrrJJebYM7JtzFt7mHd/D0dreHhAGyb1a4GPp2tdFL4UktyFcGVamzXszz8Pt9wCS5ZU6MSutWbVrhP8d9Vejp1J57oO9Zl2QxCNfSv3FExhJLkL4arOlhSYPRvuugvmzavQyx13RScy44c9/B1xmnb1a/DphB5c3tLf6rCcliR3IVyRzWaKgH3wgVkZ89ZbFfYGpROJGcz6aT/fbI+mTlUvZg7vwK3dAit8SV5Hk+QuhKvJyoKxY80yxyefNBUeK+AdmSmZOcz94xDz1kVgy9VM6tuS+65qWWF6mFpNkrsQriQtDW6+GX780RQAe/JJqyMqtWxbLp9vPsqbqw8Ql5LF4I4NeHJQu0q/tLG0JLkL4SpOnzbr1zduNHViJkywOqJSOXuxdPbP+4mIS6V7szp8NDaITk1qWx1ahSTJXQhXcPSoaWYdHg7LlpnRewWhtWZ9eDyzftrHP9GJtKlXnY/u7Er/oLqVosCXo0hyF6Ki27MHrr0WkpLgp5/gyiutjqjEtkYlMPun/Ww8HE/DWj7MHhnC8M6N5GKpHVT45H78TDoNXbwusxAX9fvv5o7TKlVMQ+uQEKsjKpHdxxJ545cD/LovFv/qXjw3pD2jewTi7VFxl2o6mwqd3HdFJzL03T+5oWNDJvdv5TKNbYUokc8+M82sW7Y0F1CbNrU6omLtOZ7EnNUH+HnPSWpV8eTxa9sy7vJmVPOu0KnIKVXo72hj3ypM6teSRRsi+X7ncW4MaciDV0uSFy5Oa7MS5umnoV8/WL4cfH2tjqpIYccTeevXg/wUdpIa3h48PKANd13RTJY1OpDSWlsdA127dtVbtmwp89fHp2Qyb10EizdGkp5t4/oODXjg6lYENahpvyCFcAZZWTBpEixcCLfdBvPng7e31VFd1D9Hz/DOmnB+2XOSGj4e3N27OXf3bk6tqpLU7UEptVVr3bXQ51whuZ91OjWLj/88zKINUaRk5jAgqB73XdWSLoHOPaoRokROn4YRI+CPP+C55+DZZ53y5iStNX9HnOadNeGsOxhHrSqe3NW7GXf1bk6tKpLU7cl1k3tiIowfbyreNWv27+a0bBZsiGDB+kgS07Pp1cKP/1zZkj6t/WVplaiY9u6FG2+EI0fMaP32262O6AK5uZrf9sXy/h+H2BqVgH91L8b3acEdPZtSXebUHcJ1k/vff5slYB4e8NVXZv6xgNTMHD77+whz1x4mNjmT9g1qMqlfC264rAEe7hWzzoaohFatgltvNStili+Hyy+3OqJzZOXksmLHMeatO8yBkyk0ql2FSf1aMKprEynB62Cum9wBDh40I5rwcFMc6T//uWCXzBwbK7Yf54O1hzh8KpVGtatw9xXNuaVbExlRCOeltXlX+tRT0KkTfPstBAZaHVW+xLRsPtt8hAXrIziZlEm7+jWY2LcFQ0Ia4imDp3Lh2skdzPTM7bfDDz+YaZp33in0IlNurmb13pN8tC6CvyNPU8Pbg9E9ArmzV1OpBy2cS2oq3H23udt05EhYsACqVbM6KgCi4lNZsD6SZVuOkpZlo3crPyb2bUlfmfYsdw5L7kqpSCAZsAE5WuuuSqk6wBdAMyASGKW1TijqdexyQdVmg+nTYeZM6NYNvv4amjS56O47jp7ho3WH+XH3CbTWDOpQn3GXN6dbM1/5ARXWCg+H4cPNnacvvQSPP275hdOzJQIWrI/gt/2xeLgpbgxpxD1XNKd9Q1mVZhVHJ/euWuu4AtteBU5rrV9WSj0F+GqtiyxNZ6/VMoB563rnneDjY27y6N+/yN2Pn0lnyaYoPv3rCInp2QQ1qMnYXk0Z2qkRVbxkvlCUsxUrTLleNzdTsveaaywNJyUzh2+2RbNkYxQHY1Pwq+bF7T0CuaNnU+rW9LE0NlH+yX0/cKXWOkYp1QD4XWvdtqjXsWtyB9i/34x89u0z7cWmTSu2UUF6lo0VO46xcEMk+04kU9PHg5tDm3B7z0BaBlS3X2xCFCYnx/ycvvoqdOliFgg0b25ZOPtOJLF00xG+2RZNapaNjo1rcWevZgwJaSAlApyII5N7BJAAaOBDrfVcpdQZrXXtAvskaK0vWGiulJoITAQIDAwMjYqKKnMchUpJgXvvhaVLzYqaTz4B/+Jbcmmt2RyZwJJNUfxvdwzZNk2vFn6M7hHItcH15Adb2N+xY+aGpLVrzc/sG2+Yd57lLCPbxqpdMSz96whboxLw8nBjcMcG3NmrmZTddVKOTO4NtdbHlVJ1gV+AB4GVJUnuBdl95H6W1qau9UMPmcS+dGmpKuadSs5k2ZajfPb3EaIT0vGt6smILo25pVsT2tSTEgfCDlatMtMw6emmJd4dd5R7CGHHE/li81GWbz9GckYOLfyrcVuPQG7q0hjfal7lHo8ouXJZLaOUeg5IASZg9bTM+XbsMJ3fDx6EZ54xD4+SL4HMzdX8GR7H55uP8Muek2TbNCFNajOqa2MGd2wod92J0svKMtMws2dDx45mVUzbIn9N7CohNYtvdxzjyy3R7IlJwsvDjes71GdUtyb0auEniwoqCIckd6VUNcBNa52c9/EvwAtAfyC+wAXVOlrrJ4p6LYcndzDTNA88AIsWwRVXwJIl59zVWlLxKZks336MZVuOcuBkCl4ebgxsX4+buzSmT2t/uTlKFG/fPrN0d9s2c1/Ga6+ZG5QcLDPHxpp9p1i+PZrf9sWSbdN0aFSTkaFNGNqpIbWryii9onFUcm8BLM/71AP4VGs9UynlBywDAoEjwEit9emiXqtckvtZn3xiusEDvPuu+SUrwyhFa82uY4l8vTWalf8cJyEtG79qXgwJacjQTg3p1KS2jH7EubSGefNgyhSoWhU+/hiGDnXoIXNzNVuiEvh2xzFW7YrhTFo2/tW9GdqpITeHNpbiehWc69/EVFqRkWZuc/16GDUK3nsP/PzK/HJZObn8vj+WFTuO88vek2Tl5BJYpypDQhowuGND2tWvIYm+souJMT1Nf/gBBgww7yAbNnTIoc4OPL7fGcMPO2M4diadKp7uDGhfjxFdGtGnlbzDdBWS3AuTkwOvvGKWSvr5mQuvQ4Zc8ssmZWTzv90n+O6f42w4FI8tV9MioBo3XNaA6zo0IKiBJPpKZ9kyM/2SlmaWOt5/f7FLc0tLa83O6ERW7Y7hx10nOHI6DQ83RZ/W/gzt1IiB7etJQwwXJMm9KDt2mNUKO3eaf994w26ND+JSMvlx9wl+3BXDpsPx5Gpo6leVa4Prc21wPTo38cVNekW6rpMnTSL/+mvo3h0WL7brRdMcWy5bohL4OewkP4Wd4NiZdDzcFL1a+jGkY0OuCa4n8+guTpJ7cTIz4YUXzEg+IMBM0wwfbtdDxKVk5v8SbjgUR7ZN41/di/7t6tE/qC5XtPanqpeMrFyC1ubu6IceguRkU3v98cdLtULrYpIzsll3MI7Ve0+yZl8sCWnZeHm40aeVP4M61Gdge0nolYkk95Latg3uuceM5keOhDffhAYN7H6YpIxs1uyLZfXeWH7fF0tyZg5eHm70bOHH1W0DuLJtXZr5O0eRKFFKUVFw331m/XqPHqb2evv2ZX45rTWHTqXy+/5Yft9/ir8i4sm2aWpV8eSqtgFcG1yfvm0CZMqlkpLkXhrZ2abM6gsvmMqSL79s2prZeY70rKycXDZHnua3fbGs2RfL4bhUwEzf9G0dQJ/W/vRs6Se9Jp1dTo4pOf3MM2b11YsvwoMPgnvp72g+k5bFhkPxrDt4inUH44hOSAegdd3qXNWuLv3b1SW0qa9cFBWS3MskPNzcCv7rr2YE9u67EBrq8MNGxqWy9uAp1h44xYZD8aRl2XBTENKkNpe39KNXC39Cm/pKUTNnsn69Ga3v3AmDB5uflVLUXU/NzGFLVAIbDsWxITye3ccT0RpqeHtweSs/+rQO4Mq2AVKWWlxAkntZaW1KFjz2GMTGmhH8zJlQp065HD4rJ5ftRxJYHx7HuvA4dkYnYsvVeLorQhrXplvzOnRvVocuTX3lLlkrnDwJTz5pljU2aQJz5phrNcWshkpIzWJrVAKbI0/zV8Rpdh379/+1cxNfLm/lxxWt/OnUpLaMzkWRJLlfqsREc1Hs7behVi0zZTNpkl0ukJVGSmYOmyNPs+lQPH9FnGb3sURycjVKQZu6NejS1JcugbXpHFibFv7VZSWOo2RmmimYGTMgIwMefRT+7/8KbaZhy9WEx6aw/UgC244ksO3IGcJjUwDy/0j3bOFHjxZ1CG3qKxfVRalIcreXXbtg8mRYs8ZcJHv9dVNx0iJpWTnsOHKGrVEJbD2SwLaoBJIycgDzlv6yxrW4rFEtOjQy/wbWqSoJ/1JobfoFPPGEmbYbPNiUDmjTBjCJPCIulbDjieyKTmRndCK7jyeSlmUDwLeqJ50DfQlt6ku3ZnXo2LiW9BgVl0SSuz1pDStXmtHaoUPmbsNXXjE1uC2Wm6s5HJfCjqOJ7DiawM7oRPbFJJNlywWgmpc7QQ1qEtSgJm3r16Bt/Rq0qVdDpnRKYuNGs5xx/XoICiL5v6+wp+PlHDiZzJ6YZPadSGL/ieT8RO7l4Ub7BjUJaVyLjo3Nu6nm/tXkBjZhV5LcHSEz05RonTED4uNNPe7nn4dWrayO7BxZObkcOJlM2PFE9hxPYk9MEntjkknJzMnfp24Nb1rVrU7LgOq0CKhGc3/zaFS7SqWe883KyeXkxq34vPAcAatXkVzbj89vuIePWl/NyfR/v3+1qnjSrn4NghrUJLhhTTo0qkWrutWlSbRwOEnujpSYaEbuc+aYMq533WWWwzlRl/rzaa05npjBgRPJ7D+ZTHhsCuGxKRyKTSG5QNJ3d1M0ql2FwDpVaexbhUa1q9CwdhUa1PahQa0q1K/pU6FX7aRm5nAiKYMTiRkcP5NOTGIG0QlpHD2dTu7Bg9z64wKGhv1OqpcP87qPYFmfkTRs7E/LgOq0qVeD1vXMvw1q+ciIXFhCknt5iIkxzYw//NBM3dx1F0ydWqaywlbRWhOXkkVkfCoRp1I5cjqNqNNpHIlP5diZDOJSMi/4mho+HgTU8KZuDW/8qnvjV80L36pe1KnmRa0qntSq4knNKp7U8PGgurcH1bw8qOLljpeH/Ua12bZcUjNzSMnMIS3LRlJ6NkkZ2SSmZ3MmLZuE1CwS0rKJT80kLjmLuJRMYpMzz3n3claXzFM8uPEL+m7+hVxPTw7fchfpDz9C41ZN8KvmJUlcOBVJ7uXpyBGT5OfPh9xcGDPGXIBr187qyC5ZRraN42fSOZGYQUxiBieSMjiVnMmp5ExikzOIT80iPiWLxPTsYl/Lw01RxdMdb083vNzd8PJww91N4eHmhpubQvHvikJbriZXa2y5mpxcTXZOLlm2XDKyc0nPtmHLLfpnWCkzdeJXzQv/6t741/AmoLo39Wv5UK+mN/Vq+tDs2CHqvf8m7su+MDev3XuvmWN3wB3KQtiLJHcrREebO13nzjXz80OHmiTfq5fVkTlcti2XxPTs/EdSejYpmTmkZJjRdXqWjfRs88jKyTUPWy45uRqbzSRw05bXvAlyc1O4K4W7m8LTXeHl4YaHuxtVPN3x8XTDx8Odqt4eVPd2p5q3BzV9zLuFmj4e+Fb1omYVT9wLWyWkNfzxh6nU+OOPUL26SeqPPQb16pXvN02IMpDkbqXYWHjnHfNISICePc1yyptuAk9ZpWKJzEz44gtTAXTHDlMsbsoUU5bXThVBhSgPRSV3uZzvaHXrmpuejhwxN77ExcHo0WYu/oUX4NgxqyOsPKKiTN/SwEBT3jk723RGioqCp5+WxC5ciiT38lK9uikktX8/fP89dOgA06dD06YwYoSZFrDZrI7S9WRnmxuPhgyBFi1MIbiePeGnn8xNaePHl0v/UiHKm0zLWOnQITMnP3++GdE3bAh33mkeQUFWR1dxaW2KeC1ZYh6xsebC6F13mbIRTrxMVYjSkDl3Z5eVZUbzCxb8O4Lv3NncGDVqlCSjkgoPhy+/NMXewsJM7Z/Bg02N/kGDyr0WkBCOJsm9Ijlxwlzs+/RT+Ptvs61bN7j5ZlNxsHVra+NzJlqbJL5ypUnqO3aY7b17w+23m4Yr/v7WxiiEA0lyr6jCw03/za++grPfnzZtzGj0+utNEvPxsTbG8paaCmvXwv/+B999BxERZnvPniaZ33yzvNMRlYYkd1cQFWWmbr77zlSlzMoyFwL79DHFy/r1M1M5rra8MisLNm+G3383571undnm42POe8gQ88euYUOrIxWi3ElydzUpKSbZ/fKLeezda7ZXq2ZGsD16mEf37lC/vqWhlorWZmnoli2waZOpxLh5M6SbNnNcdhkMHGjKLPfpI6tcRKUnyd3VnThhRrRr15qStDt3/russn596NQJQkJMDfqgIFMKoUYNa2NOTDTLQsPCzGPXLti+HU6dMs97epp3Ir16mXclffrI/LkQ55HkXtmkpcG2bWbUu2OHeezZY5o4n1W3rln33by5aRHXqJF51K1rkqifn+k65e1d8uNqbToTJSWZMsinTplliDExcPSoeURGmiWgcXH/fp23t/mj06WLSehduphHZbueIEQpSXIX5maeQ4fMFM6+fXD48L+P48fNPHZhvLzMKN/Hx3zs5WUqcWltHtnZJqFnZJjpouyLFA3z9jZ/RAIDTc37s4/gYPNHRpYpClFqRSV3+Y2qLDw9zXRMYdUpc3PNSPvYMTOiPvtISoLkZPNvZqZ5ZGWZpK6UeXh6msTv7W3uwq1Vyzx8fU3xrYAAMzXk719s42ghhP1Ichfg5maScECA1ZEIIexEassIIYQLkuQuhBAuyGHJXSk1SCm1XykVrpR6ylHHEUIIcSGHJHellDvwLnAd0B4YrZRq74hjCSGEuJCjRu7dgXCt9WGtdRbwOTDUQccSQghxHkcl90bA0QKfR+dty6eUmqiU2qKU2nLq7F2JQggh7MJRyb2wBc3n3C2ltZ6rte6qte4aIEvwhBDCrhyV3KOBJgU+bwwcd9CxhBBCnMch5QeUUh7AAaA/cAzYDNymtQ67yP6ngCi7B+J4/kBcsXu5FjnnyqGynXNFPd+mWutCpz4ccoeq1jpHKfUA8BPgDsy/WGLP279CzssopbZcrK6Dq5Jzrhwq2zm74vk6rPyA1noVsMpRry+EEOLi5A5VIYRwQZLcL81cqwOwgJxz5VDZztnlztcp6rkLIYSwLxm5CyGEC5LkLoQQLkiSu50opR5TSmmllMt3cVZKzVJK7VNK7VRKLVdK1bY6JkeobJVNlVJNlFJrlFJ7lVJhSqnJVsdUXpRS7kqp7Uqp762OxV4kuduBUqoJMBA4YnUs5eQXoIPWuiPmZrWpFsdjd5W0smkO8KjWOgjoCdxfCc75rMnAXquDsCdJ7vbxBvAE59XPcVVa65+11jl5n27ClJdwNZWusqnWOkZrvS3v42RMsmtU9FdVfEqpxsANwEdWx2JPktwvkVLqRuCY1vofq2OxyN3Aj1YH4QDFVjZ1ZUqpZkBn4C9rIykXczCDs1yrA7EnaZBdAkqp1UD9Qp6aBjwNXFO+ETleUeestV6Rt880zFv5peUZWzkptrKpq1JKVQe+BqZorZOsjseRlFKDgVit9Val1JVWx2NPktxLQGs9oLDtSqnLgObAP0opMNMT25RS3bXWJ8oxRLu72DmfpZQaCwwG+mvXvFmiUlY2VUp5YhL7Uq31N1bHUw56Azcqpa4HfICaSqlPtNZ3WBzXJZObmOxIKRUJdNVaV8TqciWmlBoEvA7001q7ZKeV0lY2dQXKjFAWAae11lOsjqe85Y3cH9NaD7Y6FnuQOXdRFu8ANYBflFI7lFIfWB2QveVdMD5b2XQvsMyVE3ue3sAY4Oq8/9cdeSNaUQHJyF0IIVyQjNyFEMIFSXIXQggXJMldCCFckCR3IYRwQZLchRDCBUlyF0IIFyTJXQghXJCUHxCiEEqpYOBNIBBYAtQFFmutN1samBAlJDcxCXEepZQPsA0YCRwG9gFbtdYjLA1MiFKQkbsQFxoAbD9bbkAp5QW8Zm1IQpSOzLkLcaHOmJE7SqmGQIrWer21IQlROpLchbhQJv92l3oJ8LIwFiHKRJK7EBf6FOirlNoP/ANsVErNsTgmIUpFLqgKIYQLkpG7EEK4IEnuQgjhgiS5CyGEC5LkLoQQLkiSuxBCuCBJ7kII4YIkuQshhAv6f33CNyMH+dt5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# function decreases fastest in the -∇f(w) direction\n",
    "w = randn(d)\n",
    "v = randn(d)\n",
    "v_normalized = v/norm(v)\n",
    "g = gradf(w)\n",
    "g_normalized = g/norm(g)\n",
    "alphas = np.linspace(-5,5,100)\n",
    "plt.plot(alphas, [f(w + alpha*v_normalized) for alpha in alphas], \n",
    "     label = r\"random direction $f(w + \\alpha v)$\")\n",
    "plt.plot(alphas, [f(w + alpha*g_normalized) for alpha in alphas], color='r', \n",
    "         label= r\"gradient direction $f(w + \\alpha \\nabla f(w))$\")\n",
    "plt.xlabel(r\"$\\alpha$\")\n",
    "plt.legend()\n",
    "# notice: compared to blue line, red line always goes down faster to the left and up faster to the right"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# gradient descent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fa3f7f71a60>]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD7CAYAAABzGc+QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAb/klEQVR4nO3de5ScdZ3n8fe3qrr6fku6kw7pzg1CQghJwFYjIHLVwCA4Z5QBxWXOupsddxSc0ePiuntm5sysxx0dR86qw2YAYdSFYRAVOV5ALiqKQHNPSAIhAdK5defWSTp9q6rv/lHVSafTSbq7qvvp56nP65ymnlvV8/0R+OTXv/o9z2PujoiIhE8s6AJERGR8FOAiIiGlABcRCSkFuIhISCnARURCSgEuIhJSpwxwM7vLzDrMbO2w7Z81s41mts7M/mHiShQRkZGMpgd+N7Bq6AYzuwS4Fljm7mcDXy98aSIicjKJUx3g7r8xs3nDNn8a+Kq79+WO6RjNyRoaGnzevOEfJSIiJ/P888/vdvfG4dtPGeAncCbwfjP7X0Av8AV3f+5Ub5o3bx5tbW3jPKWISHEys7dH2j7eAE8A9cBK4N3A/Wa2wEe4Lt/MVgOrAebMmTPO04mIyHDjnYXSDjzoWc8CGaBhpAPdfY27t7p7a2Pjcb8BiIjIOI03wH8MXApgZmcCSWB3oYoSEZFTO+UQipndC1wMNJhZO/DXwF3AXbmphf3ATSMNn4iIyMQZzSyUG06w68YC1yIiImOgKzFFREJKAS4iElKhCPDH1u/iO09uCroMEZEpJRQB/tSm3Xz7cQW4iMhQoQjwppoyuvvTHOpLBV2KiMiUEYoAn1lTBsDOrt6AKxERmTpCFeC7DijARUQGhSLAm2rVAxcRGS4UAT6zphSAXQcV4CIig0IR4BXJBNVlCXapBy4ickQoAhyyM1F2agxcROSI0AT4zJoydh3oC7oMEZEpI2QBrh64iMig0AR4U20pHQf7SGd011oREQhRgM+sKSOdcfZ0axhFRARCFuAAu7oU4CIiEKIAbxq8nF7j4CIiwCgC3MzuMrOO3OPThu/7gpm5mY34QONCmqkAFxE5xmh64HcDq4ZvNLMW4ArgnQLXNKKGqiQxgw4FuIgIMIoAd/ffAHtH2PVPwBeBSZkWkojHaKwu1f1QRERyxjUGbmbXANvc/eUC13NSuhpTROSoUz6VfjgzqwC+DHxwlMevBlYDzJkzZ6ynO8aMmjLe2XM4r88QEYmK8fTATwfmAy+b2VtAM/CCmTWNdLC7r3H3VndvbWxsHH+lqAcuIjLUmHvg7v4qMGNwPRfire6+u4B1jaiptoyungF6B9KUlcQn+nQiIlPaaKYR3gs8DSwys3Yz+9TElzWyGdW5+4KrFy4icuoeuLvfcIr98wpWzSkMfTLP3OmVk3VaEZEpKTRXYoKuxhQRGSpUAT5DDzcWETkiVAFeU5agvCSuBzuIiBCyADczmmo1lVBEBEIW4JCdiaKHG4uIhDDAm2rL2HVQAS4iEs4A7+ojo0eriUiRC12AN9eV05/O0HlIX2SKSHELX4DXVwDQvk83tRKR4ha6AG+ZVg7A1r09AVciIhKs0AX47Dr1wEVEIIQBXp6M01CVpH2feuAiUtxCF+CQHQffqh64iBS5kAZ4uXrgIlL0QhrgFWzf30Nac8FFpIiFMsBbppUzkHbdlVBEilooA/zoXHANo4hI8RrNI9XuMrMOM1s7ZNvXzGyDmb1iZj8ys7qJLfNYzfXZueCaSigixWw0PfC7gVXDtj0KLHX3ZcDrwJcKXNdJza7TxTwiIqcMcHf/DbB32LZH3D2VW/0D0DwBtZ1QWUmcGdWl6oGLSFErxBj4fwR+XoDPGRNNJRSRYpdXgJvZl4EU8IOTHLPazNrMrK2zszOf0x2jZZou5hGR4jbuADezm4CrgU+4+wknZLv7GndvdffWxsbG8Z7uOM315ezo6iWVzhTsM0VEwmRcAW5mq4D/Blzj7oF0g5vrK0hnXM/HFJGiNZpphPcCTwOLzKzdzD4FfAuoBh41s5fM7PYJrvM4Lbm54JqJIiLFKnGqA9z9hhE23zkBtYzJsXPBpwdbjIhIAEJ5JSbArLoyzHQ1pogUr9AGeGkiTlNNmWaiiEjRCm2Ag+aCi0hxC3mAV7BNAS4iRSrUAd5SX86Orh76U5oLLiLFJ9QBPq+hkozDO3u7gy5FRGTShTrAFzRWAfBmpwJcRIpPyAO8EoAtuxXgIlJ8Qh3gNWUlNFSVsrnzUNCliIhMulAHOMCChko2awhFRIpQ+AO8sZLNGkIRkSIUiQDf293P/sP9QZciIjKpQh/g8xuyM1HUCxeRYhP6AB+ciaJxcBEpNqEP8DnTKkjEjC27NRNFRIpL6AO8JB5jzrQK9cBFpOiEPsAB5msqoYgUodE8Uu0uM+sws7VDtk0zs0fN7I3ca/3ElnlyCxor2bKnm3TmhM9WFhGJnNH0wO8GVg3bdivwmLsvBB7LrQdmQWMV/akM2/fr1rIiUjxOGeDu/htg77DN1wL35JbvAT5S4LrGZEFDbiaKphKKSBEZ7xj4THffAZB7nVG4ksZu/pGphJqJIiLFY8K/xDSz1WbWZmZtnZ2dE3KOxqpSqksTuiuhiBSV8Qb4LjObBZB77TjRge6+xt1b3b21sbFxnKc7OTPL3hNFM1FEpIiMN8AfAm7KLd8E/KQw5YzfgsYqDaGISFEZzTTCe4GngUVm1m5mnwK+ClxhZm8AV+TWA7WgoZLtXb0c7k8FXYqIyKRInOoAd7/hBLsuK3AteTmzqRqA13cdYkVLXcDViIhMvEhciQlwVlMNABt2HAi4EhGRyRGZAG+uL6cyGWfDzoNBlyIiMikiE+CxmLGoqZr16oGLSJGITIADLJ5Vw4adB3HXPVFEJPoiFeBnNVXT1TPAjq7eoEsREZlwkQrwxbNyX2Tu1DCKiERfpAJ8UW4q4fod+iJTRKIvUgFeU1bC7LpyzUQRkaIQqQAHOGtWteaCi0hRiFyAL26qYfPubnoH0kGXIiIyoSIX4GfNqiGdcTZ16MZWIhJtkQvwxbOyX2RqHFxEoi5yAT5veiWliZjGwUUk8iIX4PHcJfXqgYtI1EUuwAEWN1XrYh4RibyIBngNuw/103FAl9SLSHRFMsCXNdcC8NLW/QFXIiIycfIKcDP7SzNbZ2ZrzexeMysrVGH5WDq7lnjMeLldAS4i0TXuADez2cDNQKu7LwXiwPWFKiwfZSVxFjdVqwcuIpGW7xBKAig3swRQAWzPv6TCWNFSxytbu8hkdG9wEYmmcQe4u28Dvg68A+wAutz9kUIVlq/lLXUc7EuxebeuyBSRaMpnCKUeuBaYD5wGVJrZjSMct9rM2sysrbOzc/yVjtG5uSfTv7S1a9LOKSIymfIZQrkc2OLune4+ADwInD/8IHdf4+6t7t7a2NiYx+nGZkFjFVWlCV7WOLiIRFQ+Af4OsNLMKszMgMuA9YUpK3/xmLGsuVZfZIpIZOUzBv4M8ADwAvBq7rPWFKiugljeUsf6HQd0a1kRiaS8ZqG4+1+7+2J3X+run3T3vkIVVggrWupIZZx123VZvYhETySvxBy0IvdFpsbBRSSKIh3gM2vKmFVbpnFwEYmkSAc4wPLmOl1SLyKRFPkAXzGnjrf3HGZvd3/QpYiIFFTkA7x1bj0Az721N+BKREQKK/IBvqy5jrKSGH/YvCfoUkRECiryAZ5MxGidO40/bFYPXESiJfIBDrBywTQ27DzA/sMaBxeR6CiSAJ+OOzyzRb1wEYmOoghwjYOLSBQVRYBrHFxEoqgoAhw0Di4i0VNEAa5xcBGJlqIJcI2Di0jUFE2AaxxcRKKmaAIcjo6D79N9UUQkAooqwM8/owF3eGrT7qBLERHJW14BbmZ1ZvaAmW0ws/Vm9r5CFTYRljfXUV9RwhMbOoIuRUQkb4k8338b8At3/6iZJYGKAtQ0YeIx4wNnNvLk652kM048ZkGXJCIybuPugZtZDXARcCeAu/e7+5R/csIli2ewt7ufV/SQBxEJuXyGUBYAncB3zexFM7vDzCoLVNeE+cCZjcQMDaOISOjlE+AJ4Dzgn939XKAbuHX4QWa22szazKyts7Mzj9MVRl1FknfNrefxjQpwEQm3fAK8HWh392dy6w+QDfRjuPsad29199bGxsY8Tlc4Fy+awdptB+g40Bt0KSIi4zbuAHf3ncBWM1uU23QZ8FpBqppgly6eAcAT6oWLSIjlOw/8s8APzOwVYAXwlfxLmniLm6qZVVvG4xoHF5EQy2saobu/BLQWqJZJY2ZcsngGP3lxG/2pDMlEUV3PJCIRUbTJdemiGXT3p3laN7cSkZAq2gC/cGED1aUJHn55e9CliIiMS9EGeFlJnCuWzOSX63bSn8oEXY6IyJgVbYADfHj5aRzoTfHbN4Kfny4iMlZFHeAXnNFAbXkJD7+yI+hSRETGrKgDPJmIsersJh5Zt5PegXTQ5YiIjElRBzjA1ctn0d2f5kld1CMiIVP0Af6+BdOZXpnkpxpGEZGQKfoAT8RjXHlOE4+t30V3XyrockRERq3oAxzgw8tOo3cgwy/X7Qy6FBGRUVOAA++ZP4150yu479mtQZciIjJqCnCy90b503fP4dm39rKp41DQ5YiIjIoCPOej72omETP+7bl3gi5FRGRUFOA5jdWlXLFkJg88305fSnPCRWTqU4APcf175rDv8ACPrNsVdCkiIqekAB/i/Wc0MLuunPs0jCIiIaAAHyIWM65/dwu/27SHt/d0B12OiMhJ5R3gZhY3sxfN7OFCFBS0697dQknc+O7v3gq6FBGRkypED/wWYH0BPmdKmFlTxjXLZ/Nvz21lX3d/0OWIiJxQXgFuZs3AHwF3FKacqWH1RQvoGUjz/T+8HXQpIiInlG8P/JvAF4ETPtLGzFabWZuZtXV2huPBCYuaqrlkUSN3//4t3WZWRKascQe4mV0NdLj78yc7zt3XuHuru7c2NjaO93STbvVFp7Onu58fvtAedCkiIiPKpwd+AXCNmb0F3AdcambfL0hVU8DKBdNY3lzLHb/dQjrjQZcjInKccQe4u3/J3ZvdfR5wPfC4u99YsMoCZmasvuh0tuzu5hdrdZdCEZl6NA/8JFYtbWLhjCr+8dGNpNJ6cr2ITC0FCXB3f9Ldry7EZ00l8ZjxhQ8tYnNnt8bCRWTKUQ/8FD64ZCYrWur45q/e0IwUEZlSFOCnYGZ8cdUidnT1al64iEwpCvBROP/0Bt6/sIFvP7GJg70DQZcjIgIowEftix9azL7DA3zr8U1BlyIiAijAR+2c5lqua23mzqe2sHHnwaDLERFRgI/FrVeeRVVZgv/547W46+IeEQmWAnwMplUmuXXVYp59ay8/fGFb0OWISJFTgI/Rda0tnDenjq/8bD37D+t2syISHAX4GMVixt9/5By6egb425++FnQ5IlLEFODjsOS0Gj576Rn86MVtPPTy9qDLEZEipQAfp89ccgbnzqnjf/zoVbbv7wm6HBEpQgrwcUrEY3zzT1eQyjifv/9lMrrlrIhMMgV4HuZOr+RvPnw2T2/ew3ee1AU+IjK5FOB5+lhrM9csP41/fPR1Ht+wK+hyRKSIKMDzZGb87z9ZxpJZNdxy70u82Xko6JJEpEgowAugPBlnzX9oJZmI8Z//tY0DuuGViEyCfB5q3GJmT5jZejNbZ2a3FLKwsJldV853PnEe7+w5zJ9/73ndO1xEJlw+PfAU8Hl3PwtYCfyFmS0pTFnh9N4F0/nax5bx+zf3cPO9L+oxbCIyofJ5qPEOd38ht3wQWA/MLlRhYfXH5zbzNx9ewiOv7eLWB1/V9EIRmTCJQnyImc0DzgWeKcTnhd2fXTCfrp4U//Sr1ylNxPi7a5cSi1nQZYlIxOQd4GZWBfwQ+Jy7Hxhh/2pgNcCcOXPyPV1o3HzZGfQMpLn912/S3Zfiax9bTklc3xmLSOHkFeBmVkI2vH/g7g+OdIy7rwHWALS2thbNeIKZceuVi6kpT/APv9jIob4U3/r4eZSVxIMuTUQiIp9ZKAbcCax3928UrqRo+a8Xn8HffWQpj23o4BN3PEPHwd6gSxKRiMjnd/oLgE8Cl5rZS7mfqwpUV6R8cuVcvv3x81i3vYtrv/U7Xm3vCrokEYmAfGahPOXu5u7L3H1F7udnhSwuSq46ZxYP/Pn5GPDR23/Pgy+0B12SiIScvlWbREtn1/LQZy9keUsdf3X/y9xy34u6alNExk0BPskaqkr5f//pvfzVFWfy8Cs7uOq23/Lslr1BlyUiIaQAD0AiHuPmyxZy/395H2Zw3f99mi89+IqesSkiY6IAD9C75tbzy89dxOqLFnB/WzuXf+PXPPB8u67eFJFRUYAHrCKZ4L9fdRYPfeYCZtdX8IV/f5mr/89T/PaNzqBLE5EpTgE+RZx9Wi0/+vT53Hb9Cg70DvDJO5/l4//yB37/5m7c1SMXkePZZIZDa2urt7W1Tdr5wqovleZ7T7/N7b/ezO5DfZw3p45PX3wGly6eQVz3VBEpOmb2vLu3HrddAT519Q6k+fe2rdz+681s299Dc305n1w5l+taW6ivTAZdnohMEgV4iA2kMzz62i7u+f1bPLNlL8l4jMvOmsGfnNfMBxY16iZZIhF3ogAvyO1kZWKVxGNcdc4srjpnFht2HuD+59r5yUvb+PnandRXlPDBJU2sOqeJC05vIJlQmIsUC/XAQ2ogneHXGzv56SvbeWx9B4f6UlSXJrhwYQMXL2rkA2fOoKm2LOgyRaQA1AOPmJJ4jMuXzOTyJTPpS6V56o3d/Gr9Lp7Y0MnP1+4EYEFDJStPn87KBdN519x6TqstI3sTSRGJAvXAI8bd2bjrIL99fTdPb97Ds1v2cqgvBUBTTRnnza1j6exals2uY+nsGuoq9GWoyFSnLzGLVCqd4bUdB3jh7X288M5+Xty6j617e47sb6opY1FTNYubqjl9RhVnzKji9IYqaitKAqxaRIbSEEqRSsRjLGuuY1lzHX92QXbb/sP9rN12gLXbu9i48yAbdh7k6Tf30J/OHHlffUUJc6dXMm96BS3TKmiuL2d2XQWn1ZUxq7ac8qSeLCQSNAV4EaqrSHLhwgYuXNhwZFsqnWHrvh7e7DjEm52HeHvvYd7e081zb+3joZe3M/z2LPUVJcysKWNGTRkzqktprC6lsaqUhupSGiqTTKtKMq0ySX1FUtMcRSaIAlyAbE99fkMl8xsquZyZx+wbSGfY2dXL1n2H2dnVy46uXrbv72HXgT46D/by+s6D7D7UR+oEN+GqLk1QV1lCbXkJdeVJastLqClPUFNWQk15CdVlCapKj/5UliaoKktQmUxQURqnoiROQn8JiBwn34carwJuA+LAHe7+1YJUJVNKSTxGy7TsUMqJZDJOV88Auw/1sae7nz2H+tnb3cfe7gH2He5n/+F+unoG6OoZYEdXDwd6UxzoGaAvlTnhZw6VTMSoSMapTCYoT8YpL8n+lCXjlCVilJXEKSsZfM1uKy2JU5qIUZqIkcz9lCbiJONH15OJ2JH1kniMkriRjOeWEzESsex6TLcwkClo3AFuZnHg28AVQDvwnJk95O6vFao4CY9YzKivTFJfmWThGN7Xl0pzqDfFwdzPob4U3X0puvtTdPelOdyf3dbTn+Zw7qdnILveM5Cmq2eAjoE0vQPZ9b5Uht6BNL0Do/uLYdTts+xvKSUxywV7NuzjMaMkHiMeMxLDlhNxIxE7uh7PbYvZ4HqMeIwjr4lYLLsvd0w8BnEzYjE7+jp02bL/3rPHGjFjyLLl9mW3HX01YrEhy8etgw073nLLNuSz7MhxhnF0vw17nxkYx36uYViMY94Xy01vHTw++zrkPZr+OqJ8euDvATa5+2YAM7sPuBZQgMuolSbilFbFmV5VWtDPdXcG0k5fKhvm/ekM/akMfak0/anB5QwD6QwDaac/tzx4XGpwezq7PZV2BjK519y+9OB6Jrs8kHZS6Qxpz36nkEo7Pek0qUxue8aP/riTSh9dTueOyTjHHJPWveGPMRjm2XDPhj1Hwv748DcYsn/oXxLHfg5D3scJPosh+052jB35x9DPM77yx+fwnvnTCvrvI58Anw1sHbLeDrw3v3JECsPMSCaMZCJGdcgvSM0MCfPM4GsG0p5dz2Q8G/xHlgePzf5Flt1O9ljPbh/6vsFlhyPvdXLvzWRfB49zB+foZ/uQ7cP3H1n3IcfDkWUY+t7sMhz93GwNuWOyO7Lv5dj9TnbFR3jv4LGDBmvwET4nu5TbN2w/Q8+VO9SHfd7gJwy27chpcwuVpYWfuZVPgI/0O81x3QUzWw2sBpgzZ04epxMpTrGYEcMo0cxNGSafr/bbgZYh683A9uEHufsad29199bGxsY8TiciIkPlE+DPAQvNbL6ZJYHrgYcKU5aIiJzKuIdQ3D1lZp8Bfkl2GuFd7r6uYJWJiMhJ5TUP3N1/BvysQLWIiMgY6PI2EZGQUoCLiISUAlxEJKQU4CIiITWpD3Qws07g7XG+vQHYXcBywkBtLg5qc3HIp81z3f24C2kmNcDzYWZtIz2RIsrU5uKgNheHiWizhlBEREJKAS4iElJhCvA1QRcQALW5OKjNxaHgbQ7NGLiIiBwrTD1wEREZIhQBbmarzGyjmW0ys1uDrqfQzKzFzJ4ws/Vmts7Mbsltn2Zmj5rZG7nX+qBrLTQzi5vZi2b2cG490m02szoze8DMNuT+vN9XBG3+y9x/12vN7F4zK4tam83sLjPrMLO1Q7adsI1m9qVcnm00sw+N97xTPsCHPHvzSmAJcIOZLQm2qoJLAZ9397OAlcBf5Np4K/CYuy8EHsutR80twPoh61Fv823AL9x9MbCcbNsj22Yzmw3cDLS6+1Kydy69nui1+W5g1bBtI7Yx9//29cDZufd8J5dzYzblA5whz950935g8NmbkeHuO9z9hdzyQbL/U88m2857cofdA3wkmAonhpk1A38E3DFkc2TbbGY1wEXAnQDu3u/u+4lwm3MSQLmZJYAKsg9+iVSb3f03wN5hm0/UxmuB+9y9z923AJvI5tyYhSHAR3r25uyAaplwZjYPOBd4Bpjp7jsgG/LAjOAqmxDfBL4IDH2EfJTbvADoBL6bGza6w8wqiXCb3X0b8HXgHWAH0OXujxDhNg9xojYWLNPCEOCjevZmFJhZFfBD4HPufiDoeiaSmV0NdLj780HXMokSwHnAP7v7uUA34R86OKncuO+1wHzgNKDSzG4MtqrAFSzTwhDgo3r2ZtiZWQnZ8P6Buz+Y27zLzGbl9s8COoKqbwJcAFxjZm+RHRa71My+T7Tb3A60u/szufUHyAZ6lNt8ObDF3TvdfQB4EDifaLd50InaWLBMC0OAR/7Zm2ZmZMdF17v7N4bsegi4Kbd8E/CTya5torj7l9y92d3nkf0zfdzdbyTabd4JbDWzRblNlwGvEeE2kx06WWlmFbn/zi8j+x1PlNs86ERtfAi43sxKzWw+sBB4dlxncPcp/wNcBbwOvAl8Oeh6JqB9F5L9FeoV4KXcz1XAdLLfXr+Re50WdK0T1P6LgYdzy5FuM7ACaMv9Wf8YqC+CNv8tsAFYC3wPKI1am4F7yY7xD5DtYX/qZG0EvpzLs43AleM9r67EFBEJqTAMoYiIyAgU4CIiIaUAFxEJKQW4iEhIKcBFREJKAS4iElIKcBGRkFKAi4iE1P8H73abcLXc+oYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\"\"\"This function implements gradient descent for the function f with gradient ∇f\"\"\"\n",
    "def gradient_descent(f, gradf,\n",
    "                     alpha = .01,      # small constant step size - try other values!\n",
    "                     w = np.zeros(d),     # start at a random w\n",
    "                     maxiters = 100    # maximum number of iterations \n",
    "                     ):\n",
    "    fks = []  # a list to record all the values f(w) we see\n",
    "    fks.append(f(w)) # record the initial value\n",
    "\n",
    "    # start descending!\n",
    "    for k in range(maxiters):\n",
    "        w = w - alpha*gradf(w)       # take a gradient step\n",
    "        fks.append(f(w))             # record its value\n",
    "    return w, fks\n",
    "\n",
    "w, fks = gradient_descent(f, gradf, alpha = .01)\n",
    "\n",
    "plt.plot(fks, label=\"sum of square errors\")\n",
    "#plt.semilogy(fks, label=\"sum of square errors\")\n",
    "#plt.loglog(fks, label=\"sum of square errors\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Armijo rule"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, '$\\\\alpha$')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEGCAYAAACevtWaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3hUZdrH8e+dXoAUSCAhCQkQejc0USzgCiLqWtm1sCKLorvour6W1V13XV1Rt7kWXFQUK2JBsXekCgSUFjpJILQEAoGEkPq8f8zEHcKkZ+ZMuT/XNddMTplzZzL5zZnnPOc5YoxBKaWUbwmwugCllFKtT8NdKaV8kIa7Ukr5IA13pZTyQRruSinlg4KsLgCgQ4cOJjU11eoylFLKq6xZs+aQMSbO2TyPCPfU1FQyMzOtLkMppbyKiOTWNU+bZZRSygdpuCullA/ScFdKKR+k4a6UUj5Iw10ppXyQhrtSSvkgDXellPJBXh3uB4+d5K8fZXGkpNzqUpRSyqN4dbgfPVHBi0uzeX1lnf34lVLKL3l1uPfs1JZzesTx8vJcTlZUWV2OUkp5DK8Od4Bpo7tyqLiMD37ca3UpSinlMbw+3M/s1p7eCe14fkk21dV6yUCllAIfCHcRYdroNHbkF/PdtgKry1FKKY/g9eEOcPGARDq1C2P24l1Wl6KUUh7BJ8I9ODCAKWelsmLXYTbkFVldjlJKWa7BcBeROSKSLyIbncy7S0SMiHRwmHafiOwQka0icmFrF1yXScNSaBMaxOwluveulFKN2XN/GRhXe6KIJAMXALsdpvUBJgF97es8KyKBrVJpA9qFBXPt8BQ+Xr+PPYUn3LFJpZTyWA2GuzFmMVDoZNa/gLsBxy4qlwLzjDFlxphsYAcwrDUKbYwbR6URGCC8oHvvSik/16w2dxG5BNhrjFlXa1ZnYI/Dz3n2ac6eY5qIZIpIZkFB6/Ry6RQVxmWDOvNW5h4KdUgCpZQfa3K4i0gEcD/wJ2eznUxz2vncGDPbGJNhjMmIi3N6fddmmTa6Kycrqnl1hQ5JoJTyX83Zc+8GpAHrRCQHSALWikgnbHvqyQ7LJgH7WlpkU6R3bMuYXvHMXZFDabkOSaCU8k9NDndjzAZjTLwxJtUYk4ot0IcYYw4AC4FJIhIqImlAOrCqVStuhJvP6UZhSTlvr9nT8MJKKeWDGtMV8k1gBdBTRPJE5Ka6ljXGbALmA1nAZ8Btxhi37z4PTY1hSEo0sxfvorKq2t2bV0opyzWmt8wvjDEJxphgY0ySMebFWvNTjTGHHH5+xBjTzRjT0xjzqSuKboiIMP3c7uQdKeXjDfutKEEppSzlE2eoOjOmVzzp8W2YtWgnxuiAYkop/+Kz4R4QIEw/txtbDhzn2635VpejlFJu5bPhDjBxYCKdo8OZtWin1aUopZRb+XS4BwcG8Ouz01idc4TVOc5OslVKKd/k0+EOcM3QFGIjQ3TvXSnlV3w+3MNDApkyKpVvtuSzaZ8OB6yU8g8+H+4A149MpU1okO69K6X8hl+Ee1R4MNeP7MLHG/azq6DY6nKUUsrl/CLcAaaMSiMkMIDnvtO9d6WU7/ObcI9rG8qkocm8t3Yve4+WWl2OUkq5lN+EO8CvR3cF4Hm9kLZSysf5VbgnxUTw88GdeXPVbvKPn7S6HKWUchm/CneAW8/rTkVVNS8syba6FKWUchm/C/e0DpFMHJjIa9/n6qX4lFI+y+/CHeA353WntKKKOUt1710p5Zv8MtzTO7ZlfL9OzF2eQ1FphdXlKKVUq/PLcAe47bzuHC+rZO7yHKtLUUqpVue34d43MYqxveN5cWk2x0/q3rtSyrf4bbgD3D6mB0WlFbyyItfqUpRSqlX5dbj3T4ri/F7xPL9kF8VllVaXo5RSrcavwx1gxph0jp6o4JUVOVaXopRSrcbvw31QcjTn9Ijj+cW7KNG9d6WUj2gw3EVkjojki8hGh2lPiMgWEVkvIgtEJNph3n0iskNEtorIha4qvDXdPjadIycqePV7bXtXSvmGxuy5vwyMqzXtS6CfMWYAsA24D0BE+gCTgL72dZ4VkcBWq9ZFhqTEMLpHHLN1710p5SMaDHdjzGKgsNa0L4wxNSn4PZBkf3wpMM8YU2aMyQZ2AMNasV6XuWNsOoUl5dpzRinlE1qjzX0K8Kn9cWdgj8O8PPu004jINBHJFJHMgoKCViijZYakxHBuzzhmL96pPWeUUl6vReEuIvcDlcDrNZOcLGacrWuMmW2MyTDGZMTFxbWkjFZzx9geHDlRoWetKqW8XrPDXUQmAxcD1xpjagI8D0h2WCwJ2Nf88txrUHI05/eKZ/biXXrWqlLKqzUr3EVkHHAPcIkx5oTDrIXAJBEJFZE0IB1Y1fIy3ed3Y21nrb60LMfqUpRSqtka0xXyTWAF0FNE8kTkJuBpoC3wpYj8KCLPARhjNgHzgSzgM+A2Y0yVy6p3gf5JUYzt3ZHnl+yi6ITuvSulvJP8r0XFOhkZGSYzM9PqMn6Ste8YF/1nCb89vzu//1lPq8tRSimnRGSNMSbD2Ty/P0PVmT6J7ZgwIIE5S7M5XFxmdTlKKdVkGu51+N3YdEorqvjv4l1Wl6KUUk2m4V6H7vFtuWxwZ+YuzyH/2Emry1FKqSbRcK/H7WPSqao2PP3tDqtLUUqpJtFwr0eX9pFclZHMm6t2s6fwRMMrKKWUh9Bwb8CMMd0REf791XarS1FKqUbTcG9AQlQ4k0d2YcEPeWw/eNzqcpRSqlE03Bth+rndiQgJ4h9fbLO6FKWUahQN90aIjQxh6tlpfLbpAOv2HLW6HKWUapCGeyNNPbsrsZEh/P2LrVaXopRSDdJwb6Q2oUHcem43lmw/xLIdh6wuRyml6qXh3gTXjehC5+hwHvtsC54wJo9SStVFw70JwoID+d0FPVifV8QnGw5YXY5SStVJw72Jfj64Mz07tuWJz7dQUVVtdTlKKeWUhnsTBQYId4/rSc7hE7y1ek/DKyillAU03Jvh/F7xDE2N4cmvt3OiXC+mrZTyPBruzSAi3Du+NwXHy3hhSbbV5Sil1Gk03JvpjC4xjO/Xif9+t5OC43pBD6WUZ9Fwb4G7x/WirLKaJ7/WYQmUUp5Fw70F0jpEcu3wFN5ctYedBcVWl6OUUj/RcG+hGWPSCQ8O5LFPt1hdilJK/aTBcBeROSKSLyIbHabFisiXIrLdfh/jMO8+EdkhIltF5EJXFe4p2rcJZfq53fgi6yCrsgutLkcppYDG7bm/DIyrNe1e4GtjTDrwtf1nRKQPMAnoa1/nWREJbLVqPdSUUWkkRIXx8MdZVFfrsARKKes1GO7GmMVA7V3SS4G59sdzgcscps8zxpQZY7KBHcCwVqrVY4WHBPJ/F/ZkfV4RH6zba3U5SinV7Db3jsaY/QD2+3j79M6A42mbefZppxGRaSKSKSKZBQUFzSzDc1w2qDP9O0fx+GdbKS2vsrocpZSfa+0DquJkmtN2CmPMbGNMhjEmIy4urpXLcL+AAOH+Cb3ZX3SSF5fusrocpZSfa264HxSRBAD7fb59eh6Q7LBcErCv+eV5lxFd2/OzPh2ZtWgn+cdPWl2OUsqPNTfcFwKT7Y8nAx84TJ8kIqEikgakA6taVqJ3ue+i3pRXVfP3z/WKTUop6zSmK+SbwAqgp4jkichNwEzgAhHZDlxg/xljzCZgPpAFfAbcZozxqwbotA6RTB6Zyttr8ti4t8jqcpRSfko84YpCGRkZJjMz0+oyWk1RaQXn/X0R3ePb8Na0EYg4OxShlFItIyJrjDEZzubpGaouEBUezJ0X9GBVdiGfbtQrNiml3E/D3UUmDU2mV6e2/O2TzZys8KuWKaWUB9Bwd5GgwAD+eHEf8o6U8sIS7RqplHIvDXcXGtW9A+P6duKZb3eyv6jU6nKUUn5Ew93F7p/Qm2pjePQTHTVSKeU+Gu4ulhwbwc3ndGPhun06aqRSym003N1g+jndSIwK48GFm6jSUSOVUm6g4e4G4SGB/GFCbzbvP8YbK3OtLkcp5Qc03N1kQv8ERnZtzxOfb+VwsV5QWynlWhrubiIiPHRpX06UV/H4ZzrujFLKtTTc3Si9Y1umnJXGW5l7WLv7iNXlKKV8mIa7m80Yk07HdqH86YONenBVKeUyGu5u1iY0iPsn9GHjXj24qpRyHQ13C0wckMCZ3drz+OdbKTiuB1eVUq1Pw90CtoOr/ThZUcXfPtlsdTlKKR+k4W6R7vFtuHl0Nxb8sJflOw9ZXY5SysdouFvoN+d3Jzk2nD++v5Hyymqry1FK+RANdwuFBQfy0CX92FlQwuzFO60uRynlQzTcLXZer3gu6t+J/3yzg5xDJVaXo5Ryo4Xr9pG175hLnlvD3QM8OLEvoYEB3P/+BjzhmrZKKdfbU3iCe95Zz7+/2uaS59dw9wAd24Vx9/heLNtxmAU/7LW6HKWUixljeOD9jQQI/PmSvi7Zhoa7h7h2WAqDU6J5+OPNFJaUW12OUsqFFq7bx3fbCrjrwp4kRoe7ZBstCncR+Z2IbBKRjSLypoiEiUisiHwpItvt9zGtVawvCwgQHr28P8dKK3jkY+37rpSvOnqinIc+zGJgUhQ3jEx12XaaHe4i0hmYAWQYY/oBgcAk4F7ga2NMOvC1/WfVCL06tePmc7ry7to8Fm8rsLocpZQLPPLxZo6WVvDo5QMIDBCXbaelzTJBQLiIBAERwD7gUmCuff5c4LIWbsOv/Pb8dLrGRfKHBRs4UV5pdTlKqVa0bMch3l6Tx9Sz0+iT2M6l22p2uBtj9gJ/B3YD+4EiY8wXQEdjzH77MvuBeGfri8g0EckUkcyCAt1LrREWHMjMyweQd6SUf3zhmqPoSin3O1Feyb3vrSetQyS/G9vD5dtrSbNMDLa99DQgEYgUkesau74xZrYxJsMYkxEXF9fcMnzSsLRYrh2ewkvLsvlBx31Xyif844tt7CksZebl/QkLDnT59lrSLDMWyDbGFBhjKoD3gDOBgyKSAGC/z295mf7n3vG9iG8bxj3vrqesssrqcpRSLbB29xHmLMvmuhEpDO/a3i3bbEm47wZGiEiEiAgwBtgMLAQm25eZDHzQshL9U9uwYP52eT+2HSzm6W92WF2OUqqZyiqruPfd9SS0C+Oecb3ctt2WtLmvBN4B1gIb7M81G5gJXCAi24EL7D+rZji/V0cuH9yZZxftZOPeIqvLUUo1w3++3s62g8U8cnl/2oYFu227LeotY4x50BjTyxjTzxhzvTGmzBhz2BgzxhiTbr8vbK1i/dGfJvYhNjKE/3tnvY4cqZSXWZ93lOe+28VVZyRxXk+nfUtcRs9Q9XDRESE8fFk/Nu8/xqxFOnKkUt6irLKKu95eR1ybUB64uI/bt6/h7gUu7NuJSwYm8tQ329m0T5tnlPIGT329g20Hi3n08v5EhbuvOaaGhruXeOjSvsREhvD7+eu094xSHm7dnqPM+m4nV56RxHm93NscU0PD3UtER4Qw8/L+bDlwnP98vd3qcpRSdThZUcWd838kvm0of7SgOaaGhrsXGdO7I1dnJDFr0U49uUkpD/XE51vZWVDC41cOsKQ5poaGu5d54OI+dGoXxu/nr6O0XJtnlPIk3+86zJxl2Vw/ogtnp1t75r2Gu5dpFxbM368ayK5DJTz6qQ4NrJSnOH6ygrveXkeX2Ajuu8h9JyvVRcPdC53ZvQM3nZXGKytyWbRVR3dQyhP85cMs9h0t5R9XDyQiJMjqcjTcvdX/XdiTHh3bcPc76zmiV25SylKfbtjPO2vyuO287pzRJdbqcgANd68VFhzIv64ZxJET5XphbaUsdPDYSe5bsIEBSVHMGJNudTk/0XD3Yn0To7jzgp58suEAb6/Js7ocpfxOdbXhrrfXUVZRzb+vGURwoOdEqudUoprl5tFdGdm1PX9euInsQyVWl6OUX3lpeQ5Lth/i/gm96RrXxupyTqHh7uUCAoR/XjOQ4MAA7pj3AxVVOriYUu6wcW8Rj326hQv6dOTa4SlWl3MaDXcfkBAVzmNX9GddXhH/+lIvzaeUq50or2TGvB+IiQzmsSsGYLukhWfx7nDXg4g/GdcvgUlDk5n13U6W7ThkdTlK+bSHPswi+1AJ/7p6ELGRIVaX45T1nTFbYuNGGD0a0tIgNdV27/g4NRUiIy0u0n3+NLEPmblHuOOtH/n09rPp0CbU6pKU8jkfr9/PvNV7mH5uN87s3sHqcurk3eEeEQG//CVkZ8OWLfDpp3Dy5KnLxMWdGvyO4Z+SAmFhVlTuEhEhQTz9y8Fc8vQy7py/jpd/NZSAAM/7uqiUt8o9XMK9765nUHI0d17Qw+py6uXd4d6tGzzzzP9+NgYOHoScHFvg19xnZ8PatbBgAVRUnPociYmn7/XX/JycDMHWDfzTHL06teNPF/fhgfc38vySXdx8TjerS1LKJ5RVVvGbN35ABJ76xWCP6vbojHeHe20i0KmT7TZixOnzq6pg3z5b6Nf+AFi6FObNsy1TIyAAkpJOb+qp+SBITITAQLf8ak1x7fAUlu04xBOfbyUjNcZjzphTyps99ulWNuwt4rnrziA5NsLqchoknnBmY0ZGhsnMzLS6DNte/d69/9vbd9zzz8mxfTA4vl7Bwbamnbra/Dt2tH3gWKCotIKLn1pCZZXh4xlne+xBH6W8wRebDjDt1TX86sxU/nxJX6vL+YmIrDHGZDidp+HeBGVlkJt7eujX/JxfaxCvsLBTm3lqt/3Hxro0/DfkFXHFrOWc2b09cyZr+7tSzbH78AkmPLWE1PaRvDN9JKFBnvNtvb5wb1GzjIhEAy8A/QADTAG2Am8BqUAOcLUxxjeuLBEaCj162G7OlJSc2tzj+HjlSjhS62Vo2/b04He8j4pqUbn9k6L448Q+/PH9jcz6bie3nde9Rc+nlL85WVHF9NfXIMCz1w7xqGBvSEvb3J8EPjPGXCkiIUAE8Afga2PMTBG5F7gXuKeF2/EOkZHQt6/t5kxRkfP2/uxs+PZbOH781OVjYuoP/0Z087xueAqrsgv5xxdbGZwc7dFdt5TyNA99lMWmfcd44YYMr2hnd9TsZhkRaQesA7oahycRka3AucaY/SKSACwyxvSs77m8plnGlYyBwkLnbf01t9LSU9eJi3Pe3JOaCl262L5pAMVllVz2zDKOlJTz0YyzSIgKd+uvppQ3em9tHnfOX8ct53Tj3vHWX3zDGZe0uYvIIGA2kAUMBNYAtwN7jTHRDssdMcbEOFl/GjANICUl5Yzc3Nxm1eE3jLG16dc+2Ftzn5vrvJunPewL4zvzr+3lSNc0HrjtIkLSunhdN0+l3GXjXtvxqkHJ0bw+dThBHtrt0VXhngF8D4wyxqwUkSeBY8BvGxPujnTPvRXUdPN01t6fnQ179kC1w6BiAQG2fvx1ndnrod08lXK1oyfKmfj0UioqDR/NOMujz/R21QHVPCDPGLPS/vM72NrXD4pIgkOzjF4Hzh0CA21hnZxsG5KhtooKyMvj9XmLWLf4R36VCH3K7c1AX3xh+2BwVF83z9RU27kEHjhYklItUVVtmDHvRw4WlfHWzSM8Otgb0uxwN8YcEJE9ItLTGLMVGIOtiSYLmAzMtN9/0CqVqpYJDoa0NK65uwsfd+jNZblHePvmkQxMtn/JOnnyf908a+/1f/ABFBSc+nw13Txrn9Xrpm6eSrnCv77cxuJtBTzy834MTqm3wcHjtaifu73d/QUgBNgF3IhtpMn5QAqwG7jKGFNY3/Nos4x7FZaUM/GppVRVGz787VnEtW3E3kntbp41bf+5ubZ7N3fzVKq1fbx+P7e9sZZrMpKZeUV/jxzGtzY9iUmdZtM+2wGj/p2jeH3qCEKCWnjAqKabZ11n9xYXn7q8s26ejo8jvKvbmfJuWfuOccWs5fRJbMcbvx7uNf3ZNdyVUwvX7WPGmz8waWgyj17uwj0Vx26etQ/41vxcezTP+Pi6w9+hm6dSLeX4TXbhb0cR39Z7Rop12RmqyrtdMjCRLfuP8eyinfTo2JYpZ6W5ZkMi0L697Zbh5H1YM5pnTdg7fgBkZsJ77zU8mqfjh0ByMgTpW1s1rLyymlteW8Oh4jLevmWkVwV7Q/Q/wM/d9bOebM8v5uGPs+gW34ZzesS5vwjH0TxHjjx9fu1uno4fAEuXwptvntrNMzCw7tE8tZunsjPG8MD7G1iVXciTkwYxICm64ZW8iDbLKErKKrli1nL2Hi1lwa1n0j2+rdUlNY29m6fTZp+GRvN09gFg4Wieyn1mL97J3z7Zwozzu3Pnz+o9id5jaZu7alDekRNc9swywkMCWXDrKK/u33uaukbzrHlcXzdPZ+Gv3Ty93ldZB/n1q5lc1C+Bp34x2GtHTNVwV43yw+4jTJr9PX0T2/HGr0cQFuwnTRd1dfOs+Qag3Tx9yoa8Iq7+7wrSO7bhrWkjCQ/x3ve5hrtqtE837Gf662u5eEAC/5nkvXs0rap2N8/aTT/OunnWtdffpYtfXbTd0+w9WsplzywjJDCABbed6fUHULW3jGq08f0TuHd8L2Z+uoWkmAiPHQ3PraKiYOBA2622+rp5btoEH3/cuG6ejuGv3Txdoqi0ghtfWsXJiipenzrc64O9IRru6jQ3j+7KnsITPPfdThKjw7hhZKrVJXmuhrp5VlefetH2xnbzrGsoZ+3m2SxllVVMf20NuwpKeGXKMHp09LJOA82g7xJ1GhHhoUv7cfBYGQ8u3ER82zDG9etkdVneKSAAEhJst8Z083T8EFiyxHk3T8fRPGsf+E1MtG1T/aS62nDX2+tZvvMw/7x6oN9csEbb3FWdSsur+OUL35O17xivTR3O0NRYq0vyPxUVtuGa6zrg62w0zy5dnPf2SUuzNQn5WU+fhz/K4oWl2dwzrhfTz+1mdTmtSg+oqmYrLCnnylnLOVRcxvxbRtKrUzurS1KOHEfzdNbuX7ubZ3j4qaN4+ng3z+cX7+KRTzbzqzNTeXBiH68YDKwpNNxVi+QdOcEVs5ZjDLw7/Uyvu5akX3Ps5unsBK/6unk66/HTzns+3Odn7uHud9YzYUACT/lozy8Nd9Vi2w4e56rnVhAdEcw7t5zZuGGClec7etT5QG51dfOMja2/p4+HdPP8bOMBbn19DaO6d+CFyRleM8pjU2m4q1axJvcI172wki7tI5g3bQTRESFWl6RcyRg4fLjuvf6GRvOs3d6fkuKWbp7LdhzixpdW069zO16bOpyIEN/tN6LhrlrNku0F3PRyJr0T2/H61OG0CfXdfxzVAMduns6Gdti9+9RuniKnd/N0vG+Fbp6ZOYXcMGcVKbERvDVtJFERvn0ReA131aq+zDrILa+tIaNLDC/fOMyrT99WLuTYzdNZT5+8vIZH82xCN8/1eUe59vmVdGgbyls3j/D5k5RAw125wMJ1+7h93g+c1b0Dz9+Q4T/j0KjWU7ubZ+0PgP37T10+JOTUbp4O99siOnDVeztpGx7M/JtHkhgd7v7fxwIa7sol3s7cw93vruecHnH89/ozfPaglbLIyZO2ph3HMX0c2/trdfM8GRxKYFoawd38ZzRPHVtGucRVGclUG8M9725g+mtrmXXdEA141XrCwqBHD9vNmeJistdk8eTzn5N09CBTk4XofHsz0IoVtp5Ajtq2rfuavV7WzbMxNNxVi1wzNIWqavjDgg3c+tpanrl2iDbRKLfYWmz45eIiAnuPtPXeimtz6gLOunlmZ8OuXfD117ZzABzV180zNdXrLtre4mYZEQkEMoG9xpiLRSQWeAtIBXKAq40xR+p+Bm2W8QWvfZ/LA+9v5Ox0bYNXrpe17xjXv7iSwABh3rQRdK0d7A1x7ObpbGiH3NzTu3l27Fj3mD5u6uZZm0vb3EXkTiADaGcP98eBQmPMTBG5F4gxxtxT33NouPuG+av3cM976xmR1p4Xf5Xh0/2LlXV+3HOUyXNWERESyOtThzc92Bujuhry808P/Zqfc3OhsvJ/y7uhm6czLgt3EUkC5gKPAHfaw30rcK4xZr+IJACLjDH1XqBQw913LPghj9/PX8eQlBjm3DiUdmG+3c9YudfqnEJufGk1MZHBvDF1hHVDYVRVwd69dQ/oVl83z9rt/d272y4O3wyuPKD6b+BuwHFw5I7GmP0A9oCPb+E2lBf5+eAkQgIDueOtH/jF7O+ZO2WYb12PVVnmu20F3PLqGhKiw3hj6gg6RVnYjz0w0NYUk5ICo0efPr++bp6ffXZqN88rroB33mn1Epsd7iJyMZBvjFkjIuc2Y/1pwDSAlJSU5pahPNCEAQlEhgZyy2truPq5Fbw6dTid/aTfsXKNhev28fv5P5Ie35a5U4Z5/thGwcHQtavt5kxpqa2bZ06Oy3rpNLtZRkQeBa4HKoEwoB3wHjAUbZZR2L5CT3l5NW1Cg5jrJ1e/Ua3v1RU5/GnhJoamxvLC5Axt6nNQX7NMsy/ZYoy5zxiTZIxJBSYB3xhjrgMWApPti00GPmjuNpR3G5oay1vTRlJVbbhy1nJW5xRaXZLyIsYY/vHFVv74wSbG9IrnlSnDNNibwBXX45oJXCAi24EL7D8rP9UnsR3vTj+TDm1DufaFlXy2cX/DKym/V1FVzV1vr+epb3ZwTUYys647Q7vXNpEOP6DcorCknJvmrubHPUe5/6Le3HRWms9dFUe1juMnK7j19bUs2X6I343twYwx3fW9UgeXNMso1RSxkSG8MXUE4/p24uGPN/PA+xuprKpueEXlV/YU2q76tXznYR6/cgC3j03XYG8mDXflNuEhgTzzyyHcfE5XXl+5mylzMykqrWh4ReUX1uQe4bJnlnGg6CSvTBnG1RnJVpfk1TTclVsFBAj3je/NzMv7s3zHIX7+7DJ2FRQ3vKLyae+tzeMXz39Pm7Ag3rt1FKO6d7C6JK+n4a4sMWlYCq9PHc7RExVc+swyFm3Nt7okZYHKqmoe+jCLO+evY0hKNO/fOoru8S4YTsAPabgrywzv2p6FvxlFUkwEU15ezTPf7qC62voD/Mo9CkvKuWHOKuYsy+bGUam8etNwYiL1urytRcNdWSopJoJ3p49k4sBEnvh8K9Ne1XZ4f7B29xEm/GcJmbbySTUAAA2/SURBVLlHeOLKATw4sS/BgRpHrUlfTWW5iJAg/n3NIP48sQ+LthZwydNL2bi3yOqylAsYY3h5WTbX/HcFQYHCu7ecyVV64NQlNNyVRxARfjUqjXnTRlBWUc3lzy5n7vIcPOE8DNU6ik7Y+q//+cMszukRx0e/OZv+SVFWl+WzNNyVR8lIjeWT28/mrPQOPLhwE7e8toajJ8qtLku10OqcQsY/uZgvsw5y3/hezL4+g6gIHUrAlTTclceJjQzhhRsyeGBCb77Zks+F/17M0u2HrC5LNUNFVTX/+nKbvRkmgHenn8nN53QjIEBPTHI1DXflkQIChKlnd2XBraNoExrEdS+u5KEPszhZUWV1aaqRduQXc8Ws5Tz59XYuHdSZj2ecxcDkaKvL8ht6HTTl0fp1juKj357NzE83M2dZNou25vPEVQM4o0us1aWpOlRVG15ensPjn20hIiSQZ68dwkX9E6wuy+/onrvyeOEhgfzl0n68dtNwyiqrufK5FTz8URal5boX72l25B/nyueW89ePshjVvQOf3zFag90iOiqk8irFZZXM/HQzr32/m6SYcB6+rB/n9tQrOVqtrLKK/363i6e/2UFkaCAPTuzLpYMSddAvF3PZBbJbi4a7aqrvdx3m/gUb2FlQwsSBifxxQm/i21l4TU0/tnT7If74wUayD5UwYUACf7mkr14310003JVP+mlv8dsdBAcIvx2Tzo2jUgkN0os6uMPeo6U8+slmPlq/n9T2Efzl0n6c0yPO6rL8ioa78mm5h0v460dZfLU5n7QOkfzhot6M7R2vTQIucqK8kucW7eS/i3cBMP3cbtxyTje9UpIFNNyVX/h2az5//TCLXYdKGJYWyx8u6s0g7XrXaiqqqnk7M48nv97GwWNlXDIwkXvG96JzdLjVpfktDXflNyqqqpm3eg9PfrWNQ8XljOvbid9d0IOendpaXZrXqq42fLJxP//8Yhu7DpVwRpcY7hvfi4xU7Y5qNQ135XeKyyp5fvEuXlyaTUl5JRcPSGTG+d1J76gh31hV1YZPNuznqW+2s+1gMenxbbh7XC9t8vIgGu7Kbx09Uc7sxbt4eXkOJ8qr+Fmfjkw/txuDU2KsLs1jlVVW8cGP+5i9eBc78ovpFhfJjDHpXDwgkUAdNsCjaLgrv1dYUs7Ly3OYuzyHotIKhqXFMmVUKmN7dyRIxxEH4EhJOfNW7+GlZdnkHy+jV6e23Hpedyb0T9BQ91AuCXcRSQZeAToB1cBsY8yTIhILvAWkAjnA1caYI/U9l4a7cpfiskrmrdrNy8tzyDtSSufocK4b0YUrz0girq3/9c02xrBhbxGvrMjlw3X7KKus5qzuHZg2uitnp3fQ5hcP56pwTwASjDFrRaQtsAa4DPgVUGiMmSki9wIxxph76nsuDXflblXVhq82H+SlZdl8v6uQoABhTO94rhmazNnpcT5/VaDDxWW8/+M+3s7cw5YDx4kICeTngztz/cgu9OrUzuryVCO5pVlGRD4AnrbfzjXG7Ld/ACwyxvSsb10Nd2WlnQXFvLV6D++uyeNwSTmxkSFM6J/AJYMSOSMlxmeGpz1+soIvsw7y0fr9LNleQEWVYUBSFFedkcSlgzvTLkzHV/c2Lg93EUkFFgP9gN3GmGiHeUeMMacdvRKRacA0gJSUlDNyc3NbXIdSLVFeWc132wp4/8e9fJV1kLLKajq0CeWCPvH8rE8nRnRtT3iId52oc6DoJF9vOchXWQdZtvMw5ZXVJEaFMWFAAleckaR76V7OpeEuIm2A74BHjDHvicjRxoS7I91zV57m+MkKvt1awOebDrBoSz4l5VWEBAUwLDWWs9I7MDwtln6dozyu+aaotII1uYUs3X6Y5TsPseXAcQCSY8MZ27sjFw9IYHCy73wb8XcuC3cRCQY+Aj43xvzTPm0r2iyjfEhZZRUrdxWyeFsBi7cXsO1gMQDhwYEMSo5mQHIU/TtH0S8xipTYCLcFZ2l5FVsPHidr3zE27D3KmtwjbM8vxhgIDQpgaGoso7p3YEzveNLj2+jBUR9UX7g3+2IdYnunvAhsrgl2u4XAZGCm/f6D5m5DKU8QGhTI6B5xjLYPipV//CSrs4+wOqeQzNxC5izNpqLKtpMUFhxA1w5t6B7fhuTYcDpHR9A5Jpy4NqG0bxNCbGRIo/b2jTGUlFdx9EQ5h4rLOVBUyv6ik+wpLCX7UDHZh0rYXXiCavu+WbuwIIZ0iWHigETO6BLDkC4xOtaLn2tJb5mzgCXABmxdIQH+AKwE5gMpwG7gKmNMYX3PpXvuypuVVVax7UAxm/YVsSO/mO35xewsKGZ/0Umqqk///woLDiAyJIjwkECCAwMQgQARKquqOVlRTVllFcVllT99YDiKCAkktX0kaXGRdItrQ5+EdvRNbEdSTLjumfshl+y5G2OWAnW9m8Y093mV8jahQYH0T4qif1LUKdMrq6o5eLyMfUdLOXS8jEMl5RQWl1NSXklJWSUnyquorDYYYzAGggKF0KAAQoMCaRMWRHR4MNERwbSPDCUhOoyEqHBiIoI1xFWj6DVUlXKRoMAAOkeH66iJyhKedahfKaVUq9BwV0opH6ThrpRSPkjDXSmlfJCGu1JK+SANd6WU8kEa7kop5YM03JVSygd5xGX2RKQAaMmYvx2AQ61UTmvSuppG62oaratpfLGuLsaYOGczPCLcW0pEMusaX8FKWlfTaF1No3U1jb/Vpc0ySinlgzTclVLKB/lKuM+2uoA6aF1No3U1jdbVNH5Vl0+0uSullDqVr+y5K6WUcqDhrpRSPshjw11EYkXkSxHZbr+PqWO5cSKyVUR2iMi9DtP/LCJ7ReRH++0ih3n32ZffKiIXurmuJ0Rki4isF5EFIhJtn54qIqUO9T7XyHqcbsdhvojIf+zz14vIkEbU2Kjf0RV1iUiyiHwrIptFZJOI3O6wTp1/U1fXZZ+XIyIb7NvOdJhu5evV0+H1+FFEjonIHfZ57ni9eonIChEpE5G7GrOum14vp3V5wPurvterdd9ftkt8ed4NeBy41/74XuAxJ8sEAjuBrkAIsA7oY5/3Z+AuJ+v0sS8XCqTZ1w90Y10/A4Lsjx+rWR9IBTY28TWqczsOy1wEfIrtkogjgJWNqLHB39GFdSUAQ+yP2wLbGvqbuqMu+7wcoENz3hOurKvW8xzAdmKLu16veGAo8Ijjtjzg/VVXXVa/v5zW5Yr3l8fuuQOXAnPtj+cClzlZZhiwwxizyxhTDsyzr9fQ884zxpQZY7KBHfbncUtdxpgvjDGV9uW+B5KasO1Gb6dWva8Ym++BaBFJaGDdxvyOLqnLGLPfGLMWwBhzHNgMdG7i9lu9rgae17LXq9YyY4CdxpiWnO3dpLqMMfnGmNVARRPWdfnrVVddVr+/6nm96tOs18uTw72jMWY/2P4g2D7xausM7HH4OY9T/1C/sX+FnePwVaahddxRV40p2PbGaqSJyA8i8p2InN2IWhqznbqWqW/dxvyOrqrrJyKSCgwGVjpMdvY3dVddBvhCRNaIyDSHZTzi9QImAW/Wmubq16s567rj9WqQRe+v+rTq+8vScBeRr0Rko5NbQ3vfPz2Fk2k1fTtnAd2AQcB+4B+NWMcdddVs436gEnjdPmk/kGKMGQzcCbwhIu1aup16lmnMus3VkrpsM0XaAO8Cdxhjjtkn1/U3dVddo4wxQ4DxwG0iMrqJ23dVXYhICHAJ8LbDfHe8Xq5Y1+XPbeH7qz6t+v4KasnKLWWMGVvXPBE5WPM13f71M9/JYnlAssPPScA++3MfdHiu54GPGlrHHXXZn2MycDEwxtgb0owxZUCZ/fEaEdkJ9AAyqVuDv0s9y4TUs25jfsf6tKQuRCQY2z/e68aY92oWqOdv6pa6jDE19/kisgDb1/DFWPx62Y0H1jq+Rm56vZqzrjterzpZ/P6qU2u/vzy5WWYhMNn+eDLwgZNlVgPpIpJm33OZZF+PWu2RPwc2OjzvJBEJFZE0IB1Y5ca6xgH3AJcYY07UrCAicSISaH/c1V7XrgZqqXM7teq9QWxGAEX2r3b1rduY39EldYmIAC8Cm40x/3RcoZ6/qTvqihSRtvY6IrEdGHd8T1nyejnM/wW1mmTc9Ho1Z113vF5OecD7q666Wv/91ZijrlbcgPbA18B2+32sfXoi8InDchdhO+K9E7jfYfqrwAZgvf3FSXCYd799+a3AeDfXtQNbu9yP9ttz9ulXAJuwHWFfC0xsZD2nbQe4BbjF/liAZ+zzNwAZjajR6e/YxNepWXUBZ2H7Krve4TW6qKG/qRvq6mr/26yz/5084vWyz4sADgNRtZ7THa9XJ2x7rMeAo/bH7Tzg/eW0Lg94f9VVV6u/v3T4AaWU8kGe3CyjlFKqmTTclVLKB2m4K6WUD9JwV0opH6ThrpRSPkjDXSmlfJCGu1JK+SBLhx9QylOJSF/gSSAF28kt8dhGZVxtaWFKNZKexKRULSIShu0s4auwDQGxBVhjjLnc0sKUagLdc1fqdGOBH4wxm+CnERebOkKgUpbSNnelTjcY2547IpIIFBtjlllbklJNo+Gu1OnK+N8Vsh7FNjyyUl5Fw12p070BjBaRrdhG6VshIv+2uCalmkQPqCqllA/SPXellPJBGu5KKeWDNNyVUsoHabgrpZQP0nBXSikfpOGulFI+SMNdKaV80P8Di1WqZb4r/lIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# function decreases fastest in the -∇f(w) direction\n",
    "random.seed(1)\n",
    "\n",
    "w = randn(d)\n",
    "g = gradf(w)\n",
    "c = .1\n",
    "alphas = .05*np.linspace(-1,3,100)\n",
    "plt.plot(alphas, [f(w - alpha*g) for alpha in alphas], \n",
    "     label = r\"$f(w - \\alpha \\nabla f(w))$\")\n",
    "plt.plot(alphas, [f(w) - c*alpha*norm(g)**2 for alpha in alphas],color='r', \n",
    "         label=r\"$f(w) - c\\alpha \\|\\nabla f(w)\\|^2$\")\n",
    "plt.xlabel(r\"$\\alpha$\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gradient descent with linesearch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Poll: gradient descent with linesearch will converge in ____ iterations than without\n",
    "* A. more \n",
    "* B. fewer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fa3f84cbb80>]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAUqElEQVR4nO3da4xcd3nH8e8zl72EJNhJNqlxMKZRSrkU7HSb0ga1QEgbAmqSSkhBgvoFknkBakCoKLQvCm8qKq6tVJAMpFiUpkWFNmmUXlI3aYpUhW4SExwMTdqG4GDsDZAmIcna3n36Ys6sJ86udzw7s+tzzvcjrWbmzJmd5y/HP//zzP/MPzITSVL5NNa7AEnSYAxwSSopA1ySSsoAl6SSMsAlqaRaa/lm5513Xm7dunUt31KSSu+ee+55LDOnTjy+pgG+detWZmZm1vItJan0IuJ7Sx23hSJJJWWAS1JJGeCSVFIGuCSVVN8BHhHNiLgvIm4tHp8TEbdHxIPF7cbRlSlJOtGpzMCvB/b3PL4B2JOZFwN7iseSpDXSV4BHxIXAW4DP9xy+Gthd3N8NXDPc0iRJJ9PvDPzTwAeBhZ5jF2TmQYDi9vylXhgROyNiJiJmZmdnBypyz/5DfObOhwZ6rSRV1YoBHhFvBQ5n5j2DvEFm7srM6cycnpp63oVEffm3/5pl113/M9BrJamq+rkS8zLgtyLiKmACODsi/gI4FBGbMvNgRGwCDo+qyMl2k2eOzI/q10tSKa04A8/MD2XmhZm5FbgO+NfMfAdwC7CjOG0HcPOoihxvN5k7toC7B0nScatZB/5R4IqIeBC4ong8EhPtTplzxxZWOFOS6uOUvswqM+8E7izu/wi4fPglPd9kuwnAM0fmmSjuS1LdleJKzG5oP3vMPrgkdZUkwDtlPnvUFookdZUjwFvHWyiSpI5yBPiYLRRJOlE5AryYgT971ACXpK5yBHh3GaE9cElaVIoAnyxaKM84A5ekRaUIcFsokvR85Qjw7jpwWyiStKgUAb54JaYzcElaVIoAH1+8kMcAl6SucgR4q0EEzBngkrSoFAEeEUy0mrZQJKlHKQIcOmvB/RBTko4rUYA37YFLUo9SBbgtFEk6rp9NjSci4hsR8c2IeCAiPlIc/3BEPBoRe4ufq0ZZaGcGbgtFkrr62ZFnDnhjZj4VEW3g6xHxD8Vzn8rMj4+uvOMm2g3m/DZCSVq0YoBnZyfhp4qH7eJnzXcXnmjZA5ekXn31wCOiGRF7gcPA7Zl5d/HUeyPi/oi4MSI2LvPanRExExEzs7OzAxc6OWYPXJJ69RXgmTmfmduAC4FLI+JVwGeBi4BtwEHgE8u8dldmTmfm9NTU1MCFuoxQkp7rlFahZObjdHalvzIzDxXBvgB8Drh0BPUtsoUiSc/VzyqUqYjYUNyfBN4EfCciNvWcdi2wbzQldkyMGeCS1KufVSibgN0R0aQT+F/JzFsj4ksRsY3OB5oPA+8eXZndGbgtFEnq6mcVyv3A9iWOv3MkFS2j0wN3Bi5JXaW5EnOy3eTYQnJ03lm4JEGJAvz4rjzOwiUJShXg3U0dnIFLEpQowMedgUvSc5QmwCcNcEl6jtIEuDvTS9JzlSjAix6430goSUCJArzbQnnmiAEuSVCiAHcZoSQ9V4kCvNtCsQcuSVCqAC9m4LZQJAkoY4D7IaYkAWUMcHvgkgSUKcBbnVKfOWIPXJKgRAHeajZoN8MWiiQVShPg4LZqktSrny3VJiLiGxHxzYh4ICI+Uhw/JyJuj4gHi9sld6UfpvG2u/JIUlc/M/A54I2Z+Ro6O9BfGRGvBW4A9mTmxcCe4vFITY65K48kda0Y4NnxVPGwXfwkcDWwuzi+G7hmJBX2sIUiScf11QOPiGZE7AUOA7dn5t3ABZl5EKC4PX+Z1+6MiJmImJmdnV1VsRNtA1ySuvoK8Mycz8xtwIXApRHxqn7fIDN3ZeZ0Zk5PTU0NWifQ+UKrZwxwSQJOcRVKZj4O3AlcCRyKiE0Axe3hoVd3gvF2ww8xJanQzyqUqYjYUNyfBN4EfAe4BdhRnLYDuHlURXbZQpGk41p9nLMJ2B0RTTqB/5XMvDUi/gP4SkS8C3gEeNsI6wQ6LRQDXJI6VgzwzLwf2L7E8R8Bl4+iqOVM2EKRpEXluhKz3fRSekkqlC/AbaFIElDKAF8gM9e7FEladyUL8E65c26rJkklC/CWmzpIUlepAnxyrBPgXo0pSSUL8MWd6V1KKEklC3BbKJK0qFwBbgtFkhaVK8CdgUvSonIFeHcZoT1wSSpXgLsKRZKOK1WA20KRpOPKFeDtboDbQpGkkgV4dx24M3BJKlmA2wOXpK5+tlR7cUTcERH7I+KBiLi+OP7hiHg0IvYWP1eNutjxVoMImDPAJamvLdWOAR/IzHsj4izgnoi4vXjuU5n58dGV91wRwXirwbN+G6Ek9bWl2kHgYHH/yYjYD2wedWHLmWw3eeaIM3BJOqUeeERspbM/5t3FofdGxP0RcWNEbFzmNTsjYiYiZmZnZ1dVLLgrjyR19R3gEXEm8FXgfZn5BPBZ4CJgG50Z+ieWel1m7srM6cycnpqaWnXBnX0xbaFIUl8BHhFtOuH95cz8GkBmHsrM+cxcAD4HXDq6Mo+bsIUiSUB/q1AC+AKwPzM/2XN8U89p1wL7hl/e8020G8y5M70k9bUK5TLgncC3ImJvcez3gbdHxDYggYeBd4+kwhNMtOyBSxL0twrl60As8dRtwy9nZZNjTQ4/eXQ93lqSTiuluhITOi0UvwtFksoY4LZQJAkoYYCPt5vOwCWJ/j7EPK1Mtpv8+KdzXPUn/w7AReefyZ9et43OYhlJqo/SzcDf8upNvOnlF/CiDZPMHZvn77/5A44t5HqXJUlrrnQB/osv2ciu35nm8zumue6XtgAw55WZkmqodAHea3xxk2M/1JRUP+UO8FYR4M7AJdVQqQN8zACXVGOlDvDxYpf6Iwa4pBoqeYB3Z+D2wCXVT8kDvDMDt4UiqY7KHeCLq1AMcEn1U+4At4UiqcZKHuC2UCTVVz878rw4Iu6IiP0R8UBEXF8cPycibo+IB4vbJTc1HqUxZ+CSaqyfGfgx4AOZ+XLgtcB7IuIVwA3Ansy8GNhTPF5T3RaKywgl1dGKAZ6ZBzPz3uL+k8B+YDNwNbC7OG03cM2oilyOV2JKqrNT6oFHxFZgO3A3cEFmHoROyAPnD7u4lYy3ix64q1Ak1VDfAR4RZwJfBd6XmU+cwut2RsRMRMzMzs4OUuOyXIUiqc76CvCIaNMJ7y9n5teKw4ciYlPx/Cbg8FKvzcxdmTmdmdNTU1PDqHlRqxE0whaKpHrqZxVKAF8A9mfmJ3ueugXYUdzfAdw8/PJWrI3xVtMAl1RL/WypdhnwTuBbEbG3OPb7wEeBr0TEu4BHgLeNpsSTG2s1/D5wSbW0YoBn5teB5TacvHy45Zy68VbDGbikWir1lZjQ+T4U14FLqqPyB7g9cEk1VYEAb7iMUFItVSTAnYFLqp8KBHjTKzEl1VLpA3zMFoqkmip9gNtCkVRX5Q/wdtNlhJJqqfwB7gxcUk1VJMDtgUuqnwoEuKtQJNVT+QO8bQtFUj2VPsDHmg2OzC+wsJDrXYokranSB/h4u9jYeN5ZuKR6KX+At4p9MW2jSKqZCgS4+2JKqqd+tlS7MSIOR8S+nmMfjohHI2Jv8XPVaMtc3mKAuxJFUs30MwP/InDlEsc/lZnbip/bhltW/8bbtlAk1dOKAZ6ZdwE/XoNaBmILRVJdraYH/t6IuL9osWwcWkWn6HiAOwOXVC+DBvhngYuAbcBB4BPLnRgROyNiJiJmZmdnB3y75Y3ZA5dUUwMFeGYeysz5zFwAPgdcepJzd2XmdGZOT01NDVrnsrrLCF0HLqluBgrwiNjU8/BaYN9y547a8VUo9sAl1UtrpRMi4ibg9cB5EXEA+EPg9RGxDUjgYeDdI6zxpCba9sAl1dOKAZ6Zb1/i8BdGUMtAvBJTUl15JaYklVQFAryYgbsKRVLNlD7Ax1wHLqmmKhTgtlAk1UvpA7zZCNrNcGd6SbVT+gCHYl9MA1xSzVQkwN2ZXlL9VCfAXYUiqWaqEeBtWyiS6qcSAT7WtIUiqX4qEeDj7YYzcEm1U40AbzVcRiipdioS4PbAJdVPRQLcHrik+qlGgLddRiipfqoR4LZQJNXQigFe7Dp/OCL29Rw7JyJuj4gHi9t125UeXEYoqZ76mYF/EbjyhGM3AHsy82JgT/F43biMUFIdrRjgmXkX8OMTDl8N7C7u7wauGXJdp8RlhJLqaNAe+AWZeRCguD1/uRMjYmdEzETEzOzs7IBvd3L2wCXV0cg/xMzMXZk5nZnTU1NTI3mP8VaD+YXk2LwhLqk+Bg3wQxGxCaC4PTy8kk7deNtt1STVz6ABfguwo7i/A7h5OOUMZnFjYwNcUo30s4zwJuA/gJdFxIGIeBfwUeCKiHgQuKJ4vG7G3RdTUg21VjohM9++zFOXD7mWgS1ubOzVmJJqpDJXYgIc8UNMSTVSkQB3Bi6pfqoR4G174JLqpxoB7ioUSTVUkQB3Bi6pfqoR4G174JLqpxIBPtb0SkxJ9VOJAB9vF8sIDXBJNVKNALcHLqmGKhbgzsAl1UdFAtxlhJLqpxIB3m4GETB31BaKpPqoRIBHBOMt98WUVC+VCHDo7kxvgEuqj8oE+Hi76SoUSbVSnQC3hSKpZlbc0OFkIuJh4ElgHjiWmdPDKGoQBrikullVgBfekJmPDeH3rMp4q+l3oUiqleq0UNoNe+CSamW1AZ7AP0fEPRGxc6kTImJnRMxExMzs7Owq3255tlAk1c1qA/yyzLwEeDPwnoj4tRNPyMxdmTmdmdNTU1OrfLvljbWaBrikWllVgGfmD4rbw8DfApcOo6hBjLcaXokpqVYGDvCIeEFEnNW9D/wGsG9YhZ2q8VbDXekl1cpqVqFcAPxtRHR/z19m5j8OpaoBuApFUt0MHOCZ+T/Aa4ZYy6p0VqEY4JLqozrLCFsuI5RULxUKcFehSKqXCgV4gyPHFsjM9S5FktbEMC6lPy2MFduqvf+v99JoBOefNcEHf/NlNBqxzpVJ0mhUJsAv2bKRreeewcz3fsKzRxd47Kk53vrqTbxq8wvXuzRJGonKBPivXHQud/7eGwA48JOned0f38F9j/zEAJdUWZXpgffavGGSqbPGufeRx9e7FEkamUoGeERwyZYN3PfIT9a7FEkamUoGOMD2LRt5+EdP86On5ta7FEkaicoG+CVbNgKw9/u2USRVU2UD/Bc2v5BWI7jXNoqkiqpsgE+ONXn5prO593vOwCVVU2UDHGD7lg1888DjzC94daak6ql0gF+yZSNPH5nnuz98cr1LkaShq3SAb9+yAYD7vm8fXFL1VDrAt5xzBue+YMw+uKRKWlWAR8SVEfHdiHgoIm4YVlHDEhFs94IeSRU18HehREQT+DPgCuAA8J8RcUtmfntYxQ3D9i0b+Zf9h/n1j90BQLvZ4OcuOJNXvuiF/PzPnMV4qwlAI+DcM8f5mbMnOHuyRbFVnCSdtlbzZVaXAg8VW6sREX8FXA2cVgF+7fbNPPzYTzlabHj89JF59j36BLd964fLvmas1aDd8zW0jQgajaDVCHpzvRFBu9mg1QwaywT+845G59ha/APhP0HS6eOPfvsX+KWt5wz1d64mwDcD3+95fAD45RNPioidwE6ALVu2rOLtBvOiDZN87G3P37rz/54+ykOzT7FQbABxbD557Kk5Dj3xLLNPzi0uPUxgIZOFheRYz3LEBBYWkqPzybGFBZbaR+LEQ5nZObYGqxpzLd5EUt8m282h/87VBPhSE7znpUZm7gJ2AUxPT582qfLCM9r84ks2rncZkjSw1XyIeQB4cc/jC4EfrK4cSVK/VhPg/wlcHBEvjYgx4DrgluGUJUlaycAtlMw8FhHvBf4JaAI3ZuYDQ6tMknRSq9pSLTNvA24bUi2SpFNQ6SsxJanKDHBJKikDXJJKygCXpJKKXOoSwlG9WcQs8L0BX34e8NgQyykDx1wPjrkeVjPml2Tm1IkH1zTAVyMiZjJzer3rWEuOuR4ccz2MYsy2UCSppAxwSSqpMgX4rvUuYB045npwzPUw9DGXpgcuSXquMs3AJUk9DHBJKqlSBPjpvnnyakXEiyPijojYHxEPRMT1xfFzIuL2iHiwuK3cDhQR0YyI+yLi1uJxpcccERsi4m8i4jvFn/ev1GDM7y/+u94XETdFxETVxhwRN0bE4YjY13Ns2TFGxIeKPPtuRPzmoO972gd4z+bJbwZeAbw9Il6xvlUN3THgA5n5cuC1wHuKMd4A7MnMi4E9xeOquR7Y3/O46mP+E+AfM/PngdfQGXtlxxwRm4HfBaYz81V0vnr6Oqo35i8CV55wbMkxFn+3rwNeWbzmM0XOnbLTPsDp2Tw5M48A3c2TKyMzD2bmvcX9J+n8pd5MZ5y7i9N2A9esT4WjEREXAm8BPt9zuLJjjoizgV8DvgCQmUcy83EqPOZCC5iMiBZwBp2duyo15sy8C/jxCYeXG+PVwF9l5lxm/i/wEJ2cO2VlCPClNk/evE61jFxEbAW2A3cDF2TmQeiEPHD++lU2Ep8GPggs9Byr8ph/FpgF/rxoG30+Il5AhcecmY8CHwceAQ4C/5eZ/0yFx9xjuTEOLdPKEOB9bZ5cBRFxJvBV4H2Z+cR61zNKEfFW4HBm3rPetayhFnAJ8NnM3A78lPK3Dk6q6PteDbwUeBHwgoh4x/pWte6GlmllCPBabJ4cEW064f3lzPxacfhQRGwqnt8EHF6v+kbgMuC3IuJhOm2xN0bEX1DtMR8ADmTm3cXjv6ET6FUe85uA/83M2cw8CnwN+FWqPeau5cY4tEwrQ4BXfvPkiAg6fdH9mfnJnqduAXYU93cAN691baOSmR/KzAszcyudP9N/zcx3UO0x/xD4fkS8rDh0OfBtKjxmOq2T10bEGcV/55fT+YynymPuWm6MtwDXRcR4RLwUuBj4xkDvkJmn/Q9wFfBfwH8Df7De9YxgfK+j879Q9wN7i5+rgHPpfHr9YHF7znrXOqLxvx64tbhf6TED24CZ4s/674CNNRjzR4DvAPuALwHjVRszcBOdHv9ROjPsd51sjMAfFHn2XeDNg76vl9JLUkmVoYUiSVqCAS5JJWWAS1JJGeCSVFIGuCSVlAEuSSVlgEtSSf0/YYaLHyD89u0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\"\"\"This function implements gradient descent with linesearch for the function f with gradient ∇f\"\"\"\n",
    "def gradient_descent_linesearch(f, gradf,\n",
    "                                alpha = 1,      # initial step size - try other values!\n",
    "                                w = randn(d),     # start at a random w\n",
    "                                maxiters = 100,   # maximum number of iterations \n",
    "                                c = .01            # Armijo sufficient decrease\n",
    "                                ):\n",
    "\n",
    "    fks = []  # a list to record all the values f(w) we see\n",
    "    fw = f(w)\n",
    "    fks.append(fw) # record the initial value\n",
    "    g = gradf(w)\n",
    "    alpha0 = alpha\n",
    "    \n",
    "    # start descending!\n",
    "    for k in range(maxiters):\n",
    "        w_new = w - alpha*g      # take a gradient step\n",
    "        fw_new = f(w_new)\n",
    "        while fw_new > fw - c*alpha*norm(g)**2:\n",
    "            alpha /= 2\n",
    "            w_new = w - alpha*g      # take a gradient step\n",
    "            fw_new = f(w_new)\n",
    "        w, fw, g = w_new, fw_new, gradf(w)\n",
    "        alpha = alpha0\n",
    "        fks.append(fw)\n",
    "    return w, fks\n",
    "\n",
    "w, fks = gradient_descent_linesearch(f, gradf, alpha = 1, c = .01)\n",
    "\n",
    "plt.plot(fks, label=\"sum of square errors\")\n",
    "# plt.semilogy(fks, label=\"sum of square errors\")\n",
    "# plt.loglog(fks, label=\"sum of square errors\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stochastic gradient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "This function implements the stochastic gradient method\n",
    "for the function sum_i f(x_i, y_i; w), \n",
    "where f has gradient wrt w gradf(X,y,w)\n",
    "\"\"\"\n",
    "def sgd(gradf, \n",
    "        data,                # an iterable [(X_1, y_1), ...]\n",
    "        w,                   # initial guess for w\n",
    "        alpha = .01,         # small constant step size - try other values!\n",
    "        maxiters = 100,      # maximum number of iterations \n",
    "        f_validation = None, # evaluate w on a validation set\n",
    "        validation_freq = 1  # compute value on validation set every validation_freq iterations\n",
    "        ):\n",
    "    fks = []  # a list to record the validation values f(w) we see\n",
    "\n",
    "    # start descending!\n",
    "    k = 0\n",
    "    for X,y in data:\n",
    "        \n",
    "        # stopping condition\n",
    "        if k > maxiters:\n",
    "            break\n",
    "        \n",
    "        # logging: record value on validation set\n",
    "        if (f_validation is not None) & k%validation_freq==0:\n",
    "            fks.append(f_validation(w))\n",
    "\n",
    "        # iterate\n",
    "        w = w - alpha*gradf(X,y,w)       # take a gradient step\n",
    "        k += 1                           # increment k\n",
    "    return w, fks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Apply SGD to a random synthetic dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fa3f8a10790>]"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD4CAYAAADsKpHdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de5RV9X338ffnnLkww3WAGcSZQVBHjRBFIQRjtVqMEttG08e02DbaxDw0Vle1aVcbV9Z6kj59XKtJmtq6Em2NGjVNROOlkDQarZrYpEYdlSgg6CAqAwgjIKBcZ+b7/HH2oYdhLjAXznD257XWWbPPd+99zu/H5Xzmt39776OIwMzMLFPsBpiZ2fDgQDAzM8CBYGZmCQeCmZkBDgQzM0uUFbsB/TVx4sSYOnVqsZthZnZUeeGFF96NiNru1h21gTB16lSam5uL3Qwzs6OKpLd6WudDRmZmBjgQzMws4UAwMzPgEAJB0p2SNklaVlC7T9LS5PGmpKVJfaqkXQXr/qVgn1mSXpHUIulmSUrqlcnrtUh6VtLUwe+mmZn15VBGCHcB8wsLEfEHETEzImYCDwIPFaxenV8XEV8oqN8KLASakkf+Na8CtkbEicBNwNf61RMzMxuQPgMhIp4GtnS3Lvkt//eBe3t7DUmTgTER8Uzk7qZ3D3BpsvoS4O5k+QFgXn70YGZmR85A5xDOATZGxOsFtWmSXpL0c0nnJLV6oLVgm9akll+3FiAi2oFtwITu3kzSQknNkprb2toG2HQzMys00EC4nANHBxuAKRFxBvBF4AeSxgDd/cafv+92b+sOLEbcFhGzI2J2bW2311X06fk3t/D1R1fS2enbfpuZFep3IEgqA34PuC9fi4g9EbE5WX4BWA2cRG5E0FCwewOwPlluBRoLXnMsPRyiGgy/Xvset/xsNTv2tA/VW5iZHZUGMkK4AFgZEfsPBUmqlZRNlo8nN3n8RkRsAHZImpvMD1wBLE52WwJcmSxfBjwZQ/itPeOqKwB4b+feoXoLM7Oj0qGcdnov8AxwsqRWSVclqxZw8GTyucDLkn5NboL4CxGR/23/auB2oIXcyOGRpH4HMEFSC7nDTF8aQH/6VFNdDsDWnfuG8m3MzI46fd7LKCIu76H+J93UHiR3Gmp32zcDM7qp7wY+3Vc7Bkt+hLDVIwQzswOk7krl/AjBh4zMzA6UwkBIRggf+JCRmVmh1AXCmKpyJI8QzMy6Sl0gZDNibFW5J5XNzLpIXSBA7rCRJ5XNzA6UykAYV13Oex4hmJkdIJWB4BGCmdnB0hsIHzgQzMwKpTQQPKlsZtZVOgNhZAW79nWwe19HsZtiZjZspDIQxu2/WtmjBDOzvFQGQo3vZ2RmdpBUBsK4/Xc8dSCYmeWlMhBq9n8ngg8ZmZnlpToQPEIwM/sfqQwETyqbmR0slYEwojxLVXnWF6eZmRVIZSCAL04zM+sqtYEwrrrC34lgZlYgtYFQM7Lck8pmZgX6DARJd0raJGlZQe2rktZJWpo8Li5Yd4OkFkmrJF1UUJ8l6ZVk3c2SlNQrJd2X1J+VNHVwu9i93AjBh4zMzPIOZYRwFzC/m/pNETEzefwEQNKpwAJgerLPLZKyyfa3AguBpuSRf82rgK0RcSJwE/C1fvblsOTmEDxCMDPL6zMQIuJpYMshvt4lwKKI2BMRa4AWYI6kycCYiHgmIgK4B7i0YJ+7k+UHgHn50cNQqqmuYNuufXR2xlC/lZnZUWEgcwjXSno5OaRUk9TqgbUF27QmtfpkuWv9gH0ioh3YBkzo7g0lLZTULKm5ra1tAE3PHTLqDNi+24eNzMyg/4FwK3ACMBPYAHwzqXf3m330Uu9tn4OLEbdFxOyImF1bW3t4Le6iZv/9jBwIZmbQz0CIiI0R0RERncB3gDnJqlagsWDTBmB9Um/opn7APpLKgLEc+iGqfvPtK8zMDtSvQEjmBPI+BeTPQFoCLEjOHJpGbvL4uYjYAOyQNDeZH7gCWFywz5XJ8mXAk8k8w5D6n9tXOBDMzADK+tpA0r3AecBESa3AV4DzJM0kd2jnTeBPASJiuaT7gRVAO3BNROS/luxqcmcsVQGPJA+AO4DvSWohNzJYMBgd68v+EcIHPmRkZgaHEAgRcXk35Tt62f5G4MZu6s3AjG7qu4FP99WOweZDRmZmB0rtlcqjR5SRke94amaWl9pAyGTEuOoKjxDMzBKpDQTITSx7hGBmlpPqQKjxCMHMbD8HgkcIZmZA6gOh3NchmJkl0h0II33IyMwsL9WBMK66nN37Otm1t6Pvjc3MSlyqA2F8cnHa5g/2FLklZmbFl+pAOHZcFQDrtu4qckvMzIov1YEwZXw1AGsdCGZm6Q6EY8dVkRG8vWVnsZtiZlZ0qQ6EirIMk8dWsdaBYGaW7kAAaBzvQDAzAwcCU8ZX+5CRmRkOBBprqtm0Y4+vRTCz1Et9IEyZkDvTqHWrRwlmlm6pD4TG/aeeOhDMLN0cCDW5QHh7swPBzNKtz0CQdKekTZKWFdS+IWmlpJclPSxpXFKfKmmXpKXJ418K9pkl6RVJLZJulqSkXinpvqT+rKSpg9/Nnk0cVUFVeZa3t/jiNDNLt0MZIdwFzO9SexyYERGnAa8BNxSsWx0RM5PHFwrqtwILgabkkX/Nq4CtEXEicBPwtcPuxQBIYsr4ah8yMrPU6zMQIuJpYEuX2mMR0Z48/RXQ0NtrSJoMjImIZyIigHuAS5PVlwB3J8sPAPPyo4cjpXF8ta9FMLPUG4w5hM8BjxQ8nybpJUk/l3ROUqsHWgu2aU1q+XVrAZKQ2QZM6O6NJC2U1Cypua2tbRCantM4voq3t+wkl1VmZuk0oECQ9GWgHfh+UtoATImIM4AvAj+QNAbo7jf+/Kdvb+sOLEbcFhGzI2J2bW3tQJp+gCnjq9m5t4PNH/jLcswsvfodCJKuBH4H+KPkMBARsSciNifLLwCrgZPIjQgKDys1AOuT5VagMXnNMmAsXQ5RDbX9dz31YSMzS7F+BYKk+cDfAJ+MiJ0F9VpJ2WT5eHKTx29ExAZgh6S5yfzAFcDiZLclwJXJ8mXAk3GEj93kr0XwLSzMLM3K+tpA0r3AecBESa3AV8idVVQJPJ7M//4qOaPoXOD/SmoHOoAvRET+t/2ryZ2xVEVuziE/73AH8D1JLeRGBgsGpWeHIX8tgkcIZpZmfQZCRFzeTfmOHrZ9EHiwh3XNwIxu6ruBT/fVjqFUVZGldnQla30tgpmlWOqvVM7zXU/NLO0cCInGmioHgpmlmgMhMWV8NRu27WJfR2exm2JmVhQOhETj+Go6A9Zt9TyCmaWTAyHRUJP/XgQHgpmlkwMh0VBTBcC69zyPYGbp5EBIHDN2BBl5hGBm6eVASJRnM0weW+U5BDNLLQdCgfqaKo8QzCy1HAgFGsZVse49B4KZpZMDoUBDTZWvRTCz1HIgFKivqaIz4J1tu4vdFDOzI86BUMDXIphZmjkQCtSPy12L0LrV1yKYWfo4EApMHjcCCU8sm1kqORAKVJZlqRtd6UNGZpZKDoQuGmqqfXGamaWSA6GL+nFVtPp+RmaWQg6ELhpqqtjw3m46OqPYTTEzO6IcCF001FTT3hls3O5rEcwsXfoMBEl3StokaVlBbbykxyW9nvysKVh3g6QWSaskXVRQnyXplWTdzZKU1Csl3ZfUn5U0dXC7eHjq998G2/MIZpYuhzJCuAuY36X2JeCJiGgCnkieI+lUYAEwPdnnFknZZJ9bgYVAU/LIv+ZVwNaIOBG4CfhafzszGPLfi+BrEcwsbfoMhIh4GtjSpXwJcHeyfDdwaUF9UUTsiYg1QAswR9JkYExEPBMRAdzTZZ/8az0AzMuPHoohf3GazzQys7Tp7xzCpIjYAJD8rEvq9cDagu1ak1p9sty1fsA+EdEObAMmdPemkhZKapbU3NbW1s+m925EeZaJo3wtgpmlz2BPKnf3m330Uu9tn4OLEbdFxOyImF1bW9vPJvbN34tgZmnU30DYmBwGIvm5Kam3Ao0F2zUA65N6Qzf1A/aRVAaM5eBDVEdUQ42/F8HM0qe/gbAEuDJZvhJYXFBfkJw5NI3c5PFzyWGlHZLmJvMDV3TZJ/9alwFPJvMMRdMwLvdVmp2+FsHMUqSsrw0k3QucB0yU1Ap8Bfh74H5JVwFvA58GiIjlku4HVgDtwDUR0ZG81NXkzliqAh5JHgB3AN+T1EJuZLBgUHo2AA01Vezt6OTd9/dQN2ZEsZtjZnZE9BkIEXF5D6vm9bD9jcCN3dSbgRnd1HeTBMpwkf9ehLVbdzoQzCw1fKVyN46bkAuENe/6WgQzSw8HQjemjK+mIpvh9U07it0UM7MjxoHQjbJshuNrR9Ky8f1iN8XM7IhxIPTgxLpRvL7JgWBm6eFA6EFT3WjWbt3Jrr0dfW9sZlYCHAg9aJo0ighY3eZRgpmlgwOhB011owAHgpmlhwOhB8dNGEk2I173xLKZpYQDoQcVZRmmTqj2qadmlhoOhF401Y32mUZmlhoOhF40TRrFW5t3sqfdZxqZWelzIPTixLpRdHQGb/oWFmaWAg6EXjTVjQagxYeNzCwFHAi9OL52JBnhiWUzSwUHQi9GlGdpHF/tiWUzSwUHQh+a6kb5JndmlgoOhD6cWDeaN959n/aOzmI3xcxsSDkQ+tBUN4p9HcHbW3ymkZmVNgdCH5om5e5p9OoGTyybWWnrdyBIOlnS0oLHdknXS/qqpHUF9YsL9rlBUoukVZIuKqjPkvRKsu5mSRpoxwbLhyaPYWxVOU+u3FTsppiZDal+B0JErIqImRExE5gF7AQeTlbflF8XET8BkHQqsACYDswHbpGUTba/FVgINCWP+f1t12Arz2aYd0odT6zc6HkEMytpg3XIaB6wOiLe6mWbS4BFEbEnItYALcAcSZOBMRHxTEQEcA9w6SC1a1BcOH0S7+3cx3NrthS7KWZmQ2awAmEBcG/B82slvSzpTkk1Sa0eWFuwTWtSq0+Wu9YPImmhpGZJzW1tbYPU9L6de1ItlWUZHlux8Yi9p5nZkTbgQJBUAXwS+GFSuhU4AZgJbAC+md+0m92jl/rBxYjbImJ2RMyura0dULsPR3VFGeeeVMtjy98hN4gxMys9gzFC+ATwYkRsBIiIjRHRERGdwHeAOcl2rUBjwX4NwPqk3tBNfVi58NRJrN+2m1fWbSt2U8zMhsRgBMLlFBwuSuYE8j4FLEuWlwALJFVKmkZu8vi5iNgA7JA0Nzm76Apg8SC0a1Bd8KFJZDPiseU+bGRmpWlAgSCpGvg48FBB+evJKaQvA+cDfwEQEcuB+4EVwKPANRGR/6KBq4HbyU00rwYeGUi7hkLNyArmTB3PT5e/U+ymmJkNibKB7BwRO4EJXWqf6WX7G4Ebu6k3AzMG0pYj4cLpk/jbH63gjbb3Ob52VLGbY2Y2qHyl8mG4cPoxADzus43MrAQ5EA5D/bgqTjlmNE+t8lXLZlZ6HAiH6fxT6mh+cyvbd+8rdlPMzAaVA+EwnX9yHe2dwS9ff7fYTTEzG1QOhMN05pRxjB5R5sNGZlZyHAiHqSyb4dyTanlqVZuvWjazkuJA6IfzT66jbccelq/fXuymmJkNGgdCP/zmSbn7KP3Mh43MrIQ4EPqhdnQlpzeM5alVR+6Oq2ZmQ82B0E/nnVzHS29vZesHe4vdFDOzQeFA6KfzT6mjM+DHr2wodlPMzAaFA6GfTqsfy+kNY/nK4mXc/l9v+IwjMzvqORD6KZMRP/jfc7lo+jH8v/94lb/64cvs3tfR945mZsOUA2EARlaW8e0/PJPrL2jiwRdb+aPbn2Xz+3uK3Swzs35xIAxQJiOuv+AkbvmjM1m2bhufuuW/adn0frGbZWZ22BwIg+TiD09m0cK57Nzbzu/d8ktfo2BmRx0HwiA6Y0oND//Z2UweW8WffPd5/vqBX7Ntl++KamZHBwfCIGscX83ia8/m6vNO4MEX13HhTT/n6dd8AZuZDX8OhCEwojzL38w/hX//s7MZV1XB5+56nkd8vYKZDXMOhCH04Yax/PDqszi9cRzX3vsSS369vthNMjPr0YACQdKbkl6RtFRSc1IbL+lxSa8nP2sKtr9BUoukVZIuKqjPSl6nRdLNkjSQdg0nY0aUc/fn5jDruBquX/QSi5euK3aTzMy6NRgjhPMjYmZEzE6efwl4IiKagCeS50g6FVgATAfmA7dIyib73AosBJqSx/xBaNewMaqyjLs++xFObxzH3/14Be0dncVukpnZQYbikNElwN3J8t3ApQX1RRGxJyLWAC3AHEmTgTER8Uzk7v9wT8E+JaO6ooyF5xzPu+/v5dk1W4rdHDOzgww0EAJ4TNILkhYmtUkRsQEg+VmX1OuBtQX7tia1+mS5a/0gkhZKapbU3NZ29J25c/4pdYysyLJkqecSzGz4GWggnB0RZwKfAK6RdG4v23Y3LxC91A8uRtwWEbMjYnZtbe3ht7bIRpRnuXD6MTyybAN7233YyMyGlwEFQkSsT35uAh4G5gAbk8NAJD/zl+y2Ao0FuzcA65N6Qzf1kvTJ049l++52/uv1o2+EY2alrd+BIGmkpNH5ZeBCYBmwBLgy2exKYHGyvARYIKlS0jRyk8fPJYeVdkiam5xddEXBPiXn7BMnMq663KegmtmwUzaAfScBDydniJYBP4iIRyU9D9wv6SrgbeDTABGxXNL9wAqgHbgmIvL3i74auAuoAh5JHiWpoizDJ2Ycw+Kl69m1t4OqimzfO5mZHQH9DoSIeAM4vZv6ZmBeD/vcCNzYTb0ZmNHfthxtfve0Y7n3ubU8uXITv33a5GI3x8wM8JXKRfHR4ydQO7qSh15s9Tetmdmw4UAogmxG/OGcKTyxchPffOw1h4KZDQsDmUOwAbhuXhObduzmW0+1EAR/deHJlNAdO8zsKORAKJJMRtx46YcB8e2nVpOV+OKFJxe7WWaWYj5kVES5UJjBp86o51tPtfDOtt3FbpKZpZgDocgyGfHn85roDHj4Jd8J1cyKx4EwDEybOJLZx9XwoM86MrMiciAME5fNaqBl0/v8unVbsZtiZinlQBgmLj5tMiPKMzzwwtq+NzYzGwIOhGFizIhy5k8/hiVL17N7X0ffO5iZDTIHwjBy2axGtu9u5z9f3VjspphZCjkQhpGzTpjAsWNH8MALrX1vbGY2yBwIw0g2Iy6b3cjPX2vjP1d4lGBmR5YDYZi5+jdPYMaxY7lu0Uu8tnFHsZtjZiniQBhmqiqy3HbFLKoqyvj83c1s/WBvsZtkZinhQBiGJo+t4l8/M4t3tu3mqrufZ/HSdby9eacvWjOzIeWb2w1Ts46r4RufPo0vPfgK1y1aCsDUCdXc/4WzqBs9ositM7NS5EAYxi6ZWc/FH57Mqnd28PybW/jbH61g0XNr+fN5TcVumpmVIB8yGubKsxlm1I/ls2dP45ymiSx67m06On3oyMwGX78DQVKjpKckvSppuaTrkvpXJa2TtDR5XFywzw2SWiStknRRQX2WpFeSdTfL3xTTrT+cM4X123bzs1Wbit0UMytBAxkhtAN/GREfAuYC10g6NVl3U0TMTB4/AUjWLQCmA/OBWyRlk+1vBRYCTclj/gDaVbIuOHUStaMr+f6zbxe7KWZWgvodCBGxISJeTJZ3AK8C9b3scgmwKCL2RMQaoAWYI2kyMCYinoncaTT3AJf2t12lrDybYcFHGnlq1SZat+4sdnPMrMQMyhyCpKnAGcCzSelaSS9LulNSTVKrBwpv5dma1OqT5a717t5noaRmSc1tbW2D0fSjzh98pBGA+573XVHNbHANOBAkjQIeBK6PiO3kDv+cAMwENgDfzG/aze7RS/3gYsRtETE7ImbX1tYOtOlHpYaaas4/uY5Fz69lX0dnsZtjZiVkQIEgqZxcGHw/Ih4CiIiNEdEREZ3Ad4A5yeatQGPB7g3A+qTe0E3devDHc6fQtmOPRwlmNqgGcpaRgDuAVyPiHwvqkws2+xSwLFleAiyQVClpGrnJ4+ciYgOwQ9Lc5DWvABb3t11pcP7JdXx02ni+8dNVbPGtLcxskAxkhHA28Bngt7qcYvr15BTSl4Hzgb8AiIjlwP3ACuBR4JqIyH8TzNXA7eQmmlcDjwygXSVPEn936Qze39PON366stjNMbMSoaP1/jizZ8+O5ubmYjejqG78jxXc/os1PHT1xzhjSk3fO5hZ6kl6ISJmd7fOVyofxa674CTqRlfyfxYv99XLZjZgDoSj2KjKMr7826fyyrpt/Ml3n2P9e7uK3SQzO4o5EI5yv3vaZP7u0hk0v7mVi/7paR54odW3yTazfnEgHOUk8Zm5x/Ho9edwyjGj+asf/pob/+NVh4KZHTYHQok4bsJIFi08iyvPOo7bf7GGb/x0lUPBzA6Lvw+hhGQz4qufnM7ejuCWn61mRHnW351gZofMgVBiJHHjpTPY297JPz7+Gvc9v5ZJYyo5ZuwIPn/O8Zzp01PNrAcOhBKUyYivX3YaJx8zipUbdrBxx26efWML/716Mz+69jdoHF9d7Caa2TDkQChR2YxYeO4J+5+/+e4HfPJbv+BPv/cCD179Maoqsr3sbWZp5EnllJg6cST/vOAMXn1nOzc89LInnM3sIB4hpMj5p9Txlx8/iX947DW27drHmVNqOPXYMXz0+AmMqvQ/BbO086dAyvzZeSeyfXc7jy1/h6dW5b5kaNrEkfzb5z9K/biqIrfOzIrJN7dLsff3tPPM6s188b6ljKkq5/uf/yhTJ44sdrPMbAj55nbWrVGVZXz81Encu3AuO/e28/v/+gyvbdxR7GaZWZE4EIwZ9WNZtPAsArjkW7/krl+uodN3TzVLHQeCAXDyMaNZcu3ZfGTaeL76oxVc/p1fsebdD4rdLDM7ghwItt/ksVXc/dmP8LX/9WFWrN/Ob33zZ3z2u8/x2PJ3aO/oLHbzzGyIeVLZuvXOtt18/9m3uO/5tWzasYeq8iyN46torKlmev1YPvuxqdSMrCh2M83sMPU2qexAsF61d3Ty5MpNPPPGZtZu2UXr1p28tnEHIyvLuPq8E/jc2dMYUe6rns2OFg4EG1Sr3tnB1x9dyRMrN1FVnqWqIktGYlRlljnTxnPuSbXMPX4CY6vKKcsIScVuspkljopAkDQf+GcgC9weEX/f2/YOhOL71RubeXTZO+zr6KQzgs3v7+VXb2xm++72/dtIUFmWobqijOqKLKMqy/jQ5DGcOWUcZ0yp4cS6UR5hmB1BvQXCsLhSWVIW+DbwcaAVeF7SkohYUdyWWW/mHj+BucdPOKDW3tHJy+u28eJbW9m1t4O9HZ3s3tfBrn0d7NzTwXu79vGLlnd5+KV1QC4wjh1bxXETqhlZWYbIh0iWkZVZqsrLqCjLUJYR2YwYUZ5l1Igyxowoo7IsiwQCMhKZTO723xmJrERGuTu/ZiSymdw2ZZkM2Ywoy+bqGeVuBFiezVCezVCRzZDN5vbPZv5nvUc5lgbDIhCAOUBLRLwBIGkRcAngQDjKlGUznDmlptfvXYgIWrfu4qW17/FG2/u8+e4HvLVlJ1s+2Jushz3tHXywt4Ode9rZ1xG0d3ZS7Esj8uGQSUInH0aFYXFQbBRsU5gphdt1t7+Uq+eD7YBtVLBNske+LQe9fS9BVrgmoMcbHkrJuxS8xwFbHurfS7c7H/rL9NiTbv5ce3utgUT7UP1i0NuffXeum9fE755+7KC3Y7gEQj2wtuB5K/DRrhtJWggsBJgyZcqRaZkNOkk0jq8+7O9liAh27+tkx+59bN/dzp72DvL/jzojiMj9zD2gozPo7EyWI7fc3hl0dHbSntQjgvYkcPZ2BHvbOw/YrjP/ml1eJ/9+cUD7urSX2F8r/A/f0z5RsCbfts7O3Hvufx32L+zfOiK6/QDMt6+7j5Tuts+Pzrp7je7eo6dQ607XD7yetu/pVXr6gO/pz7Wn1zooyHra6FD/0Mj9negQYqan7fbXu67qJdXGVpX3+X79MVwC4ZD++CPiNuA2yM0hDHWjbHiRRFVFbhK7bkyxW2NWeobLhWmtQGPB8wZgfZHaYmaWSsMlEJ4HmiRNk1QBLACWFLlNZmapMiwOGUVEu6RrgZ+SO+30zohYXuRmmZmlyrAIBICI+Anwk2K3w8wsrYbLISMzMysyB4KZmQEOBDMzSzgQzMwMGEY3tztcktqAt/q5+0Tg3UFsztEijf1OY58hnf1OY5/h8Pt9XETUdrfiqA2EgZDU3NPd/kpZGvudxj5DOvudxj7D4Pbbh4zMzAxwIJiZWSKtgXBbsRtQJGnsdxr7DOnsdxr7DIPY71TOIZiZ2cHSOkIwM7MuHAhmZgakMBAkzZe0SlKLpC8Vuz1DQVKjpKckvSppuaTrkvp4SY9Lej352fP3XB6lJGUlvSTpx8nzNPR5nKQHJK1M/s7PKvV+S/qL5N/2Mkn3ShpRin2WdKekTZKWFdR67KekG5LPtlWSLjrc90tVIEjKAt8GPgGcClwu6dTitmpItAN/GREfAuYC1yT9/BLwREQ0AU8kz0vNdcCrBc/T0Od/Bh6NiFOA08n1v2T7Lake+HNgdkTMIHfL/AWUZp/vAuZ3qXXbz+T/+AJgerLPLcln3iFLVSAAc4CWiHgjIvYCi4BLitymQRcRGyLixWR5B7kPiHpyfb072exu4NLitHBoSGoAfhu4vaBc6n0eA5wL3AEQEXsj4j1KvN/kbt1fJakMqCb3DYsl1+eIeBrY0qXcUz8vARZFxJ6IWAO0kPvMO2RpC4R6YG3B89akVrIkTQXOAJ4FJkXEBsiFBlBXvJYNiX8C/hroLKiVep+PB9qA7yaHym6XNJIS7ndErAP+AXgb2ABsi4jHKOE+d9FTPwf8+Za2QFA3tZI971bSKOBB4PqI2F7s9gwlSb8DbIqIF4rdliOsDDgTuDUizgA+oDQOlfQoOWZ+CTANOBYYKemPi9uqYWHAn29pC4RWoLHgeQO5oWbJkVROLgy+HxEPJeWNkiYn6ycDm4rVviFwNvBJSW+SOxT4W5L+jdLuM+T+TbdGxLPJ8wfIBUQp9/sCYE1EtEXEPuAh4GOUdp8L9dTPAX++pS0QngeaJMNd5M4AAAD+SURBVE2TVEFuAmZJkds06CSJ3DHlVyPiHwtWLQGuTJavBBYf6bYNlYi4ISIaImIqub/XJyPijynhPgNExDvAWkknJ6V5wApKu99vA3MlVSf/1ueRmycr5T4X6qmfS4AFkiolTQOagOcO65UjIlUP4GLgNWA18OVit2eI+vgb5IaKLwNLk8fFwARyZyW8nvwcX+y2DlH/zwN+nCyXfJ+BmUBz8vf970BNqfcb+FtgJbAM+B5QWYp9Bu4lN0+yj9wI4Kre+gl8OflsWwV84nDfz7euMDMzIH2HjMzMrAcOBDMzAxwIZmaWcCCYmRngQDAzs4QDwczMAAeCmZkl/j9lryCmErZMjwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n,d = 500,20\n",
    "X = randn(n,d)\n",
    "wtrue = randn(d)\n",
    "y = X.dot(wtrue) + .1*randn(n)\n",
    "\n",
    "def random_chunk(X,y,size):\n",
    "    n = len(y)\n",
    "    s = randint(n, size=size)\n",
    "    return X[s,:], y[s]\n",
    "\n",
    "\"\"\"f_validation evaluates the objective on the whole dataset\"\"\"\n",
    "def f_validation(w):\n",
    "    return norm(y - X.dot(w))**2\n",
    "\n",
    "\"\"\"gradf evaluates the gradient on the given chunk\"\"\"\n",
    "def gradf(X,y,w):\n",
    "    return 2*X.T @ (X @ w) - 2*X.T @ y\n",
    "\n",
    "maxiters = 100\n",
    "data = [random_chunk(X,y,5) for i in range(maxiters)]\n",
    "\n",
    "w, fks = sgd(gradf, data, w = np.zeros(d), f_validation=f_validation, \n",
    "             validation_freq=1, maxiters=maxiters)\n",
    "\n",
    "plt.plot(fks)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## apply SGD to a large csv\n",
    "Here we read sequential (not random) chunks of the file to compute each gradient step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate a (random) large csv \n",
    "file_name = \"random.csv\"\n",
    "n,d = 10000, 10\n",
    "X = randn(n,d)\n",
    "y = X@randn(d,1)\n",
    "Z = np.hstack([y,X])\n",
    "df = pd.DataFrame(Z)\n",
    "df.to_csv(file_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   Unnamed: 0         0         1         2         3         4         5  \\\n",
      "0           0 -2.898365  0.343570  1.487808 -0.947529 -1.711902 -0.831632   \n",
      "1           1  3.516984 -0.275772  0.338983 -1.518109  0.828724 -0.894747   \n",
      "2           2  0.695540 -1.285947  0.726503 -0.585661  1.033073 -0.895078   \n",
      "3           3  1.664532 -0.922606 -0.550442  1.891081  0.199615 -0.107702   \n",
      "4           4 -0.375763  0.373702 -1.096575 -0.398709  0.270665 -0.083821   \n",
      "\n",
      "          6         7         8         9        10  \n",
      "0 -0.131644 -0.503776  1.049698  1.351236 -0.685529  \n",
      "1  0.275944 -1.424921 -0.900341  1.388891  2.069720  \n",
      "2 -1.197002  0.742137  0.193761  0.655190 -0.256913  \n",
      "3  0.146150  0.185635 -1.813377  0.973152 -0.219059  \n",
      "4  1.108432 -1.793973  0.242327 -0.632688  2.011036  \n"
     ]
    }
   ],
   "source": [
    "for chunk in pd.read_csv(file_name, chunksize=5):\n",
    "    print(chunk)\n",
    "    break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fa3f9cabd60>]"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD4CAYAAADsKpHdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAYoUlEQVR4nO3df2xd533f8feH9/LXpSzLpClbISVKQYS2shbbNeFq8zZkVTtrrRtpQA0oWGphMyDM8NZk6NBZ3R/GBghI0aHtDMQGhNiz3HpxBDeOhKBuoykJshWuHTp2KkuyYjaKJVqqxFj+IesHJYrf/XGfSx1Rl5LIy8tL8XxeAHHP/d5zLp8Dw/zoeZ5zzqOIwMzMrKnRDTAzs7nBgWBmZoADwczMEgeCmZkBDgQzM0uKjW7AdN16662xfPnyRjfDzOyG8vrrr/88IrqrfXbDBsLy5csZGBhodDPMzG4okt6d7DMPGZmZGeBAMDOzxIFgZmaAA8HMzBIHgpmZAQ4EMzNLHAhmZgbkMBB++LOT/NFfv83FMT/228wsK3eB8OMjH/LV7/09Z86PNropZmZzSu4Cob2lAMDZ8xcb3BIzs7kld4HQ0VJ+WsdpB4KZ2WVyFwiVHsLpEQ8ZmZll5S4QKj2EsxfcQzAzy8pdIJRa3UMwM6smf4HgSWUzs6pyFwieVDYzqy53gVCZVPZ9CGZml8tdIFR6CGfcQzAzu8w1A0HSM5JOSHorU/sjSW9L+jtJL0lalPlsi6RBSQcl3Z+p3yNpb/rsCUlK9VZJ30j1VyUtn9lTvFxbcxMSnPGkspnZZa6nh/AssG5CbTewOiI+C/wE2AIgaRWwEbgjHfOkpEI65ilgM7Ay/VS+82Hgg4j4DPAnwB9O92SuhyRKzQX3EMzMJrhmIETED4CTE2rfiYjKP7H/FuhN2+uBFyJiJCIOAYPAvZKWAAsj4pWICOA5YEPmmO1p+0VgbaX3UC+l1qInlc3MJpiJOYR/B7yctnuAI5nPhlKtJ21PrF92TAqZj4Cuar9I0mZJA5IGhoeHp93gUkvBk8pmZhPUFAiS/iswCjxfKVXZLa5Sv9oxVxYjtkVEf0T0d3d3T7W540otRQ8ZmZlNMO1AkLQJeAD4N2kYCMr/8l+a2a0XOJrqvVXqlx0jqQjczIQhqpnW4R6CmdkVphUIktYB/wX4fEScyXy0C9iYrhxaQXny+LWIOAackrQmzQ88BOzMHLMpbf828N1MwNRFe4snlc3MJipeawdJXwc+B9wqaQh4nPJVRa3A7jT/+7cR8e8jYp+kHcB+ykNJj0ZE5S/vI5SvWGqnPOdQmXd4GvgzSYOUewYbZ+bUJtfRUuTExyP1/jVmZjeUawZCRHyhSvnpq+y/FdhapT4ArK5SPwc8eK12zKRSS4HTHjIyM7tM7u5UhvITT/1wOzOzy+UyEDpaiu4hmJlNkMtAaG8pcO7CGBfH6jp3bWZ2Q8llIHjVNDOzK+UyEMYfge0H3JmZjctlIHS0VtZEcA/BzKwil4FQGl81zT0EM7OKnAaC11U2M5sop4HgdZXNzCbKaSB4UtnMbKJcBoLXVTYzu1IuA6E0fpWRewhmZhX5DIQ0ZOQ5BDOzS3IZCG3FApKHjMzMsnIZCE1Nor254EllM7OMXAYCpHWV/SwjM7NxuQ2Ejlb3EMzMsnIbCO3NBU8qm5ll5DYQOlqLfnSFmVlGbgPB6yqbmV0u14HgHoKZ2SXXDARJz0g6IemtTK1T0m5J76TXWzKfbZE0KOmgpPsz9Xsk7U2fPSFJqd4q6Rup/qqk5TN7itV5XWUzs8tdTw/hWWDdhNpjwJ6IWAnsSe+RtArYCNyRjnlSUiEd8xSwGViZfirf+TDwQUR8BvgT4A+nezJT0d5S4MyIewhmZhXXDISI+AFwckJ5PbA9bW8HNmTqL0TESEQcAgaBeyUtARZGxCsREcBzE46pfNeLwNpK76GeOlqLvlPZzCxjunMIt0XEMYD0ujjVe4Ajmf2GUq0nbU+sX3ZMRIwCHwFd02zXdWtvLnD2wkXGxqLev8rM7IYw05PK1f5lH1epX+2YK79c2ixpQNLA8PDwNJtYVllX+azvVjYzA6YfCMfTMBDp9USqDwFLM/v1AkdTvbdK/bJjJBWBm7lyiAqAiNgWEf0R0d/d3T3Nppd5XWUzs8tNNxB2AZvS9iZgZ6a+MV05tILy5PFraVjplKQ1aX7goQnHVL7rt4HvpnmGurq0app7CGZmAMVr7SDp68DngFslDQGPA18Bdkh6GDgMPAgQEfsk7QD2A6PAoxFR+Yv7COUrltqBl9MPwNPAn0kapNwz2DgjZ3YNJa+aZmZ2mWsGQkR8YZKP1k6y/1Zga5X6ALC6Sv0cKVBm03gPwUNGZmZAju9U7hhfRtM9BDMzyHEgXBoycg/BzAxyHQhpXWVPKpuZAbkOhNRD8H0IZmZArgOhctmph4zMzCDHgdDe7EllM7Os3AZCU5MotRQ8qWxmluQ2EKCyapp7CGZmkPtA8LrKZmYVOQ+EAqc9qWxmBjgQ/PhrM7Mk14HQ0Vp0D8HMLMl1ILQ3F3zZqZlZkutA8LrKZmaX5DoQfB+CmdklDgT3EMzMgNwHQnnIaGys7it2mpnNeTkPhPLzjHzpqZlZ3gOh1esqm5lV5DoQOryuspnZuFwHgldNMzO7JOeBUB4yOnvBPQQzs5oCQdJ/krRP0luSvi6pTVKnpN2S3kmvt2T23yJpUNJBSfdn6vdI2ps+e0KSamnX9XIPwczskmkHgqQe4HeB/ohYDRSAjcBjwJ6IWAnsSe+RtCp9fgewDnhSUiF93VPAZmBl+lk33XZNxfi6yp5UNjOrecioCLRLKgIl4CiwHtiePt8ObEjb64EXImIkIg4Bg8C9kpYACyPilYgI4LnMMXXV0epJZTOzimkHQkS8B/wP4DBwDPgoIr4D3BYRx9I+x4DF6ZAe4EjmK4ZSrSdtT6xfQdJmSQOSBoaHh6fb9HHtlSEj9xDMzGoaMrqF8r/6VwCfAjokffFqh1SpxVXqVxYjtkVEf0T0d3d3T7XJV+ioTCq7h2BmVtOQ0a8BhyJiOCIuAN8E/glwPA0DkV5PpP2HgKWZ43spDzENpe2J9bprb/aksplZRS2BcBhYI6mUrgpaCxwAdgGb0j6bgJ1pexewUVKrpBWUJ49fS8NKpyStSd/zUOaYumpqEu3NXjXNzAzKk8LTEhGvSnoR+BEwCrwBbAMWADskPUw5NB5M+++TtAPYn/Z/NCIqf4kfAZ4F2oGX08+s6Gj1uspmZlBDIABExOPA4xPKI5R7C9X23wpsrVIfAFbX0pbpavcjsM3MgJzfqQzliWVfdmpm5kBwD8HMLMl9IHS0eF1lMzNwIFBq8aSymRk4ELyusplZ4kBo9ZCRmRk4ECg1F3yVkZkZDgRKrUXOXrjI2FjVxyeZmeVG7gOho6VABJwb9bCRmeVb7gPBq6aZmZU5EMYfge1AMLN8cyCML5LjiWUzyzcHQqvXVTYzAwcCHS1eV9nMDBwIl9ZV9qSymeVc7gNhfF3lC+4hmFm+5T4QfNmpmVmZA6HVl52amYEDgfZmX3ZqZgYOBApNoq25yZedmlnu5T4QwOsqm5mBAwFI6yp7UtnMcq6mQJC0SNKLkt6WdEDSP5bUKWm3pHfS6y2Z/bdIGpR0UNL9mfo9kvamz56QpFraNVVeV9nMrPYewv8E/ioifhG4EzgAPAbsiYiVwJ70HkmrgI3AHcA64ElJhfQ9TwGbgZXpZ12N7ZqSUmvBk8pmlnvTDgRJC4F/DjwNEBHnI+JDYD2wPe22HdiQttcDL0TESEQcAgaBeyUtARZGxCsREcBzmWNmhddVNjOrrYfwaWAY+F+S3pD0NUkdwG0RcQwgvS5O+/cARzLHD6VaT9qeWL+CpM2SBiQNDA8P19D0y5U8ZGRmVlMgFIFfBp6KiLuB06ThoUlUmxeIq9SvLEZsi4j+iOjv7u6eansnVe4heMjIzPKtlkAYAoYi4tX0/kXKAXE8DQORXk9k9l+aOb4XOJrqvVXqs8Y9BDOzGgIhIv4BOCLpF1JpLbAf2AVsSrVNwM60vQvYKKlV0grKk8evpWGlU5LWpKuLHsocMys6WgqcGXEPwczyrVjj8f8ReF5SC/BT4N9SDpkdkh4GDgMPAkTEPkk7KIfGKPBoRFT+Wf4I8CzQDrycfmZNqaXAmQsXiQhm+YpXM7M5o6ZAiIg3gf4qH62dZP+twNYq9QFgdS1tqUWptUgEnLswNr4+gplZ3vhOZbyuspkZOBCA8qQy4MdXmFmuORDIrKvsVdPMLMccCHhdZTMzcCAA0OFV08zMHAjgVdPMzMCBAFzqIfjxFWaWZw4EMpPKHjIysxxzIHBpUtmXnZpZnjkQyNyH4B6CmeWYAwEoNInWYpPnEMws1xwISUdr0VcZmVmuORASL6NpZnnnQEhKLQVPKptZrjkQklJLkTMXHAhmll8OhKTkVdPMLOccCEmppchpzyGYWY45EJKO1gJnfZWRmeWYAyEptRTcQzCzXHMgJKWWoh9/bWa55kBIyj2EUSKi0U0xM2sIB0JSaikSAecujDW6KWZmDVFzIEgqSHpD0rfT+05JuyW9k15vyey7RdKgpIOS7s/U75G0N332hCTV2q6p6mitPALbE8tmlk8z0UP4EnAg8/4xYE9ErAT2pPdIWgVsBO4A1gFPSiqkY54CNgMr08+6GWjXlHSkJ56eOudAMLN8qikQJPUCvwl8LVNeD2xP29uBDZn6CxExEhGHgEHgXklLgIUR8UqUB/Cfyxwza5YsagNg6IOzs/2rzczmhFp7CH8K/D6QHXi/LSKOAaTXxaneAxzJ7DeUaj1pe2L9CpI2SxqQNDA8PFxj0y/X19UBwLsnT8/o95qZ3SimHQiSHgBORMTr13tIlVpcpX5lMWJbRPRHRH93d/d1/trrc/vCNloKTRx+/8yMfq+Z2Y2iWMOx9wGfl/QbQBuwUNKfA8clLYmIY2k46ETafwhYmjm+Fzia6r1V6rOq0CR6O9t514FgZjk17R5CRGyJiN6IWE55svi7EfFFYBewKe22CdiZtncBGyW1SlpBefL4tTSsdErSmnR10UOZY2ZVX2eJd086EMwsn2rpIUzmK8AOSQ8Dh4EHASJin6QdwH5gFHg0Iiq3Bj8CPAu0Ay+nn1nX19XBa4dOEhE04MpXM7OGmpFAiIjvA99P2+8DayfZbyuwtUp9AFg9E22pxbLOEqfPX+T90+e5dUFro5tjZjarfKdyRl9XCYDDHjYysxxyIGSMB4Inls0shxwIGb23lAPBVxqZWR45EDLamgvcvrDNN6eZWS45ECZY1lXykJGZ5ZIDYQLfi2BmeeVAmKCvq8TwqRE/BtvMcseBMMGy9JA7X3pqZnnjQJigr9NXGplZPjkQJqjci3DEPQQzyxkHwgSLSi0sbCu6h2BmueNAqGJZl680MrP8cSBU0dfZweH3fXOameWLA6GKZV0lhj44y+jFsWvvbGY2TzgQqujrLDE6Fhz76Fyjm2JmNmscCFUs6/Klp2aWPw6EKvrSzWl+yJ2Z5YkDoYrbF7bRUmjyQ+7MLFccCFUUmkRvZ7uHjMwsVxwIk+jrLPl5RmaWKw6ESfR1dXD45BkiotFNMTObFQ6ESSztLPHJyCgnT59vdFPMzGbFtANB0lJJ35N0QNI+SV9K9U5JuyW9k15vyRyzRdKgpIOS7s/U75G0N332hCTVdlq1G3/qqYeNzCwnaukhjAK/FxG/BKwBHpW0CngM2BMRK4E96T3ps43AHcA64ElJhfRdTwGbgZXpZ10N7ZoRlaee+kojM8uLaQdCRByLiB+l7VPAAaAHWA9sT7ttBzak7fXACxExEhGHgEHgXklLgIUR8UqUB+yfyxzTMEu9LoKZ5cyMzCFIWg7cDbwK3BYRx6AcGsDitFsPcCRz2FCq9aTtifVqv2ezpAFJA8PDwzPR9Em1NRe4fWGbb04zs9yoORAkLQD+AvhyRHx8tV2r1OIq9SuLEdsioj8i+ru7u6fe2Cla1lXykJGZ5UZNgSCpmXIYPB8R30zl42kYiPR6ItWHgKWZw3uBo6neW6XecH2dXhfBzPKjlquMBDwNHIiIP858tAvYlLY3ATsz9Y2SWiWtoDx5/FoaVjolaU36zocyxzRUX1eJ4VMjnD1/sdFNMTOru2INx94H/A6wV9KbqfYHwFeAHZIeBg4DDwJExD5JO4D9lK9QejQiKn9pHwGeBdqBl9NPwy1LD7k7fPIMv3D7TQ1ujZlZfU07ECLi/1F9/B9g7STHbAW2VqkPAKun25Z6WTZ+pdFpB4KZzXu+U/kqKjen+ZlGZpYHDoSrWFRq5qa2ou9FMLNccCBchST6unylkZnlgwPhGvo6Ozj8vm9OM7P5z4FwDcu6Sgx9cJbRi2ONboqZWV05EK6hr7PE6Fhw7KNzjW6KmVldORCu4ReXLATg1UMnG9wSM7P6ciBcw529N7O0s52db77X6KaYmdWVA+EaJLHhrh7+ZvDnnPjYw0ZmNn85EK7D+rt6GAvY9eM58cw9M7O6cCBch88sXsA/6rmZnW86EMxs/nIgXKcNd/ew972PGDzxSaObYmZWFw6E6/Rbdy6hSXhy2czmLQfCdVp8Uxv3feZWvvXme5SXfjYzm18cCFOw4a4ejpw8y48Of9DoppiZzTgHwhTcv/p22pqbeOkNDxuZ2fzjQJiCBa1Ffn3V7Xz7745xftTPNjKz+cWBMEX/+u5P8eGZC/zgJ8ONboqZ2YxyIEzRP1vZTWdHCy/5aiMzm2ccCFPUXGjigc8u4f/sP86pcxca3RwzsxnjQJiG9Xf1MDI6xl+99Q+NboqZ2YxxIEzDLy9bxLLOkh9lYWbzypwJBEnrJB2UNCjpsUa352rKT0D9FH/z9z/nW2+8x5GTZ3yzmpnd8IqNbgCApALwVeDXgSHgh5J2RcT+xrZscg/2L+X5Vw/z5W+8CUBXRwt3Ll3Enb2L+OzSm1l8Uys3tTazoK3IgtYiLcU5k71mZlXNiUAA7gUGI+KnAJJeANYDczYQlnaWeGXLWn5y/BRvHPmQH6ef7x08QbXOQmuxiQWtRZoLTTQJmppEkzS+rUl+jzTZJ5Ob+hFmdiP53bUr+a07PzXj3ztXAqEHOJJ5PwT8ysSdJG0GNgMsW7Zsdlp2FS3FJlb33Mzqnpv5nTV9AJw6d4H9Rz/mgzMX+GRklE/OlV9PnRvl1MgoFy8GFyMYiyACLo6V31c1jVGomM5BZnZDubm9uS7fO1cCodo/aq/4yxYR24BtAP39/XPyL99Nbc38yqe7Gt0MM7MpmysD20PA0sz7XsCX8JiZzaK5Egg/BFZKWiGpBdgI7Gpwm8zMcmVODBlFxKik/wD8NVAAnomIfQ1ulplZrsyJQACIiL8E/rLR7TAzy6u5MmRkZmYN5kAwMzPAgWBmZokDwczMANCN+lA2ScPAu9M8/Fbg5zPYnBuBzzkffM75UMs590VEd7UPbthAqIWkgYjob3Q7ZpPPOR98zvlQr3P2kJGZmQEOBDMzS/IaCNsa3YAG8Dnng885H+pyzrmcQzAzsyvltYdgZmYTOBDMzAzIYSBIWifpoKRBSY81uj31IOkZSSckvZWpdUraLemd9HpLI9s4kyQtlfQ9SQck7ZP0pVSfz+fcJuk1ST9O5/zfUn3ennOFpIKkNyR9O72f1+cs6WeS9kp6U9JAqtXlnHMVCJIKwFeBfwWsAr4gaVVjW1UXzwLrJtQeA/ZExEpgT3o/X4wCvxcRvwSsAR5N/13n8zmPAL8aEXcCdwHrJK1hfp9zxZeAA5n3eTjnfxERd2XuPajLOecqEIB7gcGI+GlEnAdeANY3uE0zLiJ+AJycUF4PbE/b24ENs9qoOoqIYxHxo7R9ivIfix7m9zlHRHyS3jann2AenzOApF7gN4GvZcrz+pwnUZdzzlsg9ABHMu+HUi0PbouIY1D+AwosbnB76kLScuBu4FXm+TmnoZM3gRPA7oiY9+cM/Cnw+8BYpjbfzzmA70h6XdLmVKvLOc+ZBXJmiarUfN3tPCFpAfAXwJcj4mOp2n/u+SMiLgJ3SVoEvCRpdaPbVE+SHgBORMTrkj7X6PbMovsi4qikxcBuSW/X6xflrYcwBCzNvO8FjjaoLbPtuKQlAOn1RIPbM6MkNVMOg+cj4pupPK/PuSIiPgS+T3neaD6f833A5yX9jPJw769K+nPm9zkTEUfT6wngJcpD33U557wFwg+BlZJWSGoBNgK7Gtym2bIL2JS2NwE7G9iWGaVyV+Bp4EBE/HHmo/l8zt2pZ4CkduDXgLeZx+ccEVsiojcillP+f/e7EfFF5vE5S+qQdFNlG/iXwFvU6Zxzd6eypN+gPA5ZAJ6JiK0NbtKMk/R14HOUH5F7HHgc+BawA1gGHAYejIiJE883JEn/FPi/wF4ujS3/AeV5hPl6zp+lPJlYoPwPux0R8d8ldTFPzzkrDRn954h4YD6fs6RPU+4VQHmI/39HxNZ6nXPuAsHMzKrL25CRmZlNwoFgZmaAA8HMzBIHgpmZAQ4EMzNLHAhmZgY4EMzMLPn/QAlAmC5Y1Y0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "data_validation = pd.read_csv(file_name, nrows=1000)\n",
    "X_validation, y_validation = data_validation.iloc[:,2:], data_validation.iloc[:,1]\n",
    "\n",
    "\"\"\"f_validation evaluates the objective on the whole dataset\"\"\"\n",
    "def f_validation(w):\n",
    "    return norm(y_validation - X_validation.dot(w))**2\n",
    "\n",
    "\"\"\"gradf evaluates the gradient on the given chunk\"\"\"\n",
    "def gradf(X,y,w):\n",
    "    return 2*X.T @ (X @ w) - 2*X.T @ y\n",
    "\n",
    "maxiters = 100\n",
    "data = [(chunk.iloc[:,2:], chunk.iloc[:,1]) for chunk in pd.read_csv(file_name, chunksize=20)]\n",
    "\n",
    "w, fks = sgd(gradf, data, w=np.zeros(d), \n",
    "             f_validation=f_validation, \n",
    "             validation_freq=10, maxiters=maxiters)\n",
    "\n",
    "plt.plot(fks)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# preprocess some large datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import tensorflow.compat.v2 as tf\n",
    "# import tensorflow_datasets as tfds\n",
    "\n",
    "# # Construct a tf.data.Dataset\n",
    "# ds = tfds.load('higgs', split='train', shuffle_files=True)\n",
    "\n",
    "# # Build your input pipeline\n",
    "# ds = ds.shuffle(1024).batch(32).prefetch(tf.data.experimental.AUTOTUNE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for example in ds.take(1):\n",
    "#   print(example)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dr['df'].to_csv('yaleemmlc.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pytorch utils\n",
    "import torch\n",
    "from torch import optim\n",
    "import torch.nn as nn\n",
    "import torch.nn.functional as F"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we have to call `.float()` to reconcile the default types of NumPy and PyTorch. NumPy defaults to float64 (\"double\") on a 64-bit system. PyTorch models default to float32 (\"float\") as GPUs can do up to 64x the operations per clock cycle when using float32 instead of float64. \n",
    "\n",
    "[CUDA operations per cycle comparison](https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#arithmetic-instructions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'x' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-3-bc005f20f9b8>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     16\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mepoch\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn_epochs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     17\u001b[0m     \u001b[0;31m# Let's do a random order after converting to a torch array\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 18\u001b[0;31m     \u001b[0mshuffled_indices\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtorch\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrandperm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     19\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mbatch\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbatch_size\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     20\u001b[0m         \u001b[0msample_idxs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mshuffled_indices\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mbatch\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mbatch\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0mbatch_size\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'x' is not defined"
     ]
    }
   ],
   "source": [
    "# method 4: PyTorch via SGD\n",
    "\n",
    "def TorchLinRegModel(n : int = 2):\n",
    "    \"\"\"Create a simple linear regression model in PyTorch. Note that\n",
    "    we're manually including the bias in this demo.\"\"\"\n",
    "    return nn.Linear(in_features=n, out_features=1, bias=False)\n",
    "\n",
    "model = TorchLinRegModel()\n",
    "model.zero_grad()\n",
    "optimizer = optim.SGD(model.parameters(), lr=5e-6)\n",
    "loss_function = nn.MSELoss()\n",
    "\n",
    "n_epochs = 5000\n",
    "batch_size = 8\n",
    "\n",
    "for epoch in range(n_epochs):\n",
    "    # Let's do a random order after converting to a torch array\n",
    "    shuffled_indices = torch.randperm(x.shape[0])\n",
    "    for batch in range(0, x.shape[0], batch_size):\n",
    "        sample_idxs = shuffled_indices[batch:batch+batch_size]\n",
    "        \n",
    "        # Read minibatch data. PyTorch models default to float32  \n",
    "        # so we need to convert from double precision to float32\n",
    "        tensor = torch.from_numpy(X[sample_idxs, :]).float()\n",
    "        output = model(tensor)\n",
    "        \n",
    "        batch_labels = torch.from_numpy(y[sample_idxs]).float()\n",
    "        batch_labels = batch_labels.view(-1, 1)\n",
    "        \n",
    "        # loss.backward() calculates the gradient of the loss with\n",
    "        # respect to the weights\n",
    "        loss = loss_function(output, batch_labels)\n",
    "        loss.backward()\n",
    "        \n",
    "        # optimizer.step() moves in the direction of the gradient\n",
    "        # based on the learning rate (i.e. step size)\n",
    "        optimizer.step()\n",
    "        \n",
    "        # PyTorch is still saving our past gradient even though we \n",
    "        # changed our model parameters. Clearing them here prevents\n",
    "        # us from accidentally applying them again.\n",
    "        optimizer.zero_grad()\n",
    "        \n",
    "        \n",
    "w = list(model.parameters())[0][0].detach().numpy()\n",
    "w"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}