21. Gradient Descent in Linear Regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gradient Descent in Linear Regression\n",
    "\n",
    "$J = argmin(\\sum_{1}^n(y - \\hat{y} )^2)$\n",
    "\n",
    "$\\hat Y = X \\beta$\n",
    "\n",
    "$\\nabla J(\\beta) = \\begin{bmatrix} \\partial J/ \\partial \\beta_0 \\\\ \\partial J/ \\partial \\beta_1\\\\ ····· \\\\ \\partial J/ \\partial \\beta_n  \\end{bmatrix} =  \\begin{bmatrix} \\sum_{i = 1}^n 2(Y_i - X_i\\beta)(-1) \\\\ \\sum_{i = 1}^n 2(Y_i - X_i\\beta)(-X_{i1}) \\\\ ····· \\\\ \\sum_{i = 1}^n 2(Y_i - X_i\\beta)(-X_{i(k-1)})  \\end{bmatrix} = \n",
    "-2 \\begin{bmatrix} \\sum_{i = 1}^n (Y_i - X_i\\beta) \\\\ \\sum_{i = 1}^n (Y_i - X_i\\beta)(X_{i1}) \\\\ ····· \\\\ \\sum_{i = 1}^n (Y _i- X_i\\beta)(X_{i(k-1)})  \\end{bmatrix}$\n",
    "\n",
    "\n",
    "Modify:\n",
    "\n",
    "$J = MSE = argmin \\frac {1}{n}(\\sum_{1}^n(y - \\hat{y})^2)$\n",
    "\n",
    "$\\nabla J(\\beta) = \\begin{bmatrix} \\partial J/ \\partial \\beta_0 \\\\ \\partial J/ \\partial \\beta_1\\\\ ····· \\\\ \\partial J/ \\partial \\beta_n  \\end{bmatrix} = -\n",
    "\\frac {2}{n} \\begin{bmatrix} \\sum_{i = 1}^n (Y_i - X_i\\beta) \\\\ \\sum_{i = 1}^n (Y_i - X_i\\beta)(X_{i1}) \\\\ ····· \\\\ \\sum_{i = 1}^n (Y _i- X_i\\beta)(X_{i(k-1)})  \\end{bmatrix} $\n",
    "\n",
    "Sometimes:\n",
    "\n",
    "$J = argmin \\frac {1}{2n}(\\sum_{1}^n(y - \\hat{y})^2)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(666)\n",
    "x = 2 * np.random.random(size = 100)\n",
    "y = x * 3. + 4. + np.random.normal(size = 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = x.reshape(-1, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(100, 1)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(100,)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYJklEQVR4nO3dfYwdV3nH8d9jexF2oHFQlpYsuA4SMiJKicMKhbhCJBSchpdYgSpFILWIyuKPVoSqrhy1AvpHFUuu1LRSVcmiFBARGJLURQTqIEyFGppUa9YhocHl3WSTNqbgUmCrGHP6x9713r07c+/M3JnzNt+PZHn3vuycO3fuc8885zlnzDknAEB6NoVuAACgGQI4ACSKAA4AiSKAA0CiCOAAkKgtPjd2+eWXu507d/rcJAAk7+TJkz9wzs2O3u41gO/cuVMLCws+NwkAyTOz7xXdTgoFABJFAAeARBHAASBRBHAASBQBHAAS5bUKBUC8ji0u6fDx03ry3LKu2L5VB/bu0r7dc6GbhTEI4AB0bHFJd9z3qJbPX5AkLZ1b1h33PSpJBPGIkUIBoMPHT18M3quWz1/Q4eOnA7UIVRDAAejJc8u1bkccCOAAdMX2rbVuRxwI4AB0YO8ubZ3ZvO62rTObdWDvrkAtQhUMYgK4OFBJFUpaCOAAJK0EcQJ2WkihAECiCOAAkCgCOAAkigAOAIkigANAogjgAJAoAjgAJIoADgCJIoADQKKYiQkgW7lfpIIADiBLfbhIBSkUAFnqw0UqCOAAstSHi1QQwAFkqQ8XqSCAA8hSHy5SwSAmgGTUqSrpw0UqCOAAktCkqsTXRSpClSsSwAEkYVxVSZfBclJwDlmuSA4cQBJCVJWsBuelc8tyWgvOxxaXLj4mZLkiARxAEkJUlVQJziHLFQngAJIQoqqkSnAOWa5IAAeQhH2753TnrVdrbvtWmaS57Vt1561Xd5pnrhKcQ5YrMogJIBm+qkpWHdi7a90ApbQxOE8qV+yyQoUADgAlqtaSl32xdF2hQgAH0AtNe8LT9Pq7Ln0kBw4ge0XlgLcfPaVr/uyBdSWBbeu6QmViADezD5nZ02b22NBtzzOzz5vZNwb/X9ZKawCgA0U9YUk6t3x+Q113m7quUKnSA/+wpJtGbjso6QvOuZdI+sLgdwCI0rgeb5eTbrquUJkYwJ1zX5L0w5Gbb5H0kcHPH5G0r5XWAEAHJvV4u5p003XpY9NBzF92zj0lSc65p8zs+WUPNLP9kvZL0o4dOxpuDgCaKyoHHNblpJsuSx87H8R0zh1xzs075+ZnZ2e73hwAbLDaE75s28yG+1JeI7xpAP8vM3uBJA3+f7q9JgFA+/btntPi+16vu267xutszi41TaF8WtLvSDo0+P8fW2sRAHTI92zOLlUpI/y4pH+VtMvMnjCzd2klcL/OzL4h6XWD3wEAHk3sgTvn3lZy12tbbgsAoAam0gNAA6EuozaMAA40EMOHF5N19T6FvIzaMAI4UFMsH16MD9Bdvk+hrs85isWsgJpCXgMRayZdr7LL9ynkZdSGEcCBmmL58PbdpADd1ft0bHFJm8wK7/NxGbVhBHCgppDXQMSaSQG6i/dptdd/wbkN94WY0UkAB2oKeQ1ErJkUoLt4nz7w6a8Vrqey2SzIjE4COFBTiIvrYqNJAbrt9+nY4pLOLZ8vvO8XzgV5/6lCARrIaTp2qqpcr7LN92nc4Geo9BkBHECyfH6Rjhv8DJU+I4UCABWU9bIv2zYT7GyMAA6gFccWl7Tn0AldefB+7Tl0otOLBYdQlnN//5uuCtQiUihAcnxM46+7jT7MTq2Sc/fNXEE9Y1fm5+fdwsKCt+0BuRkNlNJKL7DNKpgm29hz6ISWCnLEl22b0eL7Xt9Ku/rMzE465+ZHbyeFAiTExzT+JtsoG+D70c/OZ5dKiQkpFCAhbUwPL0qPSGupgbJz8nHbuGL71sIeuCTdfvSUDh8/3Tjd0CSdc/j4aS2dW9ZmM11wTnMRpDu6QAAHElIWKEcrJMqCXlGu+sCnHpFMOn9hfDp1XK3zgb27dPvRU6X3N82J182tjz5+dcp7jjl5iRQKEK2iqo4q08PHrdJXlB45/ws3MXhPmoK+b/ectm/deMX3YU1SPXXTOUWPn2b7sSOAAxEqC8KSJk4PHxf0mqzEZ5Le8orJE2Y+8OarNny5jKq7/bopo0l/v+lKhLGWSJJCASI0Lgg/ePDGscF0XHAbl6su4yR98etnJz5uuMyubBt1p5xXTRlNenzT7Utxl0jSAwciNM1g5bjgVpSCmdlkmtlcvL51ne1KKwHtwYM36q7brmllJcC6KwoWPX6a7UtxX8CDAA5E6NKSfHKVHuS4oFe0Qt/h33q5Dr/15Zob87fr9lzbWgmw7t8Zfry0ssyrpti+FPcFPEihAJE5triknz7z8w23z2yySj3ISTMGyxaAKqpSkZr3XNtaaKru32l7gau6aRyfCOBAZA4fP11YFfKcZ2+pHJiaBrEYp4uHdmDvrta+1NpGAAciU3Zqfu5nxRcTaBtrna8X85caARyITKyn7FVmRPpYaCuEWL/UGMQEIhPjNTf/9Nijeu/RU4WTg1aNm0CEbhDAgcjEds3NY4tLuvuhMxvWSBktpYu53C5XpFAQpVxPxauK6ZT98PHTlRa4irncLlcEcEQn5plvfTRpFcLhn5vk7nP8svb1mkihIDqciselLACb1l/Mt0nuPse8uc/XRABHdDgVj0tRYDZJb79ux7peZZPcfdmX9e1HT0W1aFQdPjsgpFAQnVjL6PqqTh103dz9uC/lVFNnPjsgBHBEJ+aZb33V1aDqpNUDV3uuKQVwnx0QUiiITixldLGuAR2zuvts3OqBq1JLnfms45+qB25m75X0e1pZMvhRSe90zv1fGw1Dv4Uuo8u5EqarCokm+6yLNcRD8zn13pwbfyml0ieazUn6F0kvc84tm9knJX3WOffhsufMz8+7hYWFRtsDfNpz6ERhQJnbvlUPHrwxQIvaUbbaYBtnONPusy7bljozO+mcmx+9fdoUyhZJW81si6Rtkp6c8u8BUci1EqbLColp91ksqbOUNE6hOOeWzOwvJJ2RtCzpAefcA6OPM7P9kvZL0o4dO5puDvAq10qYLr+Y2thnoVNnqWncAzezyyTdIulKSVdIusTM3jH6OOfcEefcvHNufnZ2tnlLAY9iXFCqDeMutzatXPdZzKZJofyGpO845846585Luk/S9e00Cwgr19P5prMlq1SW5LrPYjZNFcoZSdeZ2TatpFBeK4kRSmQjx9P5uhUSdStLctxnMZsmB/6wmd0j6SuSfi5pUdKRthoGoBt1guy4QU8CdXhT1YE7594v6f0ttQVYJ8dV6mJWtL9zrcbJBVPpEaWcJ9LEqGx/X7p1RueWN16LM/VqnFwwlR5RYklZv8r2t5moLIkYARxR4tTdr7L9eu5n56ksiRgpFEQp14k0sSrb35vMJCnp5QNyRg88U6mvpMekEL/KVgW84FzyV8jJGQE8QzlcpopJIX6t7u/Ngx73MMYe4kUKJUO51O4yKcSvfbvn9N6jpwrvY+whTgTwDDEAWA115hsx9pAWUigZ6nLBolzkkGbqAmMPaSGAZ8jXhzDlgVLqzIsx9pAWUigZ8nFJp9RnSpJmKsfYQzoI4Jnq+kOY+kApuV7kgACeEZ+DcjH3YKvshwN7dxVef3GaNBODovCNAJ4J3ymNWHuwVfdD22mm1FNKSBMBPBNtpDTq9CC76MG20b46+6HNNFNqKSXOFvJAAI9I0YdKqtZLnDal0eTKK1Xb1oaq7QuV2okppTQpOHO2kA8CeCSKPlQHPvWIZNL5C+7ibWUftGlTGk16kD6rFaq2L1RqJ5aUUpXgnNrZAspRBx6Jog/V+V+4i8F7VVmt8rS13zH1IKWNNeZFwVHa2L5QE1FimQBTpb49tvcazdED92TSaW2dD0/RY6dNacTSg5SKe5EmyRU89tKtM9pz6MS613znrVd7z+/6TimVqRKcY3qvMR0CuAdVTmvLPlRFyj5o06Q0QgxKlinqRTppQxCf2WT66TM/v3jJr9X9euetVwdZvzqGCTBVgnNM7zWmQwrFgyqntUWn4DObTDOb1y/v2dUHLaYp1GW9SCeta99znr2lcoqpL6qkcmJ6rzEdeuAeVDmtLTsFL7qtqw+a7x5kWVqprBc5t33rup71lQfvL/y7MeZyfZXtVU3lxHC2gOkRwD2omnMs+1Dl+EEbl1aqeoqfSi7Xd9kewbk/SKF4EEuFQkwmlbJVOcX3uV+nWXmRlQ/RFXrgHoxLj4xWUPSl5zQprVSlF+mr8mPaHjRle+gKAdyT0YDUl9lwdfPcddMfbaULxuWop534kkqqB+khhRJIH06rx131Jqa00qSr80zbg47ptSIvBPBA+nBa3UaeO3Q7pekvURfTa0VeSKEE0ofT6jby3D5MamcbE19iea3ICz3wQPpwWt3mxZW7vP7mpHbSg0as6IEHEsvaGV1qa8p21wO+VdpJDxoxIoAHlHtQaOtLquvlT/vwZYo8EcDRqTa+pHwM+Ob+ZYo8RR/AufQT+jDgCzQR9SDmpPpc9EMfBnyBJqYK4Ga23czuMbOvm9njZvaqthom9WOyCyajCgQoNm0K5a8k/ZNz7q1m9ixJ21po00V9mOyCashRAxs17oGb2S9JerWkv5Mk59wzzrlzLbVLUrt1xACQm2lSKC+WdFbS35vZopl90MwuGX2Qme03swUzWzh79mytDZD7hE9dThYCujBNAN8i6VpJf+uc2y3pp5IOjj7IOXfEOTfvnJufnZ2ttQFyn/CFAXOkaJoc+BOSnnDOPTz4/R4VBPBpkfuED11PFgK60LgH7pz7T0nfN7PVfMZrJf17K60CPGPAHCmatgrlDyTdPahA+bakd07fpHgwiag/mCyEFE0VwJ1zpyTNt9OUuPTlijk5aOOLtq2FtwCfop9KHwo50fCqBOa2vmhZ0AopIoCXICcaVtXA3OYXLQPmSE3Ua6GExCSisKouo8AXLfqMAF6CSURhVQ3MfNGizwjgJZhEFFbVwMwXLfqs1znwSYNk5ETDqVoVwuAj+qy3AZwywbjVCcx80aKvehvAcygTzH2iEYEZGC+bAF43mKVWvTD6+m546azuPbkU1RlE7l8oQGyyGMRsspJcStULRa/v7ofORHW1IlbzA/zLIoA3ufRaStULRa/PlTy2yzOIcetlc/k7wL8sUihN0iEpVS/UCcpdnUFMGvRNLSUF5CCLAN50JblUBsnKXp9pfU+8yRlE1bz1pEFfVvMD/MsihZJSOqSJstf39ut2TDXRqE7eelIPO/f3AIhRFj3wlNIhTXT1+uqUUk7qYef+HgAxMufKhsPaNz8/7xYWFrxtD+NdefD+wsFQk/SdQ29Yd9toDlxa6WGntrwApY5IkZmddM5tuPZCFj3wGIUMFFW3XSdvnUMPm9m3yE22ATx0AA0VKOpsu+5VaFIZ9C2Tw+xbYFgWg5ijQk8qKQsUtx89taF+2te2i+qx+7biIqWOyE2WPfDQPa1xAaHr3njdIJV6r7oOSh2Rmyx74EUf0nG3t21SQOhyhmJKSwT4RqkjcpNlAN9sVuv2thUFilFdnbYTpMr1LWWE/GWZQrlQUhpZdnvbhis2ynr9ZT3iaQdfc6gW6VKfUkbIX5YBfK4k1znnMY2wGijK6qeLesRtVa8QpIB+yDKFUiWNMG5lvTbVOW1nRT8AdWTZA5+URvBdp121R0yZG4A6sgzg0vigGbrMsAxlbgDqyDKFMkmsPV0qSADU0csAHmutNGVuAOpIIoXS9romddcA8YkKEgBVRR/AuxhwpFYaw1hiFqmKPoB3NeBIT3dNnwMYS8wiZdEH8FgHHLsQIpD2PYDFWpEEVBH9IGasA45tC7UEbt8nD/Wpg4D8RB/Ay0rrbnjprJeZlJKfWZuhAmnfA1hfOgjI09QB3Mw2m9mimX2mjQaNGi6tk1ZWFFw+f0F3P3TGS2/VV884VCDtewCj9h4pa6MH/h5Jj7fwd0rt2z138YO2uqLg6LqCXfVWffWMQwXSvgcwau+RsqkGMc3shZLeIOnPJf1hKy0qURRIR3XRWy37m0vnlnVscam1D3qo2nRKKqlIQrqmrUK5S9IfS3pu2QPMbL+k/ZK0Y8eOxhuqEpy76K2WrU8iqdVqjZCBlAAGpKlxADezN0p62jl30sxeU/Y459wRSUckaX5+vvEVFcYFUqm73mpRz3hV2+VmXQbSPtd6A7maJge+R9Kbzey7kj4h6UYz+1grrSpQlKtdvUBal3nL1RxpmRSqNUKVKALoVuMeuHPuDkl3SNKgB/5Hzrl3tNOsjUKnGMouj5ZCtQaTVYA8RT8Tc5jPXO1oyuGGl87q3pNLUS6ANUnfa72BXLUykcc598/OuTe28bdiUJRyuPfkkt7yirkky836XusN5Cr6mZghlKUcPvbQGUnSX952jR48eGMSwVui1hvIVVIpFF/GpRZSXOyJWm8gTwTwApNKFrsYAOy6zI9abyA/pFAKFKUcRrU5AEiZH4AmCOAFRhfQKtLmAGDfl3QF0AwBvMS+3XN68OCNuuu2azofAKTMD0ATBPAJfKxWR5kfgCYYxKyg6wHAUCsRAkgbATwClPkBaIIUCgAkih54BPp+ZXgAzdADjwBlhACaIIBHgDJCAE0QwCNAGSGAJpIL4McWl7Tn0AldefB+7Tl0Iovp5qwWCKCJpAYxcx3so4wQQBNJBfCcLw3GaoEA6koqhcJgHwCsSSqAM9gHAGuSCuAM9gHAmqRy4Az2AcCapAK4xGAfAKxKKoUCAFiTXA+864v/AkAqkgrguU7kAYAmkkqhsGofAKxJKoAzkQcA1iQVwJnIAwBrkgrgTOQBgDVJDWIykQcA1iQVwCUm8gDAqqRSKACANQRwAEgUARwAEkUAB4BEEcABIFHmnPO3MbOzkr7X4KmXS/pBy81pA+2qL9a20a56Ym2XFG/bpmnXrzrnZkdv9BrAmzKzBefcfOh2jKJd9cXaNtpVT6ztkuJtWxftIoUCAIkigANAolIJ4EdCN6AE7aov1rbRrnpibZcUb9tab1cSOXAAwEap9MABACMI4ACQqKAB3MxuMrPTZvZNMztYcL+Z2V8P7v+qmV1b9bke2vb2QZu+amZfNrOXD933XTN71MxOmdmC53a9xsz+Z7DtU2b2vqrP7bhdB4ba9JiZXTCz5w3u63J/fcjMnjazx0ruD3KMVWhXqONrUruCHF8V2+b9GDOzF5nZF83scTP7mpm9p+Ax3R1jzrkg/yRtlvQtSS+W9CxJj0h62chjbpb0OUkm6TpJD1d9roe2XS/pssHPv7natsHv35V0eaB99hpJn2ny3C7bNfL4N0k60fX+GvztV0u6VtJjJfeHOsYmtcv78VWxXd6Pr6ptC3GMSXqBpGsHPz9X0n/4jGMhe+CvlPRN59y3nXPPSPqEpFtGHnOLpI+6FQ9J2m5mL6j43E7b5pz7snPuR4NfH5L0wha337hdHT237b/9Nkkfb2nbYznnviTph2MeEuQYm9SuQMdXlf1VpuvPZN22eTnGnHNPOee+Mvj5fyU9Lmn0ggWdHWMhA/icpO8P/f6ENr7wssdUeW7XbRv2Lq18w65ykh4ws5Nmtj9Au15lZo+Y2efM7Kqaz+2yXTKzbZJuknTv0M1d7a8qQh1jdfg6vqryfXzVEuoYM7OdknZLenjkrs6OsZBX5LGC20ZrGsseU+W506j8983sBq18wH596OY9zrknzez5kj5vZl8f9B58tOsrWlk34SdmdrOkY5JeUvG5XbZr1ZskPeicG+5JdbW/qgh1jFXi+fiqIsTxVZf3Y8zMnqOVL4zbnXM/Hr274CmtHGMhe+BPSHrR0O8vlPRkxcdUeW7XbZOZ/ZqkD0q6xTn336u3O+eeHPz/tKR/0Mqpkpd2Oed+7Jz7yeDnz0qaMbPLqzy3y3YN+W2NnNp2uL+qCHWMTRTg+Joo0PFVl9djzMxmtBK873bO3VfwkO6OsbaT+jWS/1skfVvSlVpL4F818pg3aH3y/9+qPtdD23ZI+qak60duv0TSc4d+/rKkmzy261e0NkHrlZLODPZfZ/us6t+WdKlWcpiX+NhfQ9vYqfJBuSDHWIV2eT++KrbL+/FVtW0hjrHBa/+opLvGPKazY6zVndvgxd+slVHbb0n6k8Ft75b07qGd8zeD+x+VND/uuZ7b9kFJP5J0avBvYXD7iwdvxCOSvtZ22yq06/cH231EK4Nf1497rq92DX7/XUmfGHle1/vr45KeknReKz2ed8VwjFVoV6jja1K7ghxfVdoW4hjTSmrLSfrq0Ht1s69jjKn0AJAoZmICQKII4ACQKAI4ACSKAA4AiSKAA0CiCOAAkCgCOAAk6v8BwvYkezFD44kAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(x, y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gradient Descent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def J(beta, X_b, y):\n",
    "    try:\n",
    "        return np.sum((y - X_b.dot(beta)) ** 2) / len(X_b)\n",
    "    except:\n",
    "        return float('inf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dJ(beta, X_b, y):\n",
    "    res = np.empty(len(beta))\n",
    "    res[0] = np.sum(X_b.dot(beta) - y)\n",
    "    for i in range(1, len(beta)):\n",
    "        res[i] = (X_b.dot(beta) - y).dot(X_b[:, i])\n",
    "    return res * 2 / len(X_b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gradient_descent(X_b, y, initial_beta, eta, n_iters = 1e4, epsilon = 1e-8):\n",
    "\n",
    "    beta = initial_beta\n",
    "    i_iter = 0\n",
    "    \n",
    "    while i_iter < n_iters:\n",
    "        gradient = dJ(beta, X_b, y)\n",
    "        last_beta = beta\n",
    "        beta = beta - eta * gradient\n",
    "    \n",
    "        if (abs(J(beta, X_b, y) - J(last_beta, X_b, y)) < epsilon):\n",
    "            break\n",
    "        \n",
    "        i_iter += 1\n",
    "    \n",
    "    return beta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_b = np.hstack([np.ones((len(X), 1)), X])\n",
    "initial_beta = np.zeros(X.shape[1] + 1)\n",
    "eta = 0.01\n",
    "\n",
    "beta = gradient_descent(X_b, y, initial_beta, eta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([4.02145786, 3.00706277])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "beta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gradient Descent Packaging\n",
    "\n",
    "Batch Gradient Descent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearRegression()"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from playML.LinearRegression import LinearRegression\n",
    "\n",
    "lin_reg = LinearRegression()\n",
    "lin_reg.fit_gd(X, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([3.00706277])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lin_reg.coef_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.021457858204859"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lin_reg.intercept_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Modify:\n",
    "\n",
    "$\\nabla J(\\beta) = \\begin{bmatrix} \\partial J/ \\partial \\beta_0 \\\\ \\partial J/ \\partial \\beta_1\\\\ ····· \\\\ \\partial J/ \\partial \\beta_n  \\end{bmatrix} = -\n",
    "\\frac {2}{n} \\begin{bmatrix} \\sum_{i = 1}^n (Y_i - X_i\\beta) \\\\ \\sum_{i = 1}^n (Y_i - X_i\\beta)(X_{i1}) \\\\ ····· \\\\ \\sum_{i = 1}^n (Y _i- X_i\\beta)(X_{i(k-1)})  \\end{bmatrix} = (\\frac {2}{n} (X\\beta - Y)' X)' = \\frac {2}{n} X'(X\\beta - Y)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import datasets\n",
    "from playML.LinearRegression import LinearRegression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "boston = datasets.load_boston()\n",
    "\n",
    "X = boston.data\n",
    "y = boston.target\n",
    "\n",
    "X = X[y < 50]\n",
    "y = y[y < 50]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.2, random_state = 666)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 4.01 ms, sys: 1.35 ms, total: 5.36 ms\n",
      "Wall time: 2.15 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.8129794056212895"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from playML.LinearRegression import LinearRegression\n",
    "\n",
    "regM = LinearRegression()\n",
    "%time regM.fit_normal(X_train, y_train)\n",
    "regM.score(X_test, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/yasmine/anaconda3/lib/python3.7/site-packages/numpy/core/fromnumeric.py:87: RuntimeWarning: overflow encountered in reduce\n",
      "  return ufunc.reduce(obj, axis, dtype, out, **passkwargs)\n",
      "/Users/yasmine/Documents/学习/language/python/machine learning/playML/LinearRegression.py:26: RuntimeWarning: overflow encountered in square\n",
      "  return np.sum((y - X_b.dot(beta)) ** 2) / len(X_b)\n",
      "/Users/yasmine/Documents/学习/language/python/machine learning/playML/LinearRegression.py:51: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  if (abs(J(beta, X_b, y) - J(last_beta, X_b, y)) < epsilon):\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "LinearRegression()"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "regGD = LinearRegression()\n",
    "regGD.fit_gd(X_train, y_train) # eta is too big"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearRegression()"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "regGD.fit_gd(X_train, y_train, eta = 0.000001)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.2758681872447726"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "regGD.score(X_test, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 2min 56s, sys: 2.07 s, total: 2min 58s\n",
      "Wall time: 46 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "LinearRegression()"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%time regGD.fit_gd(X_train, y_train, eta = 0.000001, n_iters = 1000000 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7542932581943915"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "regGD.score(X_test, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Normalization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import StandardScaler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "standardscaler = StandardScaler()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "StandardScaler()"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "standardscaler.fit(X_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train_standard = standardscaler.transform(X_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 845 ms, sys: 17.8 ms, total: 862 ms\n",
      "Wall time: 224 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "LinearRegression()"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "regGD2 = LinearRegression()\n",
    "%time regGD2.fit_gd(X_train_standard, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_test_standard = standardscaler.transform(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8129873310487505"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "regGD2.score(X_test_standard, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Strength\n",
    "\n",
    "Deal with relatively large n with less time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 500\n",
    "n = 5000\n",
    "\n",
    "big_X = np.random.normal(size = (m, n))\n",
    "\n",
    "true_theta = np.random.uniform(0.0, 100, size = n + 1)\n",
    "big_y = big_X.dot(true_theta[1:]) + true_theta[0] + np.random.uniform(0, 10, size = m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 11.9 s, sys: 401 ms, total: 12.3 s\n",
      "Wall time: 3.12 s\n"
     ]
    }
   ],
   "source": [
    "%time regM.fit_normal(big_X, big_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "standardscaler.fit(big_X)\n",
    "big_X_standard = standardscaler.transform(big_X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 6.39 s, sys: 58.9 ms, total: 6.44 s\n",
      "Wall time: 1.62 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "LinearRegression()"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%time regGD2.fit_gd(big_X_standard, big_y)"
   ]
  }
 ],
 "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.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}