Normal Equation for Linear regression/Normal equation for linear regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Derivation and Implementation of Normal equation for linear regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "np.set_printoptions(suppress=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "ex1data2 = pd.read_csv('ex1data2.txt',header = None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = ex1data2[[0,1]].values\n",
    "y= ex1data2[[2]].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "X0 = np.ones((len(X),1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.append(X0, X, axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inverse a matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-3.,  2.],\n",
       "       [ 2., -1.]])"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "A = np.array([[1,2],[2,3]])\n",
    "\n",
    "np.linalg.inv(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAC3CAYAAAAM7RnXAAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAAEnQAABJ0Ad5mH3gAABVMSURBVHhe7ZyBkdQ4E4VJgRhIgRwIgRhIgQzIgAyIgAhIgATIgBz4692tqHf9y+2WR56xZr+vquvYWdmWWt39JHn23vwBAABIQCgAACAFoQAAgBSEAgAAUhAKAABIQSgAACAFoQAAgBSEAgAAUhAKAABIQSgAACAFoQAAgBSEAgAAUhAKAABIQSgAACAFoQAAgBSEAgAAUhAKAABIQSgAACAFoQAAgBSEAgAAUhAKAABIQSgAACAFoQAAgBSEAgAAUhAKAABIQSgAACAFoQAAgBSEAgAAUhAKAABIQSgAACAFoQAAgBSEAgAAUhAKAABIQSjg0rx79+7Pmzdv/rFv3769fAoA9wShgEvTREL25cuXl08B4J4sJxRaVXrx0IpzdKX58+fPP2/fvv3nev33169fL7+Bq+FzPUMofvz48XfuR+3Tp08vd4GzUC63+VFuk5vXYDmh+PDhw/8lsGwkoFQs/Fr9DNfE52mGUOgefs8RU+zBufhRo4xd5DV4GqEYSWKEYh18nh4tFOwozgehuCZPIxSyasFHKNbB5+kMoWDuz8H9/P79+5dP9/n69StHTxdkaaH4+PHj33/LqgGJUKyDzxNCsQ5xQQdrs7RQ6N86DvCArLzYRijWwecJoVgHhOK5WF4otDX1gNR2dQ+EYh18nhCKdUAonovlhULEXcVeQUEo1sHnCaFYB4TiuXgKoYi7Cr0M+/379z+/63GLUOhZnz9//ufZ8fv4+ky/099pVPC/5/CdkD739y/6Xe8bNy6QfuSmPup3/g0Svb9R37ZeDur6+M5HP4/+jYr4/v37Xx/5/TRWfaYXltn8OH79KkLRfOnxIf/r2T5un78sZnweq/Ph18jfeygu1E799mtl6rv6uuereN2e9fql8bXf635VZualaGPRvfw6+SDmlv6tz7Zy6xl4CqEQsQAoMLY4IhRKcE/sPVPf9oph7IdQosRA99877otWRDVuv6ZnXmyUBCoEvXbNKmMRum+1WGiMlaLn11xdKCq+1LglpMLnL+uHX1/1QfWatqjw9pllsdBrn1mvX3F+9jgjL4Vfo7nRNT5fW1aJ6RV5GqHQRMYCu6XwsUDvFQsVgGoBdIurkUjshwqI/xwt4r6QQOwVKTc9OxOlaFppZmh12Ltuz/Z8722rRTIjFqK951fRPFd9KVM/fP6yfsTrKlSvGelzM887ZzRHejuKOD8ZZ+Wl8Pbq54if2kLgmXgaoRCxWG0Vt1igsySVAMUCrOBUkXUhUuApyGNA6dotYj/8WvVdAac2rahE3BdurX8N3ac3Bv9Zz2vJozHHhJVlyeXtdW/NRWyvfsQ+98bleFs941biuLK5HyH6Vz97wZAvspVv1g9vV/VB9Rpvp/6pz77iVozr+hjXe36L81wlzs8WZ+al8LZuvXmNY1U/no2nEgoRC2AvoPXZXptGDFwFiidSpBfAvZWTiP1oVt2+xgCVbfVPn22tiraeF8eeHedpjC1R94j93tr5CW+n/txKHFPF5LdszjVmb5/FiApLbx6yGPR2VR9Ur1FfJBDZHIgYq7omI85xlTg/W8R2M/NSeLtm2ZjjvbP5XJGnE4qYtL02Mei3JjUW172C0Yj331phxHayaiEQMRl7Y3V67y/2CruPf+/+VbQi8z5kSeXtRnyzRSwwVcv6qN2Yt407qUivaM32wZFr9vA+78VCjM0qcX56nJ2XwtvJ9oQxxvQsn1+FpxMKEYPUt4oiBsxWkt4y+bEPvRVb7MfWeLbwZyjo95IlPi/bITTiOGZQ9b/wdjOSLxaiiu0VIm+rYlrhbLE8cs0eHgt7sXo0buL89Dg7L4W32Xs/1/Br9oRlNZ5SKGIhiiuHaqGKK/C9laITAz6Klaj2Y4uKL5z4vEqCRR/MYGTc3m5GwYvzMurzyBGfirN9cOSaPUbizdvKqsT56XF2Xgpvc8TnlXxciacUChGPA/w8spqkR4NdVArISLHoUfVFo9KnSCVxeyh55XMltfrmdvaxS0Ycz6jPIyM7AyfOxWwfHLlGuyaNR+2VP3HeRo4h9XvvQ5U4Pz2O3ltEv2/5ptIm4v3a889qjHn5AlQnQ1vK1k7mxwfVJD07IKv92GI0MCt9ilQS11F77eD8mj3Lxu3tqgmbEccz6vPI0fvFuZjtg5FrJOpxYbVne/F2NHeiP3ucnZei0iYymo8rMeblCzAyGXGL2ia8mqRnB2S1H1uMBmalTxG18Wu2kAjHnULVsnF7u2rCZsTxjPo8cvR+cS5m+6B6TfzyR9X24u1o7kR/9jg7L0WlTWQ0H1dizMsXYGQyVLx6346oJunMgOx9Fa/ajy1GAzM+r5IAlcQV0VfaVWjMvfPjkXF7u2rCZsTxjPo8ojEeud/ZPqhcE/sg00tYHT+13bczEm8xHqpU4u3ovUUccy8vhbep+nzEP6sx5uULMDoZMfD0czVJR7/26MTn9p5R7ccWo76Iz6skQBxHD/nF26gvvULTGBm3t6v0d4/KvIxwxKfibB9UronxvfVitzESb95WViXOT4+z81J4m6rPR/yzGk8vFMLPzLWriNvtrWCJq8Wt1UePeAzTK5wjxaLHqC/i8yoJUEncavI1Rsbt7aoJmzHa1z3iuzAVsQojL8G9XfX+fs2W37xNJX5G4i0W8yqVeDs7L4W3qcbdaD6uxKsQiigM8WXrVpLGlXL8mu0WsRBuJXdsN1q0Rn0Rn1dJgErixjbZbkIcLZLVhM2IfR31eY8YT3t/5ay48iNRWdaPuNDZQ8dHfu8tv3kbvc/bw4vsXrxV4qZH5bqz81J4u2rcHalNq/AqhEJkL1qzJPXnyfb+kGakCMTAzfrRY9QX8XmVBKgkblzhZX/tPVokvV01YTPieEZ93iPeM5sLPS+OX5b1Ixb+rRW0BLoX51t+837sFdvYh714iz7ZO9ZqxOu2ODMvhbfb8l9kNB9X4tUIhYKiXRctC5i4epEpGVUMfeWsdgqoGIzZSi32KetHj1FfxOdVEqCSuL2vIsexyFc9/8iycXu7asJmxPGM+ryHxhbHpfnw4qj4iMXWLetHnDeZxtF2Lvrvlm9b2x6xP/o57gb1bI+zZnvxFneNyhn5QOgZyp/emOP8bHFmXgpvu+W/yGg+rsSrEQrh17rtFYt4dFW1vVVOLAB7/YiM+iI+r5IA1cSNX0WWaZWqfkW/K2n9OCUbt19XTdiMOJ4Ry3w8EiMaf+zH3tzHM//MFHfu8y2/SWB64tLmzOdI5ruVzBeNLeFq1utX9EvGWXkpvP2W/yLu84p/VuJVCcVWYlQKtNrExNkyPaPygi0W7ko/nFFfxOdVEmAkceMKtWfyoVZ53vds3H5tNWEz4nhGLUOFa684tpX16Nxrlew+65nHnbfN/Ka+7PVZpoWA97kSb3FXEa3Xr5F4E+rT7LwUfl017tznFf+sxHJC4cWo8gIuosTwwBqdUBUD9cFXVzLdU6s+/T5u37fwJNV/9fMI7ovKKikKZeXc2FdtGuMeSlz1xX2sZzbfNLzv2bj9Pn79UdQ/98GIVWJFPlZcenz0xu9FV6afK+geupePQc9SMfO484K7Vxx1ndp4oWv31Vg0JuHxqj5U0DUxHvRvfdbLk9F4a+g63XNGXoojcafntGuO1KYrs5xQADwDR4UC4BEgFAAPAKGAlUAoAB4AQgErgVAAPACEAlYCoQB4AAgFrARCAfAAEApYCYQC4AHc+tVogHuCUAAAQApCAQAAKQgFAACkIBQAAJCCUAAAQApCAQAAKQgFAACkIBQAAJCCUAAAQApCAQAAKQgFAACkIBQAAJCCUAAAQApCAQAAKQgFAACkIBQAAJCCUAAAQApCAQAAKQgFAACkIBQAAJCCUAAAQApCAQAAKQgFAACkIBQAAJCCUAAAQApCAQAAKQgFAACkIBQAAJCCUAAAQApC8Yr4+fPnnzdv3txk7969e7kbzODHjx//+PTDhw8vn1ybt2/fduNixH79+vVyN1gFhOIV8e3bt27ijtjHjx9f7ga3IIGQODS/riAUKvAeC0dMQgPrgVC8Ij5//vw3Yb9+/fry6b/47758+fLy6b+o7dbvYAwV20+fPv31Z7MVhOL79+9/+6sx/P79++U3//1dHIvvZFfZOcF/QShOREVBBVhHC1fYbrcVrJI80n4nU9JH2u+0Er4K6rN2OL3+Xg0VVYls86P67sc4KxTQ1v/e8aOPTTEfafHFQmNNEIoTkCh44ZVdoZipHypOvhJseF+z318F9VEFq/VL/76yYEjQ1E+JdFs0+FHgCkLRYrq3WPB4781Ddi1cH4RiMr6yUlHWz73Ce2+aeMUjJ6HkbX1+//79y6f/Rdf2diKPRsXWBUMF+Qr+jsjHsUi63+XfqyPfbsVAG4es53/tMjTGK84N7INQTERJ1JJFBbdy3KTE8STrWS85964beWm4d2wwyiPG5O9YJBw6F786s4RCsdbu07PeUZHYu67qw8pCA9YGoZiEi8TWqmsLbdX9vLqZPtOKeQsJUe86FZ2KSDXUvl2bPW+ER4zJj3J0jxEfPIJZQqGCvlX0sx2tPve5byZhibufDF9ojMY+rAFCMQFfzR5dUXmRa1ZZ3cdEP5KoXphnFtdHjMm/oaW52CqSV2Dm0ZPmrd2rWSUWvci3a0Z91t6/yGYtNOBaIBQ3opVzSxIV3FsKk5+1y5TEe/g1R0TKv7q4dURxC48YkxeuGUdpZzFTKITvaqv39GuO7sLOWmjAdUAobkCi4EXt1tWUr4ZleytpX7GrH0dEyp95ZDeyxyPGpGLlxat6jNLa32pVZguFi75M/suIu5Aj73XOXmjANUAobsC37TOSJCZutpr2Yqj/Hkly4SvK3jeibuURYxJ+HFgtwq39rVZltlAIxaH3JRNazUVrd3SR48KunRw8JwjFQZSAvmqdVWT92ESm4tnDk/yWvx/wwnJLYc6495hEFKjqruKenCEUcQe3JQAupLccz5290IBrgFAcxFdSsiNHJD3ifXvJ50leOfPfwoupRO8s7jkmxwXqiqvdM4QiCmRv3P5e7dbn3mOhAY8HoTiIEmxWsjkSnHZfWTyq8aJ76zsFLxhnFtJ7jsmJq+tZYj6LM4RC+M4sjlvFvO2E1e4Wn/i8nrnQgMeDUBwgFr5ZK+DG1lHNzCQXZ6zit7jXmBzd258pQboSZwnF1vGTfNtERD6/dQcwc2fS0D09VrRj2Zu3uIvya2EOCMUBPEFks8+/4/2V+DHJW6G9hXY/2ewxRO41pkgTIdnM3coMzhIK+bXdV9Z2i/4+QfNxKzMXGuqz9y/aXnwqduJiRKb7wu0gFAfwBDkrGL3AqZh6Et26EhRxFXZGkY6cPaYeKsDtGVdbYaq4tr7JHzOJRdN3Gfr3DPz9xC27NV8wyCf6Oe4Gm9hlxMVI5RqogVAcwIuPit8ZbK2uZhyfSBTi/c8+ehJnjmmLKOpXQEVQY3bhlGkOZgmm7u/3bqY5uBUV8ni8paJ8dMGk6zRPMQaj2O3dP4757F3yawKhOICvpGYeGThxdSSbkeRxJ+E2e1UbOWtMGb5qlz2yeHg/KnZLX1VU4/00v0eLuePxH23G/RtRjPaOy3yueT8xF4TiAB68ZwmFiAk5Mwkfxb3HdCWhuDdxBzfjvcQ9icdPcccR8Z3+mbvU1whCcQAP3rOEIiaJbPXgf8SYJAz3fN5VkABHUT5793YG3v8s1zTe1u6s4+DXDEJxgGrwHkVBH8+vZSu/nHvUmKJQ7K1KnwVfXTdbsYD6OLLjJH8/8Vrm+J4gFAdoASmbLRQqqO0bIErsmPD6/Wo8ckyvUSj8yCm+EF7t+CkeHW7RYkrxtWKOXB2E4gAeuLOFwpNcRc5XSrIVj04eOabXJhRxZR2/vLDa8VP8AoTmM+JjXDE/VgChOIAH7kyh8NVT+667Vkf+vLO/mTSbR48pFppnLiQuin6k13ZzzVZacUeh6/0NSFuI8E2n80AoDuCBO0sovKDFc/t4fKDkWYErjCkeXfRWpM+A/NfeAalguhjEr5muJpZtXLK4I3QhWe1YbSUQigP4t0lmCIW+DdSSQau/uOLz4wRZb1V1Na4yptcgFPJt2zXI5/K9E1flUbSvTnv/IIv51nYTM3f28P8gFAfwwJXdgpK8CU8vyYXa+PPOPKqZwZXG5O9HZM+I7862dgtNSJrJ/6vgYu/HS4qr9nkvxmAeCMUBZv6/nlx0sq3zSsdPVxqT98WLzLPgRVRxucXKx09+hClrtLnNxg1zQCgOEAP36HGGr3b3gj0e1Vw1Oa42Jr//akcue7j/9nZk8fjp6rtSx3cOMv3cxi7xX2l3tCoIxQHisUl8wVbBk7wS7DHRr7g6vtqYYoFR/54Fja29A5JVjl7kX/fHlXelEe+3Fmpt7NmOFeaBUBzEz3xHV2f+NUZZtYD5NbIrJckVxxSPW55l5alxuEhoF1chvltb6cgm9l32bDvEK4NQHMRXz7Lq6iyuBGWVo6u4+pYpea7AVcfkYv4sRUUi4eOSVXe0cUehOVtFPOOXElbq+zOAUBwkruoqySpxiQVVpiKWCY0KcSwOzZRAjzxCuOqYoggdfY90JTSGns8kANnRk3wRC20z3S+79ir4S3sZR073BaG4AQ/evfN1BbYHerStFa+SuNfe7cg7khlceUxeGFUMV0cLE/dPNIn1Fj0hd7vKzjTDY22lI7NnAaG4ASWvb+dX+EO418Az7iZeO20npf9y5HR/EIob8ZWOVm4E8ePxF5/P8m7iNaMdRMuvFY7JnhGEYgItkGXPcMyxMv5NJ1af6+NfGql+kw7mg1BMws/E9W+4P15UWH2uj7/L4r3EY0EoJhJfoN76zR2oE79YgEisjd4rtZfwHB8+HoRiMl6wFOj6meOP81BB8S8UqKjg77VoL6ol7po7zyGOD68BQnEC2kn4C1UJBt+8mY8LhP7Nd+vXQyLQ5jAaInEdEIoTkWDobFVFjICfj3YPMgRiXbb+Fod3EtcCoQCAhxG/gKD3fOy+rwdCAQAAKQgFAACkIBQAAJCCUAAAQApCAQAAKQgFAACkIBQAAJCCUAAAQApCAQAAKQgFAACkIBQAAJCCUAAAQApCAQAAKQgFAACkIBQAAJCCUAAAQApCAQAAKQgFAACkIBQAAJCCUAAAQApCAQAAKQgFAACkIBQAAJCCUAAAQApCAQAAKQgFAACkIBQAAJCCUAAAQApCAQAAKQgFAACkIBQAAJCCUAAAQApCAQAAKQgFAACkIBQAAJCCUAAAQApCAQAAKQgFAACkIBQAAJDw58//AIca+lVWu2QaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import Image\n",
    "Image(\"Normal Equations.PNG\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[89597.9095428 ],\n",
       "       [  139.21067402],\n",
       "       [-8738.01911233]])"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inv = np.linalg.inv(X.T@X)\n",
    "Multi = X.T@y\n",
    "theta = inv@Multi\n",
    "theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hypothesis(X,theta):\n",
    "    htheta = np.dot(X,theta)\n",
    "    return htheta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "ypred = hypothesis(X,theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "def r2_score(y,ypred):\n",
    "    ymean = np.mean(y)\n",
    "    num = np.sum((y-ypred)**2)\n",
    "    den = np.sum((y-ymean)**2)\n",
    "    output = 1-num/den\n",
    "    return output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7329450180289142"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r2_score(y,ypred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Verification using linear regression of sklearn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import datasets, linear_model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create linear regression object\n",
    "regr = linear_model.LinearRegression()\n",
    "\n",
    "# Train the model using the training sets\n",
    "regr.fit(X, y)\n",
    "# Make predictions using the testing set\n",
    "y_pred = regr.predict(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7329450180289142"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "regr.score(X,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Which one should we choose Gradient Descent or Normal equations?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gradient Descent\n",
    "\n",
    "1. We need to choose alpha.\n",
    "2. Works well when number of feature is very large.\n",
    "3. Need many iterations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Normal equations\n",
    "\n",
    "1. If n is very large, the computing (X^T X)^(-1) will be extremely expensive."
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}