{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Illustrations for Gradient Descent\n",
"Author: Joerg Kienitz (finciraptor.de, https://github.com/Lapsilago) for the workshop Machine Learning for Option Pricing, Calibration and Hedging Workshop with Nikolai Nowaczyk ( https://github.com/niknow; https://github.com/niknow/machine-learning-examples )"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"import scipy as sp\n",
"import matplotlib.pyplot as plt\n",
"import random\n",
"from scipy import stats\n",
"from scipy.optimize import fmin"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Gradient Descent"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Gradient descent, also known as steepest descent, is an optimization technique that finds a local minimum of a function. \n",
"To find it, the function \"steps\" in the direction of the steepest negative direction of the gradient.\n",
"The algorithm of gradient descent can be outlined as follows:\n",
"\n",
" 1: Choose initial guess $x_0$ \n",
" 2: for k = 0, 1, 2, ... do \n",
" 3: $grad_k$ = -$grad f(x_k)$ \n",
" 4: choose $\\alpha_k$ to minimize $f(x_k+\\alpha_k grad_k)$ \n",
" 5: $x_{k+1} = x_k + \\alpha_k grad_k$ \n",
" 6: end for"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As a simple example, let's find a local minimum for the function $f(x) = x^4-2x^2 - 0.5 x+2 + \\sin(x)^2$"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"f = lambda x: x**4-2*x**2 - 0.5*x + 2 + np.sin(x)**2"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD8CAYAAABq6S8VAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8XOV97/HPb0ajfbF2S7Zly8Y7BmwLm80GDCFAWLI0C9kX4pI0NyRtbtukadKm997cpIEmDWkCBF4NDQFaEkIS9s0YE2wjG7DxLm/yImu1dmud5/6hges6NhpJM3Nm+b5fL700so403zmWv370nOecY845REQkcfm8DiAiIhOjIhcRSXAqchGRBKciFxFJcCpyEZEEpyIXEUlwaeFsZGYHgC5gGBhyztVEM5SIiIQvrCIPudw51xK1JCIiMi6aWhERSXAWzpmdZrYfOA444E7n3F2n2WY1sBogJydn6bx58yIcdXQHW3vpHxpmTnlezJ9bRCKvb3CYPU3dVBVlU5AV8DpOVG3atKnFOVc6nq8Nt8grnXNHzawMeAb4H865tWfavqamxtXW1o4nz4Tc/vQu7nihju3fuZrMgD/mzy8ikfXUtmP8+X9s4vdfuoRFUwu8jhNVZrZpvMcfw5pacc4dDb1vAh4Blo3nyaJtzuQ8gg72Nnd7HUVEIuBQWy8AVUXZHieJb6MWuZnlmFneW4+Bq4A3ox1sPN6aUtnd2OVxEhGJhPq2XvIz0yjITu5plYkKZ9VKOfCImb21/a+cc09GNdU4VZfkEPAbO4+pyEWSQX1bL1XFGo2PZtQid87tA86NQZYJC/h9nFWWx44GFblIMqhv62XeZC1eGE3SLT+cX5HHjoZOr2OIyAQFg47DbSeYpvnxUSVdkS+oyKe5q5/mrn6vo4jIBDR09jEwHNSBzjAkZZEDGpWLJLgDLT0AVBfneJwk/iVdkc9XkYskhQOtI0U+o0RFPpqkK/LCnHQqCjLZriIXSWgHWnrISPMxOT/T6yhxL+mKHEamVzQiF0ls+1t6mVGcg89nXkeJe0lZ5PMr8tnb3EPf4LDXUURknA609jBda8jDkpRFvqAyn+GgY0+jTtUXSUTDQUd9ay/Vmh8PS1IWuQ54iiS2ho4TDAwHdaAzTElZ5NOLsslO9+uAp0iCOtAycrGsGVp6GJakLHKfz5g3OU9FLpKg9oeWHmpqJTxJWeQwMr2yo6GTcK63LiLx5WBLD5kBH2V5GV5HSQhJW+QLKvPp6hvi8PETXkcRkTE60NqjpYdjkLRFfnblyN1E3jzS4XESERmr/S09mh8fg6Qt8nkVeQT8xhYVuUhCGQ46DrWd0IqVMUjaIs9I8zN3ch5bD6vIRRLJ0fbQ0kOdDBS2pC1ygEVTJrHlcLsOeIokEF0sa+ySusjPmVpAZ98Q9aEbuIpI/Hv78rUq8rAldZEvmjJywHOLpldEEsb+ll6y0/1aejgGSV3kc8rzSPf72KoDniIJ40BrD1VF2YRu+C5hSOoiT0/zMb8ijy2H272OIiJh2tfczayyXK9jJJSkLnKARVMLePNIJ8GgDniKxLv+oWHq23qZVaoiH4ukL/Jzpkyiu3/o7SPhIhK/DrT0EnQwq1QHOsci6Yt80dSRA56aJxeJf3ubR+4hoBH52CR9kc8uyyUjzaeVKyIJYG+Tinw8kr7I0/w+Flbm64CnSALY29zNlElZZKX7vY6SUJK+yAHOmTqJrUc6GBwOeh1FRN7B3uYerVgZh5Qo8iXTC+kbDLKzocvrKCJyBs459jZ360DnOKRGkVdNAmBz/XGPk4jImRzr7KN3YFjz4+OQEkU+ZVIWZXkZKnKROLa3aWSJsIp87FKiyM2MpdMLVeQiceztpYdlmloZq5QocoAlVYUcajtBc1e/11FE5DTqmrrJy0yjNFcXyxqr1Cny6ZonF4lnIwc6c3WxrHFImSJfWFlAwG8qcpE49VaRy9ilTJFnBvwsrCxg80EVuUi86eobpLGzX/Pj4xR2kZuZ38xeM7M/RDNQNC2pKmTL4Q4GhnRikEg82desFSsTMZYR+a3AjmgFiYUl0yfRPxRkR0On11FE5CR1usbKhIRV5GY2FXgP8PPoxomuJVWFgA54isSb3U1dpPt9TC/O9jpKQgp3RP5D4K+BM85JmNlqM6s1s9rm5uaIhIu0yklZVBZkUntARS4ST3Yf62JmaQ4Bf8octouoUfeamV0HNDnnNr3Tds65u5xzNc65mtLS0ogFjLTlM4vZsL8V53THIJF4sbuxm7mT87yOkbDC+e/vYuAGMzsAPAisMrNfRjVVFC2rLqKle4C9zbpjkEg86Oob5Ej7CeaUq8jHa9Qid8593Tk31Tk3A/gI8Lxz7uNRTxYly6uLANiwv9XjJCICsCd0oHOuinzcUm5Cqrokh9K8DDbub/M6iogwMj8OaGplAtLGsrFzbg2wJipJYsTMWF5dxIZ9bTjndDqwiMd2NXaRne5nyqQsr6MkrJQbkcPIAc9jnX3Ut/V6HUUk5e1u7GJ2WS4+nwZV45WaRf7WPPk+Ta+IeG3XsW4d6JyglCzy2WW5FOWks0Hz5CKeau3up6W7X/PjE5SSRW5mLJtRpJUrIh7b3TiyYkUj8olJySIHWD6ziMPHT3Ck/YTXUURS1u5GrViJhJQt8gtnFQPwcl2Lx0lEUteuxi4KsgKU5emuQBORskU+tzyPktwM1u1RkYt4ZfexLuaW52kZ8ASlbJGbGZecVczLdS0Eg7ruikisOefY1djFnMm6dO1EpWyRA1wyu5TWngF2HNP1yUViraGjj66+IZ2aHwGpXeRnlQBoekXEA9uPjgygFlTme5wk8aV0kU8uyGR2WS7rdMBTJOa2N3RiBnMnq8gnKqWLHOCS2SVs3N9G3+Cw11FEUsr2o53MKM4hN2NMl3yS00j5Il8xu4T+oaDuGiQSY9sbOllQodF4JKR8kS+vLibgN17aE5+3pxNJRp19g9S39Wp+PEJSvshzMtI4f0YRL+xq8jqKSMrY2TByRqdG5JGR8kUOsGpeGbsbuzmky9qKxMT2ox2AVqxEioqckSIHNCoXiZHtDZ0U56Tr1PwIUZEDM0tzqS7J4fmdKnKRWNjR0MWCynydmh8hKvKQy+eW8ce9rfQODHkdRSSpDQ4H2dXYpfnxCFKRh1wxv4yBoSB/rNM1ykWiaV9zDwNDQc2PR5CKPOT8GUXkZqTxnKZXRKJqW+hA53yNyCNGRR6SnuZjxewSnt/ZqKshikTR1iMdZAX8zCzJ8TpK0lCRn+SqheU0dvbzxuF2r6OIJK2thztYWJlPml/1EynakydZNa+cgN948s1jXkcRSUpDw0G2He1k0dQCr6MkFRX5SQqyAlw0q4Qn3jyGc5peEYm0uuZuTgwOc46KPKJU5Ke45uzJ1Lf1sr1BN5sQibQth0cOdJ4zdZLHSZKLivwU71pQjs/Q9IpIFGw93EFuRhrVxTrQGUkq8lMU52awvLqYJ1TkIhG35UgHZ0/Jx+fTGZ2RpCI/jWsWTaauqZu6pi6vo4gkjYGhIDuOdmpaJQpU5Kdx9dmT8Rn87vWjXkcRSRq7G7sYGA7qQGcUqMhPoywvk4vPKuG3rx/V6hWRCHn7QOcUjcgjTUV+BjeeN4X6tl421+vkIJFI2HqknYKsANOKsryOknRU5Gfw7oXlZKT5+O1rR7yOIpIUXqtv55ypBbp0bRSoyM8gLzPAlQvKeWxrA4PDQa/jiCS07v4hdjd2sbiq0OsoSWnUIjezTDPbaGZvmNk2M/vHWASLB+87bwptPQO6MbPIBL1xqJ2ggyVVmh+PhnBG5P3AKufcucB5wNVmdkF0Y8WHlXNKKcwO8OtNml4RmYjNB48DsHiaRuTRMGqRuxHdoQ8DobeUWMqRnubjvYun8PT2Y7R293sdRyRhba4/zllluRRkB7yOkpTCmiM3M7+ZvQ40Ac845zacZpvVZlZrZrXNzckzFXHTsioGhx2/2axRuch4OOd47VC7plWiKKwid84NO+fOA6YCy8zs7NNsc5dzrsY5V1NaWhrpnJ6ZU57HkqpJPPBqvdaUi4zDvpYe2nsHWTpd0yrRMqZVK865dmANcHVU0sSpjyyrYl9zD68eOO51FJGE89b8+BKtWImacFatlJrZpNDjLOBKYGe0g8WT686pIC8jjQc31nsdRSThbK5vJz8zjVmluV5HSVrhjMgrgBfMbAvwKiNz5H+Ibqz4kp2exg3nVfLY1gbaega8jiOSUF6rP855VYW64mEUhbNqZYtzbrFz7hzn3NnOue/EIli8+dRFM+gfCvKARuUiYevqG2RXY5cOdEaZzuwM05zyPFbMLuG+Vw4wMKQzPUXCUXvwOM5BzfQir6MkNRX5GHz24moaO/t54s0Gr6OIJISN+9tI8xlLpmtEHk0q8jG4dE4pM0tzuGfdfi1FFAnDhn2tnDO1gOz0NK+jJDUV+Rj4fMZnLq5my+EOag9qKaLIOzkxMMyWwx0sqy72OkrSU5GP0QeWTKE4J507nq/zOopIXNtcf5yhoGP5TM2PR5uKfIyy09P43IpqXtzdzBuHdNMJkTPZsL8Nn0GNzuiMOhX5OHzigukUZAX4sUblIme0YV8rCysLyMvUhbKiTUU+DnmZAT5z8Qye3dHI9qOdXscRiTv9Q8O8dqidZdWaVokFFfk4feaianIz0vjRc7u9jiISd9441MHAUJDlKvKYUJGPU0F2gM+vmMlT2xrZdLDN6zgicWX9vlbM4PwZKvJYUJFPwM0rqinNy+C7j+/UunKRk6zb08LCynwKc9K9jpISVOQTkJORxleunE3tweM8vb3R6zgicaGnf4jN9ce55KzkuS9BvFORT9CHa6YxszSH7z2xU9dgEQE27G9lKOhYMbvE6ygpQ0U+QWl+H998z3z2tfTw83X7vI4j4rmX9rSQkebTHYFiSEUeAavmlXPVgnL+9bk9HGrr9TqOiKfW7WlhWXURmQG/11FShoo8Qr59w0IM4x9/v83rKCKeOdbRx56mbi45S9MqsaQij5Apk7L4ypWzeXZHE49v1WVuJTWtq2sB4BLNj8eUijyCPntJNYumFPB3j2ylqbPP6zgiMbduTzPFOenMn5zvdZSUoiKPoIDfx798+Fx6B4b5m19v0dpySSnDQcdLe1q4ZHaJ7s8ZYyryCDurLI+vXzOPF3Y188sNur+npI43DrfT2jPAqnllXkdJOSryKPjkhTNYOaeUf/r9dl7XpW4lRTy/owm/z7h0jk4EijUVeRT4fMaPPnwepXkZfOGXm2jp7vc6kkjUPbeziaXTC5mUrdPyY01FHiWFOenc+YmltPUM8Bf3b9ZZn5LUjrafYEdDJ1doWsUTKvIoOntKAd/7wDls2N/G1/7rDYJBHfyU5PT8ziYArpivIveCbm0dZe9dPIWGjj6+9+ROinPT+dZ1CzDTEX1JLs/vbKKqKJtZpbleR0lJKvIYuOXSmTR39XPvy/vJTvfztavmqswlaZwYGObluhZuWlaln2uPqMhjwMz45nvmc2JwiJ+8sJfegWGNzCVpvLi7if6hIFfOL/c6SspSkceIz2f8n/ctIiuQxr0v76fjxCDfff8iMtJ0YSFJbI9vPUZhdoALZupuQF5RkceQmfH3181nUnaA25/Zzf6WHu78+FLK8jO9jiYyLn2Dwzy3o5Hrz60kza+1E17Rno8xM+PLV8zmpx9bws6GLq6/Yx0v7Wn2OpbIuKzd3UzPwDDXLqrwOkpKU5F75JpFFfzmixeRm5HGJ+7ZyLcefZOe/iGvY4mMyRNvHmNSdoALZxV7HSWlqcg9NL8in8e+vILPXlzNfa8cZNVta3jktcNaby4JoX9omGe3N3LVgnICmlbxlPa+xzIDfr51/QJ+/YWLmJyfyVcfeoPr71jHY1saGFahSxx7ua6Frv4hrtG0iudU5HFi6fRCHvnixdz+oZHL4P7FrzZzxW1r+MkLdTR0nPA6nsifePT1oxRkBbh4lm4i4TWtWokjPp/x/iVTufG8KTy17Rj//vIB/vmpXfzg6V2cP6OIy+aWsnJ2KfMr8vFH+XrPzjk6+4Zo6xmgrWeA4z0DtPX+//edJ4boGxzmxMAwfUMj7weHg/jMRt584PcZOelp5GUGyM9KIz8zQFl+BpUFWVRMyqRyUhb5mYGovg6Jjq6+QZ7adowPLp1GeprGg14btcjNbBpwHzAZCAJ3Oed+FO1gqczvM65dVMG1iyo42NrDrzcf4ZntjXz/yV18/8ldZKf7WViZz4KKfKqKc5gyKYvKSZkUZAXIywyQm5GG32c453BA0Dn6BoJ09Q/S3T9Ed98QXX1DtHT309ozQEtX6H13Py3dA7R299PWM8DQGaZ20v0+8rMCZKf7yQz4yAr4yQj4yU5PI+hc6A0GB4O0dvfS1TdEZ9/Ic596r43SvAzmlucxd3Iec8vzOHfaJGaX5erGBHHuia3H6BsM8v4lU7yOIoCNdhcbM6sAKpxzm80sD9gEvNc5t/1MX1NTU+Nqa2sjm1Ro6uxjXV0LWw53sPVIBzsbOukZGJ7w980K+CnJS6c4J4OS3HRKcjMozk2nKCeDopwAhdnpFOWkv/0+O90/rrNSh4aDNHf3c7T9BEfb+zjSfoI9jd3sbuxiT1MXfYMjV4jMz0xj6fRCzq8u4rI5ZcyvyNNZsHHmw3e+QnNXP8/91aX6u4kQM9vknKsZ19eO9XZkZvYocIdz7pkzbaMijw3nHO29gxw+foKjHSdCI+2Rke9wEMzAZyNr1zMDfvIy0sjNTCM39L4kJ4OSvHSy072fYRsOOg609vBafTu1B9qoPXicuqZuAMrzM7h8bhlXzi9n5ZxS/SrvsUNtvaz4/gt87ao5fGnVbK/jJI2JFPmY/gWb2QxgMbDhNJ9bDawGqKqqGk8WGSMzozAnncKcdBZNLfA6zoT4fcas0lxmlebyZ0unAtDU1ceLu5pZs6uZx7Y08OCrhyjICnDtoslcf24ly6uLo36sQP7UI68dAUau7CnxIewRuZnlAi8C/9s595t32lYjcom0weEg6+pa+N3rR3lq2zF6B4apLMjko8ur+ND50yjL02UOYiEYdFx+2xoqCjJ5cPWFXsdJKlEfkZtZAPg1cP9oJS4SDQG/j8vnlnH53DJODAzz7I5GHnr1ED94ejc/fHYP7z57Mp+6cAbnzyjUnG0Uratr4WBrL3/5rjleR5GThLNqxYB7gB3OudujH0nknWWl+7n+3EquP7eSfc3d3L+hnoc3HeaxLQ0snV7IFy6dxap5ZVr5EgW/XH+Q4px0rj57stdR5CThHDW6GPgEsMrMXg+9XRvlXCJhmVmay99ft4D1X7+C79y4kGMdfdx8Xy3X/OglHn39iM6OjaCGjhM8u6ORD9ZM0+WX48yoI3Ln3DpAQxuJa1npfj554QxuWlbFH7Yc5adr9nLrg6/zby/s5WvvnsuV88s05TJBD2w8hAM+ukyLGeKN1nFJUgn4fbxv8VSevHUlP75pMQPDQT5/Xy0f+OkfeWVvq9fxElb/0DAPbKxn5exSqoqzvY4jp1CRS1Ly+Yzrz63k6a+u5LvvX8TR9j5uuns9n7+vloOtPV7HSziPvn6U5q5+bl5R7XUUOQ0VuSS1gN/HTcuqWPM/L+Ovr57Ly3UtvOv2tXzvyZ106/rvYXHOcffafcybnMclZ+kCWfFIRS4pITPg54uXncULX7uM686t4Kdr9rLqB2v4zebDjPXs5lSzZncze5q6Wb1ypo4zxCkVuaSU8vxMbv/QeTzyxYuomJTFX/7nG3zs5xvY36LpljO588W9TM7P5PpzK72OImegIpeUtLiqkEe+cBH/671ns/VIB+/+4Vp+/NweBoaCXkeLK+v3tbJ+Xxs3r6jWXYDimP5mJGX5fMbHL5jOc395Ke+aX85tz+zmPf/6Eq8eaPM6WlxwznH7M7spy8vg4xdM9zqOvAMVuaS8svxMfvKxJdz76Rp6B4b54M9e0c2wgVf2trJxfxt/cflZZAZ0AlA8U5GLhKyaV87TX13Jpy+awX+sP8jVP1qbsmvPnXP84OldVBRk8uHzp3kdR0ahIhc5SU5GGv9ww0Ie/PwF+My46e71fPvRN+kdSK3R+WNbG9hc386tV8zWaDwBqMhFTmP5zGKeuHUFn7l4BvetP8jVP3yJ9ftSY3TeNzjMdx/fyfyKfD5Yo9F4IlCRi5xBdnoa375+IQ+tvhAz+MhdqTE6v2fdfo60n+Dvr5uvG3ckCBW5yCiWVRfx5K0r/9voPFnnzg+19XLH83VctaCci2bpLM5EoSIXCUNWuv/t0bnP4Ka71yfdyhbnHN94ZCs+g2/fsNDrODIGKnKRMVhWXcQTt67ksxdXv72y5Y97W7yOFREPbzrMS3ta+Jtr5jFlUpbXcWQMVOQiY5SV7udb1y/gP//8QtJ8Pj569wa++dutCT06P9TWyz/9YTs10wv5+HKd/JNoVOQi43T+jCIe//IKbr6kmvs31PPuH67l5brEG50PDgf58oOvEXRw24fO1S3yEpCKXGQCstL9fPO6BTx8y4Wk+3187Ocb+MYjW+nqG/Q6Wthue3o3r9W3838/sIjpxTlex5FxUJGLRMDS6UU8fusKVq+cyQMb67n6hy/x0p5mr2ON6revHeFnL+7lo8uruO4cXd0wUanIRSIkM+DnG9fO5+FbLiIj4OMT92zkyw+8RlNnn9fRTqv2QBt//fAWllcX8Q/Xa5VKIlORi0TY0umFPP7lFdx6xWye3HaMVbe9yD3r9jM0HD+XyN1+tJOb76tlSmEWd35iKelpqoJEpr89kSjIDPj56rvm8PRXVrJ0eiH/9IftXPfjdWzc7/0lcrcf7eSjP19PdsDPLz6zjEnZ6V5HkglSkYtE0YySHP79M+fzs48voePEIB+68xU+f18tdU1dnuT5Y10LN909UuIPrr6QquJsT3JIZKnIRaLMzLj67Aqe+6tL+dpVc3hlbytX/ctavv6bLRxpPxGTDM457t9wkE/eu5GyvAwe+nOVeDKxaNx4tqamxtXW1kb8+4okg9bufn78fB33bziIc/DexVO45dJZnFWWG7Xn+8YjW3lqWyOXzinlxx9dTH5mICrPJeNnZpucczXj+loVuYg3jrSf4O61+3jw1Xr6h4JcMa+Mj5xfxWVzS0mLwP0xB4aC/GrDQX703B56+of52rvn8LlLZuqKhnFKRS6SwFq7+/nFHw/wq42HaOnupywvg/ctnsK7FpSzuKpwzMXb3jvAf9Ye4r5XDnL4+AkumlXMt65fwLzJ+VF6BRIJKnKRJDA4HOT5nU089Ooh1u5uZijoKMpJZ3l1EYurJrGgooBpRVlMLsgk3e/DzBgYCtLWM0BdUzdvHu1g7e5mNu5vYyjoWDajiFsum8nlc8sw0yg83qnIRZJMZ98ga3c389yOJjYdPE59W+9/+7wZBPw+Bob++9r0OeW5rJpXzo3nVTK/QiPwRDKRIk+LdBgRmbj8zADXnVP59mnzLd397Gns5vDxXho7++gfCjIwHCQ3PY3CnHRmluQwryKfohytCU9FKnKRBFCSm0FJbgZQ7HUUiUNaRy4ikuBU5CIiCU5FLiKS4EYtcjO718yazOzNWAQSEZGxCWdE/u/A1VHOISIi4zRqkTvn1gLeX3tTREROK2Jz5Ga22sxqzay2uTn+b3ElIpIsIlbkzrm7nHM1zrma0tLSSH1bEREZhVatiIgkOBW5iEiCC2f54QPAK8BcMztsZp+LfiwREQnXqNdacc7dFIsgIiIyPppaERFJcCpyEZEEpyIXEUlwKnIRkQSnIhcRSXAqchGRBKciFxFJcCpyEZEEpyIXEUlwKnIRkQSnIhcRSXAqchGRBKciFxFJcCpyEZEEpyIXEUlwKnIRkQSnIhcRSXAqchGRBKciFxFJcCpyEZEEpyIXEUlwKnIRkQSnIhcRSXAqchGRBKciFxFJcCpyEZEEpyIXEUlwKnIRkQSnIhcRSXAqchGRBKciFxFJcCpyEZEEpyIXEUlwKnIRkQSnIhcRSXBhFbmZXW1mu8yszsz+NtqhREQkfKMWuZn5gZ8A1wALgJvMbEG0g4mISHjCGZEvA+qcc/uccwPAg8CN0Y0lIiLhSgtjmynAoZM+PgwsP3UjM1sNrA592G9mb048XkSVAC1ehziFMoUnHjNBfOZSpvDEY6a54/3CcIrcTvNn7k/+wLm7gLsAzKzWOVcz3lDRoEzhUabwxWMuZQpPvGYa79eGM7VyGJh20sdTgaPjfUIREYmscIr8VWC2mVWbWTrwEeB30Y0lIiLhGnVqxTk3ZGZfAp4C/MC9zrlto3zZXZEIF2HKFB5lCl885lKm8CRVJnPuT6a7RUQkgejMThGRBKciFxFJcBEpcjP7ZzPbaWZbzOwRM5t0hu1idqq/mX3QzLaZWdDMzrjMyMwOmNlWM3t9Ist/IpwplvupyMyeMbM9ofeFZ9huOLSPXjezqBzsHu11m1mGmT0U+vwGM5sRjRxjzPRpM2s+ad/cHINM95pZ05nO1bAR/xrKvMXMlsRBpsvMrOOk/fStGGSaZmYvmNmO0L+7W0+zTUz3VZiZxr6vnHMTfgOuAtJCj78HfO802/iBvcBMIB14A1gQiec/Q6b5jCywXwPUvMN2B4CSaOUYayYP9tP3gb8NPf7b0/3dhT7XHeV9M+rrBr4I/Cz0+CPAQ3GQ6dPAHbH4+TnpOVcCS4A3z/D5a4EnGDkH5AJgQxxkugz4Q4z3UwWwJPQ4D9h9mr+/mO6rMDONeV9FZETunHvaOTcU+nA9I2vNTxXTU/2dczucc7ui9f3HI8xMsb4kwo3AL0KPfwG8N4rP9U7Ced0nZ30YuMLMTnfCWiwzxZxzbi3Q9g6b3Ajc50asByaZWYXHmWLOOdfgnNscetwF7GDkTPWTxXRfhZlpzKIxR/5ZRv6HO9XpTvWf8AuIAAc8bWabQpcZ8Fqs91O5c64BRn7IgLIzbJdpZrVmtt7MolH24bzut7cJDRw6gOIoZBlLJoAPhH4tf9jMpp3m87EWr//WLjSzN8zsCTNbGMsnDk3DLQY2nPIpz/bVO2SCMe6rcE7Rf+tJnwUmn+ZTf+ecezS0zd+ExIz/AAACe0lEQVQBQ8D9p/sWp/mzCa19DCdTGC52zh01szLgGTPbGRpdeJUppvtpDN+mKrSfZgLPm9lW59zeieQ6RTivO+L7ZhThPN/vgQecc/1mdgsjvzGsimKmcMR6P4VjMzDdOddtZtcCvwVmx+KJzSwX+DXwFedc56mfPs2XRH1fjZJpzPsq7CJ3zl05SrBPAdcBV7jQRM8pIn6q/2iZwvweR0Pvm8zsEUZ+nR53kUcgU0z3k5k1mlmFc64h9Ctl0xm+x1v7aZ+ZrWFkJBHJIg/ndb+1zWEzSwMKiO6v86Nmcs61nvTh3YwcI/Ja3F1W4+Sycs49bmb/ZmYlzrmoXrjKzAKMFOb9zrnfnGaTmO+r0TKNZ19FatXK1cDfADc453rPsFncnepvZjlmlvfWY0YO2np91cZY76ffAZ8KPf4U8Ce/NZhZoZllhB6XABcD2yOcI5zXfXLWPwOeP8OgIWaZTplPvYGROU+v/Q74ZGhFxgVAx1vTZ14xs8lvHc8ws2WMdE/rO3/VhJ/TgHuAHc6528+wWUz3VTiZxrWvInQkto6ReabXQ29vrSyoBB4/abtrGTlKu5eRqYZoHh1+HyP/2/YDjcBTp2ZiZDXCG6G3bfGQyYP9VAw8B+wJvS8K/XkN8PPQ44uAraH9tBX4XJSy/MnrBr7DyAABIBP4r9DP20ZgZjT3TZiZvhv62XkDeAGYF4NMDwANwGDo5+lzwC3ALaHPGyM3g9kb+vs646qtGGb60kn7aT1wUQwyXcLINMmWk7rpWi/3VZiZxryvdIq+iEiC05mdIiIJTkUuIpLgVOQiIglORS4ikuBU5CIiCU5FLiKS4FTkIiIJ7v8B7eo77GsZijwAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"x = np.linspace(-2,2.5,1000)\n",
"plt.plot(x,f(x))\n",
"plt.xlim([-2,2.5])\n",
"plt.ylim([0,5])\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can see from plot above that our local minimum is gonna be near around 1.0 but let's pretend that we don't know that, so we set our starting point (arbitrarily, in this case) at $x_0 = 2$"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Local minimum calculated (x value): 0.9374670913810693\n",
"Number of steps to convergence: 49\n"
]
}
],
"source": [
"# hands on simple implementation of the gradient descent algorithm\n",
"w_0 = 2 # The algorithm starts at x=2\n",
"alpha = 0.1 # learning rate / step size\n",
"precision = 0.0001 # break criteria\n",
"n_steps = 100 # break criteria\n",
"\n",
"x_list, y_list = [w_0], [f(w_0)]\n",
"\n",
"# returns the value of the derivative of our function\n",
"def f_prime(x):\n",
" return 4*x**3-4*x - 0.5 + 2 * np.sin(x) * np.cos(x)\n",
"\n",
"n = 0 # number of iterations\n",
"x_old = w_0 # set the starting value\n",
"\n",
"# first iteration to find the new value\n",
"ngrad = -f_prime(x_old) # gradient calculation\n",
"x_new = x_old + alpha * ngrad # the newly suggested value\n",
"\n",
"while (abs(x_new - x_old) > precision) and (n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.figure(figsize=[15,3])\n",
"plt.subplot(1,3,1)\n",
"plt.scatter(x_list,y_list,c=\"r\")\n",
"plt.plot(x_list,y_list,c=\"r\")\n",
"plt.plot(x,f(x), c=\"b\")\n",
"plt.xlim([-2,2.5])\n",
"plt.ylim([0,5])\n",
"plt.title(\"Gradient descent\")\n",
"plt.subplot(1,3,2)\n",
"plt.scatter(x_list,y_list,c=\"r\")\n",
"plt.plot(x_list,y_list,c=\"r\")\n",
"plt.plot(x,f(x), c=\"b\")\n",
"plt.xlim([0.9,1.2])\n",
"plt.ylim([0,5])\n",
"plt.title(\"Zoom around min\")\n",
"plt.subplot(1,3,3)\n",
"plt.scatter(x_list,y_list,c=\"r\")\n",
"plt.plot(x_list,y_list,c=\"r\")\n",
"plt.plot(x,f(x), c=\"b\")\n",
"plt.xlim([-0.9,-0.6])\n",
"plt.ylim([0,5])\n",
"plt.title(\"Zoom around pseudo min\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We fixed the size of the step size, aka learning rate to a constant value. This might not be the best idea and may effects the search for the minimum. An idea for improving this might be to take larger steps when far away from the minimum and smaller ones close to the minimum. \n",
"\n",
"Warning: The gradient descent is not guaranteed to converge! If it converges it may not find the global minimum but got stuck in a local one!\n",
"\n",
"Solution: Vary the step size using adaptive methods or calculating an optimal step size!"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Local minimum occurs at 0.937221769739723\n",
"Number of steps: 2\n"
]
}
],
"source": [
"# we setup this function to pass into the fmin algorithm\n",
"def stepsize(n,x,s):\n",
" x = x + n*s\n",
" return f(x)\n",
"\n",
"w_0 = 2 # The algorithm starts at x=2\n",
"precision = 0.00001\n",
"n_steps = 50\n",
"\n",
"x_list, y_list = [w_0], [f(w_0)]\n",
"\n",
"# returns the value of the derivative of our function\n",
"def f_prime(x):\n",
" return 4*x**3-4*x - 0.5 + 2 * np.sin(x) * np.cos(x)\n",
"\n",
"n = 0 # number of iterations\n",
"x_old = w_0 # set the starting value\n",
"\n",
"# first iteration to find the new value\n",
"s_k = -f_prime(x_old) # gradient calculation\n",
"alpha_calc = fmin(stepsize,precision,(x_old,s_k), full_output = False, disp = False)\n",
"x_new = x_old + alpha_calc * s_k # the newly suggested value\n",
"\n",
"while (abs(x_new - x_old) > precision) and (n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.figure(figsize=[15,3])\n",
"plt.subplot(1,3,1)\n",
"plt.scatter(x_list,y_list,c=\"r\")\n",
"plt.plot(x_list,y_list,c=\"r\")\n",
"plt.plot(x,f(x), c=\"b\")\n",
"plt.xlim([-2,2.5])\n",
"plt.ylim([0,5])\n",
"plt.title(\"Gradient descent\")\n",
"plt.subplot(1,3,2)\n",
"plt.scatter(x_list,y_list,c=\"r\")\n",
"plt.plot(x_list,y_list,c=\"r\")\n",
"plt.plot(x,f(x), c=\"b\")\n",
"plt.xlim([0.9,1.2])\n",
"plt.ylim([0,5])\n",
"plt.title(\"Zoom around min\")\n",
"plt.subplot(1,3,3)\n",
"plt.scatter(x_list,y_list,c=\"r\")\n",
"plt.plot(x_list,y_list,c=\"r\")\n",
"plt.plot(x,f(x), c=\"b\")\n",
"plt.xlim([1.3333,1.3335])\n",
"plt.ylim([0,3])\n",
"plt.title(\"zoomed in more\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Another approach to update the step size is choosing a decrease constant $dec$ that shrinks the step size over time:\n",
"$\\alpha(t+1) = \\alpha(t) / (1+t \\cdot dec)$."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Local minimum occurs at: 0.9874558760109118\n",
"Number of steps: 5\n"
]
}
],
"source": [
"x_old = 0\n",
"x_new = 2 # The algorithm starts at x=2\n",
"n_k = 0.17 # step size\n",
"precision = 0.1\n",
"t, d = 0, 1\n",
"\n",
"x_list, y_list = [x_new], [f(x_new)]\n",
"\n",
"# returns the value of the derivative of our function\n",
"def f_prime(x):\n",
" return 4*x**3-4*x - 0.5 + 2 * np.sin(x) * np.cos(x)\n",
" \n",
"while abs(x_new - x_old) > precision:\n",
" x_old = x_new\n",
" ngrad = -f_prime(x_old)\n",
" x_new = x_old + n_k * ngrad\n",
" x_list.append(x_new)\n",
" y_list.append(f(x_new))\n",
" n_k = n_k / (1 + t * d)\n",
" t += 1\n",
"\n",
"print(\"Local minimum occurs at:\", x_new)\n",
"print(\"Number of steps:\", len(x_list))"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA3gAAADfCAYAAACkn3GFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XeYVOX5//H3vbvAUgUEC13EigoqlmhisEU0scQSCyoqSrAkJmrUxHwViUSNJcYeW7Abf0YTY6wxEoMoCkpTULGCSAfpZZf798dzRmaX2d1ZmJkz5fO6rnNNOWfOuc/MzrPnnqeZuyMiIiIiIiKFryzuAERERERERCQzlOCJiIiIiIgUCSV4IiIiIiIiRUIJnoiIiIiISJFQgiciIiIiIlIklOCJiIiIiIgUCSV4RczMPjezQ6L7vzGz+2KKo7+ZzYzj2CIi2WBmZ5jZ6Azta6CZvZyJfYlI/jGzUWZ2dh3rYrs+k+KlBC8mZnaSmY01s+VmNje6f56ZWTaO5+6/d/eUhUtjmFkPM3Mzq8hEXHEzs5Fmdk3ccYjEJUoulqVY3MyujDu+UuDuj7r7D+KOQ0RyL1PXZyLJlODFwMwuBv4E3ABsBWwJDAX2B5rW8ZrynAUoIiUjSi5aJS/AL4A5wL0xh5e2YvnRSUQkH+i6s7ApwcsxM9sMGA6c5+5PuftSD95z94HuvjrabqSZ3WVmz5vZcuBAM/uhmb1nZkvMbIaZDau179PM7AszW2BmV9RaN8zMHkl6vK+ZjTGzxWY20cz6J60bZWa/M7M3zGypmb1sZh2i1a9Ht4ujX/m/k+Icm0fxLzKzD4C9aq3vZGZ/M7N5ZvaZmf08ad3eZjYuOsc5ZnZz0rrvJsU8w8zOiJ5vZmY3mtmX0WvuNrPm0br+ZjbTzC6Oakq/NrMzo3VDgIHApdG5/DONj1CkqJnZ7sAfgZPc/evouU5m9qyZLTSz6WZ2TtL2zczsFjObFS23mFmzaF3i+3dp0vfvGDM7wsw+ivb3m3piqbPMS2pNMNjMvgT+Ez1/lJm9H5UTo8xsp6TXuJn1Snr8bQ1+fWVFtH7z6D1YYmZvA9vWE3citjOjuBeZ2VAz28vMJkWx3Z60fY3mntFrh5rZx9Fr7zDLTusOkUJjZidazdYGq81sVLRuMzN7KLq++MLMfmtmZdG6sujxF9F3/CEL12SN/s5GrznLzKZG275kZt2T1h1qZtPM7JvodXV+fy3p+iwpjkHRNc18q3U9V+u1I83sTjN7IXov3jCzraJyeFEUw+5J2+8UlYuLo3LyqFr7qn3dWef1leQ5d9eSwwUYAFQBFQ1sNxL4hlCrVwZUAv2BXaPHuxF+YT8m2n5nYBlwANAMuDk6ziHR+mHAI9H9zsAC4IhoX4dGjztG60cBnwDbA82jx9dF63oAXl/8wHXA/4D2QFdgCjAzWlcGjAeuJNRW9gQ+BQ6L1r8JnBbdbwXsG93vBiwFTgaaAJsDfaN1twDPRsdrDfwTuDZa1z96H4ZHrzsCWAG0S3qfr4n770KLlnxYgLbRd/+yWs//F7gzKof6AvOAg6N1w4G3gC2AjsAY4HfRusT378ro+3dO9NrHou9qb2AV0LOOeOor8xJl0UNAy6is2h5YHpVpTYBLgelA0+g1DvRK2v+33/80yoongCejY+0CfAWMriPuRGx3R+/ZD6Lz/Hv0PnUG5gLfj7Y/I3lf0Wufiz6PbtF7NiDuvw8tWvJtAdoAU4GfRo8fAv4RlS89gI+AwdG6s6LyoCfh+uJp4OFoXWO/s8dE+9oJqAB+C4yJ1nUAlgDHR2XJL6Oy5ew6zmEY66/PEnHcG5VpfYDVwE51vHYkMB/YM4r7P8BnwOlAOXAN8Fq0bZMo5t8Qrr8OIlxX7ZC0r9rXnXVeX2nJ7yX2AEptAU4FZtd6bgywGFgJHBA9NxJ4qIF93QL8Mbp/JfBE0rqWwBpSJ3iXJQq1pO1fAgZF90cBv01adx7wYnQ/UfjUl+B9StLFCDCE9QnePsCXtbb/NfCX6P7rwNVAhxTbPJPiWEa4oNs26bnvAJ9F9/tH72tF0vq5rE8cR6IET4uWxHfpH9FiSc93BaqB1knPXQuMjO5/AhyRtO4w4PPofuL7Vx49bh2VH/skbT+eKGlLI8bkMi9RFvVMWv9/wJNJj8sIiVj/6HFDCV7KsoJwobQW2DFp3e9pOMHrnPTcAuDEpMd/A34R3T+DDRO87yY9fhK4PO6/ES1a8mmJvt/PAXdFj8sJydDOSdv8FBgV3X+V0HoqsW6H6HtdsRHf2ReIEsekWFYA3QnJ1VtJ6wyYSeMSvC5J698mtKhI9dqRwL1Jj38GTE16vCuwOLr/PWA2UJa0/nFgWNK+HqoVd53XV1rye1GfhdxbAHQwswp3rwJw9/0ALIw0mdxsdkbyC81sH0Lt2C6EX1+aAf8vWt0peXt3X25mC+qIoTtwgpkdmfRcE+C1pMezk+6vIPzala4asQBf1Dp2JzNbnPRcOaHGD2Aw4Rf0aWb2GXC1uz9HuMj8JMWxOgItgPFJLZgs2mfCgsR7vZHnI1IKLiOULXt69J880glY6O5Lk577AuiXtP6LWus6JT1e4O7V0f2V0e2cpPUrqeP72ECZl5Bc1tSIxd3XmdkMwq/v6airrOhIuAisq1yrS+3zTOu8I5tSBouUghGEH40S3Tw6EMqJ2uVR4vufqqyqIIyDkJDud7Y78CczuylpvUXHqn095lE51BiN+f6nG3MnYIa7r0tan/z+QM0yLp3rK8lT6oOXe28SfmE6Oo1tvdbjxwhV5V3dfTNCU4LEt+5rQhIEgJm1IDRjTGUGoQavbdLS0t2v24iYUqkRC6GJUfKxP6t17NbufgSAu3/s7icTmkRcDzxlZi2j16Xq8zKfUID1TtrfZh4GikhHOucjUtQs9MG9Ajje3RfXWj0LaG9mrZOe60aoGUus715r3awMhVZfmZeQ/B2uEUvUb61rUqwrCBcsCVulGcc8QhOruso1EckhMzuJ0GXjeHdfGz09n1AjV7s8qq+sqqJmQpSuGYRmocnXMs3dfQwbXo8ZNcuOuMwCuib6JEaS3x+oWZ5u6vWVxEgJXo5FF09XA3ea2fFm1irq+NuX0KyyPq0Jv6SvMrO9gVOS1j0F/MjCQCRNCbVgdX2+jwBHmtlhZlZuZpUWBhjoksYpzAPWEdqw1+VJ4Ndm1i7a58+S1r0NLDGzyywMxlJuZruY2V4AZnaqmXWMfmFKXGhWA48Ch5jZT8yswsKAB32j7e4F/mhmW0T76Gxmh6VxLhAK9vrORaSomdnWhP5lv3D392qvd/cZhGbk10ZlxW6EmvZHo00eB35rZh0tDMZ0JaGMyYT6yrxUngR+aGYHm1kT4GLCD2pjovUTgFOicmcA8P10gohqIJ8GhplZCzPbGRi0EecjIpsoGjTkNkLT7nmJ56Pv6ZPACDNrHQ16chHry6PHgV+a2TZm1orQzPqvtWrt03U34TqndxTTZmZ2QrTuX0BvMzvWwui+Pyf9H5OyaSyhyeWlZtYk+mHvSEL5v4EMXF9JjJTgxcDd/0AodC4l9PGYA/yZ0ERqTD0vPQ8YbmZLCRdRTybt833gfMIv3l8DiwhtvlMdfwahBvE3hIRtBvAr0vh7cPcVhGYRb0SjMO2bYrOrCdX+nwEvAw8nvb6aUKD0jdbPB+4DNos2GQC8b2bLCFNJnOTuq9z9S8KgBxcDCwkXan2i11xG6Dj8lpktAf5NaFufjvuBnaNz+XuarxEpJucQmij9yTacC+/uaJuTCX1DZgHPAFe5+yvRumuAccAkYDLwbvRcJtRZ5qXi7h8S+jnfRihbjgSOdPc10SYXRs8tJoyg25jv/AWEpk6zCX1V/tKI14pI5hwNtANGJ5VVL0TrfkZIYj4FRhOuiR6I1j1AuB55nXD9sYqaP0Cnzd2fIbQyeiK67pgCHB6tmw+cQGhevgDYDnhjY46TSVE5eBQhzvmEgbNOd/dp9bxsU66vJEZWs6uFiIiIiIiIFCrV4ImIiIiIiBSJtEbRNLPPCXNlVANV7t6v/leIiGSfyiYRySUzqyQ08WtGuIZ6yt2vqrVNM8J8bHuyfqj9z3McqoiUsMZMk3Bg1K5YRCSfqGwSkVxZDRzk7suigXxGm9kL7v5W0jaDgUXu3isa7fF64MQ4ghWR0qQmmiIiIiJp8GBZ9LBJtNQezOBo4MHo/lPAwZY0kZiISLalm+A58LKZjTezIdkMSESkEVQ2iUhORdNsTCCMgv2Ku4+ttUlnogmjoyH4v6HueWlFRDIu3Saa+7v7rGgejFfMbJq7v568QXRxNQSgZcuWe+64444ZDhUmTYI2baBHj4zvurh99RXMng09e0K7dnFHI5Evv4SFC6Fv37gjSc/48ePnu3vHuOOoJS/KJhGJT67Lpmi6n75m1hZ4xsx2cfcpSZukqq3bYMhylU0ixS+ua6dGT5NgZsOAZe5+Y13b9OvXz8eNG7eJoW3o0ENh8WJ4552M77q4rV0L++8P06fDxInQtWvcEQkwYADMnw9Z+KpkhZmNz+dBTOIsm0QkPnGWTWZ2FbA8udwxs5eAYe7+ZjTR9Wygo9dzwaWySaQ4xVU+NdhE08xamlnrxH3gB4QJHXOud2/44ANYty6OoxewJk3gscdgzRo47TSoro47IgE++wy22SbuKApXPpVNIlIazKxjVHOHmTUHDgFqTxT9LDAoun888J/6kjsRkUxLpw/eloRRoiYCbwP/cvcXsxtWar17w4oV8MUXcRy9wPXqBbffDv/9L1x/fdzRlLzq6pDgbbtt3JEUtLwpm0SkZGwNvGZmk4B3CH3wnjOz4WZ2VLTN/cDmZjYduAi4PKZYRaRENdgHz90/BfrkIJYG7bJLuJ08WTUfG2XQIHjxRbjySjj4YNhnn7gjKlkzZoSWs716xR1J4cqnsklESoO7TwJ2T/H8lUn3VwEn5DIuEZFkBTVNwq67glnoRiYbwQzuvhs6d4ZTToGlS+OOqGRNnx5uVYMnIiIiIplUUAleq1bhglgJ3iZo2xYefRQ+/xwuuCDuaEpWIsFTDZ6IiIiIZFJBJXgQhpSfMCHuKArcd78Lv/0tPPQQPP543NGUpE8+gWbNQmWqiIiIiEimFFyC16dPuDhW68JN9H//B9/5DgwdGmrzJKemTw+10WUF9w0UERERkXxWcJeXfaIhFSZPjjeOgldREZpqAgwcCFVV8cZTYqZPV/NMEREREcm8gkvw+vYNt2qmmQHbbAN33QVjxsCIEXFHUzLWrQu10ErwRERERCTTCi7B69IF2rXTQCsZc8opYfLz4cPhjTfijqYkfP01rFypETRFREREJPMKLsEz00ArGXf77dCjR2iquXhx3NEUvU8+CbeqwRMRERGRTCu4BA9CP7zJk6G6Ou5IikSbNqE/3syZcO654B53REVNUySIiIiISLYUbIK3cuX6C2XJgH33hauvhieegIcfjjuaojZ9ehjjplu3uCMRERERkWJTkAmeBlrJkssvhwMOgPPPV/acRdOnh/FtKirijkREREREik1BJng77wxNm8K778YdSZEpLw+1dxUVoT/e2rVxR1SUNEWCiIiIiGRLQSZ4TZvCbrvBuHFxR1KEunWDe++Ft9+GYcPijqbouIdBVjSCpoiIiIhkQ0EmeAD9+sH48WFOMcmw44+HwYPh2mth1Ki4oykq8+bBkiVK8EREREQkOwo6wfvmG3UVy5pbbgntCE87DRYujDuaovHhh+F2xx3jjUNEREREilPBJnh77RVu1UwzS1q1gscfhzlz4JxzNHVChiQSvB12iDcOERERESlOBZvg7bwzVFbCO+/EHUkR23NPGDECnn4a7r8/7miKwrRp4e9WUySIiIiISDYUbIJXUQG7764avKy7+GI4+GC48ML11U+y0T78ELbbLgxYKiIiIiKSaQWb4EFopvnuu1BdHXckRaysDB58EJo3h5NPhtWr446ooE2bpv53IiKFysy6mtlrZjbVzN43swtTbNPfzL4xswnRcmUcsYpI6SroBK9fP1ixAqZOjTuSIte5c2ii+d578Nvfxh1NwVq9Gj77TP3vREQKWBVwsbvvBOwLnG9mO6fY7n/u3jdahuc2RBEpdQWf4IGaaebE0UfD0KFw443wyitxR1OQPvkk1DarBk9EpDC5+9fu/m50fykwFegcb1QiIjUVdIK3ww5hsEcNtJIjN90EO+0EgwaFCd2kUTSCpohI8TCzHsDuwNgUq79jZhPN7AUz653TwESk5BV0gldWFvrhvfVW3JGUiBYtwtQJCxaEidA1dUKjTJsWbpXgiYgUNjNrBfwN+IW7L6m1+l2gu7v3AW4D/l7HPoaY2TgzGzdPP5qKSAYVdIIHsP/+MHEiLFsWdyQlok8fuP56+Oc/4a674o6moHz4IXTqBK1bxx2JiIhsLDNrQkjuHnX3p2uvd/cl7r4suv880MTMOqTY7h537+fu/Tp27Jj1uEWkdKSd4JlZuZm9Z2bPZTOgxtpvv9CvSc00c+jnP4cBA8IUCu+/H3c0BUMjaGZHvpZNIlJ8zMyA+4Gp7n5zHdtsFW2Hme1NuNZakLsoRaTUNaYG70JCZ+K8su++4faNN+KNo6SUlcHIkdCmTZg6YdWquCPKe+6hBk/NM7MiL8smESlK+wOnAQclTYNwhJkNNbOh0TbHA1PMbCJwK3CSu/o0iEjupJXgmVkX4IfAfdkNp/HatYPevWHMmLgjKTFbbgl/+QtMngyXXRZ3NHlv7lxYvFg1eJmWz2WTiBQfdx/t7ubuuyVNg/C8u9/t7ndH29zu7r3dvY+77+vuukIRkZxKtwbvFuBSYF1dG8TZWXi//eDNN2FdndFJVhxxRGiueeut8PzzcUeT1xIjaCrBy7i8LptEREREcq3BBM/MfgTMdffx9W0XZ2fh/fYLtSOJUQolh66/HnbdFc48E+bMiTuavPXBB+FWCV7mFELZJCIiIpJr6dTg7Q8cZWafA08Q2p0/ktWoGmn//cOt+uHFoLIyTJ2wZAmccYaqUeswZUrosti1a9yRFJW8L5tEREREcq3BBM/df+3uXdy9B3AS8B93PzXrkTVCr17QoYP64cWmd+8wCfqLL8Jtt8UdTV6aMgV22QXCuGqSCYVQNomIiIjkWsHPgwfhonm//VSDF6tzz4Ujj4RLLw0TE8q33EOC17t33JGIiIiISLFrVILn7qPc/UfZCmZTHHAAfPwxzJoVdyQlygweeAA23zxMnbBiRdwR5Y05c2DBglCDJ9mRz2WTiIiISC4VRQ0ewEEHhdvXXos3jpLWoQM8+CBMnQqXXBJ3NHljypRwqwRPRERERLKtaBK8Pn3CnHj/+U/ckZS4Qw8Nyd1dd8E//hF3NHlBCZ6IiIiI5ErRJHhlZfD976sGLy+MGAF77AGDB6vNLCHB69gRttgi7khEREREpNgVTYIHcOCB8Nln8MUXcUdS4po2hcceg5Ur4fTTS37qhPffV+2diIiIiORG0SV4oFq8vLDDDvCnP8Grr4YpFEpUYgRNJXgiIiIikgtFleD17h3G+VCClycGD4bjjoPf/AbGj487mlh8+SUsW6YpEkREREQkN4oqwSsrg/79w0Ar7nFHI5jBPffAVluFqROWLYs7opzTACsiIiIikktFleBBmC5h5swwJ57kgfbt4eGHYfp0+MUv4o4m5xIJnmrwRERERCQXii7BO+ywcPvCC/HGIUn694df/xruvx+eeiruaHJq4kTo2hXato07EhEREREpBUWX4PXsCdtvrwQv7wwbBnvvDeecAzNmxB1Nzrz3Huy+e9xRiIiIiEipKLoED+Dww2HUKFixIu5I5FtNmsCjj0JVFZx6KlRXxx1R1i1fDh9+qARPRERERHKnaBO81atDkid5pFcvuOMOeP11uO66uKPJukmTwmA/SvBEREREJFeKMsH7/veheXM108xLp50WRtS86ip46624o8mqCRPCrRI8EREREcmVokzwKivDpOdK8PKQGdx1F3TpAgMHwpIlcUeUNe+9B+3ahUFWRESk8JlZVzN7zcymmtn7ZnZhim3MzG41s+lmNsnM9ogjVhEpXUWZ4EFopvnJJ5ouIS9ttlnoj/f553DBBXFHkzWJAVbM4o5EREQypAq42N13AvYFzjeznWttcziwXbQMAe7KbYgiUuqKNsH74Q/D7bPPxhuH1GH//eHKK8MceY89Fnc0Gbd2LUyerOaZIiLFxN2/dvd3o/tLgalA51qbHQ085MFbQFsz2zrHoYpICSvaBG+bbaBvX3j66bgjkTpdcUVI9M49Fz77LO5oMmratDDQT9++cUciIiLZYGY9gN2BsbVWdQaS5wOayYZJIGY2xMzGmdm4efPmZStMESlBRZvgARx7LIwZA19/HXckklJFBTzySLg/cGCYQqFIaIAVEZHiZWatgL8Bv3D32p3JUzXM9w2ecL/H3fu5e7+OHTtmI0wRKVFFneAdd1y4feaZeOOQevToAX/+M7z5Jvzud3FHkzHvvRcG+9lhh7gjERGRTDKzJoTk7lF3T9VOaCaQPLxWF2BWLmITEYEiT/B22ilcYKuZZp476SQ4/XS45hoYPTruaDLi3Xdht91CJaWIiBQHMzPgfmCqu99cx2bPAqdHo2nuC3zj7mpLJCI5U9QJnllopjlqFCxYEHc0Uq/bbw+1eQMHwuLFcUezSaqrYfx42GuvuCMREZEM2x84DTjIzCZEyxFmNtTMhkbbPA98CkwH7gXOiylWESlRRV+/cNxxcO218I9/wFlnxR2N1Kl1a3j88TDoytCh4X6Bzi8wbRosWwZ77x13JCIikknuPprUfeySt3Hg/Mbsd/Hi7Iz67Rv0/Mvv/WZz34o5N/tWzPmh6BO8PfYII2o+8YQSvLy3995w9dVhdM3DD4dBg+KOaKO8/Xa4VYInIiLp+OQTOProuKMQkWLRYIJnZpXA60CzaPun3P2qbAeWKWah1d/vfw+zZkGnTnFHJPW67DJ4+WU4//xQm9erV9wRNdrbb4e53LffPu5Iiluhl00iIgk77bR+UOlMy1ZjmGw2slHM2d9vNvetmNfbaafs7Lch6dTgrQYOcvdl0chRo83shWjyzoJw2mlh/I7HH4eLL447GqlXeXmY/LxPHzjlFHjjDWjSJO6oGuXtt0P/u7Ki7uGaFwq+bBIRAWjRIrQ4EhHJhAYvQT1YFj1sEi0F1Vp1++1Dc7mHH447EklL165w773wzjtwVWFVyKxcCZMmqXlmLhRD2SQiIiKSaWnVMZhZuZlNAOYCr7j72OyGlXmnngoTJ8LkyXFHImk57jg4+2y47jp47bW4o0nbhAlhvnYleLlRDGWTiIiISCalleC5e7W79yVM1rm3me1SexszG2Jm48xs3Lx58zId5yY76aQwJ9nIkXFHImm75RbYbrvQxrZA5rnQACu5VQxlk4iIiEgmNaqXkLsvBkYBA1Ksu8fd+7l7v44dO2YovMzp2BGOOSYkeKtWxR2NpKVly9Bxcu5cGDKkIMaxHTs2tDDdeuu4IykthVw2iYiIiGRSgwmemXU0s7bR/ebAIcC0bAeWDUOHwsKF8NRTcUciadtjjzAE6tNPw333xR1Ng8aMgX32iTuK0lBMZZOIiIhIpqRTg7c18JqZTQLeIfRzeS67YWXHgQeGFn933RV3JNIoF10EhxwCF14YZhHPUzNmwBdfwPe+F3ckJaNoyiYRERGRTGlwmgR3nwTsnoNYsq6sDH76U7jkkjDS4W67xR2RpKWsDB56KHxgp5wCb74JzZrFHdUGRo8Ot9/9brxxlIpiKptEREREMqXkZuo64wyorITbbos7EmmUrbeG+++H996DK66IO5qURo+GVq30w4GIiIiIxKfkErzNN4dBg0KF0OzZcUcjjXLUUXDeeXDTTfDyy3FHs4HRo2G//cJorSIiIiIicSi5BA/g4oth7VrV4hWkG2+EnXcOWXoeDXm/eHGYY1HNM0VEREQkTiWZ4G23Hfz4x3DnnbBsWdzRSKM0bx6mTli0CM46K2+mTnjzzRCKEjwRERERiVNJJngAl14aal3uvjvuSKTRdtsN/vAHeO65kKXngdGjQ9NMTXAuIiIiInEq2QRvn33g0EPh+uth6dK4o5FG+9nP4PDDQ3vbKVPijob//jdM2deyZdyRiIiIiEgpK9kED2DECJg/H265Je5IpNHMYORI2GwzOPlkWLkytlCWLYOxY+Ggg2ILQUREREQEKPEEb6+9Ql+8G2+EhQvjjkYabYst4MEHQw3eZZfFFsbrr0NVFRx8cGwhiIiIiIgAJZ7gAfzud6GJ5tVXxx2JbJQBA+AXvwhDov7rX7GE8OqrYd71/feP5fDSGHkyKI+IFC4ze8DM5ppZyv4BZtbfzL4xswnRcmWuYxSR0lbyM3b17g1Dh8Idd8DZZ8Ouu8YdkTTatdfCf/4DZ54JkybBVlvl9PCvvhrmv2vePKeHlY3x7rvQoUP4G9lqK9hyy/X3az/XoQOUlfxvYCKyoZHA7cBD9WzzP3f/UW7CERGpqeQTPIBrroEnn4Tzzw+DZZjFHZE0SmVlmDphzz3D/HgvvJCzC/N582DixPA3JAVg663hmGNg9uywvPlmuE3Vh7O8PDQDbigR3Gqr0BdUBYdISXD3182sR9xxiIjURQke0L59qAQaMiR06TrjjLgjkkbbeWf44x/h3HPhT3+CX/4yJ4d97bVwq/53BaJTpw2n1nAPI+Ukkr7Zs2HOnJqPZ88OfT1nzw4dLmtr1qzhJDDxWEOtipSC75jZRGAWcIm7v197AzMbAgwB6NatW47DE5FipgQvMngwPPwwXHhhuFjv2jXuiKTRfvpTePFFuPxyOPBA6Ns364d89VVo0wb69cv6oSRbzKB167Bst139265bB4sW1Z0IzpkDn38Ob70VqndT9flr1ar+JDDx3JZbQtOmWTllEcmqd4Hu7r7MzI4A/g5sULi4+z3APQD9+vVTB2ERyRgleJGysjDqfp8+oSvXyy+r+03BMYP77gsToZ98MowfDy1aZO1w7vD88+EHgQp9k0pDWRlsvnlYeveuf9uqqpDkpaoNTDw3ZQr8+9+weHHqfbRvn35/wfKQ73O/AAAgAElEQVTyzJ+viDSauy9Juv+8md1pZh3cfX6ccYlI6dBlaZKePeHmm0NTzRtuiHXkfdlYHTqEqthDD4WLLoK7787aoSZPhpkzYdiwrB1CCllFRejzt/XWDW+7ahXMnVt/M9GxY8PtihUbvr6sLPQXTKeZaNu26i8okkVmthUwx93dzPYmjFi+IOawRKSEKMGr5eyzww/qv/kN7L47/OAHcUckjXbwwfCrX8Ef/hCmUTjmmKwcJjErwxFHZGX3UkoqK6Fbt7A0pKH+gnPmwAcfhPtr1274+qZN0+8v2KpV5s9VpMCZ2eNAf6CDmc0ErgKaALj73cDxwLlmVgWsBE5y1xwtIpI7lo0yp1+/fj5u3LiM7zdXli+H73wn1M6MHdtwtxzJQ2vWhLkLPvssTJ3QuXPGD/Hd74bBF8ePz/iu85KZjXf3gu5tWOhlU6O4r+8vWF8z0dmzQ+1hqv8FLVum31+wWbPcn6MIKptEJH/FVT6pBi+Fli3hmWdgn31CS7/Ro6FLl7ijkkZp2hQeeyxUw55+OrzySkY7VS5YEEbY/+1vM7ZLkcwyC3342rcPo8zWp6oK5s+vPxH84IMw3+SiRan30a5d/bWBifsdO6q/oIiISBYpwavDttvCSy+FwRgPPRRGjQrXKFJAtt8ebrstDJF6441w6aUZ2/VLL4UBFX/4w4ztUiQ+FRXrE7A+ferfdvXqhvsLvvNOeG7Zsg1fX1YWkrx0mom2a6f+giIiIo2kBK8ee+4Jzz0Hhx8eWvu99BL06hV3VNIoZ54ZJj6/4go46KCMzWfwz3+GMS00PYKUnGbNwjwy6cwls2zZ+uSvrtrBadPC/TVrNnx9kyb11wbW7i+oZFBEREQJXkMOOCC0SvrhD2H//eGpp+B734s7KkmbGdxzT+hMecop8O67mzxwxMqVIcEbOFBTaYjUq1WrsGy7bf3buYepIupLBGfOhHHjQu3hunUb7qNFi4YHjUncVlZm53xFRETygBK8NOyzD7zxBhx5ZGiyOXx4mEJB3UgKRLt28Mgj0L9/mMn+/vs3aXcvvRQG4jnhhMyEJ1LyzML3tF072Gmn+retrg6dYOsbNGbatNCufuHC1Pto27bhuQUT/QU1yaWIiBQY/edK0w47hB+Pf/rT0Nrv6afhrrtgr73ijkzScsABYe6LESPgsMPgJz/Z6F39v/8X5rnu3z9z4YlImsrLQ/voLbaA3Xarf9s1a2r2F0xVO/juu+F26dINX2+Wur9gqsSwXTtV6YuISF5QgtcIbdqEgRmPOQZ++ctQs3fssWEkxb59445OGnTVVWGSwyFDYN9905tzrJZVq0LzzBNP1A/7InmvadMwBHI6wyAvX95wf8GPPgr3V6/e8PUVFen3F2zdWv0FRUQkaxq8RDWzrsBDwFbAOuAed/9TtgPLV2bh4v7ww8M82rfdBn/7W8gXBg2Co4+GrbfOzrHdw5gFCxeG5ZtvwjzGa9eGUc6rq6F58zDNQ8uWoRVSp07qbvKtJk1Cht63L5x6Krz2WqPb2b70UvihX80z46eySTKqZUvo2TMs9XEPhW99U0rMmhVqBufODQVzbc2bN5wEJm6bN8/O+YqISNFKpw6iCrjY3d81s9bAeDN7xd0/yHJsea1NG7jmGrjkktCla+RIOPfcsPTuHUbd3GWXMP1U4v/0ZpuFfKKsLFwjLF8eEralS8P4AnPmpF7mzl2f1K1d2/hYN988/IC9004hpl12gb33zl4imtd69oQ77ghz4117baMnsvvrX8O0YgcemKX4pDFUNknumYVfz9q2DW3367NuXcP9BT/+OEy2On9+6n1stlnDcwsm+gs2aZL58xURkYLTYILn7l8DX0f3l5rZVKAzoIsowv/4iy+Giy6CyZNDDc8rr4RavXvv3fj9tmkT/n9vuWVIzDbfPCyJeYs33zxs07Rp+J/epEm47li1KiSOy5eH+Yi/+ioMPvfll/DWW/DEE+uPsd12oWvawQfDEUeE64iScOqp8OKLMGxYOPnvfCetly1eDM88E6bV03VU/FQ2Sd5LzPnXsSPsumv9265d23B/wQkTwu2SJRu+3iz8Y2hobsGttgr/RNRfUESkaDWqF5GZ9QB2B8ZmI5hCZhb6+++2G/zqV6GGbvZs+PDD8D97zpzwP7m6en2LndatwwjirVuH5CqR0G2xRfaaVS5ZAlOmwJtvwuuvh8Fi7r8/JCwHHhj6FJ54Ykhci5YZ3HknjBkT5jqYMCFkyw148smQQJ9xRvZDlMZR2SQFr0kT6Nw5LA1ZsWJ9E4+6agdHjw63q1Zt+Pry8pqJX33NRNu0UX9BEZECY+6e3oZmrYD/AiPc/ekU64cAQwC6deu25xdffJHJOCVL1q0LU8T9/e+hdurjj0Ny+eMfw1lnhQquov3f/uabYVLDk0+Ghx9ucPP99gsJ8uTJRfye1MPMxrt73k3trrJJpA7uodCqLxFMLHPnhs7ctVVWNjy3YOJ+ixa5P0fyt2xqjH79+vm4cePiDkNEMiyu8imtBM/MmgDPAS+5+80Nba+CqjC5w3vvwQMPhLFIFi0K/QkvuihUdDVrFneEWTB8eBhd8+GHQ9PNOnz4Iey4I9xwQ+h3WYry8SJKZZNIhqxbFzp615cEJp6bNy/1PhJ9CxpqJrrFFhlt556PZVNjqWwSKU55m+CZmQEPAgvd/Rfp7FQFVeFbtSo0SbzpJpg0Kfx/vvhiOP/82H6kzY6qqtA2deLE0FSzjhH0Lr8cbrwRZswo0cFpyL+LKJVNIjFZuzYkeXUlgsmPv/km9T7S6S+45ZbQoUOD/QXzrWzaGCqbRIpTXOVTOn3w9gdOAyab2YToud+4+/PZC0viVlkZBpo87TR49dUwJcSll8LNN4eJ3s85p0hq9Coq4JFHoE+fUE35v/9tMMHdypVw331w5JGlm9zlKZVNInFo0iTMwdOpU8PbrlzZcH/BMWPC7cqVG74+MbF9ff0FRUSkhnRG0RwNlGCPI4HQ1+yQQ8Lyv/+FWQV+9rPQVPH668OALAXfH617d/jzn+Gkk0KTzeHDa6x+4okw0vnPfx5TfJKSyiaRAtC8OfToEZb6uIc5g+qrDZw9OzQpmTMndX9BEREBGjmKppS2730PRo0K00BcfnkYm+TOO+HWW8Pc4QXtxBPD1AkjRsChh4aTJVxz3HZbmDuwf/94QxQRKVpmoQ9fmzZhDp/6rFsXOoknkr5DDslNjCIiBUIT4UijmMEPfgDvvAP33ANTp8Kee4YJ3hctiju6TXTrraEP3sCB357M6NFh4JkLLiiCmkoRkWJQVhb68PXuHYZ6FhGRGpTgyUYpLw/98D76KDTZvPfeMCH7U0+FWq+C1Lp1GD7066/h8MOhe3d+f8ALdCybz6kVTzT8ehERKXpm9oCZzTWzKXWsNzO71cymm9kkM9sj1zGKSGlTgiebpF07uOWWUKPXuTOccEKYQ++rr+KObCPttVeY7X3sWMZ92ZEXOZyL1t1Iy58PhkcfjTs6ERGJ30hgQD3rDwe2i5YhwF05iElE5FtK8CQjdt89TJh+ww3w8suw886hCWdB1ua99RYAv+c3tGUR53EnrFgRhg8VEZGS5u6vAwvr2eRo4CEP3gLampnGYBaRnFGCJxlTUREmAZ88Gfr1g5/+NLR0LLjavBkzGMeePMOxXMifaMPS8PyXX8Ybl4iIFILOwIykxzOj52owsyFmNs7Mxs2ra/J4EZGNoARPMm7bbeHf/4Y77ghTK+yyS2jdWCi1ed66DZdwIx2Zy0XcvH5Ft27xBSUiIoUi1ZBcG/wHdPd73L2fu/fr2LFjDsISkVKhBE+ywgzOOw8mTAiDr5x6KvzkJzB/ftyRNWD4cJ5b8j3+S3+GMWx97V2LFmEKBRERkfrNBLomPe4CzIopFhEpQUrwJKu22y7U4l17LfzjH6E275//jDuqOlx9NSuvupaLWt/HDlt/wzndXg6ZavfuoUPhwIFxRygiUtqqq2HBApg+PYzu9fLLcUeUyrPA6dFomvsC37j713EHJSKlQxOdS9aVl4eJ0Y84Ak47DY46Cs48M4y+2aZN3NFFhg2Dq69m+K7PMX3ylrz6d2hy0PS4oxIRKT5VVbB4cZhvtLHLkiVxR4+ZPQ70BzqY2UzgKqAJgLvfDTwPHAFMB1YAZ8YTqYiUKiV4kjO77RZ+cL36arjuutBP7777wsTpsYqSu/eOuoob/nUEZ50FBx0Uc0wiIvmsqmrjErRFi2Dp0vr3XVkZ5uBJLF26wK671nwuefne93JzzhF3P7mB9Q6cn6NwREQ2oARPcqpp09CV7cgjQy3eYYfB2WfDTTfFVJsXJXfLTh3Kye9cRceOxg03xBCHiEiurV27viZt4cLGJWnLltW/7+bNayZh3bpBnz51J2nJS2VljV2tWxcOOXfu+mXePJj7KcyencX3R0SkQCnBk1jsuy+89x5cdRXceCO89BLce29I+HLCPSR3w4fjZ5zJuWvv5OOPjX//G9q3z1EMIiKbau3aja9JayhJa9GiZuLVo0eY9DSdJK1Zszp36x4q8Wokax/UepyUzM2fH7repaLBJ0VENqQET2JTWQnXXw/HHgtnnAEDBsDgwaE2b7PNsnjgpOSOs87i5p3u5ZFfGcOGwYEHZvG4IiKprFmz8Una8uX177uxSVr79uG2bdt6k7SEREvNBQui5eNwu3Bh0nMpltWrU++vTRvYYouw9OwZfgzcYouQyCWeTzzu0CHMv2qpJiUQESlhSvAkdvvsU7M277nn4A9/CAOyZPwft3s40O9+B4MH88RB93DJwDJOOAH+7/8yfCwRKR3JSVpjmzuuWFH/vlu2rJmE9eyZXi1au3ahXXya4S9eHC1fwOKJqZO02s8tXlz3PisqYPPN1y/bbgt77x3u107WErdp5JQiItIAJXiSFxK1eT/5CZx/PgwaFGYmuOOO0G0jI9zhyivhmmtg8GAe/f49DDq9jAMOgIcegjJNGiJS2lav3viatIaStFataiZevXqll6C1bZtWkpboTvftMqPW41rLokU1H69cWf/+27Spmaz16hVu27ev+Xzy0rq1atdEROKgBE/yyp57wpgxMHIkXHYZ7LFHmCT9qqvCj9YbzT1U0Y0YgQ8+mxu3/zOXnl7GgQfC3/++QZ9+ESlUq1ZtfJLWUJbT2CQtubljkyYb7K6qKvRFW7o0jP6fuL/0K1gyNfW6GttFj7/5puH8srx8fSiJpXPnmo9rL+3arU/iUoQvIiJ5Sgme5J2yMjjrLPjxj+H3v4fbb4fHHgvPXXQR7LBDmjt69FG44gr44ovw8/OSJcweeDHnzL6B5+43jj8eHn5YyZ1IWhLfpy+/DCMijhgBAwdm51jJSVpjmzuuWlX/vlu3rpmEbb99vUmat23Hysp2LG/SluVrmrB8eUimli9fv9R4PDfFc9GybFnNBK2hfDKhWbMQduvWoShr3To0Z9x223B/s81SJ2fJj1u0UG2aiEipUIIneatdO7jhBvjlL8O15H33hWabP/gBDB0aBmVp3ryOF593Htx9d6i5AxYtKeNu+w3X/u0qVq8zbrstNAXVBY8IDSdvjz4KQ4asryb64ovwGOpO8lau3PiatBRJ2jqMNTQNS6vNWbNZR1a22ZJVrbdn5ZYdWLVNe1Y2b8+qyrasbNaWVU3bsKpJa1aWt2JVeUtWWgtWUcnK1WWsWhUOsXJldPslrPpo/ePayVtjNW8eus21bBkSq8T9zp1hxx1rJmqJ2/rup9mNTkREBFCCJwWgU6fQF+/KK8NUCnfdFUbebNkyJHnf/z7stVdoLdW+PZQ9/ihr7rqfz9iOcfTjRQbwN45jpbfgSHuFm6YcynbbxX1WEoe1a2HOnNBcrays5m15eUj4Ews0fJtt7uuXdes2/ba6OjQLTNxWVUH1P5+n6po/U7W6E1V0o/qLcqrOfIyqZ9dS3Xs3qtY6Vbe+SvWKH1FFBVVUUE05VSsqqBo8nupbqqlavpqqFWuoXrGaqpVrWbOymjXVZesTsqRlNV1YQ0/WlLdgTUUL1lQ0Z01ZZVhoxppmTVnTrAlr1lWEpaqMNVVlVFUndZJdFi1fNe79bNo0JF+Vletvk++3aRNuaydmqZK1urZp3lz9eUVEJF7mUQ1HJvXr18/HjRuX8f2KQLhIHzUKnn4ann0WZs1av668HMqq17CW9T95b858fswznMed7G4Tw5WuNJqZjXf3fnHHsSnM+jlkvmxqTFKYattUSVwh/5ka62hWUU3TinU0beI0bQpNmxpNK42mlWU0rSynaTOLnt/4pXnz1IlaqiSuslKJV7EqhrJJ100ixSmu8kk1eFJwmjSBQw8Ny513wldfwbhxMGMGzJ4N/vubqGQl2/AZvXmfvkygjOiHjG7d4w1eYtWtWxi8J1GblbhN3E8kWNC4203d1iwkH4nb5PuZvK2oCEt5edL9Y4+mgrVUUEU51VEdXVjK77uHiqZlVPzyZ5QvmFNzHdVUdNmaionja+6vHCoqygBlUyIiInFQgicFzQy6dAnLtx79c+gjlGrjESNyFpvkn44dQ/dMSdJ9YurvS/fuMDiao6TsnJp98CC0T7zuBmifmzBFREQkPfqJVYrPiBHh4jOZWRiZJVuj/okUqlTflxYtav4YMnBgGOGoe/fwXerePTzW90lERCTvNJjgmdkDZjbXzKbkIiCRTZbqYvThh0N7TikqKp8yIN3kbeBA+Pzz0I7188+V3ImIiOSpdJpojgRuBx7KbigiGTRwoC5AS8NIVD5tOn1fREREikaDNXju/jqwMAexiIg0isonERERkZoy1gfPzIaY2TgzGzdv3rxM7VZEZJOobBIREZFSkrEEz93vcfd+7t6vY8eOmdqtiMgmUdkkIiIipUSjaIqIiIikycwGmNmHZjbdzC5Psf4MM5tnZhOi5ew44hSR0qV58ERERETSYGblwB3AocBM4B0ze9bdP6i16V/d/YKcBygiQnrTJDwOvAnsYGYzzWxw9sMSEWmYyicRybG9genu/qm7rwGeAI6OOSYRkRoarMFz95NzEYiISGOpfBKRHOsMzEh6PBPYJ8V2x5nZAcBHwC/dfUaKbUREskJ98ERERETSYyme81qP/wn0cPfdgH8DD6bckUb4FZEsUYInIiIikp6ZQNekx12AWckbuPsCd18dPbwX2DPVjjTCr4hkixI8ERERkfS8A2xnZtuYWVPgJODZ5A3MbOukh0cBU3MYn4iIRtEUERERSYe7V5nZBcBLQDnwgLu/b2bDgXHu/izwczM7CqgCFgJnxBawiJQkJXgiIiIiaXL354Hnaz13ZdL9XwO/znVcIiIJaqIpIiIiIiJSJJTgiYiIiIiIFAkleCIiIiIiIkVCCZ6IiIiIiEiRUIInIiIiIiJSJJTgiYiIiIiIFAkleCIiIiIiIkVCCZ6IiIiIiEiRUIInIiIiIiJSJJTgiYiIiIiIFAkleCIiIiIiIkVCCZ6IiIiIiEiRUIInIiIiIiJSJJTgiYiIiIiIFAkleCIiIiIiIkVCCZ6IiIiIiEiRUIInIiIiIiJSJJTgiYiIiIiIFIm0EjwzG2BmH5rZdDO7PNtBiYikQ2WTiORaQ+WOmTUzs79G68eaWY/cRykipazBBM/MyoE7gMOBnYGTzWznbAcmIlIflU0ikmtpljuDgUXu3gv4I3B9bqMUkVKXTg3e3sB0d//U3dcATwBHZzcsEZEGqWwSkVxLp9w5Gngwuv8UcLCZWQ5jFJESl06C1xmYkfR4ZvSciEicVDaJSK6lU+58u427VwHfAJvnJDoREaAijW1S/erkG2xkNgQYEj1cbWZTNiWwPNUBmB93EFmicys8uT6v7jk8VjqKtWwqhL9XxZgZijEzdsjhsdIpdwqtbIrzMy7FY5fiOZfysXNZPn0rnQRvJtA16XEXYFbtjdz9HuAeADMb5+79MhJhHinW8wKdWyEq1vNqhKIsmxRjZijGzCiUGHN4uHTKncQ2M82sAtgMWFh7R/lSNunYpXFcHTu+Y8dx3HSaaL4DbGdm25hZU+Ak4NnshiUi0iCVTSKSa+mUO88Cg6L7xwP/cfcNavBERLKlwRo8d68yswuAl4By4AF3fz/rkYmI1ENlk4jkWl3ljpkNB8a5+7PA/cDDZjadUHN3UnwRi0gpSqeJJu7+PPB8I/Z7z8aFk/eK9bxA51aIivW80lakZZNizAzFmBmKsZZU5Y67X5l0fxVwQiN3G+f7rGOXxnF17BI6tqnVgIiIiIiISHFIpw+eiIiIiIiIFICsJHhmdoOZTTOzSWb2jJm1zcZx4mBmJ5jZ+2a2zszyemSxdJjZADP70Mymm9nlcceTSWb2gJnNLYBh8RvFzLqa2WtmNjX6W7ww7pjyQUN/y2bW3cxejcqlUWbWJWndIDP7OFoG1X5tnsRYbWYToiUrg8k09J2x4NYo/klmtkfSuly9h5sSY9bfwzRj3NHM3jSz1WZ2Sa11OSmTNzHGz81scvQ+Zm2EuDRiHBh9xpPMbIyZ9Ulal7X3MY3vcTMz+2u0fqyZ9Uha9+vo+Q/N7LCG9mlm95vZxOgc/2dmH0XbXFH7GEnHXWBmq8zsk8QxNvG4T5nZMdF2881sTvT8q2bWPemc10V/L8vN7JUMnXPyseeY2bKk7/DZSccea2ZrovO+OAvHXmxms6LjfmRmizN93knrbzOzlUnbZPWzrnXcZUnbZf2zrufYWf+s6zl21j9rMxtpZp8lnV/f6HmzTP2PdfeML8APgIro/vXA9dk4ThwLsBNhTotRQL+449nEcykHPgF6Ak2BicDOcceVwfM7ANgDmBJ3LBk+r62BPaL7rYGPiulz28j3pMG/ZeD/AYOi+wcBD0f32wOfRrftovvt8inG6PGyHLyP9X5ngCOAFwjzfO0LjM3le7gpMebqPUwzxi2AvYARwCWN+RuJO8Zo3edAhzx4H/dL/J0Bhyf9PWbtfUzze3wecHd0/yTgr9H9naPtmwHbRPspr2+fQJuk4y4G/hBtMzNpvycBf42O+9fo9acCz0X73WVjjxvd/yOwINruUGBSdC7nJsVwHrA2k+ec4tiDCXOZ1X6/LwGWEMqfwcAyQjmUyWN/ux3wM8LgOhk97+h1/YBHgHVJ22T1s0467sPRe5fYLuufdT3HzvpnXc+xs/5ZAyOB41OULxn7H5uVGjx3f9ndq6KHbxHmiSkK7j7V3T+MO44M2RuY7u6fuvsa4Ang6Jhjyhh3f50Ucw8VOnf/2t3fje4vBaYCneONKnbp/C3vDLwa3X8taf1hwCvuvtDdFwGvAAPyLMacSOM7czTwkAdvAW3NbGty9x5uSow501CM7j7X3d8B1tZalbMyeRNizJk0YhwT/b1BzWuNbL6P6ez7aODB6P5TwMFmZtHzT7j7anf/DJge7a/Ofbr7kqTjLgEWRNusYP178xRwcPSahdHrnyBcIE4Hzt/Y40Zxd4+O+6m7vwI8Hm2X/J4fzfq/lYycc+1jA9XAxyne79OAV919YfS+lxPKn4wdu9Z2J0fvQUbP28zKgRui59YlbZPVzzrpuJcSWvZNz9VnXdexc/FZ13PeWf+s65Gx/7G56IN3FiEblfzTGZiR9HgmShQKioWmP7sDY+ONJHbp/C1PBI6L7v8YaG1mm6f52rhjBKg0s3Fm9paZHZOF+NJR1znkU1lSXyz58B7WJ5/ex/o48LKZjTezIXEHExnM+muNbL6P6ez7222iH7u/Aeora+rdp5n9hXBuFcBt0dPNgcpax+gWPTcj6bl5QI9NOO5sYHvgjRTb1X7Pm1posjuacCG8qeec6tg7ARdbaD6ZmHB+C0IykHgvVgLbZeHYM4EdCbU0/8nCeV9AmEOxkvAdS8j2Z30B8Ky7f02oOUq1XbY+6/qOne3PuqHzzuZnDTAiaob5RzNrlnSMjPyPTWuahFTM7N/AVilWXeHu/4i2uQKoAh7d2OPEIZ1zKxKW4jkNq1ogzKwV8DfgF0m/8paqdP6WLwFuN7MzgNeBrwjlU66+B5sSI0A3d59lZj2B/5jZZHf/JAtx1qeuc8insqS+WPLhPaxPPr2P9dk/eh+3AF4xs2lRbVsszOxAwgXodxNPpdgsU+9jOvtu7Pck1Y/t3+7T3c80sxeBq4ATgb9E+0p13FTHSCXd45YTpoToWWu77YE2wPeTjr2Pu4+Pvl8fEmrANuWcax/7n9H+9gDeJ9TgHFTPMTJ57IRewFPuXh09zsh5m1knwtQa/Qk/7iXL2mdd67h1ycpn3cCxs/pZp3neWfmso9tfE35EaEqYRuEyYHg9+2p0mbbRNXjufoi775JiSSR3g4AfAQPdPR//QdWpoXMrIjOBrkmPuwCzYopFGsHMmhCSu0fd/em448kDDf4tu/ssdz/W3XcHroie+yad1+ZBjLj7rOj2U0If4N2zEGND6jqHfCpL6owlT97D+uTT+1inpPdxLvAMoSlSLMxsN+A+4Gh3XxA9nc33MZ19f7uNmVUAmxGa023K9+dLQj+hRA3/CmB1rWN8QajR6Jr0XEdCn8mNOm50cft3Qm1KwveBvsBR7r466ZybJsVaTaj92Ohzrn3s6PPdKtruXmDPaNM5hFqcxHvRnNBMLmPHTtpuG9Y32cvkee9OSCimA7cCFWY2Pdomm5/1t8c1s88J/chOSdoum591ncfOwWfd0Hln87NOdLXx6D39C+vL0Mz9j/XsdIweAHwAdMzG/vNhoTgGWakgdNTchvUdQHvHHVeGz7EHxTfIigEPAbfEHUu+LOn8LQMdgLLo/ghgeHS/PfAZoeNyu+h++zyLsR3QLGmbj8ne4Bt1fmeAH1KzA/jbuXwPNzHGnL2HDcWYtM0wag6yktMyeSNjbAm0Tro/BhgQ02fdjXBxt1+t57P2Pqb5PT6fmoOsPBnd703NwRg+JfQhSrnP6G+4V9JxFxN+7U818MaT0TyFddsAAAQlSURBVHGTB974V7TfXTfhuAbcFB17G8LAO6uBw2ud8yXAPdH9c4ClhD7Fm3LOtY/dLWm7HwNvRdv9itA/sR3rB95on+FjNwWmEVpVWKbPO8Xf2LqkbbL2Waf4+16WtF1WP+sGjp3Vz7qBY2f9swa2Tvo7uwW4Lnqcsf+x2SqQpxPaik6IlruzVfjneon+0GZGf/RzgJfijmkTz+cIwiiMnxCaoMYeUwbP7XHga0I76ZnA4LhjytB5fZdQNT8p6Tt2RNxxxb2k+lsmNHk4Krp/POGi/iPCL/7Nkl57VlRuTQfOzLcYCaMFTo7+QUzO1t9yqu8MMBQYGq034I4o/skk/ciVw/dwo2LM1XuYZoxbRc8vIVxAzmT9aIk5KZM3NkZCk7WJ0fJ+zDHeByxifTk4rr7vWgbjauh7XEkYEXc68DbQM+m1V0Sv+5CkC+c69llG6AM2GZhCGHhperTNVdExFkav65l03AXAKsLF5eEZOO6jhJrDjwi1Scui93s260f560/oB7Y62uayDJ1z8rEXEa67JhJqqs5Ner/fAdZE531pFo79CWEArOtqfdYZOe8Uf2Mrk7bJ2med4rjLkrbL6mfdwLGz+lk3cOysf9aEfn2Jv7NHgFbR8xn7H2vRi0RERERERKTA5WIUTREREREREckBJXgiIiIiIiJFQgmeiIiIiIhIkVCCJyIiIiIiUiSU4ImIiIiISFEzsxPM7H0zW2dm/erYptLM3jazidG2Vyetuz96fpKZPWVmraLnh5rZZDObYGajzWznXJ1TXZTgiYiIiIhI0TCz/mY2stbTU4Bjgdfreelq4CB370OY5H2Ame0brfulu/dx990Ik51fED3/mLvv6u59gT8AN2fqPDZWRdwBiIiIiIiIZJO7TwUws/q2ccK8eABNosWjdUui1xvQvPbzkZaJ5+OkGjwRERERERHAzMrNbAIwF3jF3ccmrfsLYeL3HYHbkp4/38w+IdTg/TzHIW9ACZ6IiIiIiBQ8MxsbJWf3AUdF/eImmNlh6e7D3auj5pZdgL3NbJekdWcCnYCpwIlJz9/h7tsClwG/zdDpbDQleCIiIiIiUvDcfZ8oOTsbeNbd+0bLSxuxr8XAKGBAreergb8Cx6V42RPAMY0OPMOU4ImIiIiISMkzs45m1ja63xw4BJhmQa/oeQOOBKZFj7dL2sUPgY9zG/WGNMiKiIiIiIgUNTP7MaHfXEfgX2Y2wd0PM7NOwH3ufgSwNfCgmZUTKsKedPfnzKwser4NYMBE4Nxo1xeY2SHAWmARMCi3Z7YhC4PFiIiIiIiISKFTE00REREREZEioQRPRERERESkSCjBExERERERKRJK8ERERERERIqEEjwREREREZEioQRPRERERESkSCjBExERERERKRJK8ERERERERIrE/wcIG76ZYoOMJAAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.figure(figsize=[15,3])\n",
"plt.subplot(1,3,1)\n",
"plt.scatter(x_list,y_list,c=\"r\")\n",
"plt.plot(x_list,y_list,c=\"r\")\n",
"plt.plot(x,f(x), c=\"b\")\n",
"plt.xlim([-2,2.5])\n",
"plt.ylim([0,5])\n",
"plt.title(\"Gradient descent\")\n",
"plt.subplot(1,3,2)\n",
"plt.scatter(x_list,y_list,c=\"r\")\n",
"plt.plot(x_list,y_list,c=\"r\")\n",
"plt.plot(x,f(x), c=\"b\")\n",
"plt.xlim([0.9,1.2])\n",
"plt.ylim([0,5])\n",
"plt.title(\"Zoom around min\")\n",
"plt.subplot(1,3,3)\n",
"plt.scatter(x_list,y_list,c=\"r\")\n",
"plt.plot(x_list,y_list,c=\"r\")\n",
"plt.plot(x,f(x), c=\"b\")\n",
"plt.xlim([1.3333,1.3335])\n",
"plt.ylim([0,3])\n",
"plt.title(\"zoomed in more\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We consider linear regression where efficient solutions exist and we do not need a numerical solver using gradient descent.\n",
"\n",
"We generate some data points and apply the method"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEWCAYAAABi5jCmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHaNJREFUeJzt3Xuc3HV97/HXe5NgSLiTcBOy6wXhAKUIq0KxlIq2iAqeI1U8Kxq8xCIVqD68tPHaGh9g7Tlq7cUocIQsIqIiUAEvGKi2IBtAFJBLIQsINJtACIiFhHzOH7/fkJnJzO7M7MzvsvN+Ph7z2JnfZX6f+e3sfn7f608RgZmZWcVA3gGYmVmxODGYmVkNJwYzM6vhxGBmZjWcGMzMrIYTg5mZ1XBisEKTtFjST7v8noskPSlpVjfft1+l5/KFecdh3ePE0MckrZb0u/QPu/L4ct5xdYuklZLeXb88Iu6PiO0i4tk84qqWJr5n03O/QdIvJL0+77jakZ7Le/OOw7pndt4BWO7eEBE/yjuIfiBpdkRsarDqPyLilZIGgPcAF0naOyLWZ3R8sxouMVhDkv5Z0iVVr8+W9GMldpZ0haQJSY+lz/eu2nalpM9I+vf0SvhySbtKGk2vim+UNFS1fUg6XdK9ktZK+rv0n2SjuPaX9ENJj0q6U9KbO/hsQ+kxZ1fF+7eSfibpCUk/kLSgavvD08+yPr2iP7pq3SmS7kj3u1fSe6vWHS3pQUkfkfQIcN5kcUXEZuACYD6wb4vHf4Gk69Lj/0jSP0paUfc53yXpfuCaFt5vcfo5npB0n6SRdPmLJV0r6fH0d/TNqn1C0ovT5ztKOj/9boxL+ljld1mpFpT0+fR7c5+k17b+m7PMRIQfffoAVgOvbrJuHnAXsBj4Q2AtsHe6blfgTek22wPfAi6t2nclcA/wImBH4Pb0vV5NUko9HzivavsAfgLsAixKt313um4x8NP0+XzgAeCU9H0OTeM6sMlnWFl5n7rlQ+kxZ1dt95/AS4Bt09dnpeueD6wDjiO5kHpN+nphuv516ecU8EfAU8Ch6bqjgU3A2cDzgG0bxFL9+WYBpwHPALu1ePz/AD4PbAO8EtgArKj7nOen527byd4v3WYDsF+6/56Vcwt8A1ia7jMXeGXd7+/F6fPzge+RfC+G0t/lu6o+60aSUtEs4FTgIUB5/y34Ufe9zDsAP3L85SeJ4UlgfdXjPVXrXw48CowDb53kfQ4BHqt6vRJYWvX674Erq16/Abil6nUAx1a9fh/w4/R59T/OtwD/VnfsrwCfbBLXSlpPDB+rO/5V6fOPABfU7X818I4mx7wUOCN9fjTJP/m5k5y7xSTJY336T/N3wJur1jc9PkkS3QTMq1q3gq0TwwtbfL/5aRxvoi6Jpf/wl5NeHNStC+DF6T/7p4EDqta9F1hZ9VnvqVo3L913j7z/Fvyofbgqyd4YETtVPb5aWRERPwfuJbkavriyXNI8SV9Jqwo2ANcBO6m2l89/VT3/XYPX29XF8UDV83FgrwaxDgKvSKtA1ktaD4wAe7T8aZt7pOr5U1XxDQJ/VnfMV5JcTSPptZKuT6u21pNciS+oeq+JiPjvKY59fUTsBOwMXEZSQquY7Ph7AY9GxFNV21efx0bLmr5fRPyWJPn+OfCwpH+VtH+634dJvgc/l3SbpHc2OM4CkpLLeNWycZJSSsVz57kq7vrvguXMicGaknQaSRXIQyT/GCo+COwHvCIidgCOquwyjcPtU/V8UXrMeg8A19Ylsu0i4tRpHHcqD5BcYVcfc35EnCXpecC3Sapydk//uX+f2vPQ8vTFEfEkSWnlZEkvner4wMPALpLmVb3NPmytOobJ3o+IuDoiXkOSeH4NfDVd/khEvCci9iIpBfxTpV2hylqSUs9g1bJFwG9aPQdWDE4M1pCklwCfAd4GnAx8WNIh6ertSa7610vaBfhkFw75ISWN2vsAZwDfbLDNFcBLJJ0saU76eJmk/zHJ+86WNLfqMafNuFYAb5D0p5Jmpe9xtJLG9m1IEucEsCltSP2TNt+/RkSsA74GfGKq40fEODAGfErSNpKOIKmm6+jzSNpd0vGS5pNUCT0JPAsg6c+0pYPBYyTJpqa7byTdfy8GlknaXtIg8IH0mFYiTgx2uWrHMXxXSW+dFcDZEfGLiLgb+GvggvQq+QskDZlrgeuBq7oQx/eAVcAtwL8C59RvEBFPkPzjPYmkRPEIWxp2m/lnkiRWeUzaM6jBMR8ATiD5/BMkV9wfAgbSeE4n+Wf4GPC/SaqCpusLwHGSDp7s+Om2I8ARJA3InyFJqE938nnSxwdJzu2jJI3p70t3fRlwg6Qn0894RkTc1+AQ7wd+S1IF+VPgQuDcts+A5UoRvlGP5UtSAPtGxD15x1J2aTfSX0dEN0px1qdcYjArsbQq7UWSBiQdS1IauDTvuKzcepYYJJ0raY2kX1Ut20XJ4KS705879+r4Zn1iD5Lutk8CXwJOjYibc43ISq9nVUmSjiL5sp4fEQelyz5H0r3uLEkfBXaOiI/0JAAzM+tIT9sYlEx7cEVVYrgTODoiHpa0J8nAl/16FoCZmbUt60n0do+IhwHS5LBbsw0lLQGWAMyfP/+w/fffv9mmZmbWwKpVq9ZGxMJ29yvs7KoRsZxkCD7Dw8MxNjaWc0RmZuUiaXzqrbaWda+k/0qrkEh/rpnuG46OwtAQDAwkP0dHp/uOZmb9LevEcBnJZF2kP783nTcbHYUlS2B8HCKSn0uWODmYmU1HL7urfoNkSuD9lMxJ/y7gLOA1ku4mme73rOkcY+lSeOqp2mVPPZUsNzOzzvSsjSEi3tpk1THdOsb997e33MzMplbqkc+LFrW33MzMplbqxLBsGcybV7ts3rxkuZmZdabUiWFkBJYvh8FBkJKfy5cny83MrDOlTgyQJIHVq2Hz5uRnN5KCu8CaWT8r7AC3vFS6wFZ6O1W6wIJLImbWH0pfYug2d4E1s37nxFDHXWDNrN85MdRxF1gzK6Is2z6dGOq4C6yZFU3W0/84MdRxF1gzK5qs2z57eqOebvG022bWzwYGkpJCPSnpqt+MpFURMdz28drdwcysGzxeqHVZt306MZhZ5jxlfnuybvt0YjCzzHm8UHuybvt0G4OZZa7TOnNrj9sYzKw0PF6o2JwYzCxzHi9UbE4MZpY5jxcqNs+uama5GBlxIigqlxjMzKyGE4OZ1fDAM3NVkpk9xzeqMnCJwcyqlGXgmUs1veUSg5k9pww3qnKppvdcYjCz55Rh4FlZSjVl5sRgZs8pw8CzMpRqys6JwcyeU4aBZ2Uo1ZSdE4OZ1RgZgdWrk8nsVq8uVlKAcpRqys6JwcxKpQylmrJzryQzKx1Pp9FbLjGYmVkNJwYzM6vhxGBmTXmEcX9yG4OZNeQRxv3LJQYza8gjjPtXLolB0l9Kuk3SryR9Q9LcPOIws+Y8wrh/ZZ4YJD0fOB0YjoiDgFnASVnHUQSjo7BgQdIXW0qeuw7XisIjjPtXXlVJs4FtJc0G5gEP5RRHbkZH4Z3vhHXrtixbtw5OOcXJwYrBI4z7V+aJISJ+A3weuB94GHg8In5Qv52kJZLGJI1NTExkHWbPLV0Kzzyz9fKNG12Ha8XgEcb9SxGR7QGlnYFvA28B1gPfAi6JiBXN9hkeHo6xsbGMIszGwAA0O/VSMk+Nmdl0SFoVEcPt7pdHVdKrgfsiYiIiNgLfAf4ghzhyNVk9retwzSxPeSSG+4HDJc2TJOAY4I4c4sjVsmWwzTZbL58zx3W4ZpavPNoYbgAuAW4CfpnGsDzrOPI2MgLnngu77rpl2a67wnnnuQ7XzPKVS6+kiPhkROwfEQdFxMkR8XQeceRtZATWrk3aGiKS5/2YFDztgjXj70Y+PCWG5crTLlgz/m7kJ/NeSZ2Yib2SLDE0lPzB1xscTO4eZv3L343pK1OvJLPneNoFa8bfjfw4MViuPO2CNePvRn6cGCxXnnbBmvF3Iz9ODJYrT7tgzfi7kR83PpuZzVBufDYzs65wYjCzvuUBdI05MUyDv1Rm5VUZQDc+nsw8UBlA579jJ4aO+UtlVm6+p3VzTgwd8pfKrNw8gK45J4YO+UtlVm4eQNecE0MbqtsUBpqcOX+pzMrBA+iac2JoUX2bwrPPbr2Nv1Rm5eEBdM15gFuLms30OGtWcn/mRYuSpOAvlZkVRacD3Hw/hhY1azvYvDl5mJnNFK5KapEbqsysXzgxtMgNVWbWL5wYWuSGKjPrF25jaMPIiBOBmc18LjGYmVkNJwYzM6vhxGBmpeJZjXvPbQxmVhqVGQgqE1hWZjUGt/91k0sMZlYantU4G04MZlYantU4G04MZlYanoEgG04MZlYYUzUsewaCbDgxmFkhtHK7XM9AkA1Pu21mhdBsavvBQVi9OutoZoZOp912icFKwX3XZz43LBeHE4MVXitVDFZ+bljuhQW7dLKXE4MVnvuu9wc3LHdXcuG0aLCTfXNJDJJ2knSJpF9LukPSEXnEYeXgKob+4Ibl7kounNTR//i8SgxfBK6KiP2B3wfuyCkOKwFXMfSPkZGkoXnz5uRnt5JCP7ZRTefCKfPEIGkH4CjgHICIeCYi1mcdh5WHqxhsOvq1jWo6F055lBheCEwA50m6WdLXJM2v30jSEkljksYmJiayj9IKw1UMNh392kaVXDjF5k72zXwcg6Rh4HrgyIi4QdIXgQ0R8fFm+3gcg5l1amAgKSnUk5Iqq5lMWnhfxMQL290vjxLDg8CDEXFD+voS4NAc4jCzPtDfbVRrH+1kr7YSg6SBtI2gYxHxCPCApP3SRccAt0/nPc3MmnEbVfumTAySLpS0Q9oOcDtwp6QPTfO47wdGJd0KHAJ8dprvZ2bWkNuo2jdlG4OkWyLiEEkjwGHAR4BVEXFwFgGC2xjMzDrRy7mS5kiaA7wR+F5EbASKP/PeDNeP/bLNLBut3PP5K8Bq4BfAdZIGgQ29DMom5/vemlkvddRdVdLsiNjUg3gaclVSLU9PbGat6FlVkqTdJZ0j6cr09QHAOzqI0brEcweZWS+10sbw/4Crgb3S13cBZ/YqIJtaf/fLNrNeayUxLIiIi4HNAGkV0rM9jcom5X7ZZtZLrSSG30ralbQnkqTDgcd7GpVNyv2yrYjcU27maKVX0geAy4AXSfoZsBA4sadR2ZRGRpwIrDjcU25mmbLEEBE3AX8E/AHwXuDAiLi114GZWXn06wym7SpLqWrKEoOkt9ctOlQSEXF+j2Iys5JxT7mplalU1Uobw8uqHn8IfAo4vocxmVnJzOSect26yi9TqWrKEkNEvL/6taQdgQt6FpGZlc6yZbVXwzAzesp18yq/TKWqTu7H8BSwb7cDMbPymqk95bp5lV+mUlUrbQyXs2XSvAHgAODiXgZlZuUzE3vKdfMqv0ylqla6q36+6vkmYDwiHuxRPGZmhbFoUeN5yTq5yq8kzaVLk8SyaFGSFIqYTFtpY7g2i0DMzIqm21f5ZSlVNU0Mkp6g8X0XBERETOsWn2ZmRVemq/xuapoYImL7LAMxMyuislzld1MrbQwASNoNmFt5HREF7GRlZmbT1cr9GI6XdDdwH3Atyd3cruxxXGZmlpNWxjH8LXA4cFdEvAA4BvhZT6MyM7PctJIYNkbEOmBA0kBE/AQ4pMdxmZlZTlppY1gvaTvgOmBU0hqS8QxmZjYDtVJiOIFkGoy/BK4C/hN4Qy+DMjOz/LRSYlgCfCsd7fz1HsdjZmY5a6XEsANwtaR/k3SapN17HZSZmeWnlTu4fToiDgROA/YCrpX0o55HZmZmuWhn2u01wCPAOmC33oRjZmZ5a2WA26mSVgI/BhYA74mIg3sdmJmZ5aOVxudB4MyIuKXXwZiZWf5amXb7o1kEYmZmxdDJrT3NzGwGc2IwM7MarTQ+/4WknbMIxoppdBSGhmBgIPk5Opp3RGbWS62UGPYAbpR0saRjJanXQVlxjI4mtzYcH4eI5OeSJU4OZjNZKwPcPgbsC5wDLAbulvRZSS+azoElzZJ0s6QrpvM+1ltLl9be7xaS10uX5hOPmfVeS20MEREkg9seIZlZdWfgEkmfm8axzwDumMb+loH7m9ynr9lyMyu/VtoYTpe0CvgcyQ16fi8iTgUOA97UyUEl7Q28DvhaJ/tbdhYtam+5mZVfKyWGBcD/iog/jYhvRcRGgIjYDLy+w+N+AfgwsLnZBpKWSBqTNDYxMdHhYWy6li2DefNql82blyw3s5mplTaGT0TEeJN1bVcFSXo9sCYiVk1x3OURMRwRwwsXLmz3MNYlIyOwfDkMDoKU/Fy+PFluZjNTK1NidNuRwPGSjgPmAjtIWhERb8shFmvByIgTgVk/yXyAW0T8VUTsHRFDwEnANU4KZmbF4ZHPZn3CAxWtVXlUJT0nIlYCK/OMwawfVAYqVsakVAYqgqsJbWsuMfSYr9KsCDxQ0dqRa4lhpvNVmhWFBypaO1xi6CFfpVlReKCitcOJoYd8lTazlama0AMVrR1ODD3kq7SZq2yzznqgorVDyfx4xTY8PBxjY2N5h9G2+jYGSK7S/AdZfkNDSTKoNzgIq1dnHY1ZY5JWRcRwu/u5xNBDvkqbuVxNaDOZeyX1mKeTmJkWLWpcYnA1oc0ELjGYdcCNuTaTOTGYdcDVhDaTuSrJrEOuJrSZyiUG60tlGoNgljWXGKzveKoSs8m5xGB9x1OVmE3OicH6jscgmE3OicH6jqcqMZucE4P1HY9BMJucE4MVSha9hTwGwWxy7pVkhZFlbyGPQTBrziUGKwz3FjIrBicGKwz3FjLrtgW7dLKXE4MVhnsLmXXbXs/vZC8nBisM9xYy67Y523SylxODFYZ7C5l128ZnOtnLvZKsUNxbyKybHvpNJ3u5xGBmNmOtfbSTvZwYzMyshhNDj3i+fzMrK7cx9IDn+zezMnOJoQc8gtfMysyJoQc8gtfMysyJoctGR5N2hUY8gtfMysCJoYsqbQvPPrv1Oo/gNbOycGLookZtCwCzZnkEb9mUoVdZGWK0csq8V5KkfYDzgT2AzcDyiPhi1nH0QrM2hM2bnRTKpAy9ysoQo5WXIiLbA0p7AntGxE2StgdWAW+MiNub7TM8PBxjY2OZxdipoaHkD7Te4CCsXp11NNapMvweyxCj5U/SqogYbne/zKuSIuLhiLgpff4EcAfQ0dSwRePZQWeGMvQqK0OMVl65tjFIGgJeCtzQYN0SSWOSxiYmJrIOrSOeHXRmKMN9IcoQo5VXbolB0nbAt4EzI2JD/fqIWB4RwxExvHDhwuwD7NDISFKU37w5+emkUD5lKPmVIUYrr1wSg6Q5JElhNCK+k0cMZs2UoeRXhhitvPJofBbwdeDRiDizlX3K0vhsZlYkpWl8Bo4ETgZeJemW9HFcDnGYmVkDmY9jiIifAsr6uGZm1hqPfDYzsxpODGYl5mkxrBd8ox6zkvK0GNYrLjGYlZRvCGW94sRgVlKeFsN6xYnBrKQ8LYb1ihNDCbiB0RrxtBjWK04MBVdpYBwfh4gtDYxODuZpMaxXMp8SoxP9PCWG5903s06VaUoMa4MbGM0sa04MBecGRjPLmhNDwbmB0cyy5sRQcG5gNLOseUqMEhgZcSIws+y4xGBmZjWcGMzMrIYTg5mZ1XBiMDOzGk4MZmZWw4nBzArFk0bmz91VzawwfFe6YnCJwcwKw3elKwYnBjMrDE8aWQxODGZWGJ40shicGMys6zptQPakkcXgxGBmXTWduw560shi8B3czKyrfNfB4vAd3MysENyAXH5ODGbWVW5ALj8nBjPrKjcgl58Tg5l1lRuQy89TYphZ1/mug+XmEoOZmdVwYjDDM3qaVXNisL43nQFZM5GTpOWSGCQdK+lOSfdI+mgeMZhVeEbPLZwkDXJIDJJmAf8IvBY4AHirpAOyjsOswgOytnCSNMinxPBy4J6IuDcingEuAk7IIQ4zwAOyqjlJGuTTXfX5wANVrx8EXlG/kaQlQHrvJp6W9KsMYpuuBcDavINoQRnizDDGBbvAokFQ1YVSbB4fv39cWvvoVDtT/HMJLcd58O/BnG3ql0ZsfEa69Zc9iKveDDufuduvk53ySAxqsGyrmfwiYjmwHEDSWCcTQWXNcXZPGWIEx9ltjrO7JHU0+2geVUkPAvtUvd4beCiHOMzMrIE8EsONwL6SXiBpG+Ak4LIc4jAzswYyr0qKiE2S/gK4GpgFnBsRt02x2/LeR9YVjrN7yhAjOM5uc5zd1VGcpbhRj5mZZccjn83MrIYTg5mZ1ShMYpB0rqQ1zcYrKPGldBqNWyUdmnWMaRxTxXm0pMcl3ZI+PpFDjPtI+omkOyTdJumMBtvkfj5bjLMI53OupJ9L+kUa56cbbPM8Sd9Mz+cNkoYKGudiSRNV5/PdWceZxjFL0s2SrmiwLvdzWRXLZHEW5VyulvTLNIatuqd29LceEYV4AEcBhwK/arL+OOBKknEQhwM3FDTOo4Ercj6XewKHps+3B+4CDija+WwxziKcTwHbpc/nADcAh9dt8z7gX9LnJwHfLGici4Ev53k+0zg+AFzY6HdbhHPZYpxFOZergQWTrG/7b70wJYaIuA6YbJTpCcD5kbge2EnSntlEt0ULceYuIh6OiJvS508Ad5CMOK+W+/lsMc7cpefoyfTlnPRR32vjBODr6fNLgGMkNRrM2TMtxpk7SXsDrwO+1mST3M8ltBRnWbT9t16YxNCCRlNpFO6fSOqItDh/paQD8wwkLYa/lOTqsVqhzuckcUIBzmdapXALsAb4YUQ0PZ8RsQl4HNg12yhbihPgTWmVwiWS9mmwvte+AHwY2NxkfSHOJVPHCfmfS0iS/w8krVIylVC9tv/Wy5QYWppKowBuAgYj4veBfwAuzSsQSdsB3wbOjIgN9asb7JLL+ZwizkKcz4h4NiIOIRmp/3JJB9VtUojz2UKclwNDEXEw8CO2XJlnQtLrgTURsWqyzRosy/RcthhnrueyypERcSjJjNWnSTqqbn3b57NMiaEUU2lExIZKcT4ivg/MkbQg6zgkzSH5ZzsaEd9psEkhzudUcRblfFbFsx5YCRxbt+q58ylpNrAjOVY5NoszItZFxNPpy68Ch2Uc2pHA8ZJWk8ys/CpJK+q2KcK5nDLOApzLShwPpT/XAN8lmcG6Wtt/62VKDJcBb09b2A8HHo+Ih/MOqp6kPSr1oZJeTnKO12Ucg4BzgDsi4v802Sz389lKnAU5nwsl7ZQ+3xZ4NfDrus0uA96RPj8RuCbSlr+stBJnXd3y8STtOpmJiL+KiL0jYoikYfmaiHhb3Wa5n8tW4sz7XKYxzJe0feU58CdAfY/Jtv/W85hdtSFJ3yDpgbJA0oPAJ0kaz4iIfwG+T9K6fg/wFHBKQeM8EThV0ibgd8BJWX+pSa52TgZ+mdY3A/w1sKgqziKcz1biLML53BP4upKbTA0AF0fEFZL+BhiLiMtIEtwFku4hubo9KeMYW43zdEnHA5vSOBfnEOdWCnguGyrgudwd+G567TQbuDAirpL059D537qnxDAzsxplqkoyM7MMODGYmVkNJwYzM6vhxGBmZjWcGMzMrIYTg9k0SXpy6q3MysOJwczMajgxWN+Q9LJ0wrO56YjR2+rnEpJ0tqT3Vb3+lKQPStpO0o8l3aRk7vsTGrz/0aqat1/SlyUtTp8fJunadKKzqyujZiWdLun2NK6LevbhzdpQmJHPZr0WETdKugz4DLAtsCIi6qcPuIhkVs1/Sl+/mWS+of8G/mdEbEjnarpe0mWtjMJO54P6B+CEiJiQ9BZgGfBO4KPACyLi6cp0FmZ5c2KwfvM3wI0k/+hPr18ZETdL2k3SXsBC4LGIuD/95/7ZdObKzSTTFu8OPNLCMfcDDgJ+mE5dMAuozFVzKzAq6VJynInXrJoTg/WbXYDtSOa3mgv8tsE2l5DM0bQHSQkCYIQkURwWERvTWTfn1u23idrq2cp6AbdFxBENjvU6krsCHg98XNKB6T0IzHLjNgbrN8uBjwOjwNlNtrmIZOK2E0mSBCRTP69Jk8IfA4MN9hsHDlByz+IdgWPS5XcCCyUdAUnVkqQDJQ0A+0TET0huCLMTSdIyy5VLDNY3JL0d2BQRF6YzkP67pFdFxDXV20XEbelUxr+pmp54FLhcyc3Wb2HrabeJiAckXUxSPXQ3cHO6/BlJJwJfShPGbJJ2jLuAFekyAf83vY+CWa48u6qZmdVwVZKZmdVwYjAzsxpODGZmVsOJwczMajgxmJlZDScGMzOr4cRgZmY1/j86f4NvZrUpHgAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"#Create data set\n",
"u01 = [0, -0.5] + np.random.rand(50,2)\n",
"\n",
"data = [1, 5] + [5, 10] * u01\n",
"\n",
"#Plot the data\n",
"plt.scatter(data[:, 0], data[:, 1], marker='o', c='b')\n",
"plt.title('Example Linear Regression')\n",
"plt.xlabel('x values')\n",
"plt.ylabel('y values')\n",
"plt.xlim([1,5])\n",
"plt.ylim([0,10])\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Our goal is to find the equation of the straight line $l_\\alpha(x) = \\alpha_0 + \\mu_1 x$ that best fits our data points. The function that we are trying to minimize in this case is:\n",
"\n",
"$Loss(\\alpha_0,\\alpha_1) = {1 \\over 2m} \\sum\\limits_{i=1}^m (l_\\mu(x_i)-y_i)^2$\n",
"\n",
"In this case, our gradient will be defined in two dimensions:\n",
"\n",
"$\\frac{\\partial}{\\partial \\alpha_0} Loss(\\alpha_0,\\alpha_1) = \\frac{1}{m} \\sum\\limits_{i=1}^m (l_\\alpha(x_i)-y_i)$\n",
"\n",
"$\\frac{\\partial}{\\partial \\alpha_1} Loss(\\alpha_0,\\alpha_1) = \\frac{1}{m} \\sum\\limits_{i=1}^m ((l_\\alpha(x_i)-y_i) \\cdot x_i)$\n",
"\n",
"Below, we set up our function for l, Loss and the gradient:"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"# linear function\n",
"l = lambda alpha_0,alpha_1,x: alpha_0 + alpha_1*x\n",
"\n",
"# mean squared differences\n",
"def msqd(x,y,m,alpha_0,alpha_1):\n",
" returnValue = 0\n",
" for i in range(m):\n",
" returnValue += (l(alpha_0,alpha_1,x[i])-y[i])**2\n",
" returnValue = returnValue/(2*m)\n",
" return returnValue\n",
"\n",
"# gradients\n",
"def grad_msqd(x,y,m,alpha_0,alpha_1):\n",
" returnValue = np.array([0.,0.])\n",
" for i in range(m):\n",
" returnValue[0] += (l(alpha_0,alpha_1,x[i])-y[i])\n",
" returnValue[1] += (l(alpha_0,alpha_1,x[i])-y[i])*x[i]\n",
" returnValue = returnValue/(m)\n",
" return returnValue"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"# to speed up we take\n",
"x = data[:, 0]\n",
"y = data[:, 1]\n",
"m = len(x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We run our gradient descent algorithm (without adaptive step sizes in this example):"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Local minimum: alpha_0 = 4.601996694119921 alpha_1 = 0.17194658082243128\n",
"No of steps 10000\n"
]
}
],
"source": [
"alpha_new = np.array([4.0,0.2]) # The starting values\n",
"alpha = 0.001 # step size\n",
"precision = 0.0001 # stopping rule\n",
"n_steps = 10000 # stopping rule\n",
"\n",
"num_steps = 0 # counter for steps\n",
"ngrad = float(\"inf\") # stores the gradient\n",
"\n",
"while (np.linalg.norm(ngrad) > precision) and (num_steps < n_steps): # for stopping\n",
" num_steps += 1\n",
" alpha_old = alpha_new # store the current values\n",
" ngrad = -grad_msqd(x,y,m,alpha_old[0],alpha_old[1]) # calculate the gradient\n",
" alpha_new = alpha_old + alpha * ngrad # update due to gradient\n",
"\n",
"print(\"Local minimum: alpha_0 =\", alpha_new[0],\"alpha_1 =\", alpha_new[1])\n",
"print(\"No of steps\",num_steps)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The method is very slow indeed and depends on the chosen initial values. For comparison we may consider the values for the linear regression $\\alpha_0$ and $\\alpha_1$ that are calculated with the standard techniques..."
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Actual values for alpha are: alpha_0 = 4.822318363640957 alpha_1 = 0.11847746555299522\n"
]
}
],
"source": [
"actualvalues = sp.stats.linregress(x,y)\n",
"print(\"Actual values for alpha are:\", \"alpha_0 =\", actualvalues.intercept, \"alpha_1 =\", actualvalues.slope)"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEWCAYAAABi5jCmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xmc5Hdd5/HXu64+ZhJyzOSEmVFww4YYIJnFsLAYCa4sIPHANTw6QNB11gMQZTFoFglKRPACYcUdDQimOUJAjSyHIJeiRCYHRwiRqDO5mUlirum76rN//H41XVVT1V1dXVW/+nW/n49HPbrq96tf/T716+7v5/f7Xj9FBGZmZnWFrAMwM7PR4sRgZmZNnBjMzKyJE4OZmTVxYjAzsyZODGZm1sSJwY4i6b9IujXrODYCSTskPSqpmHUsWZB0uaSrso7D1saJYROTtF/Sc1qXR8TfRcQZWcTUKi1YFtPC9UFJ/yDp6VnH1a2IuD0itkZEtd+fLSkkHU6PTf3xK/3ej20+Tgw2MiSVOqz6UERsBbYBnwM+POT9j7Inp4mn/nhr1gFZ/jkx2FEknS/pzobX+yX9L0lfk/SQpA9JGm9Y/wJJNzWc0Z/dsO51kv5F0iOSvinpRxvWXSLpS5L+QNIDwOUrxRURS8A0cLqk7V3u/xxJN6b7/3Aa+5sav6ekSyXdC7yni8+7VNJd6efdKumCdPnTJO2T9LCk70j6/XT5rvTMvpS+Pk3StZIekHSbpJ9p+OzLJV0t6X3p598saXfXv7gGkj4u6fcaXn9I0rvT54+X9FlJ90u6T9K0pOMa3rtf0mvT3/dhSVdKOlnSJ9K4PiPp+Jbvt0fS3ZLukfSaFeI6Lz2mD0r6qqTze/l+NmAR4ccmfQD7gee0WX4+cGfL+/4JOA04AbgF+Nl03TnAQeD7gCLwsvT9Y+n6n0i3KwA/CRwGTk3XXQIsAa8ESsBEm1guB65Kn1eA3wbuA0qr7T99/wHgF4Ey8GPAAvCmhu+5BLwlff/EKp93BnAHcFq6/S7g8enzfwRekj7fCpzX8J5oiPcLwB8B48BTgEPABQ3fdQ54XrrvNwNfXuH3F8ATOqw7Jf0ezwamgH8FjknXPQH4wfQ7bQe+CLyt5ff9ZeBk4PT0c24Anppu81ngDS3f7wPAFuB70+/0nDa/v9OB+9PvV0hjuB/YnvX/gh8tfz9ZB+BHhr/8tSWGixtevxX44/T5u4DfbNn+VuD7O+zzJuDC9PklwO2rxHg5SWH+IFBNC5LzG9Z33D/wLOAuQA3r/p7mxLAAjHf5eU9IC8nnAOWW93wReCOwrWV5veAsAY9Lv8MxDevfDPxZw3f9TMO6M4HZFY5NAA+nx6b++KGG9T9GksjuA565wuf8CHBjy+97quH1R4B3Nbx+JfCXLd/viS1/H1c2fKd6YrgU+POWfX8KeFnW/wt+ND9clWTdurfh+QzJWTHATuA1adXAg5IeJCkATwOQ9NKGapkHgbNI2grq7uhi31dHxHEkZ7DfAM5tWLfS/k8D7oq0BOqwv0MRMdfN50XEbcCrSQq7g5I+KOm0dLufBv4D8C1JX5H0gjbf4zTggYh4pGHZAZIz6brW4zy+StvHORFxXMPjUw3rPkZy5XFrRPx9faGkk9LY75L0MHAVzb8TgO80PJ9t83pr89ubjusB0t9/i53AT7Qc22cCp67w/SwDTgy2XncAV7QUTpMR8QFJO4E/AV4BnJgW7t8A1LB919P7RsR9wP8ELpdUL0w67h+4h6Q9onF/j2v92G6/TxrD+yPimSSFXJBUQxER346IFwMnpcuukbSl5bPvBk6QdEzDsh0kVzWDcAVJtd+pkl7csPzNaexnR8SxwMU0/0560Xhcd5B811Z3kFwxNB7bLRHx2+vct/WZE4OVJY03PNbaM+dPgJ+V9H1KbJH0/LTw20JSAB0CkPRykiuGnkXEt0iqH+rdMlfa/z+SVN28QlJJ0oXA03r9PpLOkPRsSWMkbQGz6ecj6WJJ2yOiRlKlQ31dQ+x3AP8AvDk91meTXGlMr+eYtCPpWcDLgZemj3dIql+ZHAM8CjyYLnttH3b5ekmTkp6U7vdDbd5zFfDDkn5IUjE9BudLemwf9m995MRgHycp4OqPy9eycUTsA34GeCfw78BtJG0HRMQ3gd8jKaC/Q9Iw+aU+xPw7wB5JJ62y/wWSevafJimsLyapXpnv5fuQNLzWG7/vJbk6+LV03XOBmyU9CrwduKiliqruxST18ncDf0HSiPvpNX7/Rl9V8ziGt0k6Fngf8IqIuCutRroSeE969fRGkkb2h4D/B3x0Hfuv+wLJsfpb4Hcj4m9a35AmxgtJjtkhkiuI1+JyaOSoufrVbGOTdB1Jw/l7so5lI5C0C/g3ksb4pWyjsX5xprYNTdL3SzolrUp6GXA28Mms4zIbZQNLDJLeLemgpG80LDtB0qclfTv9efyg9m+WOgP4Kkm1yWuAF0XEPdmGZDbaBlaVlDZ+PQq8LyLOSpe9laS73m9Leh1wfERcOpAAzMysJwNtY0jrHz/WkBhuJRmcdE/a3fDzMSKTtZmZWWLYk4adXL+MT5PDSZ3eKGkPsAdgy5Yt5z7xiU8cUohmZhvD9ddff19EbF/9nc1GdjbJiNgL7AXYvXt37Nu3L+OIzMzyRdKBXrYbdq+k79RHrKY/D673A6enYdcuKBSSn9N9HypkZra5DDsxXEsyWyXpz79az4dNT8OePXDgAEQkP/fscXIwM1uPQXZX/QDJiNczlMx5/9Mko0Z/UNK3SabcXdccKZddBjMzzctmZpLlZmbWm4G1MaQTirVzQb/2cfvta1tuZmary/XI5x071rbczMxWl+vEcMUVMDnZvGxyMlluZma9yXVimJqCvXth506Qkp979ybLzcysN7lODJAkgf37oVZLfvYjKbgLrJltZiM7wC0r9S6w9d5O9S6w4CsRM9sccn/F0G/uAmtmm50TQwt3gTWzzc6JoYW7wJrZKBpm26cTQwt3gTWzUTPs6X+cGFq4C6yZjZpht30O9EY9/eJpt81sMysUkiuFVlLSVb8TSddHxO4172+tG5iZ9YPHC3Vv2G2fTgxmNnSeMn9tht326cRgZkPn8UJrM+y2T7cxmNnQ9VpnbmvjNgYzyw2PFxptTgxmNnQeLzTanBjMbOg8Xmi0eXZVM8vE1JQTwajyFYOZmTVxYjCzJh54Zq5KMrMjfKMqA18xmFmDvAw881XNYPmKwcyOyMONqnxVM3i+YjCzI/Iw8CwvVzV55sRgZkfkYeBZHq5q8s6JwcyOyMPAszxc1eSdE4OZNZmagv37k8ns9u8fraQA+biqyTsnBjPLlTxc1eSdeyWZWe54Oo3B8hWDmZk1cWIwM7MmTgxm1pFHGG9ObmMws7Y8wnjz8hWDmbXlEcabVyaJQdIvSbpZ0jckfUDSeBZxmFlnHmG8eQ09MUg6HXgVsDsizgKKwEXDjmMUTE/Dtm1JX2wpee46XBsVHmG8eWVVlVQCJiSVgEng7oziyMz0NPzUT8H99y8vu/9+ePnLnRxsNHiE8eY19MQQEXcBvwvcDtwDPBQRf9P6Pkl7JO2TtO/QoUPDDnPgLrsMFhaOXr646DpcGw0eYbx5KSKGu0PpeOAjwE8CDwIfBq6JiKs6bbN79+7Yt2/fkCIcjkIBOh16KZmnxsxsPSRdHxG717pdFlVJzwH+LSIORcQi8FHgP2cQR6ZWqqd1Ha6ZZSmLxHA7cJ6kSUkCLgBuySCOTF1xBVQqRy8vl12Ha2bZyqKN4TrgGuAG4OtpDHuHHUfWpqbg3e+GE09cXnbiifCe97gO18yylUmvpIh4Q0Q8MSLOioiXRMR8FnFkbWoK7rsvaWuISJ5vxqTgaResE/9tZMNTYlimPO2CdeK/jewMvVdSLzZiryRL7NqV/MO32rkzuXuYbV7+21i/PPVKMjvC0y5YJ/7byI4Tg2XK0y5YJ/7byI4Tg2XK0y5YJ/7byI4Tg2XK0y5YJ/7byI4bn83MNig3PpuZWV84MZjZpuUBdO05MayD/6jM8qs+gO7AgWTmgfoAOv8fOzH0zH9UZvnme1p35sTQI/9RmeWbB9B15sTQI/9RmeWbB9B15sSwBo1tCoUOR85/VGb54AF0nTkxdKm1TaFaPfo9/qMyyw8PoOvMA9y61Gmmx2IxuT/zjh1JUvAflZmNil4HuPl+DF3q1HZQqyUPM7ONwlVJXXJDlZltFk4MXXJDlZltFk4MXXJDlZltFm5jWIOpKScCM9v4fMVgZmZNnBjMzKyJE4OZ5YpnNR48tzGYWW7UZyCoT2BZn9UY3P7XT75iMLPc8KzGw+HEYGa54VmNh8OJwcxywzMQDIcTg5mNjNUalj0DwXA4MZjZSOjmdrmegWA4PO22mY2ETlPb79wJ+/cPO5qNoddpt33FYLngvusbnxuWR4cTg428bqoYLP/csDwI207oZSsnBht57ru+Obhhub+SE6cdO3vZNpPEIOk4SddI+pakWyQ9PYs4LB9cxbA5uGG5v5ITJ/VUxmd1xfB24JMR8UTgycAtGcVhOeAqhs1jaippaK7Vkp/9SgqbsY1qPSdOQ08Mko4FngVcCRARCxHx4LDjsPxwFYOtx2Zto1rPiVMWVwzfDRwC3iPpRkl/KmlL65sk7ZG0T9K+Q4cODT9KGxmuYrD12KxtVMmJU9R62Xbo4xgk7Qa+DDwjIq6T9Hbg4Yh4fadtPI7BzHpVKCRXCq2kpMpqI5O2/1vEoe9e63ZZXDHcCdwZEdelr68BzskgDjPbBDZ3G9V9D/Sy1ZoSg6RC2kbQs4i4F7hD0hnpoguAb67nM83MOnEb1dqtmhgkvV/SsWk7wDeBWyW9dp37fSUwLelrwFOA31rn55mZteU2qrVbtY1B0k0R8RRJU8C5wKXA9RFx9jACBLcxmJn1YpBzJZUllYEfAf4qIhaB0Z95b4PbjP2yzWw4urnn8/8F9gNfBb4oaSfw8CCDspX5vrdmNkg9dVeVVIqIpQHE05arkpp5emIz68bAqpIknSzpSkmfSF+fCbyshxitTzx3kJkNUjdtDH8GfAo4LX39z8CrBxWQrW5z98s2s0Hrpo1hW0RcLelXASJiSVJ1wHHZCq64ormNAdwv22yjiQjmFmvMLCwxs1BNHx2ezy8xs5j+bFjfq24Sw2FJJ5L2RJJ0HvBQz3u0das3MF92WVJ9tGNHkhTc8GxZmp7enH+T1Vows7DE7EKVw+0K7IVqS6HdXHjPzi+wuDDH0vwc1YVZaouz1BbmqC3NUYlFxrTIGA0PLRx5XmHpyOuTtcSW4hKThSUmC4tMaIk/6/E7dZMYfhm4Fni8pC8B24EX9bg/65Opqc3xT2f5MOo95SKChWrtSOE9u7DE4fmksJ6dm2V+bpaFuRnm52ZYmJ9lcX6WpYXZpLBenCMWZ6kuzBFLc8TiHFqaJ6rzFJbmKUa9oE5/ql5oL7KVRU5sKtgXGNcS40rWV1ikTJsz+wJQWcP3K46h0jiUxuDIz7Gej1dXvZIklYAzAAG3pmMZhsa9ksxGW8895WpVWJqHpTlYmqe2OMf8/AxzM7PMHymoZ1hMC+vq/CxLi3NUF+bSs+pZYjHZXkvzUJ2jUF1A1XmKtQVKtXlKsUCpttBwxp0UyPXXBa1vWFaNAtXiGNXCGFFMH6UKlMZRaRyVxymUxylWxrnjrgn+6YYxHnhojPLkOOc9c4yzntxaoDcU7Edet3lPsbL8s9C+ubjXXkmrXjFIemnLonMkERHvW+vOzGyEREB18UihTHW+oZCeayqwl38mj+pieladnllfunuO8afPMlaeZ7w0x3h5jvHSPOPFOe79nXkKtXkK1QWKtXmKtUXKsUApFo46Wy4AE+mjW/OUWaDMoiosqkJVFaqFCtXyGLVChVpxC1EaI0rjzJfGWCiNofIExco4xfRnaWyCcvqojE1Qqoyj8kSHAnscisuFdqFYogCUV4lzehr2vKGlbfDjozk9RzdTYryj4eU4yaR3N0TE0KqTfMVgG1KttlwIVxdWKIwbC+UV3tPyGfVqj9pi8jOWFlA1WV+oJg+tcxKDaog5KkcK5/koM0/9UTnyulqosFQYI4pJoU1pnEgLV5WTM+JieZxCZZxSeZzi2ATlSlpYj09SGZ+gMjbJxMQk4+ljYmKSUnm849lyv/Sr7SSL8UcDu2KIiFe27OgxwJ+vdUdmI6V+tlxtV+C2P0tuXyivpVBvWV5bf43sososUGFRaWEcyWM2ysxGKS2Yx5hny/L6tCCfbyjIFyizpApRSuuq0+qPQjkpsEuV5My6MjZBaWySsfGk0J6cGGeyUuQr/1jinW8rMvtokVgoEYtFxopF3vG2EpdcXEBSH35pw9fPtpM8jT/qpvG51QzwPf0OxDaZWi0tlFcqUBdoX1B3Pkte05n2Os+WQwVqxTFqhbH0jLhypDpjUZXkDJox5mMrc1FiNsrM1MrMUuJwocThKPFotchMtdh0dr1caNfPxEtNZ+CUxyiVJyhVxpgcqzBRKbJlrMhEucRkw/MtY8VkXaXERKXIZKXICenz+rItY0Umy8nzSqn3M+8LnwL/8dj0zPqejdMraaW7v631u+3Y0f6KYRTHH3XTxvDXLP8HFYAzgasHGZQNWATUljqfAa96ljy/hkK9w/o+nC0faXxrqP+N0thyYa0JFscew9L42JEz63oBOxel9FFmplZiplbicDV5PLJU5NFakYcXSzyyWODhpQIPLRQ5XGsooClTpdgxtIJoKHyXC+GJSonJcpHJsaSg3lIp8ZiWwnuy0r6An6yUmCgXKRZG8+x7I/aU6+dZfp7GH3VzxfC7Dc+XgAMRceeA4tkcVjxbbi1Q+3yWXC/Ue7sV7DIV2jfKlcaWG+YmT+zYq6JaqLBYqDTVS89FhdkoMVtLHjO1Eo9W08dSkUeqRR5ZLPDQYolHl8SjC7Xl/uCH0z7hi9W2t3HspFIqHCmgkzPpIhOTy4XzSZUiu44U1EnBPFkpMjnWWMCXku0aCvWxUn6rT2xZP8/y8zT+qJs2hi8MI5ChWfFsecBnyfXPqC6s/3u0OVtu6sJW2dJSMNd/tm7X+TOiNMYC5Yaz6jKHayUOVwscXhAzS7Uj/cFnF6scbhnAM7tQ5fDMUkPf8SqH03ULS90mpgAWmazUmKxEekbNkTPrbVvHWgrqUrpuuZBu+3ysyGS5SKmYxd1tLS/6fZafl6uqjolB0iO0r4QVEBGxrlt8rsnMffDld63/LLn+er1nywhau7IVx5pftz1bHmtTULd+Rqc+zQ2fURxr6olRq0VSMNdHX85XmV1MB/AsLD+vF8yzc8sFdOuAn3oBP7vwMDOLVaq17k+/SwU1F8Jp/fXxWyqcfnxj4bxKgV1Jtqs/Hy8VKYxo9YltbHk6y++nnqbdHrbdpxVj356tywtaz5Y7nT13PEtuc7a86mc0FOqFUnKPwDVarNaYma8y01BQNw6PP1KwtxTwMy2FfeN7koJ/bVNXjZcLR+qr6/XeW1oK6omW6pWmQrttwV5aV+OlmfXfwLqrNuzgJJJxDABExPA6WZ1yFlz6+eWBJQPst9x24qrD9eczzCw8nBbUS01n241znzRWm8w2LF+sdp+ECyIpvOv13mnhfcx4iVOOHV8uvMe6K+Dr1S2j3HhpZqOhm15JLwR+j2Ta7YPATuAW4EmDDa1BoQwTxzctqk9c1XGmwZVmImw8S59PGiwbZyhcU+NlsdBQeC8X1NuPGWNHZZLJckOvlIYCvt5Q2brdljE3XppZtrq5YvhN4DzgMxHxVEk/ALx4sGE1u+3go/zg73+hqVCf77rxMrF8Vt1YRVLihC2V5i6C9cJ7rJj2QFluqFwu4Je7FpbdeGlmG0w3iWExIu6XVJBUiIjPSXrLwCNrUCqIJ5y09ahGzeX+3Uc3bDYO4HHjpZlZ97pJDA9K2gp8EZiWdBDazRM7OLu2beFdF587zF2amW1a3dSDXEgyDcYvAZ8E/gX44UEGZWZm2enmimEP8OF0tPN7BxyPmZllrJsrhmOBT0n6O0m/IOnkQQdlZmbZWTUxRMQbI+JJwC+QdFn9gqTPDDwyMzPLxFr6Wh4E7gXuB04aTDhmZpa1VRODpJ+T9Hngb4FtwM9ExNmDDszMzLLRTePzTuDVEXHToIMxM7PsdTPt9uuGEYiZmY0Gz+dgZmZNnBjMzKxJN43Pr5B0/Grvs41rehp27UpmO9+1K3ltZhtXN1cMpwBfkXS1pOfKc0FvKtPTya0NDxxI7op64EDy2snBbOPqZoDb/wa+B7gSuAT4tqTfkvT49exYUlHSjZI+tp7PscG67LLm+91C8vqyy7KJx8wGr6s2hkju/3lv+lgCjgeukfTWdez7F0lu+GMj7PYO9+nrtNzM8q+bNoZXSboeeCvwJeB7I+LngHOBH+9lp5IeCzwf+NNetrfh2bFjbcvNLP+6uWLYBvxYRPxQRHw4IhYBIqIGvKDH/b4N+BWg423YJO2RtE/SvkOHDvW4G1uvK66AycnmZZOTyXIz25i6aWP49Yg40GHdmquCJL0AOBgR16+y370RsTsidm/fvn2tu7E+mZqCvXth506Qkp979ybLzWxj6mZKjH57BvBCSc8DxoFjJV0VERdnEIt1YWrKicBsMxn6ALeI+NWIeGxE7AIuAj7rpGBmNjo88tlsk/BARetWFlVJR0TE54HPZxmD2WZQH6hYH5NSH6gIria0o/mKYcB8lmajwAMVbS0yvWLY6HyWZqPCAxVtLXzFMEA+S7NR4YGKthZODAPks7SNLU/VhB6oaGvhxDBAPkvbuPI266wHKtpaKJkfb7Tt3r079u3bl3UYa9baxgDJWZr/IfNv164kGbTauRP27x92NGbtSbo+InavdTtfMQyQz9I2LlcT2kbmXkkD5ukkNqYdO9pfMbia0DYCXzGY9cCNubaROTGY9cDVhLaRuSrJrEeuJrSNylcMtinlaQyC2bD5isE2HU9VYrYyXzHYpuOpSsxW5sRgm47HIJitzInBNh1PVWK2MicG23Q8BsFsZU4MNlKG0VvIYxDMVuZeSTYyhtlbyGMQzDrzFYONDPcWMhsNTgw2MtxbyKzftp3Qy1ZODDYy3FvIrN9OO72XrZwYbGS4t5BZv5UrvWzlxGAjw72FzPptcaGXrdwryUaKewuZ9dPdd/Wyla8YzMw2rPse6GUrJwYzM2vixDAgnu/fzPLKbQwD4Pn+zSzPfMUwAB7Ba2Z55sQwAB7Ba2Z55sTQZ9PTSbtCOx7Ba2Z54MTQR/W2hWr16HUewWtmeeHE0Eft2hYAikWP4M2bPPQqy0OMlk9D75Uk6XHA+4BTgBqwNyLePuw4BqFTG0Kt5qSQJ3noVZaHGC2/FBHD3aF0KnBqRNwg6RjgeuBHIuKbnbbZvXt37Nu3b2gx9mrXruQftNXOnbB//7CjsV7l4feYhxgte5Kuj4jda91u6FVJEXFPRNyQPn8EuAXoaWrYUePZQTeGPPQqy0OMll+ZtjFI2gU8Fbiuzbo9kvZJ2nfo0KFhh9YTzw66MeThvhB5iNHyK7PEIGkr8BHg1RHxcOv6iNgbEbsjYvf27duHH2CPpqaSS/laLfnppJA/ebjyy0OMll+ZJAZJZZKkMB0RH80iBrNO8nDll4cYLb+yaHwW8F7ggYh4dTfb5KXx2cxslOSm8Rl4BvAS4NmSbkofz8sgDjMza2Po4xgi4u8BDXu/ZmbWHY98NjOzJk4MZjnmaTFsEHyjHrOc8rQYNii+YjDLKd8QygbFicEspzwthg2KE4NZTnlaDBsUJ4YccAOjteNpMWxQnBhGXL2B8cABiFhuYHRyME+LYYMy9CkxerGZp8TwvPtm1qs8TYlha+AGRjMbNieGEecGRjMbNieGEecGRjMbNieGEecGRjMbNk+JkQNTU04EZjY8vmIwM7MmTgxmZtbEicHMzJo4MZiZWRMnBjMza+LEYGYjxZNGZs/dVc1sZPiudKPBVwxmNjJ8V7rR4MRgZiPDk0aOBicGMxsZnjRyNDgxmFnf9dqA7EkjR4MTg5n11XruOuhJI0eD7+BmZn3luw6ODt/BzcxGghuQ88+Jwcz6yg3I+efEYGZ95Qbk/HNiMLO+cgNy/nlKDDPrO991MN98xWBmZk2cGMzwjJ5mjZwYbNNbz4CsjchJ0jJJDJKeK+lWSbdJel0WMZjVeUbPZU6SBhkkBklF4P8A/w04E3ixpDOHHYdZnQdkLXOSNMjmiuFpwG0R8a8RsQB8ELgwgzjMAA/IauQkaZBNd9XTgTsaXt8JfF/rmyTtAdJ7NzEv6RtDiG29tgH3ZR1EF/IQ5xBj3HYC7NgJajhRitqBA7cfkO57YLWNGf1jCV3Hefb3QrnSujRicUH62tcHEFerDXY8M3dGLxtlkRjUZtlRM/lFxF5gL4Ckfb1MBDVsjrN/8hAjOM5+c5z9Jamn2UezqEq6E3hcw+vHAndnEIeZmbWRRWL4CvA9kr5LUgW4CLg2gzjMzKyNoVclRcSSpFcAnwKKwLsj4uZVNts7+Mj6wnH2Tx5iBMfZb46zv3qKMxc36jEzs+HxyGczM2vixGBmZk1GJjFIerekg53GKyjxh+k0Gl+TdM6wY0zjWC3O8yU9JOmm9PHrGcT4OEmfk3SLpJsl/WKb92R+PLuMcxSO57ikf5L01TTON7Z5z5ikD6XH8zpJu0Y0zkskHWo4nv9j2HGmcRQl3SjpY23WZX4sG2JZKc5ROZb7JX09jeGo7qk9/a9HxEg8gGcB5wDf6LD+ecAnSMZBnAdcN6Jxng98LONjeSpwTvr8GOCfgTNH7Xh2GecoHE8BW9PnZeA64LyW9/w88Mfp84uAD41onJcA78zyeKZx/DLw/na/21E4ll3GOSrHcj+wbYX1a/5fH5krhoj4IrDSKNMLgfdF4svAcZJOHU50y7qIM3MRcU9E3JA+fwS4hWTEeaPMj2eXcWYuPUaPpi/L6aO118aFwHvT59cAF0hqN5hzYLqMM3OSHgs8H/jTDm/J/FhCV3HmxZr/10cmMXSh3VQaI1eIpJ6eXs5/QtKTsgwkvQx/KsnZY6OROp4rxAkjcDzTKoWbgIPApyOi4/EW1mAyAAAEbklEQVSMiCXgIeDE4UbZVZwAP55WKVwj6XFt1g/a24BfAWod1o/EsWT1OCH7YwlJ8v8bSdcrmUqo1Zr/1/OUGLqaSmME3ADsjIgnA+8A/jKrQCRtBT4CvDoiHm5d3WaTTI7nKnGOxPGMiGpEPIVkpP7TJJ3V8paROJ5dxPnXwK6IOBv4DMtn5kMh6QXAwYi4fqW3tVk21GPZZZyZHssGz4iIc0hmrP4FSc9qWb/m45mnxJCLqTQi4uH65XxEfBwoS9o27DgklUkK2+mI+Gibt4zE8VwtzlE5ng3xPAh8Hnhuy6ojx1NSCXgMGVY5doozIu6PiPn05Z8A5w45tGcAL5S0n2Rm5WdLuqrlPaNwLFeNcwSOZT2Ou9OfB4G/IJnButGa/9fzlBiuBV6atrCfBzwUEfdkHVQrSafU60MlPY3kGN8/5BgEXAncEhG/3+FtmR/PbuIckeO5XdJx6fMJ4DnAt1redi3wsvT5i4DPRtryNyzdxNlSt/xCknadoYmIX42Ix0bELpKG5c9GxMUtb8v8WHYTZ9bHMo1hi6Rj6s+B/wq09phc8/96FrOrtiXpAyQ9ULZJuhN4A0njGRHxx8DHSVrXbwNmgJePaJwvAn5O0hIwC1w07D9qkrOdlwBfT+ubAX4N2NEQ5ygcz27iHIXjeSrwXiU3mSoAV0fExyT9BrAvIq4lSXB/Luk2krPbi4YcY7dxvkrSC4GlNM5LMojzKCN4LNsawWN5MvAX6blTCXh/RHxS0s9C7//rnhLDzMya5KkqyczMhsCJwczMmjgxmJlZEycGMzNr4sRgZmZNnBjM1knSo6u/yyw/nBjMzKyJE4NtGpL+Uzrh2Xg6YvTm1rmEJL1F0s83vL5c0mskbZX0t5JuUDL3/YVtPv98NczbL+mdki5Jn58r6QvpRGefqo+alfQqSd9M4/rgwL682RqMzMhns0GLiK9IuhZ4EzABXBURrdMHfJBkVs0/Sl//d5L5huaAH42Ih9O5mr4s6dpuRmGn80G9A7gwIg5J+kngCuCngNcB3xUR8/XpLMyy5sRgm81vAF8hKehf1boyIm6UdJKk04DtwL9HxO1p4f5b6cyVNZJpi08G7u1in2cAZwGfTqcuKAL1uWq+BkxL+ksynInXrJETg202JwBbSea3GgcOt3nPNSRzNJ1CcgUBMEWSKM6NiMV01s3xlu2WaK6era8XcHNEPL3Nvp5PclfAFwKvl/Sk9B4EZplxG4NtNnuB1wPTwFs6vOeDJBO3vYgkSUAy9fPBNCn8ALCzzXYHgDOV3LP4McAF6fJbge2Sng5J1ZKkJ0kqAI+LiM+R3BDmOJKkZZYpXzHYpiHppcBSRLw/nYH0HyQ9OyI+2/i+iLg5ncr4robpiaeBv1Zys/WbOHrabSLiDklXk1QPfRu4MV2+IOlFwB+mCaNE0o7xz8BV6TIBf5DeR8EsU55d1czMmrgqyczMmjgxmJlZEycGMzNr4sRgZmZNnBjMzKyJE4OZmTVxYjAzsyb/H7CNGgyT6/C7AAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"xx = np.linspace(0,5,1000)\n",
"plt.scatter(data[:, 0], data[:, 1], marker='o', c='b')\n",
"plt.plot(xx,l(alpha_new[0],alpha_new[1],xx))\n",
"plt.plot(xx,l(actualvalues.intercept,actualvalues.slope,xx))\n",
"plt.xlim([1,5])\n",
"plt.ylim([0,10])\n",
"plt.title('Linear Regression Example')\n",
"plt.xlabel('x values')\n",
"plt.ylabel('y values')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We calculate the gradient for every step running through the algorithm. This may not be a big issue when the dimensionality is small but think of real life situations with a large dimensionality...\n",
"\n",
"In machine learning, the algorithm above is often called batch gradient descent to contrast it with mini-batch gradient descent (which we will not go into here) and stochastic gradient descent."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Stochastic gradient descent"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As we said above, in batch gradient descent, we must look at every example in the entire training set on every step (in cases where a training set is used for gradient descent). This can be quite slow if the training set is sufficiently large. In stochastic gradient descent, we update our values after looking at each item in the training set, so that we can start making progress right away. Recall the linear regression example above. In that example, we calculated the gradient for each of the two alpha values as follows:\n",
"\n",
"$\\frac{\\partial}{\\partial \\alpha_0} Loss(\\alpha_0,\\alpha_1) = \\frac{1}{m} \\sum\\limits_{i=1}^m (l_\\alpha(x_i)-y_i)$\n",
"\n",
"$\\frac{\\partial}{\\partial \\theta_1} Loss(\\alpha_0,\\alpha_1) = \\frac{1}{m} \\sum\\limits_{i=1}^m ((l_\\alpha(x_i)-y_i) \\cdot x_i)$\n",
"\n",
"Where $l_\\alpha(x) = \\alpha_0 + \\alpha_1 x$\n",
"\n",
"Then we followed this algorithm (where $\\alpha$ was a non-adapting stepsize):\n",
"\n",
" 1: Choose initial guess $x_0$ \n",
" 2: for k = 0, 1, 2, ... do \n",
" 3: $grad_k$ = -$grad f(x_k)$ \n",
" 4: $x_{k+1} = x_k + \\alpha grad_k$ \n",
" 5: end for\n",
"\n",
"When the sample data had 15 data points as in the example above, calculating the gradient was not very costly. But for very large data sets, this would not be the case. So instead, we consider a stochastic gradient descent algorithm for simple linear regression such as the following, where m is the size of the data set:\n",
"\n",
" 1: Randomly shuffle the data set \n",
" 2: for k = 0, 1, 2, ... do \n",
" 3: for i = 1 to m do \n",
" 4: $\\begin{bmatrix}\n",
" \\alpha_{1} \\\\ \n",
" \\alpha_2 \\\\ \n",
" \\end{bmatrix}=\\begin{bmatrix}\n",
" \\alpha_1 \\\\ \n",
" \\alpha_2 \\\\ \n",
" \\end{bmatrix}-\\alpha\\begin{bmatrix}\n",
" 2(l_\\alpha(x_i)-y_i) \\\\ \n",
" 2x_i(l_\\alpha(x_i)-y_i) \\\\ \n",
" \\end{bmatrix}$ \n",
" 5: end for \n",
" 6: end for\n",
"\n",
"With stochastic gradient descent you run through the entire data set 1 to 10 times - see value for k in the pseudocode. But it depends on how fast the data is converging and it may be chosen with regard to the number of data in the set.\n",
"\n",
"With BGD we must run through the entire data set before we make any progress. With SGD we make progress immediately and continue to make progress by stepping through the entire set. Therefore, SGD is often preferred when dealing with large data sets. However, to profit from the both (GD and SGD) a hybrid method called Minibatch Gradient Descent is often used!\n",
"\n",
"Unlike GD, SGD tends to oscillate close to a minimum value rather than continuously approxing it and may never converge to the minimum. One way around this is to slowly decrease the step size $\\alpha$ as the algorithm runs. However, this is less common than using a fixed $\\alpha$.\n",
"\n",
"Example for SGD for linear regression with 500.000 points around the line $y = 2x+17+\\epsilon$, for values of x between 0 and 100 is considered."
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [],
"source": [
"f = lambda x: x*1.35+23+np.random.randn(len(x))*10\n",
"\n",
"x = np.random.random(500000)*100\n",
"y = f(x) \n",
"m = len(y)"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [],
"source": [
"xxnew = np.arange(0,len(x))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We randomly shuffle around the dataset. \n",
"Hint: In this example this step isn't strictly necessary since the data is already in a random order but for real life data this may not always be the case!"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [],
"source": [
"from random import shuffle\n",
"\n",
"x_shuf = []\n",
"y_shuf = []\n",
"index_shuf = xxnew\n",
"shuffle(index_shuf)\n",
"for i in index_shuf:\n",
" x_shuf.append(x[i])\n",
" y_shuf.append(y[i])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we'll setup our h function and our cost function, which we will use to check how the value is improving."
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [],
"source": [
"l = lambda alpha_0,alpha_1,x: alpha_0 + alpha_1*x\n",
"Loss = lambda alpha_0,alpha_1, x_i, y_i: 0.5*(l(alpha_0,alpha_1,x_i)-y_i)**2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we'll run our stochastic gradient descent algorithm. To see it's progress, we'll take a cost measurement at every step. Every 10,000 steps, we'll get an average cost from the last 10,000 steps and then append that to our cost_list variable. We will run through the entire list 10 times here:"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Local minimum: alpha_0 = 22.961169591278534 alpha_1 = 1.34265986023661\n"
]
}
],
"source": [
"alpha_old = np.array([0.,0.])\n",
"alpha_new = np.array([1.,1.]) # The algorithm starts at [1,1]\n",
"n_k = 0.000005 # step size\n",
"\n",
"iter_num = 0\n",
"grad_k = np.array([float(\"inf\"),float(\"inf\")])\n",
"sum_loss = 0\n",
"loss_list = []\n",
"\n",
"for j in range(10):\n",
" for i in range(m):\n",
" iter_num += 1\n",
" alpha_old = alpha_new\n",
" grad_k[0] = (l(alpha_old[0],alpha_old[1],x[i])-y[i])\n",
" grad_k[1] = (l(alpha_old[0],alpha_old[1],x[i])-y[i])*x[i]\n",
" grad_k = (-1)*grad_k\n",
" alpha_new = alpha_old + n_k * grad_k\n",
" sum_loss += Loss(alpha_old[0],alpha_old[1],x[i],y[i])\n",
" if (i+1) % 10000 == 0:\n",
" loss_list.append(sum_loss/10000.0)\n",
" sum_loss = 0 \n",
" \n",
"print(\"Local minimum:\", \"alpha_0 =\", alpha_new[0], \"alpha_1 =\", alpha_new[1])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As you can see, our values for $\\alpha_0$ and $\\alpha_1$ are close to their true values of 23 and 1.35.\n",
"\n",
"Now, we plot our cost versus the number of iterations. As you can see, the cost goes down quickly at first, but starts to level off as we go through more iterations:"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEKCAYAAAAB0GKPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8lNW9x/HPL/tCCCGEHdkFFQUxLogLgnW3Lr1Wvd5q1dbqtS61VbSLWmuvWm3r0tat1qVVqnUXFQRUBEQgLCJ7WEPYEiAEkpD93D/myWRCEpjETCbJfN+vV17zzJln5vmdYZjfPOec5xxzziEiInIoUeEOQERE2gclDBERCYoShoiIBEUJQ0REgqKEISIiQVHCEBGRoChhiIhIUJQwREQkKEoYIiISlJhwB/BtdOvWzQ0YMCDcYYiItCsLFy7c6ZzLaOrz2nXCGDBgAFlZWeEOQ0SkXTGzTc15npqkREQkKEoYIiISFCUMEREJihKGiIgERQlDRESCooQhIiJBUcIQEZGgRGTCWL19H3+YsorCkopwhyIi0m5EZMLYtKuYv32+jpzdJeEORUSk3QhZwjCzf5hZnpktCyi7zMyWm1m1mWUesP89ZrbWzFab2dmhigugV2oiAFsL94fyMCIiHUoozzBeAs45oGwZcCnwRWChmR0JXAEc5T3nb2YWHarAenVJAGDbHiUMEZFghSxhOOe+AHYfULbSObe6gd0vAv7tnCtzzm0A1gInhCq29OQ44mKi2FZYGqpDiIh0OG2lD6MPsDngfq5XVo+Z3WBmWWaWlZ+f36yDmRm9UhPYqoQhIhK0tpIwrIEy19COzrnnnHOZzrnMjIwmz87r1ys1ge3qwxARCVpbSRi5QL+A+32BraE8YFpSHAUaVisiErS2kjDeB64ws3gzGwgMBeaH8oCpibEU7lfCEBEJVsgWUDKzScA4oJuZ5QL34esEfwrIAD40syXOubOdc8vN7A1gBVAJ3OycqwpVbKCEISLSVCFLGM65Kxt56J1G9v898PtQxXOgzomxlFdWU1pRRUJsyEbwioh0GG2lSarVpSbGAugsQ0QkSEoYShgiIkFRwlDCEBEJihKGhtaKiARFCUNnGCIiQYnYhJGS4BsgVlRWGeZIRETah4hNGJ2UMEREmiRiE0Z8TDSx0aaEISISpIhNGADJ8TEUK2GIiAQlohNGp/gYikqVMEREgqGEoTMMEZGgRHzCKC5XwhARCUZEJ4zk+BjmrN3FaX/4jLve/Drc4YiItGkRnTDiYnzVz9ldwhtZuWGORkSkbYvohJG3ryzcIYiItBsRnTC27am7pndZZUjXbBIRadciOmGMGZxe5/4eTUQoItKoiE4Yj3zvGGbddQb/d8nRABSUlIc5IhGRtiuiE0ZCbDT9uiYxID0JgIJinWGIiDQmohNGjS5JcQDs0RmGiEijlDCAtGTf2hgF6sMQEWmUEgZaTElEJBhKGEBibDTRUca+UiUMEZHGhCxhmNk/zCzPzJYFlHU1s2lmlu3dpnnlZmZPmtlaM1tqZqNDFVcjsZKSEMM+zVwrItKoUJ5hvAScc0DZ3cAM59xQYIZ3H+BcYKj3dwPwdAjjapAvYegMQ0SkMSFLGM65L4DdBxRfBLzsbb8MXBxQ/orz+QroYma9QhVbQ1LiY3WGISJyEK3dh9HDObcNwLvt7pX3ATYH7JfrldVjZjeYWZaZZeXn57dYYGqSEhE5uLbS6W0NlLmGdnTOPeecy3TOZWZkZLRYACkJsexVk5SISKNaO2HsqGlq8m7zvPJcoF/Afn2Bra0ZmM4wREQOrrUTxvvANd72NcB7AeVXe6OlTgIKa5quWktKgpZrFRE5mJhQvbCZTQLGAd3MLBe4D3gYeMPMrgdygMu83T8CzgPWAiXAtaGKqzE1CcM5h1lDLWQiIpEtZAnDOXdlIw9NaGBfB9wcqliCkZoYS1W1Y19ZJZ0TYsMZiohIm9RWOr3DrldqIgBbD1hUSUREfJQwPH3TfAkjd7cShohIQ5QwPP26+tbE2FxQEuZIRETaJiUMT3pyHImx0eQW6AxDRKQhShgeM6N3lwS2KGGIiDRICSNA95QEdhaVhTsMEZE2SQkjQLeUeLI2FfDQRyvDHYqISJujhBEgo1M8AM9+sR7fpSEiIlJDCSNARkq8f1vre4uI1KWEESAlofbC9827NbxWRCSQEkaA6oBmKA2vFRGpSwkjwPdG9+WSY33rNuUWlGj2WhGRAEoYAZLjY/jz5aNITYxl2oodjLhvKpOXtuqyHCIibZYSRgP6piWStakAgC/WtNwysCIi7ZkSRgP6pSX5t2Oj9RaJiIASRoNqZq4FJQwRkRr6NmxAYMKIi9FbJCICShgN6hvQJCUiIj5KGA2oWRsDoKRcQ2tFREAJo0F9ApqkSsqrwhiJiEjboYTRgE7xMaQlxQKQt7eMyqrqMEckIhJ+ShiNGN6zMwCz1+7krreWhjkaEZHwU8JoxLNXH0fv1AQA3l60JczRiIiEX1gShpndZmbLzGy5md3ulXU1s2lmlu3dpoUjthqdE2IpqfD1XwTOYisiEqlaPWGY2Qjgx8AJwEjgAjMbCtwNzHDODQVmePfDao+3JkYXrz9DRCSSheMM4wjgK+dciXOuEpgJXAJcBLzs7fMycHEYYmtQSrwShohIOBLGMuA0M0s3syTgPKAf0MM5tw3Au+3e0JPN7AYzyzKzrPz80E4M+NK1xwNQWqGhtSIirZ4wnHMrgUeAacAU4Gsg6KvjnHPPOecynXOZGRkZIYrSZ9yw7vz3iYext1TLtYqIhKXT2zn3gnNutHPuNGA3kA3sMLNeAN5tXjhiO1BqYix791fiAlbjExGJROEaJdXduz0MuBSYBLwPXOPtcg3wXjhiO1DnhFjKq6opq9TFeyIS2cI1XvQtM0sHKoCbnXMFZvYw8IaZXQ/kAJeFKbY6Oif63qK1eUWM6JMa5mhERMInLAnDOXdqA2W7gAlhCOeghvVIAeCxT1bz0rUnhDkaEZHw0ZXeh5A5oCtnHdmDTbtKwh2KiEhYKWEEoXeXRHbuKwt3GCIiYaWEEYSMlHj2lVXqegwRiWhKGEHo1ikOgJ1FOssQkcilhBGEbp3iAdhZVB7mSEREwueQCcPMZgRT1pFlpHgJQ/0YIhLBGh1Wa2YJQBLQzZtq3LyHOgO9WyG2NiPdO8PYVayEISKR62DXYfwEuB1fclhIbcLYC/w1xHG1KTXLtRaUaE4pEYlcjSYM59wTwBNmdotz7qlWjKnNSYyNJi4mioIS9WGISOQKptN7u5mlAJjZr83sbTMbHeK42hQzIy0ploJiJQwRiVzBJIzfOOf2mdkpwNn4Fjd6OrRhtT1pSXFqkhKRiBZMwqi5Wu184Gnn3HtAXOhCapu6JMWyR01SIhLBgkkYW8zsWeD7wEdmFh/k8zoUnWGISKQL5ov/+8BU4Bzn3B6gK3BnSKNqgzrFx7A2r4i563aFOxQRkbA4ZMJwzpUA64CzzeynQHfn3Cchj6yNOdyb5vypT7Mp12JKIhKBgrnS+zbgVaC79/cvM7sl1IG1NT8cO4DTD8/gy3W7OO/JWeEOR0Sk1QXTJHU9cKJz7l7n3L3AScCPQxtW2xMbHcWwnr6zjLV5RWwr3B/miEREWlcwCcOoHSmFt22N7NuhXTq6j397VvbOMEYiItL6gkkYLwLzzOx+M7sf+Ap4IaRRtVHDe3ZmzYPnArBtT2mYoxERaV2HXNPbOfcnM/scOAXfmcW1zrnFoQ6srYqLiSI1MVYTEYpIxAmm0/skINs596Q3v9RaMzsx9KG1Xemd4nhl7iYW5xSEOxQRkVYTTJPU00BRwP1iInBqkEDpyb4L3S/525fc994ylm0pDHNEIiKhF1Snt3PO1dxxzlUTRFNWR1ZV7X87eHnuJn70clYYoxERaR3BJIz1ZnarmcV6f7cB67/NQc3sZ2a23MyWmdkkM0sws4FmNs/Mss3sdTNrs/NV7T5g1trYmIgcNCYiESaYhHEjcDKwBcgFTgRuaO4BzawPcCuQ6ZwbAUQDVwCPAH92zg0FCvBd/9Em1azxXSMpNqJPuEQkQgQzNUiec+4K51x351wP59x/O+fyvuVxY4BEM4vBtwzsNmA88Kb3+MvAxd/yGCHz16tGc+xhXfz3E+OiwxiNiEjraPVZZ51zW4DHgBx8iaIQ3xKwe5xzld5uuUCfhp5vZjeYWZaZZeXn57dGyPX06JzAj08d5L+fHK+EISIdX6snDDNLAy4CBuJbLzwZOLeBXV0DZTjnnnPOZTrnMjMyMkIX6CEc3qOTfzvK1IchIh1fONa1OBPY4JzLd85VAG/j6yPp4jVRAfQFtoYhtqAN6Z7i356VvZNrX5wfxmhERELvkL21ZnZHA8WFwELn3JJmHDMHOMnMkoD9wAQgC/gM+C/g38A1wHvNeO1WNfPOcVz2zFzy9pXx2erwNI+JiLSWYM4wMvGNlOrj/d0AjAOeN7O7mnpA59w8fJ3bi4BvvBieAyYCd5jZWiCddjBfVf/0ZE4alO6/v2WPZrAVkY4rmISRDox2zv3cOfdzfAkkAzgN+GFzDuqcu885N9w5N8I59wPnXJlzbr1z7gTn3BDn3GXOuXYxWVNpRe1EvmMf/jSMkYiIhFYwCeMwIPBKtQqgv3NuP9AuvtRDqbi88tA7iYh0AMFccfYa8JWZ1fQpXAhMMrNkYEXIImsnisqq6twvr6wmLiYcYwlEREIrmOnNf2dmH1E7vfmNzrmayZOuCmVw7UF1dd3Rv0VllXSNabOzmoiINFsw05s/AcQ7555wzj0ekCwEePyKUaTE1+bdolI1UYlIxxRM28ki4NdmttbMHjWzzFAH1Z4MzujEy9ef4L+/t7QijNGIiIROMHNJveycOw84AVgDPGJm2SGPrB0ZfVgaf758JAAXPDWbrRpeKyIdUFN6Z4cAw4EBwKqQRNOODc6onSpk2oodYYxERCQ0gunDqDmjeABYDhznnLsw5JG1MykJsf7tr9bvCmMkIiKhEcyw2g3AGOfczlAH0551Cuj43rSrJIyRiIiERjDDap8xszQzOwFICCj/IqSRtTMpCbVvZX5RxF/PKCIdUDCTD/4IuA3fDLJLgJOAufgWPBJPfMDFeruKyqiqdkRHadpzEek4gun0vg04HtjknDsDOBbQ1KwHMDM2Pnw+D1x0FNWu/rrfIiLtXTAJo9Q5VwpgZvHOuVXAsNCG1X5leOt95+0rDXMkIiItK5hO71wz6wK8C0wzswLa+OJG4dS9sy9h5O9TP4aIdCzBdHpf4m3eb2afAanAlJBG1Y5ldPKNC1DCEJGOJpgzDD/n3MxQBdJRZKT4zjA27iqmoqqa2GjNXCsiHYO+zVpYYlw0neJj+Otn67jupQXsKVHnt4h0DEoYIVAznHZW9k6Oe3B6mKMREWkZShghUFRWO8V5VbWjvLI6jNGIiLQMJYwQqDpgUaXsvH1hikREpOUoYYTAEb06A3D2UT0AWL51bzjDERFpEU0aJSXB+df1J7C5YD/H9Enl6PunsnxLIWT2C3dYIiLfis4wQiC9Uzyj+nUhKsoY1jOFl+du4p9fbQp3WCIi30qrJwwzG2ZmSwL+9prZ7WbW1cymmVm2d5vW2rGFQs9U34V8j07RmlMi0r61esJwzq12zo1yzo0CjgNKgHeAu4EZzrmhwAzvfrv3y/OOACA1KfYQe4qItG3hbpKaAKxzzm0CLgJe9spfBi4OW1QtqG9aEreOH8KWgv0aXisi7Vq4E8YVwCRvu4dzbhuAd9u9oSeY2Q1mlmVmWfn57WOW9QHdkql2kLO7ONyhiIg0W9gShpnFAd8F/tOU5znnnnPOZTrnMjMyMkITXAsb1jMFgGVbNLxWRNqvcJ5hnAsscs7t8O7vMLNeAN5tXtgia2HDe3amU3wMWZt2hzsUEZFmC2fCuJLa5iiA94FrvO1rgPdaPaIQiY4yRvXrwpLNe8IdiohIs4UlYZhZEvAd4O2A4oeB75hZtvfYw+GILVT6dU0kb6/WyBCR9issCcM5V+KcS3fOFQaU7XLOTXDODfVuO1T7TZekOPL2lTHg7g/Zumd/uMMREWmycI+Sihhdk+L82yc//CmfrtpxkL1FRNoeJYxWkpYcV+f+vA0d6gRKRCKAEkYr6Zpc/0rviqpqHp++hsL9FWGISESkaZQwWklaUt0zjC0F+5m6fDuPT8/mj5+sDlNUIiLBU8JoJQcmjNyC/f4zi+KyqnCEJCLSJFoPo5XU9GF8d2RvYqOjmLkmzz/MNlppW0TaAX1VtZLUxFg+/fnpPHbZSAZlJLOzqJxV231ThewqKg9zdCIih6YzjFY0KKMTAH3TEgGYutw3tHZbYWnYYhIRCZbOMMKgb1qSf3tkvy5s36uEISJtnxJGGPTrmujfPnN4d3YXl1NaoY5vEWnb1CQVBhmd4jn7qB5ccfxh7Cr29V/s2FtK//TkMEcmItI4nWGEgZnx7A8yOWN4d3p5a36rH0NE2joljDDr6SWM7UoYItLGKWGEWc/OvoSxZc9+dhaV4Zxj/obdTFuhyQlFpG1Rwgiz5PgYBmUk8+jU1WQ+OJ1X5m7i+8/O5cevZIU7NBGROpQw2oDfXTTCvz11+Xb/9vKthQ3tLiISFkoYbcBx/dP820tza5PE+U/OZl+pZrIVkbZBCaMNSIiN9m8XlVXWeWzzbq3OJyJtgxJGGzHrrjO465xh9cpzdpeEIRoRkfqUMNqIfl2T+N9xQzh1aLc65bkFShgi0jYoYbQxNRfy1ZizdifOuTBFIyJSS1ODtDETzxlOdFQUJw7syrQVO/jwm23MXbeLk4d0O/STRURCKCxnGGbWxczeNLNVZrbSzMaYWVczm2Zm2d5t2qFfqeNJ7xTPQ5cezcXH9uGh7x0NQNamgjBHJSISviapJ4ApzrnhwEhgJXA3MMM5NxSY4d2PaJ0TYhnSvRNfb94T7lBERFo/YZhZZ+A04AUA51y5c24PcBHwsrfby8DFrR1bW3RM31SWeRfwVVU7yio1DbqIhEc4zjAGAfnAi2a22Mz+bmbJQA/n3DYA77Z7GGJrcwZndGLH3jKKyyq5ddJihv16Cn/8ZDVz1u7k4Y9XhTs8EYkg4UgYMcBo4Gnn3LFAMU1ofjKzG8wsy8yy8vPzQxVjmzGom2+NjC/X7eLDb7YB8NSna7nq7/N4ZuY6KquqNYpKRFpFOBJGLpDrnJvn3X8TXwLZYWa9ALzbvIae7Jx7zjmX6ZzLzMjIaJWAw6lmHfDGJiNcl1/MwHs+YtL8nNYMS0QiUKsnDOfcdmCzmdVc1jwBWAG8D1zjlV0DvNfasbVF/dN963/XTIN+oAc/XAHAG1mbWy0mEYlM4boO4xbgVTOLA9YD1+JLXm+Y2fVADnBZmGJrUxJio1ly73dITYylpLyKo+6bWufxWdk7AUiKi27o6SIiLSYsCcM5twTIbOChCa0dS3vQJSkO8K2d8evzjyA5PoZJ83PqzGw7Z+0uHvp4JXeeNYyYaF3ALyItT98s7cyPTh3ElSccxt+vqZ9vn525nukrG+z6ERH51pQw2qk076zjQPM27GrlSEQkUihhtFOxXrPTiD6dSfTW0zCDL9fu4vPVeUz44+d88PXWcIYoIh2Mtecx/JmZmS4rK3LXvi4oLicxLprjH5zOvrJKzjmqJ1MClngF2Pjw+WGKTkTaKjNb6JxrqB/5oHSG0Y6lJceREBvNs1cfx4Uje5M5oP58jQPu/pBdRWVhiE5EOholjA7g5MHdeOrKY+mbltTg4yu27W3liESkI1LC6ED6piX6t6fefpp/O3tHUb19X/5yo+aiEpEmUcLoQPp5ZxjnHd2TId07+csXbNwN+Fbv+2y1b9jtfe8v55mZ61o/SBFpt7TiXgeSmhTL+z8dy9DuKURHGRseOo9bJi1m8lLfqn1X/d03fVdgR3h5ZTVxMfrdICKHpm+KDuaYvl1IjKsZZms8cNEIAK58/iv/Puvza5uoduwtbd0ARaTdUsLo4Lom17/A753FW/zbW/bsb81wRKQdU8KIQE99uta/va2wNmHkFpQw4O4PGfPQDD454HoOEREljAhy3diB9OniG0l10ajexEYbq7bvo7Siir2lFcxb7+sc31ZYyg3/XBjOUEWkDVKndwT4xw8z+Xx1PvdeeCQj+nTm+Vkb+P0lR5Ozu4RnZ67nuS/W41zt6n6NqZkVwMwAX3PW9sL9HNe/a8jrICLhp4QRAcYP78H44T0AuHR0Xy4d3ReAkX27sDhnDzWzw6zfWVzneRt2FjOwWzLOOZ77Yj0PfbyK847uyYkD0xk/vDsT/jST8spqTT8iEiHUJBXBbpswlNGHdWnwsZT4GG6dtJjqasfavCIe8i7y++ib7dz3/nImvrWU8spqwDfqqma0VWVVNVOWbae0oqp1KiEirUYJI4KlJcfx5o0n89SVx9Yp//PlI7n/u0fxzZZCPl+TR24DI6kC56wc/8eZXPDUbIrKKrnrraXc+K+F/KeRJWOdc7y3ZAtllUooIu2NZqsVwHc1+PNfrGd4zxTuOGsY5ZXVHH3/VE4dmsGMVTtozsfkJ6cN4p7zjqhT9umqHVz3Uha3jB/Cz88a1sgzRSSUNFutfCvHD+jKc1dncof3JR4XE8XIfl2YvrJusvjtd4/i6D6pQb3ms1+s929/k1uIc478fb6ZczfvLmnwOcH+gNlVVEZ1dfv9sSPSHilhSKNOPzyDzgm14yIW/OpMrjl5ALHRvlFS6Q1cFHgg5xwz1+Rz4V9m8+q8HH/CONDKbXu5+bVFDLznI/aX122uyt9XxlMzsqnyEsS2wv0c9+B0npu1vs5+L3+5kXve/qbRWApLKsjbpyvbRZpLCUMaddPpg/nqlxOY/6sJzJ54Bhkp8QCM8M4w3r15LNFRxi/OOtz/nOl3nM7AgOG5P34li1Xe9OofLt3G5t2+/pBdxeV1jnXuE7P4cOk2ABZvLvCXT1m2neN/P50/TlvD24ty+Wx1HhvyfaO53lqYy9Ofr2PZlkJyC0q47/3lTJqf0+j6H6f+4VNO+P2MJr0HZZVVVFZV++8v21JYbyXDnUVlXPzXOSzOKTjw6X45u0qCPntqilnZ+fxz7sYWf91Dqaiq5vHpa8Iytcz8Dbt5a2Fuqx+3oqqaP09bQ14Y6vzV+l2N9gu2JiUMaVRUlJEUF0P3lIQ6a2388rwj+PDWU+jXNYl1/3cePx0/lH5dfRcEDkhP4oVrMvmfkw4DYPrKPP8Iq4U5Bcxd71tzfPPuEsoqq/jg6631/iN85V1ACHDjv2ovILzzzaVc++ICNnnNWdl5RTwyZRUXPDWbUx75zL/fcQ9OZ3th/f/Ue0srgYabvR76eCXzN+yuU5a/r4xhv55S5yLGC56azS2TFlNWWcXS3D1UVzv+8ulalmzew3tLGl4Sd1FOAac9+hlvNOE//MadxUx8c6l/cMDOojIG3P1hnS/KyqpqfvDCfH7z3vJGX+evn63l3veWBX1cgNfm5TB3Xe3a8C/N2cDNry6qs88bWZt5fHo2z3+x/sCnA1BcVsnt/17cpISyvbCUx6au9ifo0ooqrn9pASu21q7n4pzj+8/O5ef/+brR15mybDuvL8gJ+rgAk+bnMG99bZ0nL93KQx+vrLPP6ws288SMbJ6f1XCdSyuq+O0Hyyk44MfQwWzds59Hp66qU+ebX11E9o59/n2cc1zx3Ffc+ebSplQpJHQdhjRZQmw0R/Wu24/x5o0ns2FnMTHRUQzK6MSDFx/N2MHduMn7oomNNqqrHQUl5RzdJ5VvthQy7NdTGnz9F+dsYOzgdE4clN7g418F/MduzOod++iZmsDvJq+goLicP35/pP+xXcXldOvkO1vaXVzO957+kg07i3l25nr/NSU7i8o4/vfTAfh0VR7/+moT44d397/G+MdmsmXPfkb06cyyLb4vtJe+3IgZ3HfhUXViWbPd959/7rpdXH78YXUee2PBZu57fzmL7/0OCd7a7DUDAwDOO6YXpx+e4X+NJ2Zk88rcjXRJimNUv9oh0fe8/Q0PXjyC6Cir8/qPTl0N4J+EMtBv3l3GOSN6MnZINwCqqx2PfbKav33um/a+5r24/4MVAFzwzTb+MHU1E88Zzn+yfInr42XbuXR0X47s3bnOa3+4dBvvLtlKYlw0D116TJ3Hlm8tZPLSbdx19jD/RaALNxXwvae/BGDskG6MGZzOopwCZqzKY/veUg7rmkRibDSnHt7N/zovztnAD08e4H+NGjU/Mg58r8GXQE8bmsHRfVP9dX5iRjZPzMiuU+efvrYYgBMHduXBySv5zQVH8q43B9tnq/O54oQiBmd0qvPaH3y9lRfnbMQ5uP+7dT8Dq7fvY9qK7dx8xpAG63za0AxOHJTOok0FfPjNNrYW7ic9OY7oKOO8o3v5X+fVeZu46sT+9erVWsJyhmFmG83sGzNbYmZZXllXM5tmZtnebf31RqXN6tE5gZMO+IIf3d/3T9g3LZHpd5zOlNtPZfZd43nhmsYHZ/zolIEkxEZz06uLeGRK3QWeaq5ED/wlf+YRtV/iJw+uPf7fZ62nvLKaF2Zv4O3FW1gU0FxU0+FeVe2YND+HDQEXLNb80sstqDuU+NfvLquzfkjNpI01yaLGi3M2+o8x7tHPWLalkIKSigbresMrWdz11lL2V1SRtbE2vppkAfDY1NVc/9ICfzw5u0v4OreQmWvy/V9y4PuFvHxrIQBTl2+v12wW2C+0t7SC7z8zl39+tck/5T3Auvwif7IA+L+PVtZ5b256dREbdhZz478WsmTzHv/7cN6TswBf890/Zm+gpLyS/d51OOWVdc/mPlm+nfOfnM3Tn69j1fbaX9E1X5zg+1J/d/EW/9nJ8q17+XjZdt5evIWfvV57ZvHbD1awxlscbFFOAVkb654h1lwnBL5f7ve8vZRHp67mwr/M9pdn5xXVeR//MXtDnX6u617KYv3OYq59aQFZm3z/Rmvzijjn8S8AXzPVGws2U1pRRYn3HheXVdaJY1Z2Pmc//gWPfbKGdQEzRQfW+emZ6/hk+XZ2eMdenLOH6SvzmLp8B7f9e4l/v1+9s4y1efUXRGst4TzDOMMmXSzFAAAQ00lEQVQ5tzPg/t3ADOfcw2Z2t3d/YnhCk5bQo3MC7908lmE9U/y/nn1ieeumk9lfXsX/vFD7hfXklcdy4TG9GNmvC7dMWszT3pdXclw0b9w4htTE2DpNT+C7cn36yjw6J8Tw/NWZHHXfVABmZe/0/6cG+N7Tc/3bP31tMdFRRk4DI7WufWkBj18+im+2FNZ7LLBp5GDW5xcxfeUONu4q4bZ/L2ZUP1/i3F1SQUFxOdHRxpzsnXyyYof/OTNW7WDskPR6v5Zr4uiZmnDI427aVcIxfbvwE68JrV/X2mbEjbuKOaJXZ0orqvjTJ2uYH/DlWlpRRUJsNPkH9P0898X6Jo1E+3DpNh6YvIIFG3f7V38sq6yiutoRFWVs3l1Sp3lv5pp8jujVud7rzF67k9lrd3LTuMFB1LmYYT1TuPRvvi/fybec4n8st6CEQRmdqKiq5u+z1jNpfm2TYEVVNbHRUfX6ux6YvILdQTQpVVT53pfJS7dy11tLmZmdT2/v36isstpf522F+/nBC/P9z/t8dT5DuqfUe73PV+fz+er8oOscuEBaawrLdRhmthHIDEwYZrYaGOec22ZmvYDPnXMHHaiv6zDavynLtvHinI2cPLgb/3vGYGKjoyjcX8HJD82g2PvFduuEodzxncOpqnYM/uVHAHz2i3HMys7ngmN6M/p30/juyN48eeWxvLt4C6/Ny6nzhXjn2cN4d/EWoqOszq/aA512eAZfrMknPTmuTqf86zecxJOfZjNnbcNNYe/ePJbrXlpQ54smMTba/ys7WMf0TeXMI3rwp2lr/GVnHdmDT1bsICbKqPS+vK8e059X5m4C4J/Xn1DnC+l3F4/gN+82rc8C4OmrRnPXm0vZ5/06HtK9E2vziuiblljvbAvgkmP78P7XW/0j1y4c2ZvE2CjeyGq4M7prclyDX8QXj+rN8QO78qt3amNOiY9hX1mlPwbwTZZZc2b5xBWj6vzqfuGaTK5/uenfA8/94Dh+8Z+v/X1bA9KT2LirhGE9Uli9o/7n5IJjejHZG5gBcHlmP6KijEnzG+4v6ZWawLYG+tIuObYPow/rUqfvKSkumpLyKob3TPF/Ri8c2dt/pvj45aO4/fXaOk++5RSG9UwhNrp5jUTt7ToMB3xiZgvN7AavrIdzbhuAd9u90WdLh3HOiF68/pMx3HbmUP+HPzUxljl3j2f+Lycwdkg6lx3nm/sqOsp46NKj+eCnpzCwWzJXjxlA1+Q4Jt9yCo98z9dOfvGxfXjjxjH87MzDiY02nr86k5vPGMK0O05nyu2n8dK1x3Pa4Rkce8CUKHPuHs8r153ACQO61kkWq353DicOSifK++V/07jBnH1UD56/OpPhPX2/FI/uk8pL1x7PMX1r+3X2V1Rx2XF9ef7qTE4c2PWg16789b9HA7A0t9CfLIZ078Tce8Zz1zm+30yV1Y4Yr3/iyhN8bfNH9OrMyYO71XmtmmTxzP8c5x+I0JjkuGguObaPr16vLvIni3svOJLpd5zOmUf0ILdgP0neglynDKk91v9dcjRXnVjbR/DB11t5IyuX0w/P4LYJQ+sd68BkcdGo3gC8u2SrP1l0Tojhs1+M49UfnwhQp+nlhycPAODwHp04d0SvOq9VkyyeuGIUCbGH/kqrabq84Z8L/cni3guO5PM7z+CMYRms3rGPmhO9MQHNrI9dNpJLR/fx3389azOvL8hh/PDu3Dp+SL3jHJgszj7KN5/bO4u3+JNF1+Q4Zt45jtd+fBJAnR8014711Xl4zxTOPbpnnde64KnZPNfIgINQCleT1Fjn3FYz6w5MM7NVh3yGx0swNwAcdlj9Ti3pGLok+a7xePVHJ9Upr/myDDSigS/j284cyk/HD6nXCTxuWHfGDevOnpJyRj0wDfA1hdVM+37TGYNJnRfLNK+5qKYp7fTDM5iVvZOrx/SnV6pv3+P6p7F1z36io4xj+nbhvZvHsi6/mP8s3Mzow9I4Y1h34mKi+M6RPXDOMXX5Do7q3ZlT/1DbrPbmjWPIHNCVMYO/w9iHP/WflVw8qje9UhPrtMOv/N057CmpoFunOO48exgXHtOb6Chj48Pn8/7XW/lsVR5V1Y5Lju3DGcO7c/ZRPSgqq2TF1r1kbSqgc2JsnbOPZb89GzMjNTGWl77c6C+/7pSB3pbv7OE3FxzJCQO70qdLIsu37mXhpt0kxkVz7wVH8oOT+vPA5BVUVTsy+6fx/eP70Tcticsy++Ic/GveJsYMSueHLy7wv/4z/zOas4/qyT3nHsGYh2f4Lwy9adwQBnZLpiigD2Dhr8+koKScQd06cfMZg7l0dF/iYqJY+Osz+fvsDXy1fhd905I4/+ienDOiF2cd2ZN1+UXkFuxn+dZCUhNjefDD2tFOKx84h/iYKH77wXJe9s7SAutc88Ng4jnDGT+8O/3SklixbS+LcwpIiI3m4UuP4aoTD+P+91fQJSmWY/t14QdjBpCREs+FI3vj8A33HjukG1f/o/bM7+mrRnPOiJ5sKyzllEc+paal7yenDaJ/ejJpybX9XPN/NYG9+ysZ1C2ZW8YP4dLRfYmPiWb+Lyfw9Mx1LMrZw+BuyRzRq37TVsg558L6B9wP/AJYDfTyynoBqw/13OOOO86JNNfSzXtcUWlFg4+9Nm+Te3bmWv/9qqpqV1Bc1iLHfWzqKtd/4mS3dPOeOuX5+0rd6/NzXP+Jk91r8zb5y1+cvd4t3LT7Wx93T0m56z9xsjv38S9c9o59/vKSskq3JKfA9Z842fWfONlfviSnwN06aZErraj81sf+r6fnuP4TJ7vX5+e46upqf/m6vH3ul28vdf0nTnZ/+TTbX/6LN5a495ds+dbHzd9X6vpPnOxOfmiGW719r798f3mlW7Rpd706z1u/y/3wH/NcSdm3r/N3n5rl+k+c7N5YULfOG3cWuYlvfu36T5zs/vpZbZ1vnbTIvbMo91sfNxhAlmvG93Wr92GYWTIQ5Zzb521PAx4AJgC7XG2nd1fn3F0Hey31YUhH45xjztpdjBmcXu/sqCVs2lVMz9QE4mOi6z02f8NuSiuqOO3wjBY/7t7SCrbtKWVYz/q/iksrqvjLp2v5yemDSEmIbfFjL83dw+CMTiTH129Qmb5iB2WV1Zx/TK8GnvntFJZUkLunpN4QdICS8kr+9Mkabj1zKJ1DUOdDaW4fRjgSxiDgHe9uDPCac+73ZpYOvAEcBuQAlznndjfyMoAShohIczQ3YbR6H4Zzbj0wsoHyXfjOMkREpA3S1CAiIhIUJQwREQmKEoaIiARFCUNERIKihCEiIkFRwhARkaAoYYiISFDCMlttSzGzfGDTIXdsWDdg5yH36lhU58igOkeGb1Pn/s65Jl/S364TxrdhZlnNudKxPVOdI4PqHBnCUWc1SYmISFCUMEREJCiRnDCeC3cAYaA6RwbVOTK0ep0jtg9DRESaJpLPMEREpAkiMmGY2TlmttrM1nqLNbVJZvYPM8szs2UBZV3NbJqZZXu3aV65mdmTXp2WmtnogOdc4+2fbWbXBJQfZ2bfeM950sy3PmVzjtFC9e1nZp+Z2UozW25mt0VAnRPMbL6Zfe3V+bde+UAzm+fF87qZxXnl8d79td7jAwJe6x6vfLWZnR1Q3uDnvTnHaOG6R5vZYjObHAl1NrON3mdviZlleWXt67PdnGX62vMfEA2sAwYBccDXwJHhjquRWE8DRgPLAsr+ANztbd8NPOJtnwd8DBhwEjDPK+8KrPdu07ztNO+x+cAY7zkfA+c25xgtWN9ewGhvOwVYAxzZwetsQCdvOxaY5x3nDeAKr/wZ4CZv+3+BZ7ztK4DXve0jvc9yPDDQ+4xHH+zz3tRjhODzfQfwGjC5OfG0tzoDG4FuB5S1q8922L8UW/vPe0OnBty/B7gn3HEdJN4B1E0YDa59DjwLXHngfsCVwLMB5c96Zb2AVQHl/v2aeowQ1v094DuRUmcgCVgEnIjvgqyYAz+zwFRgjLcd4+1nB36Oa/Zr7PPuPadJx2jhuvYFZgDjgcnNiacd1nkj9RNGu/psR2KTVB9gc8D9XK+svejhnNsG4N1298obq9fBynMbKG/OMVqc1yRwLL5f3B26zl7TzBIgD98a9+uAPc65ygaO6Y/He7wQSD9InI2VpzfjGC3pceAuoNq735x42ludHfCJmS00sxu8snb12W71JVrbAGugrCMMFWusXk0tb84xWpSZdQLeAm53zu31mmKbEk+7qrNzrgoYZWZd8K13f8RBjtnUujX0o/BQ70VI62xmFwB5zrmFZjYuiGO2+zp7xjrntppZd2Cama06yL5t8rMdiWcYuUC/gPt9ga1hiqU5dphZLwDvNs8rb6xeByvv20B5c47RYswsFl+yeNU593Yz42lXda7hnNsDfI6vPbmLmdX8oAs8pj8e7/FUYPdB4mysfGczjtFSxgLfNbONwL/xNUs93ox42lOdcc5t9W7z8P0wOIF29tmOxISxABjqjZaIw9fB9X6YY2qK94GakRHX4Gvnrym/2hv5cBJQ6J1+TgXOMrM0b3TEWfjabbcB+8zsJG80xdUHvFZTjtEivDheAFY65/4UIXXO8M4sMLNE4ExgJfAZ8F+NxFMT538BnzpfA/T7wBXeaJ+BwFB8naANft695zT1GC3COXePc66vc26AF8+nzrmrOnKdzSzZzFJqtvF9JpfR3j7bLdmp017+8I0OWIOvrfhX4Y7nIHFOArYBFfh+DVyPr111BpDt3Xb19jXgr16dvgEyA17nOmCt93dtQHmm96FdB/yF2gs5m3yMFqrvKfhOiZcCS7y/8zp4nY8BFnt1Xgbc65UPwvfltxb4DxDvlSd499d6jw8KeK1feXGuxhshc7DPe3OOEYLP+DhqR0l12Dp7x/3a+1teE1N7+2zrSm8REQlKJDZJiYhIMyhhiIhIUJQwREQkKEoYIiISFCUMEREJihKGRDQz+9K7HWBm/93Cr/3Lho4l0l5pWK0I4E1R8Qvn3AVNeE60803r0djjRc65Ti0Rn0hboDMMiWhmVuRtPgycar61Cn7mTQj4qJkt8NYK+Im3/zjzrdnxGr6LnTCzd70J5ZbXTCpnZg8Did7rvRp4LO/K2kfNbJn51i+4POC1PzezN81slZm96l21i5k9bGYrvFgea833SKRGJE4+KNKQuwk4w/C++Audc8ebWTwwx8w+8fY9ARjhnNvg3b/OObfbm9pjgZm95Zy728x+6pwb1cCxLgVGASOBbt5zvvAeOxY4Ct+cPnOAsWa2ArgEGO6cczVTiYi0Np1hiDTsLHzz7CzBN8V6Or65igDmByQLgFvN7GvgK3yTuQ3l4E4BJjnnqpxzO4CZwPEBr53rnKvGNzXKAGAvUAr83cwuBUq+de1EmkEJQ6RhBtzinBvl/Q10ztWcYRT7d/L1fZyJb/GdkfjmhUoI4rUbUxawXYVvsZ9KfGc1bwEXA1OaVBORFqKEIeKzD9+ysDWmAjeZb7p1zOxwb5bRA6UCBc65EjMbjm9q8hoVNc8/wBfA5V4/SQa+pXjnNxaY+dYHSXXOfQTcjq85S6TVqQ9DxGcpUOk1Lb0EPIGvOWiR1/Gcj+/X/YGmADea2VJ8M6Z+FfDYc8BSM1vkfNN313gH3/KgX+Obnfcu59x2L+E0JAV4z8wS8J2d/Kx5VRT5djSsVkREgqImKRERCYoShoiIBEUJQ0REgqKEISIiQVHCEBGRoChhiIhIUJQwREQkKEoYIiISlP8H0JKqyf4LYBMAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"iterations = np.arange(len(cost_list))*10000\n",
"plt.plot(iterations,cost_list)\n",
"plt.xlabel(\"iterations\")\n",
"plt.ylabel(\"avg cost\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" "
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.4"
}
},
"nbformat": 4,
"nbformat_minor": 1
}