GradientDescent.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Local search\n",
    "\n",
    "In unconstrained optimization, we wish to solve problems of the form\n",
    "\\begin{align}\n",
    "\\text{minimize} & & E(w) \n",
    "\\end{align}\n",
    "\n",
    "* The local search algorithms have the form :\n",
    "\\begin{align}\n",
    "w_0 & = \\text{some initial value} \\\\\n",
    "\\text{for}\\;\\; & \\tau = 1, 2,\\dots \\\\\n",
    "& w_\\tau = w_{\\tau-1} + g_\\tau\n",
    "\\end{align}\n",
    "\n",
    "Here, $g_\\tau$ is a search direction. The loop is executed until a convergence condition is satisfied or \n",
    "the maximum number of iterations is reached. The algorithm iteratively search for solutions that achieve a lower objective value by moving in the search direction.\n",
    "\n",
    "# Gradient Descent \n",
    "* Gradient descent is a popular local search method with the search direction chosen as the negative gradient direction:\n",
    "\\begin{align}\n",
    "g_\\tau & =  - \\eta \\nabla E(w_{\\tau-1})\n",
    "\\end{align}\n",
    "\n",
    "* When the gradient vanishes, i.e., $\\nabla E(w) = 0$, the algorithm does not make any progress. Such points are also called fixed points. \n",
    "\n",
    "* The iterates, under certain conditions, converge to the minimum $w^* = \\arg\\min_{w} E(w)$. A natural question here finding the conditions for guaranteed convergence to a fixed point and the rate -- how fast convergence happens as a function of iterations\n",
    "\n",
    "* The parameter $\\eta$ is called the *learning rate*, to be chosen depending on the problem. If the learning rate is not properly chosen, the algorithm can (and will) diverge.\n",
    "\n",
    "* There is a well developed theory on how to choose $\\eta$ adaptively to speed up covergence.\n",
    "\n",
    "* Even for minimizing quadratic objective functions, or equivalently for solving linear systems, gradient descent can have a quite poor converge properties: it takes a lot of iterations to find the minimum. However, it is applicable as a practical method in many problems as it requires only the calculation of the gradient.\n",
    "\n",
    "* For maximization problems\n",
    "\\begin{align}\n",
    "\\text{maximize} & & E(w) \n",
    "\\end{align}\n",
    "we just move in the direction of the gradient so the search direction is $g_\\tau = \\eta \\nabla E(w_{\\tau-1})$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true,
    "hide_input": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pylab as plt\n",
    "\n",
    "from ipywidgets import interact, interactive, fixed\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import clear_output, display, HTML\n",
    "from matplotlib import rc\n",
    "\n",
    "mpl.rc('font',**{'size': 20, 'family':'sans-serif','sans-serif':['Helvetica']})\n",
    "mpl.rc('text', usetex=True)\n",
    "\n",
    "import time\n",
    "import numpy as np\n",
    "\n",
    "y = np.array([7.04, 7.95, 7.58, 7.81, 8.33, 7.96, 8.24, 8.26, 7.84, 6.82, 5.68])\n",
    "x = np.array(np.linspace(-1,1,11))\n",
    "N = len(x)\n",
    "\n",
    "# Design matrix\n",
    "#A = np.vstack((np.ones(N), x, x**2, x**3)).T\n",
    "degree = 9\n",
    "A = np.hstack([np.power(x.reshape(N,1),i) for i in range(degree+1)])\n",
    "\n",
    "# Learning rate\n",
    "eta = 0.001\n",
    "              \n",
    "# initial parameters\n",
    "w = np.array(np.random.randn(degree+1))\n",
    "W = []\n",
    "Err = []\n",
    "for epoch in range(50000):  \n",
    "    # Error\n",
    "    err = y-A.dot(w)\n",
    "    \n",
    "    # Total error\n",
    "    E = np.sum(err**2)/N\n",
    "    \n",
    "    # Gradient\n",
    "    dE = -2.*A.T.dot(err)/N\n",
    "    \n",
    "    if epoch%100 == 0: \n",
    "        #print(epoch,':',E)\n",
    "        # print(w)    \n",
    "        W.append(w)\n",
    "        Err.append(E)\n",
    "\n",
    "    # Perfom one descent step\n",
    "    w = w - eta*dE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell demonstrates interactively the progress of plain gradient descent \n",
    "and how its solution differs from the optimum found by solving the corresponding least squares problem.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "339343d8a19e49efaf00e0faa815ea25",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>interactive</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in Jupyter Notebook or JupyterLab, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another notebook frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "interactive(children=(IntSlider(value=0, description='iteration', max=499), Output()), _dom_classes=('widget-interact',))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "fig = plt.figure(figsize=(5,5))\n",
    "\n",
    "left = -1.5\n",
    "right = 1.5\n",
    "xx = np.linspace(left,right,50)\n",
    "AA = np.hstack((np.power(xx.reshape(len(xx),1),i) for i in range(degree+1)))\n",
    "\n",
    "# Find best\n",
    "A_orth, R = np.linalg.qr(A)\n",
    "w_orth, res, rank, s = np.linalg.lstsq(A_orth, y)\n",
    "w_star = np.linalg.solve(R, w_orth)\n",
    "yy = AA.dot(w_star)\n",
    "\n",
    "#ax.set_xlim((2,15))\n",
    "\n",
    "#dots = plt.Line2D(x,y, linestyle='', markerfacecolor='b',marker='o', alpha=0.5, markersize=5)\n",
    "#ax.add_line(dots)\n",
    "plt.plot(x,y, linestyle='', markerfacecolor='b',marker='o', alpha=0.5, markersize=5)\n",
    "plt.plot(xx, yy, linestyle=':', color='k', alpha=0.3)\n",
    "ln = plt.Line2D(xdata=[], ydata=[], linestyle='-',linewidth=2)\n",
    "\n",
    "\n",
    "ax = fig.gca()\n",
    "ax.add_line(ln)\n",
    "plt.close(fig)\n",
    "\n",
    "ax.set_xlim((left,right))\n",
    "ax.set_ylim((5,9))\n",
    "\n",
    "def plot_gd(iteration=0):\n",
    "    w = W[iteration]\n",
    "    f = AA.dot(w)\n",
    "    #print(w)\n",
    "    ln.set_ydata(f)\n",
    "    ln.set_xdata(xx)\n",
    "    \n",
    "    ax.set_title('$E = '+str(Err[iteration])+'$')\n",
    "    \n",
    "    display(fig)\n",
    "    \n",
    "res = interact(plot_gd, iteration=(0,len(W)-1))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting the Error Surface"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEICAYAAAC3Y/QeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG8JJREFUeJzt3X1sXNd55/HfMyQlMbLlsWQ5jejSshTbabxOLIraLZKg\nxSq0hW63gSPLcYBFgWzR0unutl201ouDLRCgwMqS0mzTfxrJQXa3SRd1zLhGXy1bCpykTtJYL1a8\nG0d2RFmJqdjW28iSxbchn/1j7jjD4b3kvdTM3Lkz3w9AcObycu65Gurh4XPOeY65uwAA7SWXdgMA\nAI1H8AeANkTwB4A2RPAHgDZE8AeANkTwB4A2RPAHgDZE8AeANkTwB4A21Jl2A6LccMMNvnr16rSb\nAQCZcfjw4bPuvjLOuU0b/FevXq1Dhw6l3QwAyAwzOxX3XNI+ANCGCP4A0IYI/gDQhgj+ANCG6jrg\na2Z97n6k4vlg8HCtu2+v57UBANHqFvzNbEDSXklrK54fcPdhM3vczAbc/UC9rg8AWfLk0RHt2X9c\npwujWpXv1tZNt+vedT11u17d0j5BYB+uOLRG0kDweDh4DgBt78mjI3r4iRc1UhiVSxopjOrhJ17U\nk0dH6nbNhuX83X2fu+8LnvZJYhI/AEjas/+4RienZhwbnZzSnv3H63bNhg/4mlmfpCOVYwEVXxs0\ns0NmdujMmTONbhoApOJ0YTTR8VpIY7bPQNRgb/DXQb+7969cGWuFMgBk3qp8d6LjtdDQ4G9mg+6+\nO3g8MN/5ANAOHvzV2UOg3V0d2rrp9rpds27B38y2SOoPPpeD/S4zO2FmF+p1XQDImp+evyJJevey\nxTJJPflu7dx8Z11n+9Rtqqe7D0kaqnh+QNL19boeAGTRm5fG9JXvndLmvh59/hN3Ney6rPAFgBTt\n/eawJqdcv7/x1oZel+APACl589KYvvq9U7r3rh6tvmFpQ69N8AeAlHzx2WEVp12//9H3NvzaBH8A\nSMGbb43pr/7llDav69HNKxrb65eaeCcvAGhF5Ro+I8ECrtt/4dpU2kHPHwAapLKGT9mfPv1yXWv4\nRCH4A0CDpFHDJwrBHwAaJI0aPlEI/gDQIGnU8IlC8AeABkmjhk8Ugj8ANEgaNXyiMNUTABrg7OXx\nUg2fdT36/AONq+EThZ4/ADTAo98a1kRxWv95Y+NX84Yh+ANAnZ27PK6//O4p/cYHV2ntymvSbo4k\ngj8A1N2X/vmkxopT+r0m6fVLBH8AqKsLb0/oL7/zqn79zvfovTemU8ohDAO+AFAH1TV87li1LOUW\nzUTPHwBqLKyGz58f/HEqNXyiEPwBoMaaqYZPFII/ANRYM9XwiULwB4Aaa6YaPlEI/gBQY7+3ce2s\nY2nV8IlC8AeAGjt7eUKStPLa9Gv4RGGqJwDU0OXxor70zye18X036suf2pB2cyLVtedvZn1Vz7eY\n2YCZbavndQEgLV/57ikVrkw21WreMHUL/mY2IOnxiud9kuTuByQVqn8xAEDWXZko6tFvD+tXblup\ndb3Xp92cOdUt+AdBfrji0AOSCsHjYUkD9bo2AKThr773E51/e0J/8NFb027KvBo54JuXdL7i+YoG\nXhsA6mp0Ykp7v3VCH3nvDVp/c3P3+qUmG/A1s0FJg5LU29ubcmsAYH7VNXzW9eZTblE8jez5FyQt\nDx7nJZ2rPsHd97l7v7v3r1y5soFNA4Dkwmr4fOnbJ5uqhk+URgb/xySVdy9eI+lAA68NADWXhRo+\nUeo522eLpP7gs9z9SHB8QFKh/BwAsioLNXyi1C3n7+5Dkoaqju2r1/UAoNFW5btnpHwqjzc7yjsA\nwAJtuuPds441Ww2fKE012wcAsuLyeFH/+OLres+yxbKc6WeFMa3Kd2vrptubqoZPFII/ACzA5/Yf\n1xuXxvTE736o6VfzhiHtAwAJHftpQf/7u6/qN3/55kwGfomePwDEUl7Mdbowqo6c6drFHXooA7n9\nKPT8AWAelYu5XFJx2jU2Oa1vvPRm2k1bMII/AMwjbDHXxJRnYjFXFII/AMwjy4u5ohD8AWAeWdiQ\nPSmCPwDM44/uvk1mM49lZTFXFII/AMyj6C536fp3dTXthuxJMdUTAOZwebyoPfuPa11vXk/87odk\n1X8CZBQ9fwCYwxefPaEzl8b1x//+/S0T+CWCPwBEeu3CFT367WF97IOr1JfRlbxRCP4AEGHXU6V5\n/Nt/7X0pt6T2yPkDQIXKMg4u6Z7336ieDE/pjELPHwAC1WUcJOlbr5zNxJ68SRH8ASAQVsZhbHI6\n02UcohD8ASDQimUcohD8ASCwKr8k4jg5fwBoWb9y28pZx7JexiEKs30AQKXUzt8f+5nW3LBUY8Wp\nzO3JmxTBH0Dbc3dt//oPVJx2/c//uEE3r1iadpPqjrQPgLb318//VN9+5aw+8+/e1xaBX6LnD6BN\nVS7mkqRbb1yq//Bvbk65VY3T0J6/mW0xswEzG2zkdQGgUvViLpf0k/Oj+ttjp9NuWsM0LPibWZ+k\nYXc/IGk4eA4ADRe2mGu82JqLuaLETvuY2WpJfZI2SHpe0hF3fzXh9XZJulvSmuCXAAA0XDst5ooy\nb8/fzNaZ2dckbZe0QtKB4PMOM3vMzO6KcyF3P6JSj/+CpPNX0WYAuCqtuCdvUnF6/v3u/omqYwcl\nPSpJZvY7kl6Y70XMLC+pIGmnpEfN7Ii7D1edMyhpUJJ6e3tjNA0Aktt0x7v15edenXGsVRdzRTF3\nn/8sSWa2TNJylXrtg5KGkqR9zGybpH3uXjCzLSqlfnZHnd/f3++HDh2K+/IAEMulsUkNfP6bykmy\nnLXUYi4zO+zu/XHOTTLV81FJ+1RK/wyrlL9/IHnzJHcfYsYPgDTsfuq43rw0rr/5Tx/WXb+YT7s5\nqUky2yfv7gdV6rHvkZToXy3o5Q8G0z0H3X1fku8HgKt1+NR5ffVfTulTH1rd1oFfStbzNzPbKelo\nMMi7NunF5krzAEA9VC7m6siZrlvSqT+6p31y+1GS9PwfVCnfv1Ol6Z7316VFAFAj1Yu5itOuKxPT\nOvDDN9JuWupiB393P+nue9z9ors/6u5H69kwALhaYYu5JqbaazFXlNjB38zuM7Mfm9kr5c/1bBgA\nXC0Wc0VLkvPfIWm9u1+sV2MAoJbek1+i04WxWcfbaTFXlCQ5/wsEfgBZ8sGbZs/oabfFXFGS9PwL\nZva8SuUdJEnu/nDtmwQAV+/gS2/oqf/3utb35vX6W2M63UKLuWohSfDfW7dWAEANnThzWf/1r1/Q\n+9+zTF/97V9W96KOtJvUdGIH/2CBFwA0per5/Is7c9r7m+sJ/BGSzPb5qJmdN7Nzwce/rWfDACCu\nsPn8k9OuQ69eSLtpTSvJgO8jkm5x9xUqLfJitS6AphA6n7/NNmdJakGzfYJSzPxKBdAUmM+fXNLZ\nPjslPSPpHpVq8wNAqn7wWkFmUlh1eubzR0sy4PsJM9sq6ROSTrj7jvo1CwDCVQ7s5t/VpbdGJ7Vs\nSZdGJ6c0Xpx+5zzm888tSc9fQSlnAEhFeWC3nN+/cGVSOZMe2nSbrlnc9c4vBebzz2/e4G9mD7n7\n58zsEUkz/rBikReARgob2J126S+eHdZzOzYS7BOI0/Mv77P7fNXxePs/AkCNMLBbO3GC/782sw1V\nx0yl4P9E7ZsEALO9eWlMHTlTcXp2v5OB3eTiBP/qHn8ZPX8AdVM5sLvy2sUqTk/LTFrUkdPEFAO7\nVytO8N+g2YGenj+Auqke2H3z0rgk6Q/vvlW9y5cysFsD9PwBNJ2wgV1Jeuz51xjYrZE4wf+WkNk+\n9PwB1A0Du/WXdLZPdfAHgJpb1t2li6OTs44zsFs7cYL/3WZ2d/DYKo7T8wdQE5WDu9cs6dSlsaJy\nVprDX8bAbm3FCf4PSDonaUiluj429+kAEF/14O6lsaI6zPTJDTfp2ZfPMrBbJ/MGf3dfbmbrJD2o\nUhnnZyTtdfdXk17MzPokrQledyjp9wNoPWGDu1Puevbls3pux8aUWtX6YpV0dvej7v5pd+9XaQ/f\nXWb2ygKu93AQ9NcEvwgAtDkGd9ORqLCbmW2UdL+ktZL2JfzeLQqmjbo7G8EAUHFqWku6chqdnJ71\nNQZ36ytOYbe7VEr5lHv9X3T3FxZwrQ3B6/VJGuAXANDepqZdDz1+TKOT0+qsKtvA4G79xen5H5F0\nQtJRlXr8D5uVxnzd/YGE1zvn7kfMbMDMtlTn/c1sUNKgJPX29iZ8aQDNrnJWT/eiDl2ZmNLWTber\nJ9/Nqt0GizXVs0bXOqefrxkoqPSXwIzg7+77FKST+vv7WUcAtJDqWT1XJqbUmTP15Lt177oegn2D\nxRnwvVulHvvB6g8zWxds7RjHkIKZPpLyii4bAaAFhc3qKU47m6ynJM5Uzx1mttXMdqu0aft5SSsk\nXSfpmbgburj7sJkVgoHfFeT8gfbCrJ7mEmu2T7B94x4zu06l3vuwu19MerEgrSNVpXsAtDZ319LF\nnbo8Xpz1NWb1pCPpHr4XVRr4BYDYvnDwFV0eL6ojZ5piVk9TiLXIK0wwBRQA5vQXz57Qnx14RVvW\n36TP3fcB9eS7ZZJ68t3auflOBnpTErvnb2b3qTTf31Wq73OLpFvr1C4AGVae0jkS5PP7fjGvXfd9\nQB0508fX35Ry6yAl6/kPStou6aSkXZK+XpcWAci08pTOkYqB3B++/pb+7tjpFFuFaonSPu5+tPTJ\nD6rU8weAGcKmdI5NTjOls8kkCf5mZpuDz78t6fo6tQlAhjGlMxtiB393v0elmT7bJb03+AwAM6y4\nZlHocaZ0NpekaZ+Tkq539x1BCggA3uHuunbx7HkkTOlsPrGDv5ndZ2bnJe0zs3Nm9vE6tgtABn3z\n5TM6ee6Ktqy/iSmdTS7JIq8d7r5ckswsr9KOXn9Tl1YByBx3158+/bJ68t367x+/U4s6F7yMCA2Q\n5N25UH7g7oXK5wDw9A/f0IsjF/UHA7cS+DMgzmYum4OHBTPbr1KP/24R/AEEpqdd/+OZl3XLDUu1\nmfROJsRJ+6wNPpdLMJtKO3oBgCTpH178mX70+iV94ZN3qbODXn8WxCnpvCfsOLV9ADx5dES79/9I\npwtj6syZpqfYgykrktT22axSbZ+yNaK2D9C2qnfmKk67PvPk/5XljJk9GZBkts+DknYEnx9X7bZ3\nBNDkKvfeLe+xu/upH80q4zA6OaU9+48T/DMgaT3/o2bmwRaOg/VqFIB0hAV5STN6+COFUf3h117Q\ndESGhzIO2ZAk+FPbB2hh1WmckcKotg4dU4eZxorTM86d9tLMj7D4TxmHbEia9pGkg5IeFrV9gMyq\n7uE/dM9t2vlPs9M4k1OuydAQXwr83V0dM76HMg7ZkWRO1tPuftLdL1LbB8iuynr7rnIa55jevDSe\n6HXKZRso45BNSXr+B83seVXM8Xf3h2vfJAD1FFZv3yWZSR7Syc93d2m8OB3aw793XQ/BPqOSBP9n\ngg8AGRY1IOsensb57MfukKRZA8EE/WyLHfzdfca2jRVlHwBkxA9eK0SO1PYEQT0qyBPsW0uiqZ7B\nqt4HJT2gUrmHJxZyUTPb5u67F/K9ABbm8KkL+tSXv6/r39Wlt8enNF4xg4c0TvuJU9httaRPS9oi\nqSDpOkm3uPvFhVzQzAZUWiBG8AfqrHJWj1TaZetv/8tH9P2T50njtLk4Pf9hSXslDbj7q2b29EID\nP4DGqZ63L0mXx4r6/snz9PARa6rnDkkbJD0eLO5acOUmM+tzdyqCAg0QVn5hrDitPfuPp9QiNJM4\nVT13S9ptZreolO/fENT13+vuSXP+yxfQRgDzqF60dX//TTp9cSz0XMovQEqwyCtY4FXeynGPpE8m\nuVCcXr+ZDZrZITM7dObMmSQvD7StsEVbf3bgFZmFn0/5BUgJZ/uUBUE8afpmjZmtUan3vzz4ZXCk\n6nX3SdonSf39/RQGB2IIW7QlScuWdGqi6JRfQKiGbbnj7kPuPhQ8zTfqukCri0rjvDVapPwCIi2o\n5381Knv3AK7OSGFUnR2myZAdtFblu5nVg0gND/4AFq5yYPf6pYs0OlFUTlJX1S8A0juYD8EfyIjq\nefvn356QmfSZX/slrbx2MYu2kAjBH8iI0GqcLv2v77yq53ZsJNgjkYYN+AK4OlEDu8zbx0IQ/IEM\n+Prh1yKX1jNvHwtB2gdoQpUDu0sXd+ryeFG33rhUP70wqrHJ2dU4gaQI/kCTqR7YvTxeVEfO9Olf\nXauOXI6BXdQEwR9oMmEDu1PTrs8/8woDu6gZgj+QouqCbL/1kdUaYWAXDUDwB1JSnd4ZKYzqT/7+\npcjzGdhFLTHbB0hJdEG2DnV3dcw4xsAuao3gD6QkKo1zaWyKgmyoO9I+QANU5vZvXLZYq67rnnPe\nPgXZUG8Ef6DOqnP7b7w1rjfeGte/WnWtfnzmbebtIxWkfYA627N/9l66knThSlGPbP4A6R2kgp4/\nUEdnL49rpBC9ly7pHaSF4A/UUGVuP/+uLk0UpyPPZeom0kTaB6iR6o3UL1yZ1JXJKf3GB36BqZto\nOgR/oEai6u0f+clFpm6i6ZD2AWpkrnr75PbRbOj5AzXSvagj9Di5fTQjgj9QA1/57qu6MjGlzpzN\nOE5uH82K4A9cpe/8+Kw++3c/1Mb33ajd9zFvH9lAzh9YgMopnTJp5TWL9YVP3qVrl3Rp8/qb0m4e\nMC96/kBC1VM63aWLo5M6+NKbaTcNiK2hwd/MBoOPXY28LlBLYVM6x4vT2rP/eEotApJrWPA3swFJ\nB9x9n6Q1wXMgc+aa0glkRSN7/msklQP+cPAcyJQTZy4rVzWjp4wpnciShg34Bj3+sj5JjzXq2sBC\nVQ7srrhmkd4eL6q7K6fJKdd4kVLMyK6Gz/Yxsz5JR9z9SMjXBiUNSlJvb2+jmwbMUF2H/+zlCZmk\n//brv6QV1yyesfH61k23M6UTmWLuUfsJ1emCZtvcffd85/X39/uhQ4ca0SQg1Icf+YZGQvL4Pflu\nPbdjYwotAuZmZofdvT/OuQ2f7VMO/Az4otmFBX6JgV20hoalfYJgv8vMtktaLun+Rl0bmE/1Hrvv\nXrY48lwGdtEKGjnge0DS9Y26HhBX1B67H7xpmY6/cZk9dtGSWOGLthe2aEuSzl6eZI9dtCxq+6Ct\nVKZ3VuW79VsfXj1nbp86/GhVBH+0jer0zkhhVH/yDy9Fnk9uH62M4I+WVN3D37rp9sj0zrIlnZqc\n8hlfI7ePVkfOHy2nuurmSGFU24aORaZ3Lo0V2WMXbYeePzItbg9/Yip6MeOqfDe5fbQdgj8yKyyH\nv23o2JyBvrurg/QOIII/MqIWPfyeiu+jJg/aXcNr+8RFbR+UVffwJakzZypOJ+vhk8dHq0tS24ee\nP5pKWA9/11M/mtXDnyvw08MH5kfPH6kIC/KSZvXwcybNEefp4QMV6PmjacQJ8iOFUW0dOqac2YwN\nUqRS4DdJYfGfHj6wcPT8kUhYML93XU/snnxXh6nDTGNVQX4+9PCB+SXp+RP8ESpuMO/u6tB963s0\ndPi1GdUvO3OmXM40kTDIh6GHD8RD8EdscYP8os6cunKmtydml0eolXx3l8aL0/TwgQUi549Z4ube\nH3r8mHI2e778RHFaEzVqS1SQ/+zH7pAkevhAA9DzbzFxe/KdOZOZNDnHoqi4Osw0FfJzNFdPXiLI\nA7VG2qcNJBlgNatN7j0qmN+3vkdfPzxCkAdSRtqnhcRN12wb+oE6ctLo5MwgX+rZJ/sFv5C0TP/N\nyyODPMEeaD70/JtYWFmD7q6cujpyemuseNWvT1oGaC30/DOqupd/ZaI4q6zB6OT0rN79fBY6wEqw\nB1oXwT8FcVM5SRHkAcRF2qcGarHqNay0QRTSNQDCtO1snyRBuFbHpehVr9UzYBZ15NSRs9B9ZOMi\nyAOI0rTB38y2SCpI6nP33XOdmzT4hw+ORk9DnPv4azPy6ku6crr3rh49+cLIjBIGiztzWtSZ06Ua\nDL5GyXd3aeniToI8gHk1ZfA3sz5Ja9x9yMwGJR1y9yNR5ycN/h9+5BuhefKoipDNhtIGAK5WkuCf\nq3djKjygUq9fkoYlDdTyxU9HDJCmFfg7zEKP57u71N3VMeNYeVB25+Y71ZPvlqlUzIzAD6BeGjnb\nJy/pfMXzFbV88VX57tCef4dJYRUMokoSJD2edNUrM28ANIOmmuoZpIMGJam3tzfR927ddHsNc/7x\nj7PqFUAWNTL4FyQtDx7nJZ2rPsHd90naJ5Vy/klevBxMkwThWh2vvH51mwjyAJpRowd8+919n5lt\nk3SglgO+ANDumnLAtxzozWxAUmGuwA8AqK+G5vyDtA4AIGWNnOoJAGgSBH8AaEMEfwBoQ01b2M3M\nzkg6lXY7Itwg6WzajWiQdrpXifttZe1wrze7+8o4JzZt8G9mZnYo7nSqrGune5W431bWTvcaB2kf\nAGhDBH8AaEME/4Vpp/UK7XSvEvfbytrpXudFzh8A2hA9/xBBHaLK59vMbEtQdVRm1mdmbmYngo+9\nwfEtZjYQ1C7KhPnuNThWvq+wY5m5Vyn2/c71b5Cp+20XZjYYfOyqOBb6noX8DLTle0vwrxLUHnq8\n6rncfUjSWjNbI2m5u5u7r5V0v6Rd5R8odz8gqVD9A9aM4txrcB/DwX0NB7/4MnevUuz7jfo3yOL9\nxgqIcY81q+A9OxCUj1kTtDv0PQv5Gcjke1sLBP8q5SBXcejuiucnJA0E55T1u/uw6rxTWT3Eudfg\ncTl4rAkK8mXuXqXY9xt2LHP3Gzcgxj2W0m3EtUY/f0+Gg+eh71nIz0Dm3ttaIfjP75xm7kOwtvyF\n4D/Y1yq+Vredyhpk1r0GwX7YzC7o5/fXCvcqhb+3YceyeL9xA2LcY03L3fdVFI3sk3RI8d+zLL63\nNUHwn9+Qfh7wy8Gh7G53L8z+lsyada9mllcpEOyU9GiQ9moVYe/tXO93ZiQIiHGPNb3gL5QjlIuP\np6m2cWxG7j5sZo8FP1gFzfyTsfLP4Xl3Kmt2Efc6KGmnuxfMbFjSFrXAvUrh9xvxb7BCGb3fyoBo\nZmk3p94G3H178Djuz2hL/CwvBD3/eVTsQHZEUj4YCFRID/gxlf60VvD5gDIm6l7LgucFtcC9SuH3\nG/FvkOX7nS8gxj3W1Mxs0N13B48HFP89y/J7e1UI/lXMbIuk/uBzeQey88HzvVWnv/NXQBZ3Kotz\nr8F/qMHy1McgnZC5e5Vi32/UsSzeb5yAGPdY0wrubVcw7fqCFP2eRfwMZO69rQUWeQEtqGJK43mV\nevH3u/uBYO3CsEozt/YF58Y6htZC8AeANkTaBwDaEMEfANoQwR8A2hDBHwDaEMEfANoQK3yBOQQV\nMfPu/mDFsROS/o+kSypNh9xQsZAKyASCPzAHd98eLB7qC0okbFOpWuYfS+8sGgIyh+APzO9+SY+b\n2XpJD0paHwT9DZKeT7VlwAIR/IF5BD3+A5JOSvqdoMhd2s0CrgoDvkA826V3itsBmUfwB+I7P/8p\nQDYQ/AGgDVHYDQDaED1/AGhDBH8AaEMEfwBoQwR/AGhDBH8AaEMEfwBoQwR/AGhDBH8AaEP/HzZB\nADr8PjBiAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1086f15c0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import scipy as sc\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pylab as plt\n",
    "\n",
    "df_arac = pd.read_csv(u'data/arac.csv',sep=';')\n",
    "#df_arac[['Year','Car']]\n",
    "\n",
    "BaseYear = 1995\n",
    "x = np.matrix(df_arac.Year[0:]).T-BaseYear\n",
    "y = np.matrix(df_arac.Car[0:]).T/1000000.\n",
    "\n",
    "plt.plot(x+BaseYear, y, 'o-')\n",
    "plt.xlabel('Yil')\n",
    "plt.ylabel('Araba (Millions)')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from itertools import product\n",
    "\n",
    "def Error_Surface(y, A, left=0, right=1, bottom=0, top=1, step=0.1): \n",
    "    W0 = np.arange(left,right, step)\n",
    "    W1 = np.arange(bottom,top, step)\n",
    "\n",
    "    ErrSurf = np.zeros((len(W1),len(W0)))\n",
    "\n",
    "    for i,j in product(range(len(W1)), range(len(W0))):\n",
    "        e = y - A*np.matrix([W0[j], W1[i]]).T\n",
    "        ErrSurf[i,j] = e.T*e/2\n",
    "    \n",
    "    return ErrSurf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'df_arac' 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-4a9af2fc1034>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0mBaseYear\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m1995\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmatrix\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdf_arac\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mYear\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[0mT\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mBaseYear\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m \u001b[0my\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmatrix\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdf_arac\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mCar\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[0mT\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m1000000.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0;31m# Setup the vandermonde matrix\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'df_arac' is not defined"
     ]
    }
   ],
   "source": [
    "BaseYear = 1995\n",
    "x = np.matrix(df_arac.Year[0:]).T-BaseYear\n",
    "y = np.matrix(df_arac.Car[0:]).T/1000000.\n",
    "\n",
    "# Setup the vandermonde matrix\n",
    "N = len(x)\n",
    "A = np.hstack((np.ones((N,1)), x))\n",
    "\n",
    "\n",
    "\n",
    "left = -5\n",
    "right = 15\n",
    "bottom = -4\n",
    "top = 6\n",
    "step = 0.05\n",
    "ErrSurf = Error_Surface(y, A, left=left, right=right, top=top, bottom=bottom)\n",
    "\n",
    "plt.figure(figsize=(10,10))\n",
    "#plt.imshow(ErrSurf, interpolation='nearest', \n",
    "#           vmin=0, vmax=10000,origin='lower',\n",
    "#           extent=(left,right,bottom,top), cmap='jet')\n",
    "\n",
    "plt.contour(ErrSurf, \n",
    "            vmin=0, vmax=10000,origin='lower', levels=np.linspace(100,5000,10),\n",
    "            extent=(left,right,bottom,top), cmap='jet')\n",
    "\n",
    "plt.xlabel('$w_0$')\n",
    "plt.ylabel('$w_1$')\n",
    "plt.title('Error Surface')\n",
    "#plt.colorbar(orientation='horizontal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Animation of Gradient descent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "all the input arrays must have same number of dimensions",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-4-f759bc6fe2b2>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[0;31m# Setup the Design matrix\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      9\u001b[0m \u001b[0mN\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 10\u001b[0;31m \u001b[0mA\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mhstack\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mones\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mN\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\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     11\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     12\u001b[0m \u001b[0;31m# Starting point\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/cemgil/anaconda/lib/python3.6/site-packages/numpy/core/shape_base.py\u001b[0m in \u001b[0;36mhstack\u001b[0;34m(tup)\u001b[0m\n\u001b[1;32m    291\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0m_nx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconcatenate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marrs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    292\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 293\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0m_nx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconcatenate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marrs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    294\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    295\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: all the input arrays must have same number of dimensions"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pylab as plt\n",
    "\n",
    "import time\n",
    "from IPython import display\n",
    "import numpy as np\n",
    "\n",
    "# Setup the Design matrix\n",
    "N = len(x)\n",
    "A = np.hstack((np.ones((N,1)), x))\n",
    "\n",
    "# Starting point\n",
    "w = np.matrix('[15; -6]')\n",
    "\n",
    "# Number of iterations\n",
    "EPOCH = 200\n",
    "\n",
    "# Learning rate: The following is the largest possible fixed rate for this problem\n",
    "#eta = 0.0001696\n",
    "eta = 0.0001696\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.gca()\n",
    "\n",
    "plt.plot(x+BaseYear, y, 'o-')\n",
    "\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "\n",
    "f = A.dot(w)\n",
    "ln = plt.Line2D(xdata=x+BaseYear, ydata=f, linestyle='-',linewidth=2,color='red')\n",
    "ax.add_line(ln)\n",
    "\n",
    "for epoch in range(EPOCH):\n",
    "    f = A.dot(w)\n",
    "    err = y-f\n",
    "    \n",
    "    ln.set_xdata(x)\n",
    "    ln.set_ydata(f)\n",
    "    \n",
    "    E = np.sum(err.T*err)/2\n",
    "    dE = -A.T.dot(err)\n",
    "    \n",
    "#    if epoch%1 == 0: \n",
    "#        print(epoch,':',E)\n",
    "        # print(w)    \n",
    "        \n",
    "    w = w - eta*dE\n",
    "\n",
    "    ax.set_title(E)\n",
    "    display.clear_output(wait=True)\n",
    "    display.display(plt.gcf())\n",
    "    time.sleep(0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "all the input arrays must have same number of dimensions",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-5-b067099a950a>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;31m# Setup the Design matrix\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mN\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m \u001b[0mA\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mhstack\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mones\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mN\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\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      6\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0;31m# Starting point\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/cemgil/anaconda/lib/python3.6/site-packages/numpy/core/shape_base.py\u001b[0m in \u001b[0;36mhstack\u001b[0;34m(tup)\u001b[0m\n\u001b[1;32m    291\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0m_nx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconcatenate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marrs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    292\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 293\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0m_nx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconcatenate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marrs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    294\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    295\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: all the input arrays must have same number of dimensions"
     ]
    }
   ],
   "source": [
    "# An implementation of Gradient Descent for solving linear a system\n",
    "\n",
    "# Setup the Design matrix\n",
    "N = len(x)\n",
    "A = np.hstack((np.ones((N,1)), x))\n",
    "\n",
    "# Starting point\n",
    "w = np.matrix('[15; -6]')\n",
    "\n",
    "# Number of iterations\n",
    "EPOCH = 5000\n",
    "\n",
    "# Learning rate: The following is the largest possible fixed rate for this problem\n",
    "#eta = 0.00016\n",
    "eta = 0.000161\n",
    "\n",
    "Error = np.zeros((EPOCH))\n",
    "W = np.zeros((2,EPOCH))\n",
    "\n",
    "for tau in range(EPOCH):\n",
    "    # Calculate the error\n",
    "    e = y - A*w    \n",
    "    \n",
    "    # Store the intermediate results\n",
    "    W[0,tau] = w[0]\n",
    "    W[1,tau] = w[1]\n",
    "    Error[tau] = (e.T*e)/2\n",
    "    \n",
    "    # Compute the gradient descent step\n",
    "    g = -A.T*e\n",
    "    w = w - eta*g\n",
    "    #print(w.T)\n",
    "    \n",
    "w_star = w    \n",
    "plt.figure(figsize=(8,8))\n",
    "plt.imshow(ErrSurf, interpolation='nearest', \n",
    "           vmin=0, vmax=1000,origin='lower',\n",
    "           extent=(left,right,bottom,top))\n",
    "plt.xlabel('w0')\n",
    "plt.ylabel('w1')\n",
    "\n",
    "ln = plt.Line2D(W[0,:300:1], W[1,:300:1], marker='o',markerfacecolor='w')\n",
    "plt.gca().add_line(ln)\n",
    "\n",
    "ln = plt.Line2D(w_star[0], w_star[1], marker='x',markerfacecolor='w')\n",
    "plt.gca().add_line(ln)\n",
    "plt.show()\n",
    "\n",
    "plt.figure(figsize=(8,3))\n",
    "plt.semilogy(Error)\n",
    "plt.xlabel('Iteration tau')\n",
    "plt.ylabel('Error')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "* The illustration shows the convergence of GD with learning rate near the limit where the convergence is oscillatory.\n",
    "\n",
    "* $\\eta$, Learning rate is a parameter of the algorithm\n",
    "\n",
    "* $w$, the variable are the parameters of the Model \n",
    "\n",
    "* $y$: Targets\n",
    "\n",
    "* $x$: Inputs, "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Accelerating Gradient descent\n",
    "\n",
    "## Momentum methods, a.k.a., heavy ball\n",
    "\n",
    "\\begin{align}\n",
    "p(\\tau) & =  \\nabla E(w(\\tau-1)) + \\beta p(\\tau-1) \\\\\n",
    "w(\\tau) & =  w(\\tau-1) - \\alpha p(\\tau) \n",
    "\\end{align}\n",
    "\n",
    "When $\\beta=0$, we recover gradient descent.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/cemgil/anaconda/lib/python3.6/site-packages/matplotlib/font_manager.py:1328: UserWarning: findfont: Font family ['sans-serif'] not found. Falling back to DejaVu Sans\n",
      "  (prop.get_family(), self.defaultFamily[fontext]))\n",
      "/Users/cemgil/anaconda/lib/python3.6/site-packages/matplotlib/font_manager.py:1328: UserWarning: findfont: Font family ['sans-serif'] not found. Falling back to DejaVu Sans\n",
      "  (prop.get_family(), self.defaultFamily[fontext]))\n",
      "/Users/cemgil/anaconda/lib/python3.6/site-packages/matplotlib/font_manager.py:1328: UserWarning: findfont: Font family ['sans-serif'] not found. Falling back to DejaVu Sans\n",
      "  (prop.get_family(), self.defaultFamily[fontext]))\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXt4VdWZ/7/rJIQQbrmQAAGSkABF\nBYUQ0LReUIPWqr2MUNuhdUY7jXbs9TdTW0un03kqdqCdaqdjW2mrrRXrBe+9iCQaFQlKoKCIAklI\nwk1uSRAJCUnO+v3x7sXa52SfcwLkdna+n+fZz1p7r5199tnZ57vf/a53vUtprUEIIST+CfT3CRBC\nCOkZKOiEEOITKOiEEOITKOiEEOITKOiEEOITKOiEEOITKOiEEOITKOikT1BKFSqlNiqlmpRS95/J\nft05htOeH+XYTd38rFJnu1ZKPaGUSu1OW29wGteuVClVE+H7bHTO1yxP9OY5k/6Bgk76inIA9wOY\nDKBIKVV6BvtFbFNKlTgiVhjlHH4DwC2+nsdzHgj3A1gEIA1AI4Blsdq8cMS4UCm1UClVEuXcohHz\n2jnHXgbgVgBzAKQrpdznle9sL3CWL5/huZABDAWd9DqO2DRqrVdorZsBfAciPN3erxvHiCbkUErd\nAaC2m+fUaEqnDRBBjNUW/pmpAKC13qS1XgV5CJwW3b12kO//uNa6TGtd6+y30NWe6pxHrbM0exyD\nxDmJ/X0CZFCQD2CTa70K3iIYbb+ox9BaLwcApdRnww/qWNW3AlgAK3IRj6e1blZK3Qpgo1LKtC+K\n1eZBkda6zDmHhQDOxM3R3WtXC+BOpdQyR9BvNfs5379ZKbXGdbwvU9T9By104olSKtXxyRY61u0p\ni/MMSIW1bN3bTme/7h7DiycgAuf++4jHcwRwGYA5WmsFsXZrY7V54biC7gAw14j7adKt7+28ATwO\noEYpFZ6gKdU5x2WQhxpwZg8XMsChhU4i8YTWegEgIub4bWsBhIiSI/JdrGIXZQCaAaR34zOj7dfd\nY4TgnHeV1ros7IEU7XgLIe6LTYBY/0qpGojvOVqb+3MLnc9tBlCmlLpDKZXqtop7+NpBa30rrIsq\nH0Cps30TxH9uPvfLAJq8jkHiGwo66YLjHljj2tQM4EZ4dKQ5ArUixvFqEeomKEKoG8EQbb/uHiOc\nBQAWul0xTqTLN7p7PEd4o/nJvdrSjeg7zAWwCnItAfTstXMeIDdqrb/jbCqBfasodD7PPISaXe4i\n4iMo6MSLfHR1IzR6+Vy7Y2U61nG+UqrEcTt8B8BjrmMsjLVfrGNE4cuwD6J8ABsBTHZE7ecRjrcK\n4iN/AuKzXgYrotHa3BTCeZtxdY6GXNOevHaQ/9cdSqkNzvl8x1kM5UqpK2FdL6uifC6JV7TWXLiE\nLBCf6x0QUSpxyvsBFJ7FMQsB1EBe9e8Pa9MASrqxX8Q21z5NAPKjfC/dzXNa6LRpyNtKfnfaXPuU\nOteuxKmn9sG1W+js0wTgDo/zMef8xNmcD5eBuyjnn00I6SGUUoU61N1CSJ/QLZdLtBvUNcihQFv/\nHSGDmdPuvCWkJ4gZtugMbPAMcXLayrTWKwDkn8VIOEJ8gz6z8ERCzpqYgu7cnJHibPMhPkKgaxQC\nIYSQPuSsolwcy9xQiO5FHRBCCOkFeiRs0Ylz3eTlZ3d87KUAMHz48DnTp0/viY8khJBBw8aNGw9r\nrTNj7ddTceglkTpEHSt+BQAUFRXpqqqqHvpIQggZHCil6ruz3xnlcgnPDa1tYiR2ihJCSD/RnSiX\nhZAczO5UnOVOWwmAZSapfi+dIyGEkG4Q0+WiJYvbqrBtc5yyDJLknxBCSD/D9LmEEOITKOiEEOIT\nKOiEEOITKOiEEOITKOiEEOITKOiEEOITKOiEEOITKOiEEOITKOiEEOITKOiEEOITKOiEEOITKOiE\nEOITKOiEEOITKOiEEOITKOiEEOITKOiEEOITKOiEEOITKOiEEOITKOiEEOITKOiEEOITKOiEEOIT\nKOiEEOITKOiEEOITKOiEEOITKOiEEOITKOiEEOITKOiEEOITuiXoSqnCKG0LlVIlSqk7eu60CCGE\nnC4xBV0pVQLgiQhthQCgtS4D0BxN+AkhhPQuMQXdEevaCM03Amh26rUASnrovM6YlStXIi8vD4FA\nAHl5eVi5cmV/nxIhhPQJiWf596kAGl3rGWd5vLNi5cqVKC0tRUtLCwCgvr4epaWlAIDFixf356kR\nQkiv46tO0SVLlpwSc0NLSwuWLFnST2dECCF9x9kKejOAdKeeCuBI+A5KqVKlVJVSqurQoUNn+XHR\naWhoOK3thBDiJ85I0JVSqU71MQD5Tj0fQFn4vlrrFVrrIq11UWZm5pmdZTfJyck5re2EEOInuhPl\nshBAkVMaygFAa73J2acEQLNZ7y+WLl2KoUNTQralpKRg6dKl/XRGhBDSd8TsFNVarwKwKmzbHFd9\nRS+c1xmxePFiVFUB9967BEo1ICcnB0uXLmWHKCFkUHC2US4DjhkzFgNYjLo6gJ4WQshgwldRLgDQ\n6ARRpqdH348QQvyGLwV9yBBg+PD+PhNCCOlbfCfoR44AGRmAUv19JoQQ0rf4TtAbG+luIYQMTijo\nhBDiEyjohBDiEyjohBDiEyjohBDiE3wl6G1twPHjFHRCyODEV4JuBhVl9GtWdkII6R98Kei00Akh\ngxFfCfqrr0r5/vv9ex6EENIf+CY5V2Ul8M1vSv2OO4ChQ4HDh4ErrgCKi/v33AghpC/wjaBXVADt\n7VJvbwduu03qQ4cCL79MUSeE+B/fuFzmzwcSnceTO49LWxvw0EP9ckqEENKn+EbQi4uBz39exPzS\nS/v7bAghpO/xjaADQEoKMGYMcP31dltSEnDTTf13ToQQ0lf4StDNKNFUZwrrMWPEt07/OSFkMOA7\nQc/IAPbvl/W5cynmhJDBg68E/cgRsdCNoHPEKCFkMOErQTcuFzOwiNPQEUIGE74UdGOhp6T07/kQ\nQkhf4htBb28Hjh0LFXRa6ISQwYRvBL2pScq0NGDXLqnTQieEDCZ8I+gm02JiIqC11GmhE0IGE74T\n9LY2u62iQpJ2EULIYMDXgv7ss8CVV1pRr6wEfvxjijwhxJ/4Jtuil6AHg8DJk2KpA5JK9+RJycBY\nXs5BR4QQfxHTQldKLVRKlSil7ojRXtrzp9d9jKDX1NhtgYDkcpk/X0S9tVVEvq3NijwhhPiFqIKu\nlCoEAK11GYBmsx7WXuu014a39yVG0FeutNv+4R+sJX7ypN0eDHIUKSHEf8Sy0G8E0OzUawGUeOyz\nzCnztdabeurETpfGRmDYMKCjw27Ly7Nulffes9sDAUkTQAghfiKWoKcCaHSth9i1joDXKqWawvY7\nhVKqVClVpZSqOnTo0FmdbDTMKNGA6xtdfrmtuyeOHjpU3DCEEOInzirKRSmVCrHgfwzgN0qp/PB9\ntNYrtNZFWuuizMzMs/m4qDQ2AhMmAFOm2G1XXuk+DymTk9khSgjxJ7GiXJoBGNs2FUC4o6IUwI+1\n1s1KqVoACwEs79lT7B4mde6ePbKelCSWuKG+XsphwyjmhBB/EstCfwyAsbrzAZQBpyzzELTWq2D9\n7X1OY6MM+z94UNZHjQptN4Lunm+UEEL8RFRBN52cSqkSAM2uTs9yp305gFIndLFUa72iV882Co2N\n4j83naJuQdcaqKuTemsrBxYRQvxJzIFFXiKttZ7jqveLi8VNZyfQ3Bwa4eIW9MZG4Phxqbe0iG+d\nfnRCiN/wxdD/o0fFCncL+siRtm7cLQb36FFCCPELvhB0M6jICHpiYqiFHu4/N6NH3axbB9x9N90x\nhJD4xRe5XMIFPTk5VNCN//yCC4B9+4Bnngl1t6xYAdx2m1j5w4bRHUMIiU98aaEnJYUK+vr1wJAh\nIvQ5OaFiXVkJ3H67jVNnnhdCSLzSt4Le2Ai8+qqYzO3tPXpYQAQ9IwM4ccIKemUlsGqVfNybb3b9\n2IqKUN97QgJHkRJC4pO+dbk0NwPf+x7Q0AC8/z6QlSUmc6QlLa1bgeNG0Nvb5U+OHLGCXlEhybgA\nKU20i2H+fBHxzk5Z/8Uv6G4hhMQnfSvo+fnA2rVS7+gQh3ZDg122bQNeeMGud3SECvykSV3Xk5Kw\nZYsc0iToAmyUy/z5oc+E8eNDT6m4GJg5E9i8Wda/+MVe+/aEENKr9F+naGKiFeZIfPABsHt3qOiX\nlUlZXw/s34+20Zm4+VAuLkcu6rfmoW1cLrKRi0kf5gItuZg3LwVKAZdeClRXS76XcMzoUgB45RXg\nmmt6/usSQkhvo7TpDewDioqKdFVVVc8dsKMD9y3Zh8eW1yEX9chFPWaMqEfGh3W4aHw9RjY2oHPE\nKGw6kou0Wbkor87FiPPysPh7uUCuLB8mpmLkSLHitZaO05deotuFEDJwUEpt1FoXxdovvsMWExNR\n+OkcfP2nOXjN8ZPPO1c6P1/8A7DgyiDeKTuIr11dh/+9oh77t9fj8uPvAStWS8dsfT2SkYDNyEW9\nlkdCQ1se9v8iF0h0RD8zE1AKlZXij58/n2JPCBmYxLegQ8T13HPFbXLwoM2wWF8PIBBAw8lxeAPj\ngBsvwvJfAcc/Dlz2E+ePtcafH2zCf32pHgWBOkwM1qMgUI/L9r0O3FonBzlxAi1ZuTi2OxcZOhcv\nJORizPfzMLXEEfzx46VXlRBC+pm4F3RAoltSU0XQ162TbV/7GnDeedLvCojunjhhO00BAErhr+vT\nsRnpKPjMbPz8SeDPzwIZ17r2OXYMd32hAVvqxa0zubMe7U8+B/ytTgS/sRGYOPGUCwe5uTJVkqlP\nnCiB8YQQ0sv4QtAbG21Uiwk/bG8XF0lbm/jH09Jku1vQKyuBBx6Q+nPPSXnRRaHHrtw6Esv/ch46\ncR4AYGgScPH9AIzbpbVVOm4dFw7q68UJb9b37wfGjg0Ve7NMniydwhR8QkgPEPeCrrUIenKyrJuY\ncpOv5YEHJNzdDChyC3pFhX0AuPPAuHHvAwA33xzmQ09OBqZOlcWLjg5g714RdyPy69cDf/qTrO/d\nKydoBN4t9nl5EpoZflKEEOJB3CvFhx+K4BrBvuEG4PHHgTVrRHjvugvIzhZ3CxAq6O4Y9cREOUa4\ndl58cej67NmneYKJidY6v/TSru1G8HftEoHftUtiJ3//e1k/cED8RUbgw8vsbPrwCSEAfCDo7mH/\nSUmifaNGAR/7mGzft0/c2EbQy8qA888XsZ8zRyz8K64Apk0Dfv3rroIePg3qN78pA5F6LNLFLfhe\nnDxpXTpG9FevtutHjogVH0nwx43jNE2EDBJ8JehZWTIWyZ2Ya98+YN48mxb3qaeAv/xFMiqOGyfb\nFi8WIxnoauy+/XbousmlHi7ovRbWmJQEFBTI4kVrq3XnGMF/9lm7fuyY9d17if6YMRR8QnyCbwS9\nrQ1ISRENGz1atrW3S+RLdraNfgkGu4pybq4MPgW6CvrWraJ3CQn2LSA8eVdlpWzr6JCwyT5Nv5uc\nDHzkI7J4cfy4iLtb8N98025ra+vqt3eXqakUfELiBN8I+okTsrz/vsSlA1IHRNAzMqQeCFhR3rFD\ntuXkiLYpJf2VbjF+5RX528suA558UqJhwsX67rvlIQFEtuD7jeHDJX7zvPO8248eFQvf7cN/7TUp\nd+2Si+IWeLPk50uZktKX34YQEgXfCLohGBSjE7BulOxsm+/8ttuAL3xBBLeszO63cqXs455vtLJS\nsv1qbcMai8IG31ZWigvHEAjEWfrd0aOlU+H887u2aQ00NVmh37VLnoIvvADU1sqDYPRoEXcj8O76\nhAnssCWkD/GVoBvPwKRJUppBRRMmWGv8X//VGqv19RIiXllpQxPdFnZ5uX0QmHaTitfw0kt2H6/2\nuEYpID1dlsLCru3BoMTZ19aK2NfWAi+/LLGitbXSYZuTE2rRuwXfDA4ghPQIcS/oTU22fvXVQFWV\n6AVgBT07G/j736U+fLjdv6FB/OcmJ3q4j9y4pQMBaQ8GQ2PSga7BKVr3cadpfxIIyNNywgTgkku6\ntp84IU/N2lor+pWVdj0x0Vvo8/PlwnLAFSGnRdwLuttC/+xnxf1rolzefFM0Z+dOoKVFtrldvvX1\nMs9ocTFw7bVikb/4ohVcc5x/+ReJd3/kEeCNN4DrrrPHMAOSAgER/EgTUF9xhX1gDJo5S4cNA6ZP\nlyUcrcWCd1v3mzbJ9FK1teIHGzs2sjtn7Fh21hIShq8EPT1dgjpGjRJD8JFHRGQXLABuuUX2MYKu\ntVjo119v/zYtLVRoa2qkvO46GbAEAIsWhabX/fOfRaRvvBH44x/FFx8u1t/9rvXrmzlLB4WgR0Mp\nCZkcM0biSsPp6JD4e7fgP/+8rbe0SEdtJMF3v4oRMkjwlaCbUaCjRoUO2T95Enj3XakbQT90SEK4\njcvk5Mmub/jV1bL/229bS9zkiDGdpk8/LQ+Nxx6T9vDowddfl7cGQzBoI25IFBITbUSNF8eOWXGv\nrZWn75o1Uq+rk5vAS+gLCthZS3yLrwTdpM4dPVqSbJlJK5KSZBBRcrK4RgAbd24mTPIS9Joa0YHL\nL7epAYYMsS6Vl1+2naBG8H/5SxtFAwBPPBF6TKXE00DOkpEjI0fnBIMSs+oW/FdfBR580HbW5uba\nAVvuZfLksJSchMQPvhT0UaNEUEeOBM45B7jnHuDhh0Pfwl98UcrmZinb2uzfG2pqJOdWcbEM+f/J\nT0SgjVibeHd3p+mvfy26YfzkQ4aEHtPLx/7KK7IsWEBXTI8QCEhPeHa2zQHhpqVFxL6mxlr3L74o\nZX29vEJ5iX1Bgfjm6LsnA5SYgq6UWgigGUCh1nq5R3shgHwA0Fqv6vEzjIER9GHDrItl1CgR6A8+\nkM7O4mJgxQrrbqmsBH74Q6l/5SuSxyXcQg8G5ff98Y/Lel6elHPn2n3M8W65BTh8GHjmma4jUbds\nkf0mTQK2b7dJwwz/93+Sux0AfvQj4L77gNLSHro4xJuUlMiDrTo7gT17QsX+6aelrKmRV75IYj9x\nIl05pF+JKuiOWENrXaaUyldKFWqtN4XtdqfWepFS6o4I7b1Ga6tNumXyuAAi6GaU6PjxUh4/bi30\nigrrIjHiGy7o+/fL8U0KFRNr7jbOjF/+rrtkcNEzz4SORK2slMFLWos/HhDXr6GyEvjGN+x6Rwfw\n1a/2cPIvcnokJNhkaVdcEdpmcjUboa+pkX/iww9L/fBh8eFFcuVwVC3pZWJZ6DcCWOPUawGUADgl\n2I71vgEAvKz33sYdg56ZKaPYAfGh798vdSPoLS329zR/vgivO2/6X/8aKujPPiulGdJvBN344AER\n9LQ0eZiYcTf/+I8yeKm4GPiP/+g6MOmGG2yUTEVF14FInZ2hUTC+jF+PV5QSd0xGRuirmqG11bpy\nzFJWZjtq09K8xT4/n0nSSI8QS9BTAbgH14fHZ8wFTlnyJRFcMqUASgEgx/RA9hBu/3lmZqiF/t57\nUvey0IuLxc3S3g489JCspx+rxzy1AfjFfuzefBgnf38M9yCIhH8D9r+ucfEB4H+QgJT/Hg5kpAAp\nKSh4KQVfHzMc6tkUpB5MQTFG4AtF6SjOTwNOpiErK8m5BlbY3VEyI0Z0/U6JidbHfv/94hLSWrbT\nHTPASU6WTptzzunaFgxKbL1b7J991lr7nZ0i7F5iz0lOSDfpibvkiNZ6k1KqRCm1MNyPrrVeAWAF\nABQVFWnPI5whbkEPd7l4WejukeZNTda/jnvuwaNvLcFbYxcA2yehvi4De4IT0IkAVFDhnRZgSIbC\nPrQDSSfkg/bvx7j6FmSOaEHj/xxHZrAFP8MxnLusCVjaBDQ14bZAEj6NdGB0GnYdTUOjTsMHOg0l\nb6Wj/l+z8NaKLHwCWTissnBAZ+IgshBU8hpRWQncfrt9ENAdE+cEAiLMkyZ5J/tpagoV+zfflFmt\namokZeikSSLwU6cCU6bYJS+va28+GbTEEvRmAOlOPRVAeMDdEYgrxuw7F0CfdYx6WehKiSW+f7/8\nhrKypL2lRcKPAXGjHDhgc75g4kScVMnYmzoD+H9fQsKBfNx7iRhNw5KBRUskzvznAK6aD1x1lcwx\ncdMvANUEJG8US3/ROuChZcAXvwhAa3x18YfYVN6EDS824e2VTfjDT5pw6z80IruwEeufP4SPdb6L\nTBxClj6ILBzEWBxAR1si2q/NQsuILDzVmQVpkeVQexZ2/XYsikePlwiO0aP5mu4X0tIk81t49jdA\nXDl1dSLu1dWyvPCClLt3i9XiFnmz5OczBHOQEUvQHwNg7rB8AGUAoJRK1Vo3Q8R7odOeCsef3leE\nW+h79oh1HgjYuZlN0IHb5bJ3r1i+pwR90SJc87U5+Nqhe9A+5yIUT5uMZWMWonzIx/Efj80AoPDL\nX8qun/60hCQ+/bSsay0PiA3ON3/mGfktFRcrlL85EoHRI1HZkoPJNwPP/AQYmgCMvhjAxcCXLrbp\nAjo6gGBQIzP5Q7xw/0HsfP0gHvm5lfPJ2IUL8QYuW/s+8Pp+SVTT3i7CPn68DdMzdXfJnObxTXJy\n5BQK7e0yqMIIfXW1+PSqq8Wfn5npLfYFBd4+PxLXRBV0x5VSpJQqAdDsimApBzBHa12rlGp2Okcz\n+rpjNNxC37bN5l/Zv9+6WwCJN9++XVwZZv5RI+iVlcD6g/lYj1+gdOg9eOMfX0bGnU/hD4mfQubC\nVrw1/ip8rv1KvIpLsKctFxUV6tRoz4QEEeQxY2T96aeBv/0N+OlPxaBSSlLy/uAH0v7445KK9557\n5G+nTwd++1vgW98CtmxRuOvnI1G4aCTuea4Az4V930AAuOufgTvvdDZ8+KF80f2OwJvy7bdD10+e\nDBX97Gz58hMnWjfA+PH008YjQ4ZYf/vVV4e2mRBMt9ivXy9lTY086COJvZklhsQVMX/Bjg88fNsc\nj/Z+i0EHrMvFLejZ2VKvrBQXZVWViOsdd8h200dbUWF91SfaE/F40wL814kF+PHdGt9dVIMRv16N\nT/39OSzX30ZnMBFDX7oEa3ExZidcghv+cwauKAngr3+VvzcW+8MPh66bgUxaS4z8V78qVnl1tejv\nxo2ybuYsNZki3R2q7g5TAGJhTZ0qSzSOHw8V/r175Ye+fr28su/ZI37arKxQkTd1U44bR9GPJ9wh\nmFdeGdoWDMq94Bb7xx+39ZQUb7GfMkUGV5EBSVz/OsNdLg0NEuddWSn36hznsfPyy1IacTXzixoL\n3YikUmJtT5ki6/kFCpgyBfk/nYJ7T96Ohb/QWPv7Gnws+BpG/HAtnkr4OfJ+egB4eQ7GZ83DNsxF\nFebi0JBJmD5dobLSWvCLFsl5KGVDJgEpn3wyNO/MQw/Zcw4E5LyDwbPwmgwfbn+MkWhvl+D93but\nyDc0SKpIs+3wYfFjeYn9xIkiHGPHhsZ2koFJICD/s4kTu3bSai2dTG6xf+45KXfulJs6kthnZtK9\n14/EtaC749B37xYr18w6dOKEdblcfLGURrCHDxcjw8SlFxdLfdYscZUcPCjbTV51wIwUVZj5mSnA\nqCn42vKbcc61wFP3HwKqqjChcgP+Gb/Hb5Jux4hhGtWvzkNOYC6mL56Lgs/NxdxrxuAb35DPWrxY\nYtVNHPwNN4iAm1wxra3WKnfHqXd09GKmxiFDrGUeifZ2eVLu2WNFv65Oso/t3i0PgKNH5dXHWIZm\nycuTcsIEWvkDHaXkbWzcOPvjMZi0x26xX71aYmqrq+UmjST248ZR7HuZuP5luS30t9+2ImgGAxlB\nv+ACKa+5Bvj+92UO0HDd0lrSfhQXi38bCBV0N+3t4oL89KchFsk11yBw1TX45I+A//qexg9u3o0n\nP7kBuWM24HN7fgp8rgpITcWzahaCx2fh2qzZeDRnFhpULv74sEJxsXyXO++UUMrkZPmcQEC0T6nI\nE1T3KUOGWIGOREuLCHt9vSx1dfKDN+sHD8oP2y3y7iUnx14AMvBwpz2+6KKu7Y2NodE4FRXSSWQm\nJTAuQrNMmyYlB1b1CL4Q9JQUCSX84Q/lnkhMFOvXCLpJD3DddSLYu3d3FfT2dms47tolvnivGdK0\nlvb29tCgg4QEWdo7FJCTg5/vyUFu7g2YfhdQfGEQqK3Fk3M348I9m3Fk2W/xx11/x6jEFgRvvwC4\nbDYubJqF8zELf37qHHQEkhAIAF/6EnDzzXL8uBktmpISOSIDkKftnj1W4OvrJcfwI49Ifc8eeX0K\nt+5NKt28PAr+QMZMWeg1kvboURF2s5SVAb/6ldSDwVCBdwt+amrff484xReCnpUlQjdsGDB7NvCZ\nzwDf/nZXQTchubt3Ax/9qD2O1mIBm8yItbVinbsNBnc9PC2AISFBsrT+7GfyVtrYaCadDgCYggeO\nTsHv9ELc3gi0AxjbeRBF72zBfR/djNQ3V+NPWIY8XYftnR/B25iJ7Q/OQOr0mThn0QwUf3eiPyyY\npCSbn9yLzk7x5dfVWcHfssWOqmxokLeiSFPXjRtHH/5AZfToyLH2R47IxL9G7J9/3q4nJ3e16M3C\n0MsQ4lrQja972DB5wJ84AZSU2Igr9yhRQIzH48dFaN1ZCEyHpBH0rVvlPqms7GoRv/kmsGSJ1L/x\nDWDGDDvZxcmT4k5+/XVpN52wFRV2HbBhkwd0Fl7oXIBHJizAJSuAwkuAYWjBDGzFDGzFBZ1vI/lX\na4Dlb4tjfcYMCYGZMcPW/RZxkJBg5yn1Sn3b2SlROu55Sl94wdY/+ECseiP0brGfPFlyKpOBR0aG\n/JDCf3Cmg9Yt9o89Jusm9NJL7AsKBuWgqrgV9NdeE3EGJL78pZekPnKkHfY/bpyUbgv9OSe4u7XV\nHssI7JAhIsb19TZ+3OQ1r6uTff70J7u/O02uW7TNA8JEuBi/t5l31ODOzGg4gRRswDxswDwMTQLm\nPQRMLoZMsfTOO9JZ8NZb4qLYulV6eI3Im/Lcc/2b2S8hQZ7GOTneHQrHj4uwuye3ePllK/gpKd5i\nP2WKRHzQuh9YuDtoL700tM3kx3GL/dq1sl5XJxFXXv76yZN9OwF53Ar66tW2rrUILyCCXlMjfSzm\nf2YEfdcuG4P+3/8t4zCKi0ONPTKOAAAYvUlEQVQF3Qh+uHX9q19JuXJl6ExIRlPmz7fbTcqBr3/d\n+r0rK7t6TAIB4N57pf0rX+n6HW++2WWwZGbKwdwiprX4j95+W8S9vFwOuH27vJ6YRFHnnmvrfvdH\nDh9u32DC0Vpe69zW/bp1MhlsTY289k+e7B2hkZPD6JyBhjs/TnicfUeHuOfcYr9mjazv2SMPby9/\nfW5uXOe0j9s71D3zWEKCXR81SgzZhATrMjGCvm2b9Xu7QwCNoCcmyv8ZCLWe3fnTOzulrbgYWL7c\nCm5xsUTTNDeLLowfH9qJ6R68ZDARYF4kJAA33RTjIihlrdVrr7XbOzpErLZtkxy/r7wiUym9+674\nktwCbwR/7Fh/+OijoZR8z7FjvXuXW1pE6E2ExltvAU89JfX335fr7CX2eXm+tfjilsRE+xZmZqkx\nnDwpv4+dO0Xgt22TPpodO+SBn5cXatFPmybLhAkD/jcSt4Lujpy79VY7EcXevfLWZeLRy8utoF90\nEfC733UNATSCXl5uJ3n+93+XsETzuzdziiYmykjPL32pqyaMGiWGckcHsHlzqMtm/nz5W3dHqnvk\n5003ydR1bW0i5r/85VlEtCQmWqvjU5+y27UW6+Tdd+Umfust8Udu2yZPKi+LPicnri2WU6xcKZ0f\nDQ3ynZYulQEBblJSIlv3bW0iAkbsd+yQJPrV1XJNs7MjJ8hiVM7AIilJfujhM7oDIhY1NVbsq6rE\nvbljh6TamDrV/u20abY0Q9T7mbgVdPegorlz7UxAO3aExqNXVNjp4y66SObtXLdO8q0YwVy/Xsrn\nn7cTXfz4x9adaiar+MEPRHh/85uuv/nKSvG/G/+522Vj+npuuUVynGstD3q3S6W4WFy9vRqeqJR9\nRb3qqtC2Q4dE6I3Yr14t9SNH5Glpbl73jRwvHbIrV0oiedM7Xl9vE8uHi3okhg6NniCrvj50sM1L\nL0lZVyf+N/d1M3W/PCz9xLBhkR/qzc0iMDt2iFvzueek3LlTBN3rNzJ5cteJhXuRvhX0tjZRvB64\nid2CfviwfUAWF4vgmlGh8+fbqeKGDZPfdHhO8XXrpDThi16j12fNknLbNindU8kB3rMPhedeuekm\n4A9/sNPdhbtUvDr5+4zMTFnCO56OH5cbdvt2uZHXrJFRgdu3y43qJVRTpgysHN1LllgxN7S0yPbu\nCno0hgyJnFqhs1PE3oiAEYIdO+QhWlDgfQ0zwueSIf1Oaiowb54sbtyds+Z3UlYm9X37xJ3gJfa9\n4ObsW0HfsUN8uD1wE1dV2fqSJbaz07xFXXedjLwsLgY2OTkiU1LktxUeDWdE1ORN8XoTM0EjJiTx\n2mutOwUQ4U5IsL72cAvcfE55eRwNEgKkk3HWLPtEM7jDycxN/OCDNsIgO9v7fzxhQt9HkjQ0nN72\nniQhIbIv1zwszTUsLxdfm3lYul/p3Q9LunAGFtE6Z1tbxYVj/seVlcDvfy/rJ09a/7z7NzJtms31\nfZr0raDPnCkjb4zFt327PMncFp/5QuFxpWEDCN56y9aNzxoQNxcQKqbGh56UJO7O8JHrF14o5XXX\nyWTPXkJbUxO63tYWmleluFi+xle/KkbZ0KHenZr9aoX3JNHCydrbxd9sbuItWyST344dMlrQ/G/D\n/c3jx/dOp1NOjjzJvbb3J9EelgcP2gfl9u3yard9u31YGhGYPt32ezAx1sAjORk47zxZwmlsDP0f\nP/GElNXVYty6H+bdpO996CkpEg5iEqwYzE1sfFQ7dwKPPipldbX4a10if1XrNOzFVNSgAAlJQ0/N\nRmRcISZ1LmDftpuaRPzDBd14gEaMsFEs4YQLekJC1zDo0lJ5ZsWVBd4buB/M110X2vbBB/L/Nb7m\nV1+VnurqankaFxRYgXdPt3Y2lv3SpaE+dEDuw6VLz/w79ibuaJzwh6WJYDIG0aZN0kewbZv8XXin\n9jnniOXI+PqBR3q6t4UXDNqQS/N/7iYDp1PUfRNfckloWzAo8dbm9XTnTlxVV4FPYSdyUQ+dNgEN\nq6dihpqG8U9MxdWYhknt04BOiSmtrhYB/tvf5HCRBP3RR6X82c+A668Pvc4misbwrW95C7ZvLPDe\nYtSoyMO/P/ggNLHT+vWSWL66WqwZMwAofIk1ibLxk8eKcokH3BFM7oelMYjcHdt/+YvUjx4NteSN\n0BcUMLZ+IBIISCRHXp4NXrjvvm79aXz8NwMBm6ippAQA8J+7gGf3AFPz2rFj9S6U/7+d2Pf+Dsyr\nfQ//jucwfvFO4MABtIybjEUNU3G+noa3vzYN8zEVBcnTAJ196vXUCLrp1PRKU3vkSOjAIb+Pz+kX\nRo2SZDyzZ3dtO348NEZ882Zg1SqpHzgg90Z4uGB+vvwohg0T8Y5HAe8uboMo/NXx6FEr9O++K9kP\n331XOuwKCkTczzvPjjSeMoXRN3FKfAi6ByYxV1qWvN6vz5iGisxrcaQQeLoeOFAPoKUFj9xZg7/9\n705MxQ4Udb6Bz+OPKFi0Azj+4SlLJyV3Gm7CNMhe09AyJAPz54f6IufPF3eYiVDp1zS2gxGT4mDm\nzK5tra2hMeLvvSfxp7W14jtPT4+czCs72//uiNGjJWY3PN3tiRN2YM0778jMKlu3Su6M6dO75g6a\n6JMEcT4mbgXdhC1mZUl57JgYePv2ufznKSk473Mzcdt9M09FS44cCTQdREgqT7V1BxZgDW7HfZge\n2IGUIQqJ3wztnC2eOhUVz01F+YZRg9s/PhBJTrZuhHDMVGvu4f5lZbbe1CQumHChN3U/z605bJh3\nf9aHH4rIb90qaSVefFHKwZIgLo5ROnw8ei9SVFSkq9zxhmfBxIkS+nnzzcADD4gnpqVFLOixY8V9\naLj+ehn9fs458vvesCH0WK+9ZvueEgIalc8fxtzRO6zP3p0PYvTorsOCp04VAfBrQiw/09IikSPu\nZF5G7GtrJVzJK2tjXp48CAZSvH1v404QZ8TepCa94AIbsTN7trht/P7m04copTZqrT06nkKJWwvd\nnQsdEAs9NVXusXAX7JgxMlnF+++LFR+eFveRR2y9M6jwwPOZmPurzK4B68ba2+ES+7Vr5TV/1y4J\nG/OK0JgyhWI/UElJkY7Cc8/t2qa1iJhb7N98U3rP3ZNxeM28ZLb5KV93pARxDQ0Smrp5s6SSuPNO\nuW4zZ1qBnzVLLPpBmNK2L4lLQW9ttbHlmZlSHjsmrpYDB0JDFgEJRTx5UgQ9PC3uaeGeWPeKK7p+\nyO7ddiLd6moZhbRzpwiCCbsMF/yCAn/96P2EUmIxZGXZwQpuOjvF3+yebm/LFhkJarYNGxZ5ur28\nPLE04tkvrZT9Pp/8pN3e3GxFft06O9YkP18EvqhIcnbMnk1jpweJS0F3D/s3FvoHH4jIB4Pegn7o\nkNTDc6wANjFWpCH53SIhwYYaOZE4ISewd68VehOSV10tYXppaZEte07IMHBJSLAPeK/JOIyF755u\nr7ZW8ryY9c7OriJvRh1OnCg3czxmckxNBS67TBZDW5tE12zaJEO9H35YfPXTpom4z50rw+rPO69P\n85/4ibgXdGOhNzVJcAPgLegmaZbXpBK9nhjLPSlD+NBgkwfCCP3OnTKLxs6dIvajRnlb9lOmDJgM\nbyQCbgvfa45NQDrnjXVvRH7jRnHn7N4tr5UZGVbgvcrs7PgQwKFDrZ/9lltkW2urDPvesEHeaO+9\nV1w4558vUTkXXyyLsdxIVOKyU/T11+V/DMi939oaaiA98ICdXBkAFi4EnnxS6t//PvCJT8RJlEow\nKK/0bsveXR8xwlvop071d3TGYMLMsbp7txV5U5r6gQPSURRN9MePjw/RB+R1e+NGeYt97TVx2ZgB\nhxdfLBEMkeak9Sm+7hQ1HaKAPLh/97vQ9p07Q9eNuyUtDfjRj3r33HqUQMDOrxke+K61iL1b6J98\n0q6biXW9BD8trV++DjkD3HOsRqKjw4q+W/DfeCNU9NPTxZrPzhaB9yrHju1/4R81Crj8clnuvFMe\nalu3SgDC6tWybfhwmXLs6qulP4v9UAC6YaErpRYCaAZQqLVeHmW/O6K1Az1noT/0EPBP/yT11lZx\nSX7iE7b9tdesBV9ZKfVgUPRx7do4sc7PBpMJMdyiN/WkpFA/fUGBzaHCBE/+pKNDUgPs3y+RWuGl\nqR86ZIXfS/SNCykrS/p3+uNe0VrC2VavlmXDBhH2z31OYpT7+4HUC/SIha6UKgQArXWZUipfKVWo\ntd7ksV8JgAUAogp6T+G20IcODc1FP2aMFXMgNE95MCgPA98LujsTovtiALajzsTV19RIkpuaGlna\n2qzAhy+xcqaQgUtiorXO58yJvF9npwh/uOhv3iz3yaFD0n7woGTVzMoSI8At9F5LerpEs/TEA0Ap\n8bGffz7w7W9LXo6nnhL/+7/9G7BsGfDZz57958QhsX6dNwJY49RrAZQA6CLofY3pFJ08WUoTwgiI\nwLvjzE2ectMp+uCDEsXie1GPhLujzisyo7nZintNjVg/jz4q9YMHpWPXS+zz8xlj7AcSEsQSHz8+\n9r4tLVbg3UJ/4IBY0Gb94EER3WBQ3H3hS3q6rY8eLcI/fLiU4fWhQ+2UZKZsaZGH1KRJIuY33igT\nFpxhTvF4JpagpwJw2cPoMgOFY7WXKaW+06NnFgUj6OPGSenOirp3b2iceXEx8MUvSk55wDvxFnGR\nmio/Di8rzuRMcQt+WZmUdXUSjeF24bgXDg/3HykpNtyyO7S2yo/XLI2Noeu1tdIh2tIiy/Hjtm7W\n29rEKDGWvlLSXzRmjCyzZgF33z0oxRzomU7RqL9UpVQpgFIAyOmhCQWMy8WELLotdKBrnHlpqQxg\nY2KtsyRazpTOTul8c4v9k0/adLgJCZHFfjAkyCJy/3TX+idnRCxBb4YV7FQAR9yNxjqPdgCt9QoA\nKwDpFD3D8wwhPDGX20J3zyVqiMup3+KNhARrrYWPotVaJn51i/0rr0h8aU2NuHkmT/YW/Ly8+BxY\nQ0g/EEvQHwNgelbzAZQBgFIqVWvdDCBfKZUPEf30SJ2mPU00C/3zn5dp4MJFmxNP9CNK2Umow1O4\nApLdr7bWiv0778jw+epq8aGNH2+jcMIFnyNpCTlFVEHXWm9SShU5USzNLrEuBzBHa70KOOVW6bMp\nH3btktJY6G5B/9a3vCfDIQOYESNs1EI47e0yetJt3a9bJ2JfWyt/a7Ihhi8TJnCiBjKoiOlDd1wm\n4dvmeOzTZb/e4sABKY2F7na5dLd/hsQJQ4bYePlwTLy9O+3tq69KD3htrbh5IuU6z8/naFriO+Iu\nqNg9DsrLQh8zpm/Ph/Qj7nj7j360a/uJEzYhlllef93W3bnOwxfG3JM4JO7u2OPHbX3fPimNoGdn\nc5AjcTFsmEylNn161zYzwMot9pWVwMqVUj9wQPKgRBJ8pk8gA5C4E/QyV0zNbbdJ5s1IWRYJiYh7\ngJVXR21bW1fr/o03bD0hIbLY5+T4cvg5GfjEnaC//LKtnzwpQ/kfeEDWN24EVqyQuHNCzoqhQ8Va\nmData5vWMvLRPV1dVRXw+OOybia2jST46el8lSS9QtwJujskORiUJHNmWL/WErI4cyZDFEkvopQd\nmThvXtf29nbJ6e227letsqGZWttOWvdiJkgZpKMcydkTd4JuIlwAGVw4bpy8/XZ0yLbOTg7tJ/3M\nkCE2Tt4LM8zdWPfvvAP8+c92kotRo0JF3i36OTkcaEUiEneCnpdn60OHSqKt2bPFMu/slG0c2k8G\nNGlpkfPlmNfOXbtkqasT3/2jj8r6vn3i9w+37E2dsfeDmrgTdBOqmJcHPPKIHQE6cyaH9hMfEAjY\nNLde2TA7OiRnjlvwy8vt+uHDEp0TSfDHjqX/3sfEnaA3N0u5YEGocHNoPxkUJCZaX/vll3dtb20V\n/70R+F27gGeftfXjx2X0nZf/fvJkeXug4MctcSfoBw9KOXFi/54HIQOS5OTI0TkAcOyY+Ondgr92\nrVj6JqdGJP99Xh6nehvgxJ2g79gh5dix/XsehMQlI0fKFF/uab4MWkuHrXHl7NoFbN8OvPCCrNfV\nSQSOW+Bzc+0bQ24uI3T6mbgT9O3bpaSgE9LDKCUx8unp3h22JneOsezr64EtWyQzponQGTHCirxX\nyfw5vUrcCXpdnZQUdEL6GHfuHK8OK63FJ2rEva5OhnG/8IJdHzIkstjn5dGHf5bEnaAbKOiEDDCU\nkh/m2LHAhRd2bddaJjMw7hvjy6+osNu0jm7hZ2ZS8KNAQSeE9A1KybyzGRneLh1AwtiMNW9Ef906\nu6211c6M5SX6Y8cO6ukM41bQ2fdCiA9JTZXlggu8202Ujlv0N22y4n/0qIymjWThZ2f7euBVXAm6\nOxc6IWQQEi1KB5DZbhoaQv34f/2rFf8jR2Q0bW6uCL8p3UtKSt99nx4mrgTdPTMRIYR0ISUlcg58\nQFw2u3eL6Jtl3TpJrdDQIG0jRniLvVnPyhqwfvy4EnQzShSQuQg4MpQQclokJwNTp8riRTAoE58Y\nsa+vl3LtWrv+4Ycyo5WX2OfkSFtyct9+L4e4EvSKClufP59ZFQkhPUwgYCN15s713qelRSx5I/YN\nDcArr9j1PXsk/NJL7E09I6NXrPy4EvTycls3k1tQ0AkhfUpKCvCRj8jiRTAoA7Dcgl9TA7z0kl0/\ncSKy2OfkSG6TM0iTHFeCPnRof58BIYTEIBAAxo+XxWt6Q0Ciddy+/Pp6YM0au75vn0ygYgS+m8SV\noN90E/Dgg2KdJyXJOiGExB0jRwLnniuLF52dwP79Vuwff7xbh1W6D2MBi4qKdFVV1Vkdo7KSec8J\nIYMLpdRGrXVRrP3iykIHmPecEEIiMXjHyBJCiM+goBNCiE+goBNCiE+I6UNXSi0E0AygUGu93KO9\n1KkWaK2/08PnRwghpJtEtdCVUoUAoLUuA9Bs1l3tJQDKtNYrAOQ764QQQvqBWC6XGyHWOQDUAggX\n7HzXtlpnnRBCSD8Qy+WSCqDRtZ7hbnQsc0MhgMd66LwIIYScJj3SKeq4YjZprTd5tJUqpaqUUlWH\nDh3qiY8jhBDiQSxBbwaQ7tRTARyJsF9JpA5RrfUKrXWR1rooMzPzDE+TEEJILGIJ+mOwfvF8AGUA\noJRKNTsopUpN9As7RQkhpP+IKujGheIIdbPLpVLu2r5MKVWjlGrq1TMlhBASlZhx6GEdn2bbHKcs\nA5DWC+dFCCHkNOFIUUII8QkUdEII8QkUdEII8QkUdEII8QkUdEII8QkUdEII8QkUdEII8QkUdEII\n8QkUdEII8QkUdEII8QkUdEII8QkUdEII8QkUdEII8QkUdEII8QkUdEII8QkUdEII8QkUdEII8QkU\ndEII8QkUdEII8QkUdEII8QkUdEII8QkUdEII8QkUdEII8QkUdEII8QkUdEII8QkUdEII8QkUdEII\n8QkUdEII8QmJsXZQSi0E0AygUGu9/HTbCSGE9A1RLXSlVCEAaK3LADSb9e62E0II6TtiuVxuhFjf\nAFALoOQ02wkhhPQRsVwuqQAaXesZp9kOpVQpgFJn9UOl1PbTPckzZAyAw330WQMdXgsLr4WF18Iy\n0K9Fbnd2iulDP1u01isArOjtzwlHKVWltS7q688diPBaWHgtLLwWFr9ci1gul2YA6U49FcCR02wn\nhBDSR8QS9McA5Dv1fABlAKCUSo3WTgghpO+JKuha600AoJQqAdBs1gGUx2gfCPS5m2cAw2th4bWw\n8FpYfHEtlNa6v8+BEEJID8CRooQQ4hMo6ISQQYdS6o4I2xcqpUoitQ904l7QlVKlzrIsQntc/4NO\nh25ci2Vmv749s77H+Z+X8L7o1rUYNPcFcKrPb4HH9rgf+R7Xgu78Y8qcWPd8Z93dHvf/oO4S61o4\nlCqlaiCjen2L890XOf/3wsGcsiLWtXAYFPdFN4j7ke9xLeiQUElz0WthQygNcf8POg1iXQsA+LLW\nusD5cfsWrXWZ1vpWZzXfI/pq0NwX3bgWwCC5LwB5mEf5njFHvg90en2kaG/iWKOGQkhcvJu4/wd1\nl25cC8Ba7oMiM6bjTrnVo2nQ3BeGKNcCGFz3RXrsXeKXeLfQAZx6hd40wOLg+4Vo10JrvdyxTjIi\nuGR8hSNOt7oGwg1aol2LwXJfxLDOAR+MfPeFoAMo0Vp/x2N73P+DzgDPa+F0li50Vo/A2yXjC5RS\nbl9xLWxyOMOguS9iXYvBdF9A3kQWOp2/6ea6+Gnke9wLulKq1LwmGuvCT/+g0yHGtaiC/f4Fzrpf\nKUGoYNcCg/a+iHUtBs19obVepbVe5ay631TiYeR7t4jrkaLOhX8C4g9Nh9Obr5TaqLWe4+xTCqeT\nMMzP7CtO41o0Qq6Fb32ljlh91lmdYzoFB+l90d1r4fv7YjAQ14JOCCHEEvcuF0IIIQIFnRBCfAIF\nnRBCfAIFnRBCfAIFnRBCfAIFnRBCfAIFnRBCfML/B92cmKBbrIZwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pylab as plt\n",
    "\n",
    "from notes_utilities import pnorm_ball_line\n",
    "\n",
    "\n",
    "import time\n",
    "from IPython import display\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "#y = np.array([7.46, 6.77, 12.74, 7.11, 7.81, 8.84, 6.08, 5.39, 8.15, 6.42, 5.73])\n",
    "y = np.array([8.04, 6.95, 7.58, 8.81, 8.33, 9.96, 7.24, 4.26, 10.84, 4.82, 5.68])\n",
    "#y = np.array([9.14, 8.14, 8.74, 8.77, 9.26, 8.10, 6.13, 3.10, 9.13, 7.26, 4.74])\n",
    "x = np.array([10., 8., 13., 9., 11., 14., 6., 4., 12., 7., 5.])\n",
    "N = len(x)\n",
    "\n",
    "# Design matrix\n",
    "A = np.vstack((np.ones(N), x)).T\n",
    "\n",
    "w_best, E, rank, s = np.linalg.lstsq(A, y)\n",
    "err = y-A.dot(w_best)\n",
    "E_min = np.sum(err**2)/N\n",
    "\n",
    "\n",
    "def inspect_momentum(alpha = 0.005, beta = 0.97):\n",
    "    ln = pnorm_ball_line(mu=w_best, A=np.linalg.cholesky(np.linalg.inv(A.T.dot(A))),linewidth=1)\n",
    "    ln2 = pnorm_ball_line(mu=w_best, A=4*np.linalg.cholesky(np.linalg.inv(A.T.dot(A))),linewidth=1)\n",
    "\n",
    "    # initial parameters\n",
    "    w0 = np.array([2., 1.])\n",
    "    w = w0.copy()\n",
    "    p = np.zeros(2)\n",
    "\n",
    "    EPOCHS = 100\n",
    "    W = np.zeros((2,EPOCHS))\n",
    "    for epoch in range(EPOCHS):\n",
    "        # Error\n",
    "        err = y-A.dot(w)\n",
    "        W[:,epoch] = w    \n",
    "        # Mean square error\n",
    "        E = np.sum(err**2)/N\n",
    "\n",
    "        # Gradient\n",
    "        dE = -2.*A.T.dot(err)/N\n",
    "        p = dE + beta*p\n",
    "\n",
    "#        if epoch%10 == 1: \n",
    "#            print(epoch,':',E)\n",
    "            # print(w)    \n",
    "\n",
    "        # Perfom one descent step\n",
    "        w = w - alpha*p\n",
    "\n",
    "\n",
    "#    print(E_min)\n",
    "\n",
    "    plt.plot(W[0,:],W[1,:],'.-b')\n",
    "    plt.plot(w_best[0],w_best[1],'ro')\n",
    "    plt.plot(w0[0],w0[1],'ko')\n",
    "    plt.xlim((1.8,4.3))\n",
    "    plt.ylim((0,1.2))\n",
    "    plt.title('$\\\\alpha = $'+str(alpha)+' $\\\\beta = $'+str(beta))\n",
    "    plt.gca().add_line(ln)\n",
    "    plt.gca().add_line(ln2)\n",
    "    plt.show()\n",
    "    \n",
    "inspect_momentum(alpha=0.0014088, beta=0.95)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7914f4b92752492da520347e406fac3a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>interactive</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in Jupyter Notebook or JupyterLab, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another notebook frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.005, description='alpha', max=0.02, step=0.001), FloatSlider(value=0.97, description='beta', max=0.99, step=0.001), Output()), _dom_classes=('widget-interact',))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.inspect_momentum>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "from __future__ import print_function\n",
    "from ipywidgets import interact, interactive, fixed\n",
    "import ipywidgets as widgets\n",
    "import matplotlib.pylab as plt\n",
    "from IPython.display import clear_output, display, HTML\n",
    "\n",
    "interact(inspect_momentum, alpha=(0, 0.02, 0.001), beta=(0, 0.99, 0.001))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Advanced Material\n",
    "\n",
    "Why does momentum work? by Goh. Interactive articles\n",
    "https://distill.pub/2017/momentum/\n",
    "\n",
    "Beyond vanilla Gradient Descent, Sebastian Ruder blog\n",
    "http://ruder.io/optimizing-gradient-descent/\n",
    "\n",
    "Blog post by M. Hardt\n",
    "http://blog.mrtz.org/2013/09/07/the-zen-of-gradient-descent.html\n",
    "\n",
    "A great talk by Ben Recht:\n",
    "https://simons.berkeley.edu/talks/ben-recht-2013-09-04\n",
    "\n",
    "Backpropagation, by blog post by S. Arora\n",
    "http://www.offconvex.org/2016/12/20/backprop/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis of convergence of Gradient descent for a quadratic function\n",
    "\n",
    "Recall that the error function we minimize is\n",
    "$$\n",
    "E(w) = \\frac{1}{2} (y-Aw)^T(y-Aw) = \\frac{1}{2}(y^\\top y - 2 y^\\top A w + w^\\top A^\\top A w)\n",
    "$$\n",
    "\n",
    "The gradient at the point $w$ will be denoted as $\\nabla E(w) = g(w)$ where\n",
    "$$g(w) = -A^\\top (y - Aw) = A^\\top A w - A^\\top y$$\n",
    "\n",
    "Moreover, the gradient at the minimum will vanish:\n",
    "$$\n",
    "g(w_\\star) = 0\n",
    "$$\n",
    "Indeed, we can solve \n",
    "$$0 = A^\\top (Aw_\\star - y)$$\n",
    "as\n",
    "$$w_\\star =  (A^\\top A)^{-1}A^\\top y  $$\n",
    "but this is not our point.\n",
    "\n",
    "For a constant learning rate $\\eta$, gradient descent executes the following iteration\n",
    "$$\n",
    "w_t = w_{t-1} - \\eta g(w_{t-1}) = w_{t-1} - \\eta A^\\top (Aw_{t-1} - y)\n",
    "$$\n",
    "\n",
    "$$\n",
    "w_t  =  (I - \\eta A^\\top A) w_{t-1} + \\eta A^\\top y\n",
    "$$\n",
    "\n",
    "This is a fixed point equation of form\n",
    "$$\n",
    "w_t  =  T(w_{t-1}) \n",
    "$$\n",
    "where $T$ is an affine transformation. \n",
    "\n",
    "We will assume that $T$ is a contraction, i.e. for any two different\n",
    "parameters $w$ and $w'$ in the domain we have\n",
    "$$\n",
    "\\| T(w) - T(w') \\| \\leq L_\\eta \\|w-w' \\|\n",
    "$$\n",
    "\n",
    "where $L_\\eta < 1$, then the distance shrinks. Hence the mapping converges to a fixed point (this is a consequence of a deeper result in analysis called the Brouwer fixed-point theorem (https://en.0wikipedia.org/wiki/Brouwer_fixed-point_theorem))\n",
    "\n",
    "We will consider in particular the distance between the optimum and the current point $w(t)$\n",
    "\n",
    "$$\n",
    "\\| T(w_t) - T(w_\\star) \\| \\leq L_\\eta \\|w_t - w_\\star \\|\n",
    "$$\n",
    "But we have \n",
    "$T(w_\\star) = w_\\star$ and $w_t = T(w_{t-1})$ so\n",
    "$\\|w_t - w_\\star \\| = \\|T(w_{t-1}) - T(w_\\star) \\|$. \n",
    "\n",
    "\\begin{align}\n",
    "\\| T(w_t) - T(w_\\star) \\| & \\leq L_\\eta \\|T(w_{t-1}) - T(w_\\star) \\| \\\\\n",
    "& \\leq L^2_\\eta \\|T(w_{t-2}) - T(w_\\star) \\| \\\\\n",
    "\\vdots \\\\\n",
    "& \\leq L^{t+1}_\\eta \\| w_{0} - w_\\star \\| \n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "T(w)   =  (I - \\eta A^\\top A) w + \\eta A^\\top y\n",
    "$$\n",
    "\n",
    "$$\n",
    "T(w_\\star)   =  (I - \\eta A^\\top A) w_\\star + \\eta A^\\top y\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\| T(w) - T(w') \\| = \\| (I - \\eta A^\\top A) (w-w') \\| \\leq \\| I - \\eta A^\\top A \\| \\| w-w' \\|\n",
    "$$\n",
    "\n",
    "When the norm of the matrix $\\| I - \\eta A^\\top A \\| < 1$ we have convergence. Here we take the operator norm, i.e., the magnitude of the largest eigenvalue.\n",
    "\n",
    "Below, we plot the absolute value of the maximum eigenvalues of $I - \\eta A^\\top A$ as a function of $\\eta$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEKCAYAAAA4t9PUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XucXWV97/HPdyb3ZCYh95AbAQIh3jgwXPUUKNoCVoPW\nY0GtFCgpylV7o1gP2p5aqva0UFFMSwTUglhQgwfk9pJErQES5BJIAjESSMgkIVwySch1fuePtWZm\nzziXNZm9ZvZif9+v137NXs9aa+aXwM5vnvX8nudRRGBmZnagagY6ADMzKzYnEjMz6xMnEjMz6xMn\nEjMz6xMnEjMz6xMnEjMz65PcEomkhZI2S1rRxXlJul7SGklPSTqmp3sljZX0gKTn068H5RW/mZll\nk2eP5GbgjG7OnwnMTl/zgW9kuPcq4KGImA08lB6bmdkAyi2RRMQS4NVuLpkH3BqJpcAYSVN6uHce\ncEv6/hbg7DKGbGZmB2DQAP7sqcBLJcfr07aN3dwzKSJazjcCk7q6UNJ8kp4OI0eOPHbOnDl9i9bM\n+sWLr+5k1979HDGpbqBDqXrLly9/JSIm9HTdQCaSPomIkNTl+i4RsQBYANDQ0BDLli3rt9jM7MB9\n6jvL+fWW7dz/mVMGOpSqJ2ldlusGsmprAzC95Hha2tadTS2Pv9Kvm3OKzcwGSHMENdJAh2G9MJCJ\nZBHwybR660TgjZLHVt3dc176/jzgR3kGaGb9r9nryBZObo+2JN0GnAqMl7QeuAYYDBARNwL3AGcB\na4CdwPnd3RsRNwHXAndIuhBYB3w0r/jNbGBE4B5JweSWSCLi3B7OB3BJb+6NiK3A6X2PzswqVURQ\n46nSheL/XGZWUZojEO6RFIkTiZlVlABqnEcKxYnEzCpKc4A8RlIoTiRmVlEiAueRYsmUSCTNlPTe\n9P1wSZ5yama5cNVW8fSYSCRdBPwX8M20aRrwwzyDMrPqlUxIHOgorDey9EguAd4NbAOIiOeBiXkG\nZWbVqznCYyQFkyWR7I6IPS0HkgaRFFaYmZVdBC7+LZgsiWSxpKuB4ZLeB3wfuDvfsMysWnmMpHiy\nJJKrgC3A08CfkSxt8rd5BmVm1avZM9sLp8clUiKiGfj39GVmlivPbC+eHhOJpN/QyZhIRByaS0Rm\nVtUCPI+kYLIs2thQ8n4Y8L+AsfmEY2bVrtljJIXT45PIiNha8toQEf8KvL8fYjOzKuSZ7cWT5dHW\nMSWHNSQ9lMJu0Wtmlc1VW8WTJSH8c8n7fcALeEMpM8uJZ7YXT5aqrdP6IxAzM2jZateZpEi6TCSS\nPtvdjRHxf8sfjplVu3CPpHC665F4hV8z63ceIymeLhNJRHyxPwMxMwPPbC+iLFVbw4ALgbeRzCMB\nICIuyDEuM6tSAZ7ZXjBZ8v63gcnA7wOLSfYjacozKDOrXs2eR1I4WRLJ4RHxeWBHRNxCMhnxhHzD\nMrNq5TGS4smSSPamX1+X9HZgNN7Yysxy4pntxZNlQuICSQcBnwcWAaPS92ZmZee1toonS4/kWxHx\nWkQsjohDI2JiRHyzp5skLZS0WdKKLs5L0vWS1kh6qnQpFklnSFqdnruqpP1oSUslPSFpmaTjM/0p\nzawwPEZSPFkSyW8kLZB0unq3kfLNwBndnD8TmJ2+5gPfAJBUC9yQnp8LnCtpbnrPl4EvRsTRwP9O\nj83sLSTZateZpEiyJJI5wIPAJcALkr4m6T093RQRS4BXu7lkHnBrJJYCYyRNAY4H1kTE2nSv+NvT\nayGpDKxP348GXs4Qv5kViGe2F0+WtbZ2AncAd6RjJdeRlAHX9vFnTwVeKjlen7Z11t5SJXYlcJ+k\nr5IkwZO7+uaS5pP0dJgxY0YfQzWz/uIxkuLJNH9U0imSvg4sJ5mUOFCr/34K+ExETAc+A9zU1YUR\nsSAiGiKiYcKECf0WoJn1jcdIiifLzPYXgF+R9Er+MiJ2lOlnbwCmlxxPS9sGd9EOcB5wRfr++8B/\nlCkWM6sQyVa7ziRFkqX8950RsS2Hn70IuFTS7SSPrt6IiI2StgCzJc0iSSDnAB9L73kZOAV4GPhd\n4Pkc4jKzAeQxkuLJkkiGSroaOKT0+p7W2pJ0G3AqMF7SeuAakt4GEXEjcA9wFrAG2Amcn57bJ+lS\n4D6ScZiFEfFM+m0vAq6TNAjYRToGYmZvHR4jKZ4sieRHwM9IKrf2Z/3GEXFuD+eDpBKss3P3kCSa\nju0/B47NGoOZFY9nthdPlkQyIiL+OvdIzMxwj6SIslRt/VjSWblHYmaGq7aKKEsiuYIkmeyStE1S\nk6Q8Bt/NzMAz2wsny4REb7lrZv2m2VVbhdNjjyRdXPETkj6fHk/3YolmlpfmgBpnkkLJ8mjr68BJ\ntM3l2E6yqKKZWdkF4QdbBZOlauuEiDhG0q8AIuI1SUNyjsvMqlRzeGZ70WTaITFd2j0AJE0AmnON\nysyqlme2F0+WRHI98ANgoqR/AH4OfCnXqMysaiU9koGOwnojS9XWdyUtB04HBJwdEStzj8zMqlLS\nI3EmKZIsq/+OBTYDt5W0DY6IvXkGZmbVyWMkxZPl0dbjwBbgOZLVdreQ7JT4uCSve2VmZZMswYfH\nSAomSyJ5ADgrIsZHxDiSvdR/DHyapDTYzKws0jzime0FkyWRnBgR97UcRMT9wEnpPutDc4vMzKpO\ns3skhZRlHslGSX8N3J4e/xGwKS0JdhmwmZVNc9oj8cz2YsnSI/kYyXa3P0xfM9K2WgZu73YzewsK\nYqBDsAOQpfz3FeCyLk6vKW84ZlbNWsZIXP5bLF0mEkn/GhFXSrobfvvXhIj4YK6RmVnV8RhJMXXX\nI/l2+vWr/RGImVlr1ZYTSaF0mUgiYnn6dbGk4cCMiFjdb5GZWdVp65E4kxRJlv1IPgA8AfwkPT5a\n0qK8AzOz6tPc2iNxIimSLFVbXwCOB14HiIgngFk5xmRm1ap1QqIVSaZl5CPijQ5trtEzs7LzYHsx\nZZmQ+IykjwG1kmYDlwP/nW9YZlaNWhOJM0mhZOmRXAa8DdgN/CfwBnBlTzdJWihps6QVXZyXpOsl\nrZH0lKRjSs6dIWl1eu6qDvddJmmVpGckfTlD/GZWEB4jKaYsExJ3Ap9LX71xM/A14NYuzp8JzE5f\nJwDfAE5Il165AXgfsB54TNKiiHhW0mnAPOBdEbFb0sRexmRmFaxlZrvTSLFk6ZEckIhYArzazSXz\ngFsjsRQYI2kKycD+mohYGxF7SNb4mpfe8yng2ojYnf6MzXnFb2b9zzPbiym3RJLBVOClkuP1aVtX\n7QBHAP9T0iOSFks6rqtvLmm+pGWSlm3ZsqXMoZtZHjzYXkxZ5pG8O0tbPxkEjAVOBP4SuENdPEyN\niAUR0RARDRMmTOjPGM3sAHlmezFl6ZH8W8a23toATC85npa2ddUOSe/krvRx2KMky9iPL0MsZlYB\nWnokHmwvlu4WbTwJOBmYIOmzJafqSZaQ76tFwKWSbicZbH8jIjZK2gLMljSLJIGcQ7JsPSTL2J8G\n/FTSEcAQ4JUyxGJmFcBjJMXUXdXWEGBUek1dSfs24CM9fWNJtwGnAuMlrQeuAQYDRMSNwD3AWSRL\n0e8Ezk/P7ZN0KXAfScJaGBHPpN92IbAwLSneA5wXLZs8m1nhhWe2F1J3izYuBhZLujki1gFIqgFG\nRcS2nr5xRJzbw/kALuni3D0kiaZj+x7gEz39bDMrprYJiQMciPVKlv9c/yipXtJIYAXwrKS/zDku\nM6tCXv23mLIkkrlpD+Rs4F6SBRv/ONeozKwq+Tl1MWVJJIMlDSZJJIsiYi/+721mOQj3SAopSyL5\nJvACMBJYImkmyYC7mVlZNbtqq5CyrLV1PXB9SdO6dM0rM7Oy8sz2Ysoys32SpJsk3ZsezwXOyz0y\nM6s6ntleTFkebd1MMqfj4PT4OTIsI29m1lue2V5MWRLJ+Ii4g2Q5EiJiH7A/16jMrCp5ZnsxZUkk\nOySNI63UknQiyeZWZmZl5ZntxZRlq93PkqyLdZikXwATyLBEiplZb3lmezFlqdp6XNIpwJEkvyis\nTueSmJmVlcdIiqnHRCLpwx2ajpD0BvC0dyg0s3JqmensNFIsWR5tXQicBPw0PT4VWA7MkvR3EfHt\nnGIzsyrjme3FlCWRDAKOiohNkMwrAW4l2UNkCeBEYmZl4ZntxZRlSGt6SxJJbU7bXgU8VmJmZeMJ\nicWUpUfysKQfA99Pj/8wbRsJvJ5bZGZWddoG2wc4EOuVLInkEuDDwHvS41uBO9ONqbzmlpmVjfcj\nKaZuE4mkWuDBiDgNuLN/QjKzquUxkkLqdowkIvYDzZJG91M8ZlbFmj1GUkhZHm1tB56W9ACwo6Ux\nIi7PLSozq0peRr6YsiSSu9KXmVmuPLO9mLIskXJLfwRiZuaZ7cWUZYmU2cA/AnOBYS3tEXFojnGZ\nWRXyzPZiyjIh8VvAN4B9JOW+twLfyTMoM6tOzc3JVyeSYsmSSIZHxEOAImJdRHwBeH++YZlZNWp9\ntOU8UihZEsluSTXA85IulfQhYFRPN0laKGmzpBVdnJek6yWtkfSUpGNKzp0haXV67qpO7v1zSSFp\nfIb4zawgPLO9mLIkkiuAEcDlwLHAHwPnZbjvZuCMbs6fCcxOX/NJHp+1TIK8IT0/FzhX0tyWmyRN\nB34PeDFDDGZWIB4jKaYsVVuPpW+3A+dn/cYRsUTSId1cMg+4NV1qZamkMZKmAIcAayJiLYCk29Nr\nn03v+xfgr4AfZY3FzIrBizYWU5eJRNLdtD2y/C0R8cE+/uypwEslx+vTts7aT0hjmgdsiIgne6oz\nlzSfpKfDjBkz+hiqmfUHLyNfTN31SL7ab1FkIGkEcDXJY60eRcQCYAFAQ0NDlwnRzCqHZ7YXU5eJ\nJCIW5/yzNwDTS46npW2Du2g/DJgFtPRGpgGPSzo+IhpzjtXM+kFb1ZYzSZFkGWzPyyLgk2n11onA\nGxGxEXgMmC1plqQhwDnAooh4OiImRsQhEXEIySOvY5xEzN46WgbbnUaKJctaWwdE0m0k+7uPl7Qe\nuIakt0FE3AjcA5wFrAF2kg7kR8Q+SZcC9wG1wMKIeCavOM2scng/kmLKnEgkjYiInVmvj4hzezgf\nJJtmdXbuHpJE0939h2SNxcyKITzYXkg9PtqSdLKkZ4FV6fG7JH0998jMrOp4P5JiyjJG8i/A7wNb\nASLiSeB38gzKzKqTZ7YXU6bB9oh4qUPT/hxiMbMq55ntxZRljOQlSScDIWkwyZIpK/MNy8yqkWe2\nF1OWHsnFJIPiU0nmcxxNF4PkZmZ94ZntxZRlra1XgI/3QyxmVuU8RlJMWXZInABcRLKYYuv1EXFB\nfmGZWTVq22rXmaRIsoyR/Aj4GfAgHmQ3sxyF19oqpCyJZERE/HXukZhZ1WtudtVWEWUZbP+xpLNy\nj8TMql7Loy0nkmLpbj+SJpL/rgKulrQb2JseR0TU90+IZlYtmtsGSaxAultGvq4/AzEz8xhJMWVZ\na+uhLG1mZn3lRRuLqbtHW8OAkSTLwB9EW2eznmRyoplZWXkeSTF1V7X1Z8CVwMHActoSyTbgaznH\nZWZVyDPbi6m7MZLrgOskXRYR/9aPMZlZlQrcIymiHsdInETMrL+0Ltrosq1CGcg9283M2mmbkDjA\ngVivdJlIJL07/Tq0/8Ixs2rmMZJi6q5Hcn369Zf9EYiZmcdIiqm7qq29khYAUyVd3/FkRFyeX1hm\nVo3a9mx3JimS7hLJHwDvJdmvfXn/hGNm1SwiPD5SQN2V/74C3C5pZUQ82Y8xmVmVivD4SBFlqdra\nKukHkjanrzslTcs9MjOrOs0RHh8poCyJ5FvAIpIZ7gcDd6dt3ZK0ME08K7o4L0nXS1oj6SlJx5Sc\nO0PS6vTcVSXtX5G0Kr3+B5LGZIjfzAqiOTw+UkRZEsnEiPhWROxLXzcDEzLcdzNwRjfnzwRmp6/5\nwDcAJNUCN6Tn5wLnSpqb3vMA8PaIeCfwHPA3GeIws4IIPEZSRFkSySuSPiGpNn19Atja000RsQR4\ntZtL5gG3RmIpMEbSFOB4YE1ErI2IPcDt6bVExP0RsS+9fyngR2xmbyERntVeRFkSyQXAR4FGYCPw\nEeD8MvzsqcBLJcfr07au2juL696uvrmk+ZKWSVq2ZcuWMoRrZnlrbnaPpIh63LM9ItYBH+yHWDKT\n9DlgH/Ddrq6JiAXAAoCGhobo6jozqxyBq7aKqMdEkqMNwPSS42lp2+Au2gGQ9Cckc1xOj5bt1Mzs\nLaE5ebZlBTOQizYuAj6ZVm+dCLwRERuBx4DZkmZJGgKck16LpDOAvwI+GBE7BypwM8uH55EUU249\nEkm3AaeS7LC4HriGpLdBRNwI3AOcBawBdpKOu0TEPkmXAvcBtcDCiHgm/bZfA4YCD6Qlgksj4uK8\n/gxm1r+aPbO9kHpMJJImAV8CDo6IM9NS3JMi4qbu7ouIc3s4H8AlXZy7hyTRdGw/vKd4zay4wvNI\nCinLo62bSXoHB6fHz5FswWtmVlbukRRTlkQyPiLuAJohefQE7M81KjOrSp7ZXkxZEskOSeNIKvNo\nGRjPNSozq1LukRRRlsH2z5JUTR0m6Rcky6N8JNeozKwqNTd7ZnsRZZmQ+LikU4AjSSq8V0fE3twj\nM7Oq4zGSYspStfXJDk3HSCIibs0pJjOrUoHHSIooy6Ot40reDwNOBx4HnEjMrKy8H0kxZXm0dVnp\ncboHyO25RWRmVcsz24vpQJZI2QHMKncgZmbes72YsoyR3E1a+kuSeOYCd+QZlJlVJ88jKaYsYyRf\nLXm/D1gXEetzisfMqpjHSIopyxjJ4v4IxMzM+5EUU5eJRFITbY+02p0iWXOxPreozKwqeYykmLpM\nJBFR15+BmJl5ZnsxZd6PRNJEknkkAETEi7lEZGZVy2MkxdRj+a+kD0p6HvgNsBh4Abg357jMrAp5\njKSYsswj+XvgROC5iJhFMrN9aa5RmVlVCvdICilLItkbEVuBGkk1EfFToCHnuMysCjV7ZnshZRkj\neV3SKGAJ8F1Jm0lmt5uZlZWrtoopS49kHrAT+AzwE+DXwAfyDMrMqlNzsvzvQIdhvZSlR/JnwPci\nYgNwS87xmFkV834kxZSlR1IH3C/pZ5IulTQp76DMrHp5jKR4ekwkEfHFiHgbcAkwBVgs6cHcIzOz\nqtMc4emIBdSbZeQ3A43AVmBiPuGYWTVrbnaPpIiyTEj8tKSHgYeAccBFEfHODPctlLRZ0oouzkvS\n9ZLWSHpK0jEl586QtDo9d1VJ+1hJD0h6Pv16UJY/pJkVQ+B5JEWUpUcyHbgyIt4WEV+IiGczfu+b\ngTO6OX8mMDt9zQe+ASCpFrghPT8XOFfS3PSeq4CHImI2SWK7quM3NbPi8jySYsoyRvI3EfFEb79x\nRCwBXu3mknnArZFYCoyRNAU4HlgTEWsjYg/Jtr7zSu5pqRy7BTi7t3H1xqZtu9i6fXeeP8LMSnhm\nezFlXrQxB1OBl0qO16dtnbWfkL6fFBEb0/eNQK4VZF//6Rpu+eU6JtQNZc7kuvRVz5wpdRw+cRRD\nB9Xm+ePNqo73bC+mgUwkfRIRIamz/VIAkDSf5JEZM2bMOKCf8ZFjpzN97AhWNTaxqnEbt/xyHXv2\nNQNQWyMOmzCyNbEclX6dXD/MW4WaHSCv/ltMmRKJpJnA7Ih4UNJwYFBENPXxZ28gGX9pMS1tG9xF\nO8AmSVMiYmP6GGxzV988IhYACwAaGhq6TDjdece00bxj2ujW4337m3lh6w5WbmxidZpclq97jUVP\nvtx6zejhg5kzuY6jptQzZ3IdR6avEUMKm7PN+o33bC+mHv91k3QRyW/2Y4HDSP5hv5FkFeC+WARc\nKul2kkdXb6QJYgswW9IskgRyDvCxknvOA65Nv/6ojzH0yqDaGg6fWMfhE+v4wLva2t94cy/PbWpi\n1cZtrGxMvn5/2Uvs2LMfSFZ8mDl2RGvvZc7keo6aUsf0g0ZQ42m8Zq281lYxZfk1+RKSAfBHACLi\n+XSTq25Jug04FRgvaT1wDUlvg4i4EbgHOAtYQ7KW1/npuX2SLgXuA2qBhRHxTPptrwXukHQhsA74\naLY/Zr5GDx/McYeM5bhDxra2NTcH6197k5WN21i1Mem9rGps4r5nG4m0fzRiSC1HTm5LLHMm13Pk\n5DpGDx88QH8Ss4Hl/UiKKUsi2R0Re1q6m5IG0fle7u1ExLk9nA+SJNXZuXtIEk3H9q30vSfUL2pq\nxIxxI5gxbgS//7bJre079+zj+U3bWdW4jZVpgrnn6Y3c9mjbhpNTxwxPBvan1HHk5HqOmlzHrPEj\nGVTbm/mjZsXjme3FlCWRLJZ0NTBc0vuATwN35xvWW9eIIYN41/QxvGv6mNa2iGDTtt3tey8bm1j8\n3Bb2NSc5e8igGmZPHNWu9zJnSh3jRw0dqD+KWdk1N3uMpIiyJJKrgAuBp0lWAr4nIv4916iqjCQm\njx7G5NHDOO3ItqeGu/ft59ebd7B6U5JYVjY28bPnt3Dn4+tbrxk/amiaWFyabMWXPNoa6Cist7Ik\nkssi4jqgNXlIuiJtsxwNHVTL3IPrmXtwPfyPtvat23ezurGpdWB/VWMTt/5yHbtLSpMPHT+SOWnl\nWEsPZspolyZbZfOExGLKkkjOAzomjT/ppM36ybhRQzn58KGcfPj41rakNHln62OxVY3beHzda9xd\nUppcP2wQc6YkYy5zSsqTXZpslSLZj8SZpGi6/BdE0rkkZbezJC0qOVVH90uf2ABISpNHcfjEUfxB\nyZKa23bt5bkOvZf/Wr6+29LkOZPrmDHWpcnW/zyzvZi6+1X0v4GNwHjgn0vam4Cn8gzKyqd+2GAa\nDhlLQ4fS5A2vv8nKNLG09GI6liYfMantsVjLGMzoES5Ntvx4ZnsxdZlIImIdyVyNk/ovHOsPNTVi\n+tgRTB87gt8rKU1+c8/+ZGJlOudl1cYm7l3RyG2Pti19dvDoYa2PxVoek7k02colPLO9kLLMbG+i\nbd7IEJJJhTsioj7PwKz/DR9S22lp8uam3W29l/TrktLS5NoaZk8a1dZzSXsxE+pcmmy946qtYuox\nkUREXct7Jb8qzANOzDMoqxySmFQ/jEn1wzi1pDR5z75mfr1le+tjsc5Lk4eUJJfk6+ETRzFssEuT\nrXMebC+mXpXrpLPRfyjpGrypVFUbMqiGo6bUc9SU9qXJr+7Y065ybFVjE99e6tJky8Yz24spy6Ot\nD5cc1gANwK7cIrJCGztyCCcfNp6TD2srTd7fHLywdUdrclm5sYlfvdhJaXJp5diUOo6cVMfIoS5N\nriYeIymmLJ/SD5S83we8QNuOhWY9SvZuGcVhE0bx/ndOaW1vKU0urRy76/ENbN+9rvWameNGtFaM\ntfReXJr81pWU/w50FNZbWcZIzu+PQKz6dFaaHJGsmlw6sL+ycRsPPLuJdGyf4YOTVZOPSnstLY/J\nxowYMkB/EisXl/8WU3cTEv+Nblb5jYjLc4nIqprUVpr8vrltOym/uWc/z29uSgf2k97LTzqUJk8Z\nPazdwP5RU+qZNX4kg12aXBgebC+m7noky/otCrMeDB9SyzunjeGd0367NLld72XjNn6+5hX27m8r\nTT584qh22yG7NLlyeYykmLqbkHhL6bGkUWn79ryDMsuitDT5lCMmtLbv2dfM2le2t+u9/GLNK9z1\n+IbWa1pKk5ONxZLei0uTB16y1e5AR2G9laVq6+3At0m22lW6Fe4nS3YtNKsoQwbVpPNX6jmbqa3t\nHUuTVzc28d1H1rFrb1tp8qzxI1sTS8tjsoNdmtxvvNVuMWWp2loAfDYifgog6VSSJeVPzjEus7Lr\nqjR53dYdrY/HVjY28eT61/nxUxtbr6kbNqjdY7Ej01WTR7k0uey81W4xZfkkjGxJIgAR8bCkkTnG\nZNZvamvEoRNGceiEUZz1jrbS5KZde3luU1PrdsidlSbPGDui3Zpjc6Ykpcm1/pX6gHmwvZiyJJK1\nkj5P8ngL4BPA2vxCMht4dcMGc+zMsRw7s31p8obX32ybWJn2Yh5c2b40+YjJdUlimezS5N5qbu6y\nUNQqWJZEcgHwReCu9HhJ2mZWVSQx7aARTDtoBO8tKU3etXc/z2/a3jqwv6pxG/c908jtj3Vdmjxn\ncj2HTnBpckd+tFVMWSYkvgZcDiCpluRR17a8AzMrimGDa3nHtNG8Y9ro1raIYEvT7tZeS8vWyD9f\ns7ZdafJhE0elj8XaloaZMGpo1Q7ue2Z7MWWp2vpP4GJgP/AYUC/puoj4St7BmRWVJCbWD2Nih9Lk\nvfubWbtlR+uaY6sat/Hfv97KXb9qK00eN3JIu90qq6k02TPbiynLo625EbFN0seBe0lW/V0OOJGY\n9dLg2prWqq95R7e1v7ZjT7s1x1Y1bmtXmlwjktLkloH9tPcydczwt1TvxVvtFlOWRDJY0mDgbOBr\nEbFXkkfEzMrooJFDOOmwcZx02LjWtpbS5JbHYqs2buPp9W/w/zqUJs8pSSwt5clFLU1OeiROJEWT\n5f+2b5Ks+PsksETSTCDTGImkM4DrgFrgPyLi2g7nDwIWAoeRLE1/QUSsSM9dAVwECPj3iPjXtP1o\n4EZgGMlqxJ+OiEezxGNWJKWlyWeWlCZv372P1R16Lz/81Qaalu5rvaaopcnhme2FlGWw/Xrg+pKm\ndZJO6+m+dGD+BuB9wHrgMUmLIuLZksuuBp6IiA9JmpNef3o6m/4i4HhgD/ATST+OiDXAl4EvRsS9\nks5Kj0/N8Gc1e0sYNXQQx848iGNnHtTa1llp8urGpnalycMG1ySrJZfu+zK5joNGVk5pcuCZ7UWU\nZbB9HHAN8B6S6ryfA38HbO3h1uOBNRGxNv0+t5PsY1KaSOYC1wJExCpJh0iaBBwFPBIRO9N7FwMf\nJkkaAbTsFz8aeBmzKtddafKazdtZmS5quapxGw+s3MT3lrWVJk+uH9aaWFr2fBmo0uRmj5EUUpZH\nW7eTzB35w/T448D3gPf2cN9U4KWS4/XACR2ueZIkQfxM0vHATGAasAL4hzSJvQmcRdtqxFcC90n6\nKsmOjZ3q2hxHAAAIQUlEQVQu1SJpPjAfYMaMGT2EavbWNGxwLW+fOpq3T+1Qmrx9d9t2yBuTMZhf\nlJQmD64Vh0+say1NPnJy8ohsQl2+pcneareYsiSSKRHx9yXH/0fSH5Xp518LXCfpCeBp4FfA/ohY\nKemfgPuBHcATJOXHAJ8CPhMRd0r6KHATnSS1iFhAsk4YDQ0NLg4wS0liYt0wJtYN43c6lCb/5pUd\nbb2Xjdv45dr2pcljRw5pN7h/1OR6Zk8qX2myl5EvpiyJ5H5J5wB3pMcfAe7LcN8GYHrJ8bS0rVU6\nsfF8SJYVBn5DuvxKRNxEkiSQ9CWSHg3AecAV6fvvA/+RIRYz68Hg2hqOmFTHEZPq2u2l/frOPR12\nrGzitkdf5M29ye925SpN3pV+vyGDPNu/aLrbIbGJZDxCJI+TvpOeqgG2A3/Rw/d+DJgtaRZJAjkH\n+FiHnzEG2BkRe4A/BZa0zJqXNDEiNkuaQfL468T0tpeBU4CHgd8Fns/0JzWzAzJmxBBOPHQcJx7a\nvjT5xVd3tq6YvLpxGys2dChNHjoo2e+lZPzliEl11A0b3OnPeX3nXgAO8rpkhdPdxlZ1ffnGEbFP\n0qUkvZdaYGFEPCPp4vT8jSSD6rek81KeAS4s+RZ3pmMke4FLIuL1tP0iksdhg0hKhuf3JU4z672W\nvVtmjR/5W6XJz21qajf+8qMnXuY7u15svWb62OFJYilZe2zmuJG8umMPAGNHdp5orHIpoufhg3S+\nx2ySuRsARMSSHOMqq4aGhli2zDsHmw2EiODlN3a12w55VWMTa7dsb1eaPLl+GC9s3cn35p/ICSW9\nHxs4kpZHRENP12Up//1TkjGJaSSD3icCvyR5rGRm1i1JTB0znKljhnP6Ub9dmlw6/jJy6CCOmNSn\nhyE2ALIMtl8BHAcsjYjT0omDX8o3LDN7q+usNNmKKUt5xK6I2AUgaWhErAKOzDcsMzMriiw9kvVp\nddUPgQckvQas6+EeMzOrElnW2vpQ+vYLkn5KsizJT3KNyszMCqNXa01HxOK8AjEzs2LyFFIzM+sT\nJxIzM+sTJxIzM+sTJxIzM+sTJxIzM+uTTGttFZ2kLRz43JfxwCtlDCcPlR5jpccHlR9jpccHjrEc\nKi2+mRExoaeLqiKR9IWkZVkWLRtIlR5jpccHlR9jpccHjrEcKj2+rvjRlpmZ9YkTiZmZ9YkTSc8W\nDHQAGVR6jJUeH1R+jJUeHzjGcqj0+DrlMRIzM+sT90jMzKxPnEjMzKxPqiqRSDpD0mpJayRd1cl5\nSbo+Pf+UpGN6ulfSWEkPSHo+/XpQBcb4FUmr0ut/kO4vU1Exlpz/c0khaXylxSfpsvTv8RlJXz7Q\n+PKKUdLRkpZKekLSMknHD1B8CyVtlrSiwz2V9FnpKsayfVbyiK/kfJ8/J2UVEVXxAmqBXwOHAkOA\nJ4G5Ha45C7gXEMne9I/0dC/wZeCq9P1VwD9VYIy/BwxK3/9TJcaYnp8O3EcyeXR8JcUHnAY8CAxN\njydW2t8hcD9wZsn9D/d3fOm53wGOAVZ0uKciPis9xFiWz0pe8ZXrc1LuVzX1SI4H1kTE2ojYA9wO\nzOtwzTzg1kgsBcZImtLDvfOAW9L3twBnV1qMEXF/ROxL718KTKu0GFP/AvwV0JcKkLzi+xRwbUTs\nBoiIzRUYYwD16fvRwMsDEB8RsQR4tZPvWymflS5jLONnJa+/QyjP56SsqimRTAVeKjlen7Zluaa7\neydFxMb0fSMwqQJjLHUByW9BFRWjpHnAhoh4sg+x5RYfcATwPyU9ImmxpOMqMMYrga9Iegn4KvA3\nAxBfdyrls5JVXz4rucRXxs9JWfVqh0TrXkSEpIr5LaEjSZ8D9gHfHehYSkkaAVxN8lihUg0CxpI8\ngjgOuEPSoZE+a6gQnwI+ExF3SvoocBPw3gGOqVP+rPReJX9OqqlHsoHk2WKLaWlblmu6u3dTS3c0\n/dqXRx55xYikPwH+APh4H//xyyPGw4BZwJOSXkjbH5c0uULig+S3xbvSxxCPAs0kC+wdiLxiPA+4\nK33/fZLHK/0dX3cq5bPSrTJ9VvKIr5yfk/Ia6EGa/nqR/Ea5luQ/RMvg19s6XPN+2g9+PdrTvcBX\naD+A+OUKjPEM4FlgQqX+PXa4/wUOfLA9r7/Di4G/S98fQfJIQhUW40rg1PT96cDy/o6v5Pwh/PZA\ndkV8VnqIsSyflbziK9fnpNyvAQ+gX/+wSZXEcyTVFJ9L2y4GLk7fC7ghPf800NDdvWn7OOAh4HmS\nqp6xFRjjGpJ/+J5IXzdWWowdvn+fPiA5/R0OAb4DrAAeB3630v4OgfcAy0n+0XoEOHaA4rsN2Ajs\nJenJXViBn5WuYizbZyWP+Mr5OSnny0ukmJlZn1TTGImZmeXAicTMzPrEicTMzPrEicTMzPrEicTM\nzPrEicRsgEi6eqBjMCsHl/+aDRBJ2yNi1EDHYdZXXmvLrB9I+gRwOcnExkeAbcBwSU8Az0TExyX9\nkGTJjGHAdRFRyP27rfo4kZjlTNJRwB8B746IvZK+TjKT+c2IOLrk0gsi4lVJw4HHJN0ZEVsHImaz\n3nAiMcvf6cCxJMkBYDidL1h4uaQPpe+nA7MBJxKreE4kZvkTcEtEtNsfRNJflLw/lWTJ95MiYqek\nh0kecZlVPFdtmeXvIeAjkiZC697lM4G9kgan14wGXkuTyByS1WDNCsE9ErOcRcSzkv4WuF9SDcmK\nrpcAC4CnJD1OshvfxZJWAqtJtnk1KwSX/5qZWZ/40ZaZmfWJE4mZmfWJE4mZmfWJE4mZmfWJE4mZ\nmfWJE4mZmfWJE4mZmfXJ/wdbJoKE5/A8YwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11bc50080>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "left = 0.0000\n",
    "right = 0.015\n",
    "\n",
    "N = 1000\n",
    "ETA = np.linspace(left,right,N)\n",
    "\n",
    "def compute_largest_eig(ETA, A):\n",
    "    \n",
    "    LAM = np.zeros(N)\n",
    "    D = A.shape[1]\n",
    "    n = A.shape[0]\n",
    "    \n",
    "    for i,eta in enumerate(ETA):\n",
    "        #print(eta)\n",
    "        lam,v = np.linalg.eig(np.eye(D) - 2*eta*A.T.dot(A)/n)\n",
    "        LAM[i] = np.max(np.abs(lam))\n",
    "        \n",
    "    return LAM\n",
    "\n",
    "# This number is L_\\eta\n",
    "LAM = compute_largest_eig(ETA, A)\n",
    "\n",
    "plt.plot(ETA, LAM)\n",
    "#plt.plot(ETA, np.ones((N,1)))\n",
    "#plt.gca().set_ylim([0.98, 1.02])\n",
    "plt.ylim([0.997,1.01])\n",
    "plt.xlabel('eta')\n",
    "plt.ylabel('absolute value of the largest eigenvalue')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD+CAYAAAA09s7qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGnNJREFUeJzt3Xd4VGXexvHvMymEhN5LQm/SlQ7BrmsJICggNlARKQm6\nuuu66/puebcobpOEqqCAiqJSJGt3LSRAMAhIFQFBQBCliFLEwHn/mIkvspA6yTnPmftzXbkcJjPh\nZ8i5c+Z57pkxjuMgIiL+FXB7ABERKVsKehERn1PQi4j4nIJeRMTnFPQiIj6noBcR8TkFvYiIzyno\nRUR8LjrcX9AYkwBMBk4A7zmO82y4/w4RESm6Ip3RG2NmGmP2GWPWnXH9VcaYT4wxW4wxD4auHgS8\n5DjOXUD/MM8rIiLFVNSlm6eBq06/whgTBUwCrgbaAsOMMW2BRGBn6GYnwzOmiIiUVJGWbhzH+cAY\n0+SMq7sDWxzH2QZgjHkeGADsIhj2qyngF4kxZhQwCiAhIaFLmzZtijs7e745ztfffU/1+FgaVquI\nMcX+EiJSDDv2H+XEyVO0rFPJ7VEEWLly5deO49Qu7HalWaNvyP+fuUMw4HsAE4EMY8y1wOJz3dlx\nnOnAdICuXbs6ubm5xR7g5CmHx9/ezMT/bKFu/SpMvaULjWrGF/vriEjRjJyVyxeHjvHqPX3dHkUA\nY8yOotwu7K0bx3GOOI5zu+M4Y8p6IzYqYLjvytbMHNGVXQePkpK+hHc2flmWf6VIhNOr3dqoNEG/\nG0g67c+JoevK3aVt6pKZ1pfE6vHcOSuXv7/5CSdP6QdSpCxoidQ+pQn6D4GWxpimxphY4EbglfCM\nVXyNasYzf2xvBndJJP0/Wxjx1AoOHDnh1jgiIp5R1HrlXGAZ0NoYs8sYc6fjOHlAKvAGsBGY5zjO\n+rIbtXBxMVE8NrgTjwzqQM5nB+iXnsWanYfcHEnEV/Q+RXYqUtA7jjPMcZz6juPEOI6T6DjOjND1\nrzqO08pxnOaO4/y5bEctuhu7N+Kl0b0AGDx1Gc/m7EDvpCUSHlq6sY9vXwKhY2I1MtOS6dW8Jg8t\nWMcvXvyYYydU6xcpDZ0u2cm3QQ9QPSGWmSO6cc9lLZm/aheDpixlx/4jbo8lYjWDTult4+ugh2AF\n8+dXtGLm8G58cegYKelZvL1BFUwRiRy+D/p8l7SpQ2ZaMo1rxjNydi5/e0MVTJHi0l6XnSIm6AGS\nasTz0ujeDO2aRMa7Wxg+cwX7v/ve7bFErKLNWPtEVNBDsIL56A0defT6DqzYHqxgrlYFU6RIdD5v\np4gL+nxDuzXi5dG9CQQMg6cuZc5yVTBFikIn9PaJ2KAH6JBYlcy0ZPq0qMXDC9dx/7w1qmCKiO9E\ndNADVIuPZebwbvz88lYsWL2bgZOz2f61KpgiZ6MHvXaK+KAHCAQM91zekqdGdGPv4eP0y8jiLVUw\nRc5Ou7HWUdCf5uLWdVicmkyTmgncNTuXCa9vIu/kKbfHEhEpFQX9GZJqxPPi6F4M657E5Pe2ctvM\nFXytCqYIoNaNrRT0ZxEXE8VfB3Vkwg0dyd1xkH7pWaz6/KDbY4l4ghZu7KOgL8CQrknMH9Ob6CjD\nkGnLmLNsuyqYEtH0828nBX0h2jesSmZqX5Jb1OLhReu5b94ajp7Ic3ssEddoL9Y+CvoiqBofw4zh\n3bjvilYsXL2bQZOX8pkqmCJiCQV9EQUChvGXteTp27uz9/Bx+qdn8cb6vW6PJSJSKAV9MV3UqjaZ\nack0rZ3A3XNW8shrqmBKZNHKjX0U9CWQWD1YwbypRyOmvr+VW2eogimRQXuxdlLQl1CF6Cj+MrAD\nj93QkY8+P0jKxCxW7lAFU/zPaDfWOgr6UhrcNYn5Y3sTGx3gxunLmLVUFUwR8RYFfRi0a1CVxanJ\nXNiyNr97ZT33vrBaFUzxJUfPjbWSgj5MqsbH8MRtXfnFla14Zc0XDJy0lG1ffef2WCJhp4Ub+yjo\nwygQMKRe2pLZd3Rn37fH6Z+RzevrVMEU/9CqpJ0U9GWgb8vaZI7vS/PaCYx+ZiV/fW2jKpjiG9qL\ntY+Cvow0rFaReaN7cXOPRkx7fxu3zMjhq29VwRSR8qegL0MVoqP488AO/H1wJ1Z9foiU9CWs3HHA\n7bFESkxLN3ZS0JeD67sksmBsH+Jiohg6bTlPZX+mCqZYy2g71joK+nLStkEVXklN5uLWtfnD4g3c\n8/xqjnyvCqaIlD0FfTmqWjGG6bd25Zc/a03mx19w3aRstqqCKRZRj95OCvpyFggYxl3Sgtl39GD/\nkRMMyMjmtbV73B5LpOi0cmMdBb1LklvWIjMtmRZ1KjHm2Y/4y6uqYIr3aWvJTgp6FzWoVpEX7u7J\nrT0bM/2Dbdz8ZA77vj3u9lgiBdIJvX0U9C6rEB3F/17Xnn8O7cSaXYdImZhF7nZVMEUkfBT0HjHw\n/GAFMz42ihunL2dmliqY4j36ibSTgt5DzqtfhVfSkrmkTR3+mLmBtLmrVMEUz9FLINhHQe8xVeJi\nmHZLFx64qjWvrt3DgEnZbNmnCqZ4hE7praSg96BAwDD24hbMubMHB4+cYEBGFq+qgikeoWfG2kdB\n72F9WtQic3wyrepVZuyzH/GnzA38oAqmiBSTgt7j6letyAujejG8V2OezPqMm5/IYd9hVTDFHXpm\nrJ0U9BaIjQ7whwHt+dfQzqzd/Q3Xpmex4jNVMMUd2oy1j4LeIted35AF43pTqUI0w55YzpNLtqmC\nKeVKP252UtBbpk29KixK7cPl59XhT//eSOrcVXynCqaUI53R20dBb6EqcTFMvaULD17dhtfW7mFA\nRhZb9n3r9lgi4lEKeksZYxh9UXOeGdmDb479wICMbDI//sLtscTntHJjJwW95Xo3r0VmWl9a16tM\n6nOr+ONiVTClbKlHbx8FvQ/UqxrH86N6MaJ3E2Zmf8ZNTyxXBVNEfqSg94nY6AC/79+Ox2/szLrd\nh7lmYhY52/a7PZb4jFpedlLQ+8yAzg1ZOK4PVeKiuenJHJ74QBVMCS+1buyjoPeh1vUqsyi1D1ec\nV5c/v7qRcc99pAqmhIVOGeykoPepynExTLnlAn5zTRveWP8l/TOy+PRLVTBFIpGC3seMMYy6sDnP\njuzB4WM/MGBSNovXqIIpEmkU9BGgZ7Oa/Ht8X9rWr0La3FX8YfF6VTClRLTdYycFfYSoWyWOuaN6\ncnufJjyVvZ1h05fzpSqYUgJGu7HWUdBHkJioAL/r146Jw85nw57DXDsxi+WqYEox6ITeTgr6CNS/\nUwMWjetDlYrR3PxkDtM/2KoKphSZzufto6CPUC3rVmbRuD78rF1d/vLqJsY88xHfHv/B7bFEpAwo\n6CNY5bgYJt10Ab+99jze2vglAzKy2awKphREj/yspKCPcMYYRvZtxnMje3D4eB4DMrJZtHq322OJ\nh2kv1j4KegGgR7Oa/Ht8Mu0bVuGe51fz+1fWcyJPFUz5KZ3P20lBLz+qWyWO5+7qyZ3JTXl66XaG\nPbGcvd+ogik/pRN6+yjo5SdiogI8nNKWjJvOZ+Oew6SkL2Hp1q/dHktESkFBL2eV0rEBr6T2oWrF\nGG55Moep76uCKdqLtZWCXs6pRZ3KLEpN5ur29XnktU2MfmYlh1XBjHh6Zqx9FPRSoEoVosm46Xx+\ne+15vL1xHwMysvlkryqYIjZR0Euh8iuYc+/qyXff53HdpGwWrlIFMxI56t1YSUEvRda9aQ3+nZZM\nh4ZVufeF1fxu0TpVMCOQFm7so6CXYqlTJY5n7+rBXX2bMmvZDoZOX8aeb465PZaUE23G2klBL8UW\nExXgoWvbMvnmC9i891tSJmaxdIsqmJFCe7H2UdBLiV3ToT6LUpOpnhDLLTNymPzeFlUwRTxIQS+l\n0qJOJRaN68M1Heoz4fVPGDVHFUw/0+9xOynopdQSKkSTPux8/ielLe9u2kf/9Cw27jns9lhSZrR2\nYxsFvYSFMYY7kpsyd1RPjp44ycDJ2SxYtcvtsSTMdEJvJwW9hFW3JjXIHJ9Mp8Rq/PyFNTy8cB3f\n5510eywJI23G2kdBL2FXp3Icz47swagLmzFn+Q6GTlvOF4dUwRRxi4JeykR0VIDfXHMeU26+gC37\nviMlPYtsVTCtp1aVnRT0Uqau7lCfRal9qJkQy60zcpj07hZOnVJY2EwrN/ZR0EuZa167EgvH9SGl\nYwMeeyNYwfzmmCqYIuVFQS/lIqFCNI/f2Jnf92vLe5/so39GFhu+UAXTRtqMtY+CXsqNMYYRfZry\nwt09Of7DSQZNyebllapgipQ1Bb2Uuy6Na5CZ1pfOSdW4/8U1PLRgrSqYltBerJ0U9OKK2pUr8Myd\nPbj7omY8m/M5Q6YtZ7cqmFYw2o61joJeXBMdFeDXV5/H1Fu6sHXfd6RMXMKST79yeywR31HQi+uu\nal+PV1L7UKdyHLfNXEHGfz5VBdOj9A5TdlLQiyc0q12JBeN6079TA/725mZGzcnlm6OqYHqRWjf2\nUdCLZ8THRvOvoZ35Q/92vL/5K/plZLH+i2/cHktOo81YOynoxVOMMQzv3YTnR/XiRN4pBk1eykuq\nYHqKzujto6AXT+rSuDqZ45O5oFF1fvHiGn49XxVMkZJS0Itn1apUgTl3dmfMxc2Zu+Jzhkxdpgqm\ny7RyYycFvXhadFSAX13Vhmm3dmHbV0dImbiEDzargukm9ejto6AXK/ysXT1eSUumbpU4hj+1gonv\nqILpBr1MsZ0U9GKNprUSmD+2N9d1bsg/3trMyNmqYLpCJ/TWUdCLVeJjo/nHkE7874B2LPn0K1Iy\nlrButyqYIgVR0It1jDHc2qsJL9zdi7yTDtdPWcq83J1ujxURtHBjJwW9WOuCRtXJTEuma5PqPPDS\nx/x6/scc/0EVzLKmlRv7KOjFajUrVWD2HT0Ye3Fz5q7YyeCpy9h54KjbY/mXTumtpKAX60UFDA9c\n1Ybpt3Zh+/4j9MvI4r1P9rk9lm8ZPTXWOgp68Y0r29VjcWoy9arEcfvTH/L426pgioCCXnymSa0E\nFoztw8DODfnn25u5Y9aHHDp6wu2xfEO/Nu2koBffqRgbxd+HdOJP17Une8vXpKRnqYIZRlq4sY+C\nXnzJGMMtPRvz4ujenDrlMGjKUl748HO3xxJxhYJefK1zUjUyx/ele5Ma/OrltfzqJVUwS0MvgWAn\nBb34Xo2EWGbd0Z3US1rwQu5Obpi6VBXMUlDpxj4KeokIUQHDL37Wmidv68qO/UdJSc/iXVUwi03n\n83ZS0EtEubxtXTLTkmlQrSJ3PP0h/3xrsyqYxaQTevso6CXiNK6ZwPwxvRl0fiKPv/Mptz/9IQeP\nqIIp/qWgl4hUMTaKvw3uyJ8HtmfZ1v2kpGexdpcqmIXRXqydFPQSsYwx3NyjMS+O7gXA9VOW8vwK\nVTALo5dAsI+CXiJep6RqLE5LpkezGjw4fy0PvLRGFcxzcLQdayUFvQjBCubTt3dn/KUtmJe7i+un\nLOXz/apgno3O5+2joBcJiQoY7ruyNTNHdGXngaOkpC/h3U2qYIr9FPQiZ7i0TV0y0/qSWD2e25/+\nkH+8+QknVcEEtBlrKwW9yFk0qhnP/LG9GdwlkYn/2cKIp1aogplPazfWUdCLnENcTBQTbujIXwd1\nIGfbAVLSs1iz85DbY7lKZ/R2UtCLFMAYw7DujXhpTLCCOXjqMp7L+TyiX9zL6JTeOgp6kSLomFiN\nzLRkejavyW8WrOWXehVMsYiCXqSIqifE8tSIboy/rCUvrdzFoMmqYIodFPQixRAVMNx3RSueGtGN\n3YeOkZK+hHc2fun2WOVKT4y1j4JepAQuaVOHzLRkkmrEc+esXP6uCqZ4mIJepISSasTz8pjeDO2a\nRHqognnA5xXMSN6EtpmCXqQU4mKiePSGjjx6fQdyPjtAysQlrPZ5BVMrN/ZR0IuEwdBujXh5dG8C\nAcOQqct4ZvkOX579+u//KDIo6EXCpENiVTLTkundoia/XbiO+19cw7ET/qtgajPWPgp6kTCqFh/L\nzOHduPfylixYtZuBk7PZsf+I22NJhFPQi4RZIGC49/JgBXPv4eOkpGfx9gZ/VDB9uBoVERT0ImXk\n4tZ1WJyaTJOaCYycnctjb2zyRQVTL4FgHwW9SBlKqhHPi6N7cWO3JCa9u5XhM1ew/7vv3R6rxPQO\nU3ZS0IuUsbiYKB65viMTru/Iiu0H6JeexarPD7o9VolpM9Y+CnqRcjKkWxLzx4QqmNOWMcenFUzx\nHgW9SDlq3zBYwUxuUYuHF67j/nl2VTD1e8lOCnqRclYtPpYZw7tx3xWtWLA6WMHc/rU9FUwt3dhH\nQS/igkDAMP6yljx9e3f2Hj5Ov/Qs3ly/1+2xCqUTejsp6EVcdFGr2mSmJdO0dgKj5qzk0dc3kXfy\nlNtjFUKn9LZR0Iu4LLF6PPPu7sWw7o2Y8t5Wbpu5gq8trmCK9yjoRTwgLiaKvw7qwGM3dGTljoOk\nTMziIw9WMLUZaycFvYiHDO6axPyxvYmNDjB02jJmL9vuuQqmNmPto6AX8Zh2DaqyODWZC1vW5n8W\nree+eWs4eiLP7bHEYgp6EQ+qGh/DE7d15f4rWrFw9W4GTV7KZ56oYHrr0YUUjYJexKMCAUPaZS2Z\ndXt3vjx8nP4eqWBq5cY+CnoRj7uwVW0Wn1bBfOQ19yqYHtsukCJS0ItYILF68FUwb+rRiKnvu1vB\n1GasfRT0IpaoEB3FXwZ6v4Ip3qOgF7FMfgUzJtqUewVTKzd2UtCLWKhdg6pkpvalrwsVTL3DlH0U\n9CKWqhofw5PlXMH02pO3pGgU9CIWy69g5r8KZnlUMLUZax8FvYgP2PkqmFJeFPQiPnF6BbOsXgVT\nCzd2UtCL+MiZFcxwvxF53kmH6IBiwzb6FxPxofwKZnRU6I3Il4Wngpl36hQxUVqkt42CXsSnTq9g\nPhyqYJb2jcjzTjpEBRT0tlHQi/jYmRXM0rwRueM45J1yiI5SbNhG/2IiPndmBbNfRhZvbfiy2F/n\n5Kng0k+0zuito6AXiRA/VjBrJXDX7FwmvL7px/Auirz8oNcavXUU9CIR5PQ3Ip/83lZum5nD/iJW\nMPODPkatG+voX0wkwuS/EfmEGzqSu/0gKUWsYOY/AUubsfZR0ItEqCFdk3h5TNErmD+e0WvpxjoK\nepEI1r5hsIKZ3KJWoRXMvJPBoI/S0o119C8mEuGqxscwY3g37iukgpl3Krh0o81Y+yjoRYRAwDC+\nkApm/hm96pX2UdCLyI8ualWbxanJNKn53xXM/69XKjZso38xEfmJpBrBV8Ec1j3pJxXMH5dudEZv\nHQW9iPyXYAWz408qmLnbgxVMBb19FPQick6nVzB/u3AdoM1YGynoRaRA+RXMS1rXBqBShRiXJ5Li\ninZ7ABHxvvwK5ppdh+iUWM3tcaSYFPQiUiSBgOH8RtXdHkNKQEs3IiI+p6AXEfE5Bb2IiM8p6EVE\nfE5BLyLicwp6ERGfU9CLiPicgl5ExOcU9CIiPqegFxHxOQW9iIjPKehFRHxOQS8i4nMKehERn1PQ\ni4j4nIJeRMTnFPQiIj6noBcR8TkFvYiIzynoRUR8TkEvIuJzCnoREZ9T0IuI+JyCXkTE5xT0IiI+\np6AXEfE5Bb2IiM8p6EVEfE5BLyLicwp6ERGfU9CLiPicgl5ExOcU9CIiPqegFxHxOQW9iIjPKehF\nRHxOQS8i4nMKehERn1PQi4j4nIJeRMTnFPQiIj6noBcR8TkFvYiIzynoRUR8zjiO4/YMGGO+AnaU\n8O61gK/DOE5Z8PqMXp8PNGM4eH0+8P6MXpuvseM4tQu7kSeCvjSMMbmO43R1e46CeH1Gr88HmjEc\nvD4feH9Gr893Llq6ERHxOQW9iIjP+SHop7s9QBF4fUavzweaMRy8Ph94f0avz3dW1q/Ri4hIwfxw\nRi8iIgVQ0IuI+Jyngt4Yc5Ux5hNjzBZjzINn+bwxxkwMff5jY8wFhd3XGFPDGPOWMebT0H+re3DG\nx4wxm0K3X2CMqea1GU/7/P3GGMcYU8tr8xlj0kLfx/XGmAklna+sZjTGdDbGLDfGrDbG5Bpjurs0\n30xjzD5jzLoz7uOlY+VcM3rpWDnrjKd9vtTHStg4juOJDyAK2Ao0A2KBNUDbM25zDfAaYICeQE5h\n9wUmAA+GLj8IPOrBGa8EokOXH/XijKHPJwFvEHxyWy0vzQdcArwNVAj9uY7XvofAm8DVp93/vfKe\nL/S5C4ELgHVn3McTx0ohM3riWCloxnAdK+H88NIZfXdgi+M42xzHOQE8Dww44zYDgNlO0HKgmjGm\nfiH3HQDMCl2eBVzntRkdx3nTcZy80P2XA4lemzHkn8ADQGl28MtqvjHAI47jfA/gOM4+D87oAFVC\nl6sCX7gwH47jfAAcOMvX9cqxcs4ZPXSsFPR9hPAcK2HjpaBvCOw87c+7QtcV5TYF3beu4zh7Qpf3\nAnU9OOPp7iB4BuGpGY0xA4DdjuOsKcVsZTYf0Aroa4zJMca8b4zp5sEZ7wUeM8bsBP4G/NqF+Qri\nlWOlqNw8Vs4pjMdK2ES7PUB5chzHMcZ44jfs2RhjHgLygGfdnuV0xph44DcEHzZ7VTRQg+DD627A\nPGNMMyf0ONojxgA/dxznZWPMEGAGcLnLM52VjpWS8eqx4qUz+t0E17XyJYauK8ptCrrvl/kPtUL/\nLc1D+rKaEWPMCCAFuLmU4VQWMzYHmgJrjDHbQ9d/ZIyp55H5IHimNT/0EHsFcIrgC1CVRFnNOByY\nH7r8IsGlg/KeryBeOVYK5JFj5VzCeayEj9ubBPkfBM/IthH8JuVvjLQ74zbX8tONkRWF3Rd4jJ9u\nME3w4IxXARuA2l79Pp5x/+2UfDO2rL6Ho4E/hi63Ivhw23hsxo3AxaHLlwEry3u+0z7fhP/e6PTE\nsVLIjJ44VgqaMVzHSjg/XP3Lz/JNuQbYTHAn/KHQdaOB0aHLBpgU+vxaoGtB9w1dXxN4B/iUYCuj\nhgdn3EIwmFaHPqZ6bcZw/vCW0fcwFngGWAd8BFzqte8hkAysJBgoOUAXl+abC+wBfiD4SOhODx4r\n55rRS8fKWWcM57ESrg+9BIKIiM95aY1eRETKgIJeRMTnFPQiIj6noBcR8TkFvYiIzynoRUR8TkEv\nIuJz/weWbd6vybFN0gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11b8481d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogy(ETA,LAM)\n",
    "plt.ylim([0.997,1])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If $E$ is twice differentiable, contractivity means that $E$ is convex.\n",
    "\n",
    "For $t>0$\n",
    "\\begin{align}\n",
    "\\|T(x + t \\Delta x) - T(x) \\| & \\leq \\rho \\|t \\Delta x\\| \\\\\n",
    "\\frac{1}{t} \\|T(x + t \\Delta x) - T(x) \\| &\\leq \\rho \\|\\Delta x\\| \n",
    "\\end{align}\n",
    "\n",
    "If we can show that $\\rho< 1$, then $T$ is a contraction.\n",
    "\n",
    "By definitions\n",
    "$$\n",
    "T(x)  = x - \\alpha \\nabla E(x)\n",
    "$$\n",
    "\n",
    "$$\n",
    "T(x + t \\Delta x)  = x + t \\Delta x - \\alpha \\nabla E(x + t \\Delta x)\n",
    "$$\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{1}{t} \\|T(x + t \\Delta x) - T(x) \\| & = \\frac{1}{t} \\|x + t \\Delta x - \\alpha \\nabla E(x + t \\Delta x) - x + \\alpha \\nabla E(x) \\| \\\\\n",
    "& = \\| \\Delta x - \\frac{\\alpha}{t}  (\\nabla E(x + t \\Delta x) - \\nabla E(x) ) \\|  \\\\\n",
    "\\end{align}\n",
    "\n",
    "As this relation holds for all $t$, we take the limit when $t\\rightarrow 0^+$\n",
    "\n",
    "\\begin{align}\n",
    "\\| \\Delta x - \\alpha \\nabla^2 E(x) \\Delta x  \\| & =  \\| (I - \\alpha \\nabla^2 E(x)) \\Delta x  \\| \\\\\n",
    "& \\leq \\| I - \\alpha \\nabla^2 E(x) \\| \\| \\Delta x  \\| \n",
    "\\end{align}\n",
    "\n",
    "If we can choose $\\alpha$ for all $\\xi$ in the domain such that \n",
    "$$\n",
    "\\| I - \\alpha \\nabla^2 E(\\xi) \\| \\leq \\rho < 1\n",
    "$$\n",
    "is satisfied, we have a sufficient condition for a contraction."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lemma:\n",
    "\n",
    "Assume that for $0 \\leq \\rho < 1$, $\\alpha> 0$ and $U(\\xi)$ is a symmetric matrix valued function for all $\\xi \\in \\mathcal{D}$ and we have \n",
    "$$\n",
    "\\| I - \\alpha U(\\xi) \\| \\leq \\rho \n",
    "$$\n",
    "then $U = U(\\xi)$ is positive semidefinite with $$\\frac{1 - \\rho}{\\alpha} I \\preceq U $$ for every $\\xi$.\n",
    "\n",
    "Proof:\n",
    "\n",
    "$$\n",
    "\\|I - \\alpha U \\| = \\sup_{x\\neq 0} \\frac{x^\\top(I - \\alpha U )x }{x^\\top x} \\leq \\rho\n",
    "$$\n",
    "\n",
    "$$\n",
    "x^\\top(I - \\alpha U )x \\leq \\rho x^\\top x\n",
    "$$\n",
    "\n",
    "$$\n",
    "(1- \\rho) x^\\top x  \\leq \\alpha x^\\top U x\n",
    "$$\n",
    "\n",
    "This implies that for all $x$ we have\n",
    "$$\n",
    "0 \\leq x^\\top (U - \\frac{1 - \\rho}{\\alpha} I) x\n",
    "$$\n",
    "In other words, the matrix $U - \\frac{1 - \\rho}{\\alpha} I$ is positive semidefinite, or:\n",
    "\n",
    "$$\n",
    "\\frac{1 - \\rho}{\\alpha} I \\preceq U\n",
    "$$\n",
    "We now see that $\\rho<1$ we have the guarantee that $U$ is positive semidefinite."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "T(x) = M x + b\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\|T(x) - T(x_\\star) \\| = \\|Mx + b - M x_\\star + b \\| = \\| M(x-x_\\star) \\|\n",
    "$$\n",
    "\n",
    "By Schwarz inequality\n",
    "\n",
    "$$\n",
    "\\|T(x) - T(x_\\star) \\| \\leq \\|M\\| \\|x-x_\\star\\|\n",
    "$$\n",
    "If $\\|M\\| < 1$, we have a contraction. Assume the existence of a fixed point $x_\\star$ such that $x_\\star = T(x_\\star)$. (Does a fixed point always exist for a contraction?)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Try to fit with GD to the original data\n",
    "BaseYear2 = 0\n",
    "x2 = np.matrix(df_arac.Year[31:]).T-BaseYear2\n",
    "\n",
    "# Setup the vandermonde matrix\n",
    "N = len(x2)\n",
    "A = np.hstack((np.ones((N,1)), x2))\n",
    "\n",
    "left = -8\n",
    "right = -7.55\n",
    "N = 100\n",
    "ETA = np.logspace(left,right,N)\n",
    "\n",
    "LAM = compute_largest_eig(ETA, A)\n",
    "\n",
    "plt.plot(ETA, LAM)\n",
    "plt.plot(ETA, np.ones((N,1)))\n",
    "plt.gca().set_ylim([0.98, 1.02])\n",
    "plt.xlabel('eta')\n",
    "plt.ylabel('absolute value of the largest eigenvalue')\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Analysis of Momentum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "p(\\tau) & =  \\nabla E(w(\\tau-1)) + \\beta p(\\tau-1) \\\\\n",
    "w(\\tau) & =  w(\\tau-1) - \\alpha p(\\tau) \\\\\n",
    "w(\\tau-1) & =  w(\\tau-2) - \\alpha p(\\tau-1) \\\\\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "\\left(\\begin{array}{c}\n",
    "w(\\tau) \\\\\n",
    "w(\\tau-1)\n",
    "\\end{array}\n",
    "\\right)\n",
    "& = &\\left(\\begin{array}{cc}\n",
    "\\cdot & \\cdot \\\\\n",
    "\\cdot & \\cdot\n",
    "\\end{array}\n",
    "\\right)\n",
    "\\left(\\begin{array}{c}\n",
    "w(\\tau-1) \\\\\n",
    "w(\\tau-2)\n",
    "\\end{array}\n",
    "\\right)\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{1}{\\alpha}(w(\\tau-1) - w(\\tau))  & =  p(\\tau) =  \\nabla E(w(\\tau-1)) + \\beta \\frac{1}{\\alpha}(w(\\tau-2) - w(\\tau-1)) = \\\\\n",
    "\\frac{1}{\\alpha}(w(\\tau-2) - w(\\tau-1))  & =  p(\\tau-1) \\\\\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{1}{\\alpha}(w(\\tau-1) - w(\\tau))  & = \\nabla E(w(\\tau-1)) + \\beta \\frac{1}{\\alpha}(w(\\tau-2) - w(\\tau-1)) \\\\\n",
    "  w(\\tau)  & = -\\alpha \\nabla E(w(\\tau-1)) - \\beta w(\\tau-2) + (\\beta+1)  w(\\tau-1)\n",
    "\\end{align}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Note that GD is sensetive to scaling of data\n",
    "* For example, if we would not have shifted the $x$ axis our original data, GD might not have worked. The maximum eigenvalue is very close to $1$ for all $\\eta$ upto numerical precision"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWwAAAD7CAYAAABOi672AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAEY1JREFUeJzt3d1vHOd1x/HfmRm+qZEskVrDLhxb\npoG2QBwEYFeOnYsWKGigL0BRFAp64yJXof8DGblK39DAuutN0LK9Lxzzqg2KXjAtCiRFnVAqijZB\n0No03Lito5VkWrYrSuTu6cXO7C6pedmKnN1nd78fwDBXXJHPg13+dHjmzIy5uwAA4YvGvQAAwHAI\nbACYEAQ2AEwIAhsAJgSBDQATgsAGgAlBYAPAhCCwAWBCENgAMCGS0/xiFy9e9EuXLp3mlwSAqXf9\n+vVb7t6oet6pBvalS5e0s7Nzml8SAKaemb03zPNoiQDAhCCwAWBCENgAMCEIbACYEKWBbWZrZuZm\n9k7635+NamEAgKOqpkSW3d2kbnhL2qt/SQCAPKUVtrtvDzxsuvtuzesBABQYqodtZuuSvlXXIr73\n9i298hdv6YOP9uv6FgAw8YY96Piyu+e2Q8xsw8x2zGyn1Wo90iJufXJf3337lj59cPhIfx8AZsGw\ngb1W9Al333T3prs3G43KMytzJVF3Ge0ONwQGgCKVgW1mq3UvIo5MknTYJrABoMiwFXatBxuTNLCp\nsAGgWOXFn9LJkFfrXEQcpxV2p1PntwGAiRbEmY5U2ABQLYjA7vWwCWwAKBREYDMlAgDVgghsKmwA\nqBZEYPd72Bx0BIAiQQQ2c9gAUC2IwE5ipkQAoEoYgU0PGwAqBRHYMVMiAFApiMCmwgaAakEEdsyU\nCABUCiKwqbABoFoQgR1zLREAqBREYGenpjOHDQDFggjsmDlsAKgURGDTwwaAakEENlMiAFAtjMA2\nKmwAqBJEYEeRKTJ62ABQpvKejma2JmlVktx9q7aFRBEVNgCUGKbC/loa1KtpeNcijowKGwBKlFbY\nZnZF0g8kyd2v1bqQyJjDBoASVRX2ZUkrZrZmZlfrXEgcG1MiAFBimJbIbXe/IfUq7iPMbMPMdsxs\np9VqPfJCksjoYQNAiarAvi1pN/14T92K+wh333T3prs3G43GIy+EHjYAlKsK7C2lEyKSzivtZ9eB\nKREAKFca2O6+K2kvbYWs1DnWR4UNAOUq57DdfTP9sLawluhhA0CVIM50lLIKmykRACgSVGAzhw0A\nxYIJ7CSmhw0AZYIJ7JgpEQAoFUxgJ0yJAECpYAI7jkyHHHQEgELBBDYVNgCUCyawY+awAaBUMIFN\nhQ0A5YIJ7DiKmMMGgBLBBDYVNgCUCyaw45gpEQAoE05gGxU2AJQJJrC5Wh8AlAsmsLkeNgCUCyaw\nk5gKGwDKBBPYVNgAUC6YwE6iSIdtpkQAoEgwgU2FDQDlgglspkQAoFxlYJvZ6+n/N+pcCBU2AJQb\npsLeMLN3JO3WuZCswnYntAEgTzLEc77q7lt1LySOuv92dFyKre7vBgCTZ5gKe9XM1s3sap0LSdKU\n5noiAJCvMrDd/Zq7b0taMbP14583sw0z2zGznVar9cgLiaNuYNPHBoB8pYGdhvGV9OFtSavHn+Pu\nm+7edPdmo9F45IUkUVZhE9gAkKeqwt6RtJ1+/Fz6uBa9CpubGABArtKDju5+I62y70h6x91v1LYQ\nKmwAKFU5JeLum6NYSDYlQg8bAPIFdaajxJQIABQJJrCZEgGAcsEEdn8Om8AGgDzBBDYVNgCUCyaw\nEwIbAEoFE9hMiQBAuWACmzlsACgXTGD3e9iM9QFAnmACu1dhc2o6AOQKJrCZEgGAcsEENnPYAFAu\nmMBmSgQAygUT2EyJAEC5YAKbKREAKBdMYFNhA0C5YAKbKREAKBdMYCfpQUfmsAEgXzCBHcdU2ABQ\nZujANrOrdS6EHjYAlBsqsM1sXdLLdS6EKREAKBdMS4QKGwDKVQa2ma25+3bdC2FKBADKDVNhL9e+\nCg1MiRDYAJCrNLCHqa7NbMPMdsxsp9VqPfJCqLABoFxVhb1qZlfMbEPSspmtHX+Cu2+6e9Pdm41G\n45EXwvWwAaBcaWC7+5a7b6UPz9e6kMhkxpQIABQZakokraKfc/cbdS4miYweNgAUCGasT+r2selh\nA0C+oAI7iSIqbAAoEFRgU2EDQLGgArvbw+agIwDkCSqwqbABoFhQgZ1Exhw2ABQIKrDjmAobAIoE\nFdhMiQBAsaACmx42ABQLKrCZEgGAYkEFdmRU2ABQJKjATmKuJQIARYIKbHrYAFAsqMBmDhsAigUV\n2FTYAFAsqMDuzmEzJQIAeYIKbCpsACgWVGBzxxkAKBZUYFNhA0CxoAKbOWwAKJZUPcHM1tMPX3b3\n1+pcTBxFVNgAUKC0wk7D+svuvi1pzczW6lwM1xIBgGKlFXYa1Nvpw1V3v1HnYuLI1ObEGQDINVQP\n28yuSnq15rUwJQIAJYYKbHe/JulVMzt//HNmtmFmO2a202q1TrQYpkQAoFhVD3uwb70raeP4c9x9\n092b7t5sNBonWgwVNgAUq6qw1yUtpx+fVze0a8OUCAAUqwrsTUmrZrYhSe6+VediunPYTIkAQJ6q\nKZE9dUN7JPJ62A8OO5pPgjq/BwDGIqgkPN7Dfv/D/9Xnvv63+pef7I1xVQAQhqACO45M7lInDe3/\n3tvXQdv14w/ujnllADB+QQV2Epkk9ars/YO2JKn18f2xrQkAQhFUYMdRdzlZH/teGtg3CWwACCuw\n+xV2d1Ikq7Bv3iWwASCowI7TwO5V2A+yCnt/bGsCgFAEFdhJfLSHTUsEAPqCCuyHKuyBwHbnDEgA\nsy2owH5oSiRtiTw47OjuvcOxrQsAQhBUYPemRNpHK2yJPjYABBXYx6dEjgY2fWwAsy2owH54SqR/\nISgqbACzLqjAzjvT8Ylzi5KYxQaAoAI7b0qkcXZBi3MRp6cDmHlBBfZDc9gP2lqaj/X42UV62ABm\nXlCB3b+WSP+g4+JcrMfPLtDDBjDzggrsXg+73e9hL81FevzcAhU2gJkXVGD3etje72EvzXVbIi0O\nOgKYcUEFdpJz8ael+ViNswv6+P5h72JQknR3/4DT1QHMlMrANrON9L/X615MHD188aeshy31Z7E/\nunegl/74O/qT7/xH3UsCgGCUBraZrUvadvfs7unrdS4mOXZq+n7WEslmsdM+9vX37ujTB2198+/f\n0W7rkzqXBADBqKqwVyVlIb2bPq7NYIV90O7ooO1pDzutsNM+9lu7dzQfR1pIIn39r354pDXy4LCj\ng3bn4S8OABMuKftkWlln1iS9Ueti4n4PO7vbTNbDlvotkbfevaMvfPYx/cbnn9Tv/fWP9Df/+oF+\n/fNP6M3r7+uPvv0jtTuuL66u6EvPrej8mXn99O6+PvhoXw8OO1qaj3VmPlYSRzpsd9TuuDruiswk\nk0wmV35v3GR1bh/ABPvScyv6pZ9r1Po9SgM7Y2Zrkm64+42cz21I2pCkp59++kSLiQcu/pRd+Glx\nLtbymXklkenmx/f16f1D/dt/faRXf3lVr7z4jN68/r7+4Ns/1Ld2fqJ/+PeWLl+6oJ9/4qy+9/Zt\n/d2Pb/a+9rnFRItzse4dtHXvQVuHHVcSmZLYeiHdcUmuNLiP4vAmgDKLc1EYgS1p3d1fy/tEWoVv\nSlKz2TxRrg1Oidw/6LY1luZiRZHp4mcWdPPuff3zf+7psON64dkVJXGkP/yt5/Xb3/xHfbx/R7//\nm5/T7774jKL06/zPR/e0f9DRE+cWtTQfH1+3zKiYAUyOysA2sw13v5Z+vO7u23UtZrCHfW+gJSJJ\nj59bUOuT+/r+u7cVR6ZffOaCJGnt6Qv6y6++qM8uL+mpC2eOfL0nH1sq/F6ENYBJM8yUyOtm9o6Z\nfVj3YnpTIh3vzVwvzaWBfXZBN+/u65/evaPnf/acPrPQ/7fmpedWHgprAJg2VQcdtyVdGNFacivs\nxTSwG2cX9f1372j31qf6ykvPjGpJABCMYXvYI9HrYbc7D7dEzi7o7n73vo4vPLsyngUCwBgFdWp6\nPHB51f3jLZFzC73nXb40sqIfAIIRVGAPTon0KuxeD7t7tuMvPHFW58/Mj2eBADBGQQV2bg97vrvE\n7OSZF55dHs/iAGDMggrssimRZy/+jJ66sKRfe/7Jsa0PAMYpqIOOaYHd7WEfmxJ5bGlO333tV8a1\nNAAYu6AC28yURKZ2p6N7HWkuNs3FQf0SAABjE1RgS1IUWfdqfYeuxSSu/gsAMCOCK1+TyNRudw86\nLs4T2ACQCS6w47TCzm5eAADoCi6wuz3s7pQIgQ0AfcEFdhxFvTlsWiIA0BdcYPemRA7aWpoLbnkA\nMDbBJSI9bADIF1xgJ/FAD5uWCAD0BBfYWYV976DdO8sRABBgYGdz2LREAOCo4AK7NyXCWB8AHBFc\nYCeR6bDT0f5hhx42AAwYKrDNbK3uhWTiyLR/0Fa74/SwAWBAZWCnd05/cwRrkdStsD+53713Iy0R\nAOirDOz0zum7I1iLpG6F/Ul6s11aIgDQF14PO6bCBoA8Jw5sM9swsx0z22m1WideUBxF+jitsOlh\nA0DfiQPb3TfdvenuzUajceIFJZHp/mFHEi0RABgUXEsku3O6REsEAAYNMyVyRVIz/X/tkoHAXuRq\nfQDQU3lPR3ffkrQ1grVIosIGgCLBlbBHK2wCGwAywQV2HPWXxEFHAOgLLrATWiIAkCu4wI5jWiIA\nkCe4wM4q7PkkOnIAEgBmXXCBnYU07RAAOCq4wE4IbADIFVxgZ1MiTIgAwFHBBXZWYXPAEQCOCi6w\n+z3s4JYGAGMVXCr2eti0RADgiOACO5vD5qAjABwVXGDTwwaAfMEFdm9KhMAGgCOCC2x62ACQL7jA\n5kxHAMgXXGDTwwaAfMEFdkxLBAByBRfYSTrWt5gEtzQAGKvgUpFriQBAvsqb8KZ3S9+TtObu12pf\nED1sAMhVWmGb2Zokufu2pL3scZ2YEgGAfFUtkd9Rt7qWpF1J6/UuhzlsAChSFdjnJd0ZeLxS41ok\n9Svrs4tzdX8rAJgoJz7oaGYbZrZjZjutVuvEC/ri6or+9JU1feGpx078tQBgmlQF9p6k5fTj85Ju\nH3+Cu2+6e9Pdm41G48QLiiPTrz7/pMy4AS8ADKoK7DckraYfr0rarnc5AIAipYHt7jckyczWJe1l\njwEAo1c5h+3um6NYCACgXHBnOgIA8hHYADAhCGwAmBAENgBMCAIbACaEufvpfTGzlqT3HvGvX5R0\n69QWMzlmcd+zuGdpNvc9i3uW/v/7fsbdK888PNXAPgkz23H35rjXMWqzuO9Z3LM0m/uexT1L9e2b\nlggATAgCGwAmREiBPatnVM7ivmdxz9Js7nsW9yzVtO9getiYfmZ2NbvN3KhvPQfUwczWBq+xlPe+\nPs33ehAVtpldMbN1M7s67rXULb1++IaZvT7wZ1O///QCYi+nH4/81nPjYGZr6Wt7ZeDPpvq1Htjf\nRs6fTdWe0/f0mwOPH3pfn/Z7feyBPSs/vFLvBd5OL6i1mr6JZ2b/A0Z+67kx+Zq7b6n7Wp/6D29o\n0v3spvvbnfY9Z/sc+KO89/WpvtfHHtianR9eqXtN8Wx/u+njqd9/+mvj4LXUR37ruVFLq+ofSJK7\nX0t/bZ7611pS9pvj6gztOZP3vj7V93oIgT31P7yZ9O482cGINUk7mo39L1c/ZepclrSSVplZK2Cq\nX+s0oHfN7EP19znVex61EAJ75qS/Ft6YhRtC5FTX0hC3npsStwduAnKl6smTzszOq/vafkPSn5vZ\nasVfmTZ57+tTfa9X3sBgBGblh3fQuru/ln487ftfTX9wlyUtp/9YvSEpOwtsWm89d1v9/uaeuhX3\ntL/WG5K+4e57ZrYrKZuOmOY9Dyp6X5/aez2ECnum7htpZhsD4z7rmvL9u/tWeuBN6v7Azsqt57bU\nf13Pq9vPnurXelD6mu9pivec/tbUzH57yntfn/Z7PYg57HQEaFfdAxVTO2g/MAZ0R92q48vuvj0r\n+5816et6R9Ll7DeqaX+t0379rqTlbH/TvudRCiKwAQDVQmiJAACGQGADwIQgsAFgQhDYADAhCGwA\nmBAENgBMCAIbACYEgQ0AE+L/ALoWXgJPnEGvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11ec7ceb8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "w = 7\n",
    "alpha = 0.7\n",
    "EPOCH = 100\n",
    "W = []\n",
    "for tau in range(EPOCH):\n",
    "    W.append(w)\n",
    "    w = w - alpha*(2*w - 4)\n",
    "\n",
    "plt.plot(W)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD7CAYAAAB68m/qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAHK9JREFUeJzt3XuMXOd53/Hfe87M7JXL5WVFyZQo\nipRkKZZRa00prd2qqbJygjR/xIVSI+glRQMwDRLHaYPIMWC3BdrGsYwaTQ2kBZ07WseCZSMOmqZO\naBWOlMaxKMq2FCm2pZVEmaKk5WV3xb3O5e0f57xnzjkzs2Q4l933zPcDCDOzXO2ew9398dnnPO97\njLVWAAD/Bdt9AACA3iDQAaAgCHQAKAgCHQAKgkAHgIIg0AGgIAh0ACgIAh0ACoJAB4CCKA3yk+3f\nv98ePnx4kJ8SALz31FNPnbfWzlzp/QYa6IcPH9apU6cG+SkBwHvGmFeu5v1ouQBAQRDoAFAQBDoA\nFASBDgAFQaADQEEQ6ABQEAQ6ABSEF4H+G4/P64O///R2HwYA7GheBPqLCyv62vyF7T4MANjRvAj0\nMJDqDW5mDQBb8SLQS0GgWr2x3YcBADuaF4EeBkYU6ACwNW8CvdagQgeArXgT6OQ5AGzNi0AvUaED\nwBV5EeiBiXro1tJIB4BOvAj0UmAkMboIAFvxItDDMAr0GoEOAB35EeiGCh0ArsSPQHctF3roANCR\nF4Ge9NDrBDoAdOJFoLsKnR46AHTmSaBHh9mg5QIAHXkR6CUqdAC4Ii8CPaCHDgBX5EWgl5hyAYAr\nKnX7AYwxs5KOSJK19tGuj6iNZGyR/VwAoKNeVOgfiYP8SBzuPceUCwBcWVcVujHmQUlPSpK19uGe\nHFEbIXu5AMAVdVuh3yNpnzFm1hjzULt3MMYcN8acMsacWlhYuKZPwtJ/ALiyXrRcLlhrT0tJxZ5h\nrT1hrT1mrT02MzNzTZ8gHkPnNnQAsIVuA/2CpPn4+aKiir3nTFyhsx86AHTWbaA/qnjCRdK04n56\nrwVxoFOhA0BnXQW6tXZe0mLcatnXr7HF+JooFToAbKHrOXRr7Yn4aV/CXKJCB4Cr4cVK0TjP2ZwL\nALbgRaA3K3QCHQA68SrQXZ4vrVZ1bmltG48IAHYeTwI9enQV+vv+y1f1dz7+2DYeEQDsPF4EerOH\nHj2+sbyxfQcDADuUJ4HefmHRerW+HYcDADuSF4Ge76E7F1Y2t+FoAGBn8iTQo8f8lMvl9do2HA0A\n7EyeBHr7hUU1bngBAAkvAr3TwiK20wWAJi8CPehwUbTKTaMBIOFVoOcLcip0AGjyJNCjx3zLhR46\nADR5Eej5hUVOjZYLACQ8CfT2PXRaLgDQ5EWgd1pYVCPQASDhSaBHjy099Do9dABwPAn0TguLqNAB\nwPEi0FlYBABX5kWgd15YRMsFAByvAp2FRQDQmSeBHj02rM1U6fTQAaDJi0BXamFRuuvClAsANHkR\n6K7lImtVp0IHgLa8CvSGzU660EMHgCZPAj16jHrozbe7Cv3Jly/qvb/6mN5ar27D0QHAzuBFoJtU\nhZ6uyhvx809++ds6u7imZ84ubcvxAcBO4EWgB0kL3WZbLvHzUvwO7KYLYJh5EuiuQreZ0HYVehgH\nOguNAAwzzwI9e1HUdV9cS2aZHjqAIeZFoKf3cmnXcgnjP19erw360ABgx/Aq0K1VZg7dtVzcjYs2\nqvVBHxoA7Bg9C3RjzEO9+lh56c250mOLbuLFbQdQ5ZZ0AIZYTwLdGDMn6YFefKx2Oi4sip+7Hjpb\nAQAYZl60XNILi9rNoTtMuQAYZl0HujFm1lp7shcHs8XnkNS6OZfLcxfsVbYCADDEelGh7+3Bx7ii\nwHReWOSqdlouAIZZV4F+NdW5Mea4MeaUMebUwsLCNX+uwJiOLRf3Ni6KAhhm3VboR4wxDxpjjkva\na4yZzb+DtfaEtfaYtfbYzMzMNX+iKNCzdy2qJ2OLLtCp0AEMr64C3Vr7qLX20fjldA+OpzPTesei\n1pYLFTqA4dWTKZe4Cj9qrT3di4/XTtRDzy4sck+bLRcqdADDy4uxRSlqudjc5lz1fA899fpTf/Jt\nvb60PvDjBIDt4lWgd1pYlJ9yOfXyRf3Xx17Qx7707OAPFAC2iTeBbkzr5lyNDhdFL61uSsrOrANA\n0XkT6FHLpf2USyM3tnhhJQr0XaOlwR4kAGwjjwK9TYUeP3f3Fq3FDfa1zXrmzwFgGHgU6Ca+Y1Hr\nDS7yY4sbtSjYN6pMvQAYHt4EutlqYVHu0QX6JmOMAIaIR4HeupeLe17PtV42alHLhbl0AMPEm0B3\nC4vSLZf8RVH3ejOu0Al0AMPEo0CPe+jWvU71zhv5Ct0FOhdFAQwPzwK92WYphUEyZ97ILSxyF0PZ\nThfAMPEm0N3CItcvLwempXdez/XQN6nQAQwRbwLdLSyyqQo9v31urUMPvd6w+pFfe1yf/sp3B33Y\nADAwHgW6u6do9LocmqT9kr8oupEL9JfOX9Zz55b1n//0OwM+agAYHI8CPdtDL6cq9PxKUfdYjYN9\neb026MMFgIHzJtBN7gYXpdCo3rBtxxjddIvroa9sEOgAis+jQI/2Q3ctl1IQTbnUUoGevzjqKnUC\nHcAw8CbQk4VFScslmnLJ7I9ez44vupbL5Y36gI8WAAbPo0A3md0WS0GgRsMm1bjUusCoSssFwBDx\nJtBN/qJoKVDd2iS8K6WgeVE06aE3ZK3V6mazQnfV+5tvrSfb7AJAEXgT6IHbnMuNLQbZi6IjYdAy\n7RI9t8lcutQcabz3P31F9/7KyQEdPQD0n0eBnl/6bzIXRUfKQbJ5V7oNU603Mpt0baTC/S3GGQEU\niEeBnr1jkZtDd69HSqGkKODTm3JVazazL7rbFsBZr9J2AVAM3gR6/gYXpXgvF1eNj5SiU6nnK/RG\nI9tyqTYym3YtrlYHcPQA0H/eBHqQu8FFKcxOuVTiQK81Gqo1GgoDIylquWzmWi6rqar80urmoE4B\nAPrKm0A3uXuKVsIgU6FXUhV6rWE1Xo5aMNVa/qJoPTPdsrwWVeif+/oZ/dLnvzmQcwGAfvAm0JsL\ni6LXbul/LddyqTWsanWr0UoU6Jv1RsuUSzrQXbX+y198Rp9/6ns6t7Q2iNMBgJ7zJtBNm4VF6ZWj\n2Qq9oTFXoeenXKqNzFx6fhb95fOrfT0PAOgXbwI9mnJp7tNSjiv05kXRZoDX6jYT6C0tl2pzXDG/\nivT85Y2+ngcA9ItHgR5tzmXTLZd0Dz2MWy71qA3jWi7uouhE/Hqjlq3QVzfrmVFGF+h/8eIF/cAn\n/y8BD8AbXgV6ZmFRbi+XkXJ0Km6iJbkoWo8uik6OliRFFXo+0C+nFhhdXImmXv79H/6VXr6wqsee\nf7PPZwYAveFNoOfvKVrJ7eXiLoq6G0SP5Sr0XaPl5M8zF0U3a7qcaru4qRdX4Z+52Oyps8kXgJ3M\nm0B3FbpNLSyyqZ66uyjq2icu0Ddr0UXRyRFXobe2XNJbACzFgf7G0rqkaBMvSfra/AW94999WV88\n/b1+nSIAdMWjQHebc7mWS7RwyK36dBdF112Fnmu57Mq0XGrx/xO0VOhLa1U1GjYJ8oW3oh76Y38d\ntV7+8JuvJe/72uJay1YCALBdvAl0N7ZYT+3lIjV75i0Vetnt7dLIBPpmag59/+RIpoe+a7Sk5fWa\n3tqoJfPub8aB/p033pIknb0Uzam/dH5F7/nVx/TBzz7dv5MGgL+BUrcfwBhzPH561Fr74W4/Xif5\nhUVhGFXobiTRTbmst+mhV+tWY+WSAhO1XOoNq3JoNDVW1spGXStxxX5wekxLa1Utxfu7hIFJAv31\nuAXzvUtrstbqK8+/IUn6k+feULXeUGiMfvazp7VnoqJfef87+/XXAAAddVWhG2PmJJ201p6QdCR+\n3RfGRHPn1loFRgqN26slO+XiKvTRVMtlo9ZQpRRopBQmPfSxcqjxSqi1ai3pobtAX1yLJl1uu25S\nFy5vqN6wen05CvS1al3nL2/qhTcvJ8c2v7Cip19d1B8/+7o++5dn9NevL0uSPvoHz+hHfu1xXYon\nZ9Y263r27FJyo2sA6KVuWy5HJLkQn49f90XcMle9YRUYk9l8S2pToWcWFtU1Ugo0Ug60Xo32chmr\nRIG+slFPeugH90SBfimu0G87sEsNG/XKF1ermj00LUl69dKqXly4rOnxaHLm+XPLevrMpeRYv/rt\nBb2+tK7/8bUzeu7csj735Kuy1uonf/vr+tFPP6HPPD4fvd93FvSP//tf6EvfOCspuh7wyJNnMh/r\nlQsrOrvY3I6gWm+w5S+AtrpqucSVuTMr6ZHuDqez5j1FpSAwCuIK3fXQR8ruoqir0JsLjap1G1fo\ngTbj3RbHKyWNV0K9ubyhy+s1BUY6MDWqzVpDb8TV+K0zk5KiwJakew7v1ekzi3r14qpeXFjRA3ce\n0Je+8ZqeP7esVy+t6sY9YxqvhHrihfNJy2f/5IgeefKM3n3zHn39pYuSpE9/5QXdf8cBfehzT2tx\ntapTr1zU5EhJX3z6rP7oW+ckSR/9h3fq3NK6fvOJlxQY6SfuPaR9ExX9zv97WW9t1PTeo/v1t27a\nrce/e17zCyu6/cCk3nlwt+bPr+il8yu6fmpUt143qUurmzpzcU0TlVA3TI9po1rX4mpVm/Hkz3jF\n/dZSUzkMNFoOZSRV49+GwvjvOr33vDFGgZFM/DzP/QZik9fX9jVv86EBb/3SD71d73jb7r5+jq57\n6JJkjJmVdNpae7rNnx2XdFySDh06dM2fI72wKDDNir2aTLnEFXruoqibQy+HRiOlMK7Qaxorh5qo\nlLQST7lMjpQ0NRZV3GcuRLPnR6+bkCT91WtRoL/75j2SpGe+t6SLK5t6+/W7dNuBST13blnffeOy\n7rllr/ZPVvT7Xz+j1c26js5M6KfvO6qHvvAtfehzT2tqtKTf+6nv1/t//c8196mvqhwa/cHPvlcf\n+eIz+qnfPSVJ+vkfvE3Pnl3Sf/yj5yVJP3HvTSoFgf7nX76ihpUe+L4DOjIzof/9zDk98cJ53XVw\nSu+/+2Dym8At+yc0e2iPzi2t6U+fe0PT42XdvG9Cq5s1ffPVRY2WA+0Zr2jXaEkrGzWdv7yhkVKg\nsUqozVpDy+tVWRtddDZGSZCHxiThbSUp/sfVysqoNXldGJv8G64WbSkUTK3e/+/pngS6pLlOF0Tj\nKv6EJB07duyaz8gtLGrkWi7uominhUXuImglDDVSCpIe+ngl1Fgl1Fo8h75rtKzdcaC/Ei8mOhpX\n6C7Qb9k/of2TFf3ZdxeSP7/zhik9+lQ0mz57aFqH9o7rt//8ZT31yiX93D+4VT901/X62Jee1bml\ndX3w/lv1rpum9S/fe4t+84mX9G8eeLveddO0futfHNMn/8+39Xdv269/NHujqvWG/te3XtPB6XHd\ne8teSdK/+oGjqtUbunlf9I/ML//wHbq8UUsWTAFAT6ZcrLUPx8/nrLV9ufNyYIwajagqDI1RkOuh\n5yt0d1HUre6sxD10F+i7RkuaGHEVelWTI6Uk0M9cWNHkSEnXT41Kkp49uyRJumF6TDfuGdc3Xl2U\nFAX8nTdMJcd496E9uvW6yeT1j919ULvHyvpv/3RWf/ad8/rpv39UkvSxH/0+/eL7btd4Jfrrv2H3\nmD71gXcl/185DPT+u2/MnP/B6bHMa2MMYQ4go6tAj6daPmGM+bCkvZJ+vCdH1YbbbbFhrUybKZdK\nrkJ3Ae9WhZZDo9FSmNzg4sDUiMbKodarDS2v1TQ5WtJUPKv+ysVV7R6LKvYwMHp9eV27RkqaHCnp\npr1RoJdDoxv3jOnew3uTY7zrbVMqhYG+8DPv0Xq1noT7/Xcc0P13HMicjwtzAOiVbi+KnpS0p0fH\nsqUgaO6HHgTNCn2jll0p6l6Xw0CVMEgqdDflslFtaC2+KDoxEv0/C5c3dHB6LKnQF1erOjg9piAw\n2jdR0ZtvbeiG6ahadxdK90+OqBQGeueNu/Uffuwu3XH9LpXiSRvXaweAQfJmpWi0fa6SC3SByY0t\nJhV6VJGHgVEpNEmF3jKHXgmTKvnN5fVMy0WS9oxXJEnXTY1Ikq7fHbU83veOAwqM9OC7my2Rf/a3\nb9Y9qUodALaDN7/3B263xUbUP46LYVVzK0WbFbpROQySGfNyGMQXRZtTLuPxhdPl9eyUiyTtjmfM\nb943oWfPLuttu6MK/c4bpvT0v31fstkXAOwUXlXoDetmo9VSobuVoutJhR6oHJpkIy43h75edXPo\nYaaPPTlaUjkMkpDfNxFV6LfEUyV3HWzOj7reOgDsJN6UmSYzh54aW3QXRcPslEspiCr0lY168ucj\npTCZs3YrRR1XcUc7MNZ1IJ5w+bn7b9XMrpFMiwUAdiJvAt1tn1tvRNW5m3JJ5tDL2e1zS6HroWfH\nFhfjZf3j5TC5KCop2Y3R/f/X7Yp656PlUD/5nsN9PjsA6J5nLZco1INAHefQN5IKPYh76M2Lom42\nXYrGBsfKqZZLXKG75e23HdjV5zMCgN7yKNDVHFvcYsolqdADo3IQNCv0+KKoM1bJV+jRRdDP/PNj\n+uF3XK93HuzvngsA0GvetFxMvFK0buOWi5ty6VChtx9bbAa6W/rv7I0vgt53+4zuu32m7+cDAL3m\nUYXenEOPNufqcFG02lxY5O5qJDXn0J2xSqip1NL5fZOVvp8DAPSTR4GuZGwxuzlXPdrK1RiVApNZ\nWFQOm6OFlTBIRhulaDfGdE/dVegA4Ct/Aj2ILoomN7hI7eXiwj0MTMvCIqcc5lsu2W7TPgIdgOc8\n6qErc4MLk7oo6tovpcBoZTPdQ28G+EjLlEv0/Mu/cF/88VkoBMBv3gS666En9xRN7YdeSlXoTrQ5\nV6rlkuuhu0B/+/WMJwIoBn9aLsleLrZlysXNpKcr8jAwKgXZi6JTY81/v9L7tgBAEXgU6J3vKdqu\nQi8FRuVUz7wSBpndFNP9dQAoAm9SLbuXSzO8q7XmRdFyKtiNMclrE7//9DgXPgEUlzeB7opv13IJ\n2lwUDcPmxVGpWYVXwkDGmEyFDgBF49VFUUmqNbI3uNisN7SrFJ2G65kngV6KHt22ABOVUP967nYu\nhAIoJI8CPXqsN6xKQXNhUfqiaJh7HI2nWtz8uTFGH5q7bZCHDQAD403Lxc2J1+IWi7ummR5bzLda\n3F4tXAAFMAy8SbpMyyU15eKmXqQ2FXq8kCjeERcACs2jQI8e6w2bTK04oelQoceBXifRAQwBjwI9\nO9USpJbq5ytz95i0XLj/J4Ah4E2gm1SFHgYmabNIzQB3Uy5uqsVV6CV66ACGgDdJF6R2VwxMs80i\nNVstpTDbcnE99PQ2ugBQVB4FevQY9dCNUtu0tFwUdZtyTY9HC4ncTS8AoMj8CfQgu7AoXaF3uih6\nV3xf0EN7xwd5qACwLbxZWJTMoTcaCoLclEtSoUdB7gJ9cqSkL/zMe3TL/okBHy0ADJ43gZ60XOrR\nXi6mzZRLc8l/8xePd9+8Z3AHCQDbyJ+Wi7so2mhk7ikqpSr0MNtDB4Bh4lGgR4+1NlMuLRU6Y4oA\nhpA3yZde+h8E2SkXF+4hgQ5giHmTfOmVoUF+yoUKHQA8CvT03Hl+L5fcPUUrJXroAIZP11MuxpgH\nJS1KmrXWPtz9IbWX37tlyykXKnQAQ6ir5DPGzEqStfakpEX3uh/SAW5yPfP8Y4VABzCEuk2+Dyiq\nziVpXtJclx+vo/SGie55cjE0fnRV/BT3DgUwhLoN9GlJF1Ov93X58ToK2iz1d311Wi4AMICLosaY\n48aYU8aYUwsLC9f8cdIVuslV5vmWS73BZlwAhk+3gb4oaW/8fFrShfw7WGtPWGuPWWuPzczMXPMn\nMrmxxfSj2zbXVei1BncoAjB8ug30RyQdiZ8fkXSyy4/XUXbKJX5b7sYWxw5H/7bcfYj9WwAMn67G\nFq21p40xx4wxc5IWrbWne3RcLbIXRbMtFleZ33f7jE59dE77J0f6dRgAsGN1PYdurT3RiwO5kqDN\n2GKz5dL8RYMwBzCsvBkHSeV50nJxj9xiDgA8CvT8Xi6SZJRtvQDAMPMy0F3LpWGjaRbmzgHAq0Bv\nPg+TQI9el6jQAcCfQM/OoUePNq7QS1ToAOBPoGfGFoN8y4UKHQD8CfSg9aKoa7lwURQAfAr0Nrst\nJhV64M1pAEDfeJOE7fZyad6piAodALwJ9Mwcem7JPxdFAcCrQG997kK+TA8dAHwK9M73EOWiKAB4\nFOjpvVxMbsplrBJuwxEBwM7iTaAHbRYWbdTqkqTpscp2HBIA7CheBrpb+r9Ri241Nz3OTaEBwKNA\nTz+PXqxXowp9aoxABwBvAt1kdluMHu+7PbpH6a6Rru/TAQDe8yYJM7stxi9+/Z/M6s3ljcy2AAAw\nrDwK9NaVouOVkg7v9+YUAKCvvGm5pGfNDQU5ALTwJtBNm4uiAIAmbwI9HeLcoQgAWnkZ6Cz1B4BW\nHgV68znb5QJAK28C3WQqdG8OGwAGxptkzFTotFwAoIU3gZ7um9NDB4BW3gR6+q5EVOgA0MqfQKdC\nB4AteRnoJS6KAkALb5KRHjoAbM2bQE+PLTKHDgCtvAn0NCp0AGjlZaAz5QIArbwMdCp0AGjV9d0h\njDHH46dHrbUf7vbjXQ2mXACgVVfJaIyZk3TSWntC0pH4dd+R5wDQqttoPCLJhfh8/LrvqNABoFVX\nLZe4MndmJT2Sf5+4JXNckg4dOtTNp0vQQweAVj0pdY0xs5JOW2tP5//MWnvCWnvMWntsZmamF5+O\nKRcAaOOKFXrqomfavLX2ZOr13KAuiEosLAKAdq4Y6Lm2SgtjzHFr7cPx87lc0PfFSCns96cAAO/0\nYsrlE8aYF40xl3p0TACAa9DtRdGTkvb06Fiu6I8/9PdUpt0CAG11vbBokO68YWq7DwEAdiwGugGg\nIAh0ACgIAh0ACoJAB4CCINABoCAIdAAoCAIdAAqCQAeAgjDW2sF9MmMWJL1yjf/7fknne3g4PuCc\nhwPnPBy6OeebrbVX3K52oIHeDWPMKWvtse0+jkHinIcD5zwcBnHOtFwAoCAIdAAoCJ8Cfct92QuK\ncx4OnPNw6Ps5e9NDRzEZYx5K3SDlQUmLkma3ehuwExljZtO34bza7+defo97UaEbYx40xswZYx7a\n7mPpB2PM8fi/T6Te1nLORft7iG+Q8kD8fFZK9thfNMbMtnvbth1sl+LzeTD+4XVvK/TXOHUux9u8\nrVDnHH8vfz71+qq+n3v9Pb7jA71IP9TtxN8IJ+Nb/R2Jv7ELHW4dfEBRlSJJ85LmOrzNVx+x1j6q\n6Gs8kB/u7RQfu7v38HzRz9mdZ+pNV/v93NPv8R0f6CrWD3U7R9Q8p/n4ddHDzf16mr7/7LSki6nX\n+zq8zTtxVf6kJFlrH45/LS/811iS+43zyBCds3O13889/R73IdAL8UPdibX2ROpG3LOSTqnA4Zay\nd7sPYIDukbQvrkhdW6HQX+M4wOfjew27cyr0Oe8EPgT6UIh/1TydvqhSVG2qcymq0lzIT0u60OFt\nvrrgvrbpPnpRGWOmFX39Pi7pM8aYI9t8SIN2td/PPf0e9+GeokX6od7KnLX2w/HzTudclL+HI/EP\n+F5Je+N/zB6R5FbRHZHkAr/d23xzQc3+6qKiir3oX+Pjkj5urV00xsxLcpMcRT7ntL/J93PPvsd9\nqNAfUXSikt8/1B0ZY46nRpjm1P6cC/P3YK19NL5AKEU/xEpVr3OSFq21p9u9bVsOuHuPqvm1m1bU\nTy/01zgt/lovqsDnHP/Wdcz99nW138+9/h73Yg49HnuaV3RxpVALElLjThcVVSo/bq092e6ci/z3\nUHTx1+6ipHvcb2JF/xrH1wvmJe3d6vyKdM7bzYtABwBcmQ8tFwDAVSDQAaAgCHQAKAgCHQAKgkAH\ngIIg0AGgIAh0ACgIAh0ACuL/A18OoXvAqZqEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11eca52b0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "w = 7\n",
    "alpha = 0.1\n",
    "beta = 0.95\n",
    "p = 0\n",
    "EPOCH = 1000\n",
    "W = []\n",
    "for tau in range(EPOCH):\n",
    "    W.append(w)\n",
    "    p = (2*w - 4) + beta*p\n",
    "    w = w - alpha*p\n",
    "\n",
    "plt.plot(W)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "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.6.5"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": false,
   "autocomplete": false,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": false,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": true
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}