pca_linear_regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Understanding linear regression when m < n -- when you do want model selectoin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us first generate data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAdxElEQVR4nO3df4wc533f8fd3Znb3flASSZGhJVISKZuRosq1ZTOSUrdpYcmAnQSW/jAMO2kjtCrUom7iOAYSpf0jSIECMWDEcYDWrWI5YVrHtiI7lWAYTl3ZrmGgonz6UYv6UUumSJEUKZ1J6ijd8W53Z779Y2b29o534u3tHHdn9vMCiLv9cbvP7Ox+9uH3eWYec3dERKRagkE3QEREiqdwFxGpIIW7iEgFKdxFRCpI4S4iUkHRoBsAsG3bNt+9e/egmyEiUiqPP/74z9x9+0q3DUW47969m6mpqUE3Q0SkVMzsyGq3XbAsY2ZfMrPXzOxg13Vbzew7ZvZC9nNLdr2Z2Z+Z2Ytm9mMze08xmyAiIr1YS839L4EPLrvuXuARd98LPJJdBvgQsDf7dw/whWKaKSIivbhguLv7D4DTy66+A9if/b4fuLPr+r/y1KPAZjO7oqjGiojI2qx3tswOdz+R/X4S2JH9vhM42nW/Y9l15zGze8xsysympqen19kMERFZSd9TIT09OU3PJ6hx9/vcfZ+779u+fcXBXhERWaf1hvurebkl+/ladv1x4Kqu++3KrhMRkYtoveH+MHBX9vtdwENd1/9mNmvmVmCmq3wjIiIXyVqmQn4F+D/AdWZ2zMzuBv4Y+ICZvQDcnl0G+BZwCHgR+HPg32xIqzfAtw+e4NWz84NuhohIIdYyW+bj7n6Fu9fcfZe73+/up9z9Nnff6+63u/vp7L7u7p9w97e7+zvdfaiOTJprtvn0A/+Xo6fnllz/2hvz/Ov//gRffezoKn8pIlIuI3VumR++8DO+/sQxHnpq6TDAE0fOAPDGfGsQzRIRKdxIhfuBl04v+ZmbOpyG+2yzfdHbJCKyEUYq3B89dAqAx4+coRUnneunsp777EI8kHaJiBRtZMJ95lyLZ0+c5bodlzDXjHnmlbMAzLdinnllBoDZBfXcRaQaRibcf/TSadzht257BwCPvZT24n98bIZW7AQGbyrcRaQiRibcHz10inoUcPsv7ODabZMcOJTW3aeOpD///q7NzDVVlhGRahiZcD/w0mluumozY7WQm/ds5bHDp4kT54kjZ7h2+yS7toyrLCMilTES4X52vsUzr8xw67WXA3Dznq28Md/m+ZNnefzIGfZds4VNjUhlGRGpjJEI96nDp0mcTrjfkv386mNHOTPX4r3XbGGiHqksIyKVMRLh/uih09TDgJuu3gzAzs3j7Nw8zgNT6RGp771mK5saIbPNNulJLkVEym0kwv3AoVO8++q03p67Zc9WFtoJWyZqvH37JBONCHc411LvXUTKr/Lh/sZ8i6ePz3Drnq1Lrr/l2vTye6/Zgpkx2UjXClfdXUSqoPLh/vzJN0gcbrp6y5Lrb9mT1t1/cXca8psaaa9eR6mKSBVEg27ARjuXDZJeOr50U3dvm+Sv/+UtndCfqKe3azqkiFRB5cN9oZ2eQ6Yehufd9g/esa3z+6aGwl1EqqPyZZlmFu6N2ltval5z15khRaQKKh/uC+20LNOILhDuddXcRaQ6RiDcs557dH5ZptukyjIiUiHVD/fWWnvueVlGPXcRKb/qh/uaa+55WUY9dxEpv5EJ93r41psahQGNKFC4i0gljEC4x0SBEV0g3CGtu2u2jIhUQfXDvZVcsN6em2yEmi0jIpVQ/XBvJzRqbz1TJjdZ1zndRaQaRiDc4x567hFzKsuISAWMQLj3UpaJeFNlGRGpgOqHeyu54AFMuU2NkDmVZUSkAqof7u34gnPccxP1SFMhRaQSRiDc116W0SLZIlIVIxLuayvLTNRD5pqx1lEVkdIbgXDvbbZMO/HOUa0iImVV/XBvJWuuuWvBDhGpir7C3cw+ZWbPmNlBM/uKmY2Z2R4zO2BmL5rZ18ysXlRj16PXsgzAnM4MKSIlt+5wN7OdwG8D+9z9RiAEPgZ8Bvicu78DOAPcXURD16uXskzec9egqoiUXb9lmQgYN7MImABOAO8HHsxu3w/c2edz9GWhnVDvoeYO6ChVESm9dYe7ux8HPgu8TBrqM8DjwOvunqfjMWDnSn9vZveY2ZSZTU1PT6+3GRfU7OkI1bQso6NURaTs+inLbAHuAPYAVwKTwAfX+vfufp+773P3fdu3b19vMy6ol5q7ltoTkaropyxzO/CSu0+7ewv4BvA+YHNWpgHYBRzvs43r1o4T4sTX3nOvK9xFpBr6CfeXgVvNbMLMDLgNeBb4HvCR7D53AQ/118T1W+sSezn13EWkKvqpuR8gHTh9Ang6e6z7gN8HftfMXgQuB+4voJ3r0gn3NZdlsnVUNRVSREouuvBdVufufwj84bKrDwE39/O4RVlopyG91rJMPQyIAlPPXURKr9JHqC60eivLmFm6jqrCXURKrtrh3mNZBvIzQ6osIyLlVvFw760sA/mZIdVzF5Fyq3i4995zn9Q53UWkAqod7j3W3CGdMaMTh4lI2VU73NdRlpnUUnsiUgEVD/f1Dqgq3EWk3Coe7usYUFVZRkQqoNrhvq6au3ruIlJ+1Q739ZRl6hHNdkIr1jqqIlJeFQ/39ZRlsgU7dCCTiJRYtcM9L8v0EO6b8gU7dCCTiJRYtcO9nRAGRhT2coRq3nNXuItIeVU83Ne+OHZOi2SLSBVUOtybPSyOnVtcsEM1dxEpr0qH+0IPi2PnJur5gh3quYtIeY1AuK99GiQslmVmF9q4O3/+g0P88IWfbUTzREQ2TMXDvfeae/c6qv/5+z/lP37rOR6YOroRzRMR2TB9LbM37BZaSU9Hp8LiOqrfePI4T778OpDW7kVEyqTa4b6Ossx4LSQwePLl17n12q28PteiqaNVRaRkVJZZxszYMlFn789t4r/+031MNqLOka4iImVR+Z57PkDai/3/4mZ2bh7nsoka9TBQWUZESqfa4d7qvSwDcOPOyzq/16OAuTlNixSRcql+WabHAdXl6lHQObukiEhZVDzcez+Iabl6FGhAVURKZwTCvfeyTLeGau4iUkLVDvdW77NllmvUFO4iUj7VDvd27wcxLVcPVXMXkfKpbLi344R24n2XZeqReu4iUj6VDfd8EFQDqiIyiiob7utZYm8l9TAkTpw48SKaJSJyUfSVfGa22cweNLPnzew5M/slM9tqZt8xsxeyn1uKamwv8jp5o9Z/WQZ08jARKZd+e+6fB77t7tcD7wKeA+4FHnH3vcAj2eWLLj8fTN+zZRTuIlJC604+M7sM+GXgfgB3b7r768AdwP7sbvuBO/tt5Hp0eu4FDKimj6eTh4lIefTTrd0DTAN/YWZPmtkXzWwS2OHuJ7L7nAR2rPTHZnaPmU2Z2dT09HQfzVhZs13cgCqg6ZAiUir9JF8EvAf4grvfBMyyrATj7g6sOBLp7ve5+z5337d9+/Y+mrGyvKfd6wLZy3XKMpoxIyIl0k/yHQOOufuB7PKDpGH/qpldAZD9fK2/Jq5PcbNlVHMXkfJZd/K5+0ngqJldl111G/As8DBwV3bdXcBDfbVwnTRbRkRGWb/nc/8t4MtmVgcOAf+c9AvjATO7GzgCfLTP51iX4mbLpF8OKsuISJn0Fe7u/hSwb4WbbuvncYuwUPSAakvhLiLlUf0jVIsqy8SaCiki5VHdcC+oLKMBVREpowqHu+a5i8joGoFw73MlJs2WEZESqm64t2LMoBZaX4+jg5hEpIyqG+7Z4thm/YW7ZsuISBlVPNz7K8lA92wZhbuIlEeFw73/xbFBs2VEpJyqG+6t/hfHBojCgMAU7iJSLtUN94LKMqB1VEWkfCoc7sWUZSCdTqmeu4iUSYXDPSks3OtRoJWYRKRUqhvurQLLMmGgI1RFpFSqG+7tuJABVUgPZFJZRkTKpMLhXmxZRuEuImVS2XBvthPqmi0jIiOqsuFeZM9dZRkRKZsKh3txUyHT2TIKdxEpj+qGe8GzZdRzF5EyqW64t4s5/QBoQFVEyqeS4Z4kTjMucrZMqAFVESmVSoZ7HsQqy4jIqKpkuOcLa2hAVURGVTXDPTsPTJFHqOrcMiJSJhUN92LLMprnLiJlU9Fwz3ruRZ5+IE5w9xVvP3DoFH/75LFCnktEpAiVDPf5omvuYYA7tJOVw/2//O+f8tm/+0khzyUiUoRKhnunLFMr7twysPpSe4dPzWnAVUSGSkXDvfiyDKwc7u044ejpORZaGnAVkeFR0XAvtiyTD8yu1Ds/duYc7cSZ12waERkilQz3M7NNAC4ZqxXyeG/Vc3/p1CwArdiJV6nJi4hcbJUM94PHzzJWC9izbbKQx+uEe3x+7/zIz2Y7v2suvIgMi77D3cxCM3vSzL6ZXd5jZgfM7EUz+5qZ1ftvZm8OHp/hhisuJQyskMerh+nLtFJZ5vCpuc7v+SwdEZFBK6Ln/kngua7LnwE+5+7vAM4AdxfwHGuWJM4zr8zwzp2XFfaYjbcqy3T13Oc1qCoiQ6KvcDezXcCvAl/MLhvwfuDB7C77gTv7eY5evXRqltlmzI0Fhvtb1dwPn5rt/A9B4S4iw6LfnvufAr8H5Kl3OfC6u7ezy8eAnSv9oZndY2ZTZjY1PT3dZzMWHTw+A1BouHd67stO+9uKE46dOdep7assIyLDYt3hbma/Brzm7o+v5+/d/T533+fu+7Zv377eZpzn4PEZGlHA3p/bVNhj5j33hWXhffT0HHHiXP+2S9LbNaAqIkOin577+4APm9lh4Kuk5ZjPA5vNLMrusws43lcLe/T08Rmuv+JSorC4iUD1VXruR7LB1F+44lJAPXcRGR7rTkB3/wN33+Xuu4GPAd91998Avgd8JLvbXcBDfbdyjZLEeeb4Wd6589JCHzefLbO85p4Ppl63I+2560AmERkWGzHP/feB3zWzF0lr8PdvwHOs6OXTc7yx0C50pgysPqB6+NQslzQirtw8DqBTEIjI0IgufJcLc/fvA9/Pfj8E3FzE4/bq6Www9e9duTHhvhCf33PfvW2Ssdrq8+BFRAahUkeoHnxlhnoY8PNZmaQo+bllVuq5p+Ge3q6pkCIyLKoV7sdnuO5tl3R62kXJp0J2z4ZpthOOnznHnssnusJdPXcRGQ6VCXd35+Dxs4XOb8+tNKB69MwcicPubZOd8FfPXUSGRWXC/diZc8ycaxU+mAoQBEYU2JJwP5zNlFlallHPXUSGQ2XC/amjrwNwY8HTIHP1ZYtk59Mg91w+SRgYtdB0EJOIDI3KhPtDTx1n+yWNzgFFRcsXyc4dPjXLZeM1tkymJ70ci0L13EVkaFQi3E/MnOO7z7/GR/ftolbgkandGst67kdPn+OayycWb6+FOohJRIZGJcL9gR8dI3H42C9evWHPUY+CJfPYT8822bap0bnciAINqIrI0Ch9uMeJ87Ufvcw/2ruNq7ZOXPgP1qkeLu25n55tsnlicRm/sVqgg5hEZGiUPtx/8JNpXpmZ59dv3rheO0A9CpeE95m5JlsnFheZGquFOv2AiAyN0of7Xz/2Mts2Nbj9hh0b+jzdA6rzrZi5ZtwZTIU03DWgKiLDotThfnJmfsMHUnONMKCZDZiemWsCsLUr3FVzF5FhUupwf2DqKHHiGzqQmuue535mtgXAluVlGdXcRWRIFHJWyEH59VuuZve2Sa6+fOMGUnONKODMXBbuK/Tcx2rquYvI8Ch1z33bpgYffteVF+W5unvup2fzcO+aLRNpnruIDI9Sh/vF1D2gmvfcN3eVZRoaUBWRIaJwX6Puee55z33z+GLPXQOqIjJMFO5rtHRAtcll47Uli3BrQFVEhonCfY2W1NznWksGUyEdUG22E5LEB9E8EZElFO5r1Og6QvX1uSZbuk49AHTO6a7eu4gMA4X7GuUDqu7O6dnmeT13rcYkIsNE4b5GeXg344Qzs80lBzCBeu4iMlwU7mvUvY7q6bnmkvPKQFpzB/XcRWQ4KNzXqJ713GfOtZhvJef33KNsHVUdyCQiQ0DhvkZ5uL96dh5YenQqQKPTc1dZRkQGT+G+RnnN/cRMGu6r9txVlhGRIaBwX6PFnvsCwPmzZTSgKiJDROG+RvmAal6W0YCqiAwzhfsa5T33k6uVZWoqy4jI8FC4r1En3M/OYwaXjS8bUM1uX9CAqogMAYX7GjW6ZstsHq8RBrbk9sWDmNRzF5HBW3e4m9lVZvY9M3vWzJ4xs09m1281s++Y2QvZzy3FNXdwGtlsmJMz8+fV26G7LKOeu4gMXj899zbwaXe/AbgV+ISZ3QDcCzzi7nuBR7LLpZeXZRbaCVsnVgh3nVtGRIbIusPd3U+4+xPZ728AzwE7gTuA/dnd9gN39tvIYVDvOnf7Sj33KAwIA9MRqiIyFAqpuZvZbuAm4ACww91PZDedBHas8jf3mNmUmU1NT08X0YwNlffcgRV77pD23lWWEZFh0He4m9km4OvA77j72e7b3N2BFVevcPf73H2fu+/bvn17v83YcN3hvnnZqQdy6WpM6rmLyOD1Fe5mViMN9i+7+zeyq181syuy268AXuuvicNhTT13LZItIkOin9kyBtwPPOfuf9J108PAXdnvdwEPrb95w6MRvXXNPb+PBlRFZBhEffzt+4B/BjxtZk9l1/074I+BB8zsbuAI8NH+mjgcugdUV+u5N9RzF5Ehse5wd/cfArbKzbet93GHlZlRD9Ol9lbruY/VAtXcRWQo6AjVHuR19+VnhMyNRaFOPyAiQ0Hh3oM83LdMrDZbJtA8dxEZCgr3HtTDgMDg0rGVw70RhRpQFZGhoHDvQaMWsGWiThCsPNQwVtNBTCIyHBTuPaiHwaqDqaCDmERkePQzFXLk1KOAyfrqL5kOYhKRYaFw78GH33Ulk43VXzIdxCQiw0Lh3oN/9Y/f/pa3N2ohC+0Edyc9gFdEZDBUcy9Qvkj2QlulGREZLIV7gcay1Zp0IJOIDJrCvUCdpfY0Y0ZEBkzhXqCGltoTkSGhcC+QFskWkWGhcC/Q4oCqeu4iMlgK9wKp5y4iw0LhXiDV3EVkWCjcC7TYc1e4i8hgKdwLpIOYRGRYKNwL1IjUcxeR4aBwL1Aj67nPq+cuIgOmcC9QXnNfUM9dRAZM4V6gMZVlRGRIKNwLVAuNwDSgKiKDp3AvkJllqzGp5y4ig6VwL1i6GpN67iIyWAr3gqnnLiLDQOFesLFsqT0RkUFSuBdMi2SLyDBQuBesUQt1EJOIDJzCvWBj6rmLyBBQuBdMNXcRGQYK94KN1QKdfkBEBm5Dwt3MPmhm/8/MXjSzezfiOYZVI9JUSBEZvKjoBzSzEPhPwAeAY8CPzOxhd3+26OcaRmO1jT+Iyd2B9IjY7uviJL0+DGzJbYOUJI4ZQ9Oefrg77qy4Pe5O4hCscVvdnWacvk8CM0IzHEiy54gCIwiG4zXLtxuWbnt+/bC082JZ6fO3kiTxgb42hYc7cDPworsfAjCzrwJ3ACMR7uO1kJNn57npP/xPABKHVpzQbCfE7kzUQiYaEeO1kCQL5DhxEnfaiZMkThgYtTCgFgbp38YJrXZCK1m8P6TnsgkD6zxH/gGELBzMcLIPoBm10KhFAYFZ+pyJ46Qf2DB7E7ZjpxUnJO6YpefKMawTOvnjdT0V1vUYtSAgDI127Cy0Y1rxYlvrYfofxdidJEnvP14PGYsCgsBothNaWeA1opBGLX0Nkuz1yf+unSQkvvjcDrTjhHacvoZRkG5nFBj5588d2onTbCedUA0tff26/8X5fdoJ7STpbGf3awtQDwPqUYAZNNvJknGWWmhEQZBuc5RuQxgYUfYav7nQ5uy5dqcdq4mC9O/D7hCx9PrlX+BG+noG2TYFlu7zxL3Tvmacvr6tePE9ZJa+DhP1kMlGxFgtZK7ZZm4hZq4Vd+73VuphwGQjZKKexslCO2ahlb7fu9uXPp+lHZFsX5ql41TjtbDzeub7tNPu7H3RThx3Z6IeMVEPOwcMnmvGzLdjAlt8vSH9osw/T2G2P4LsNQsCSBJoxgntOOn8bSMKiN0514yZa8a0YycI0tcyfQ8lnff0ZPaajddD3Ol8nueaMXPNNq3YqUdB2tYoxFn8/IZB+ly10PjUB36eO96984Kvc682Itx3Ake7Lh8Dbll+JzO7B7gH4Oqrr96AZgzGx29JtyV/Wxt03nCBGeda6Y4/14wJAlsSMFH2gY0Tp50kNNu+JCBqYbDYo8u+DNqJp2/M0IjCAANaidOO0wDMPyyJpwHYitMPXWjp4ywGdxr0eSgFQfpmTnuk3rmvWdo7TbfNloR9nHgWsAlRENCoBdTDoPNBbWYBGAZpby9JnHOtmPlWQpL4kg9ms50w384/XGlY5W1Ofy593dN2B0TZF0srXvyiyOWvYS1KNyBJnDgh+2JNvxzCLFDrUUAtWAwbSIPJLH0t8+1J3GnUAhpRSGhGnKQhlIdo9xd7O05fy0vGalw6HnFJI8IsfR1iT/dj/kXZ/bfd+Zr4Ykeg+wsn/59b7Itf2okvvv8a2WtbzwKl+wujnQXSmwtt5ltxJ+gn6iFhEHS+4Dv7OvvvS/4+mG8lzC60mV1oZ+dXSp8n/zLL29n9NZF/Ebl7GtCtmGY7/TLN/3fUyPdDmL6PwuzL+lwzST9DrZjxWph2EGohSeKdLzDDOv+Lyl+v/DORd1TSTlT6Rdz9JZh3OiZqIVEYdD4D+XsoCgNwZ7YZ8+Z8u/PF0v1FOdGIqIcB8+30y+dcM17Siejev5dPNlbNk35sRLivibvfB9wHsG/fvgt3D0ri+rddyh/dceOgmyEiI24jBlSPA1d1Xd6VXSciIhfJRoT7j4C9ZrbHzOrAx4CHN+B5RERkFYWXZdy9bWb/Fvg7IAS+5O7PFP08IiKyug2pubv7t4BvbcRji4jIhekIVRGRClK4i4hUkMJdRKSCFO4iIhVk7oM/fsjMpoEj6/zzbcDPCmxOWYzido/iNsNobvcobjP0vt3XuPv2lW4YinDvh5lNufu+QbfjYhvF7R7FbYbR3O5R3GYodrtVlhERqSCFu4hIBVUh3O8bdAMGZBS3exS3GUZzu0dxm6HA7S59zV1ERM5XhZ67iIgso3AXEamgUof7KCzEbWZXmdn3zOxZM3vGzD6ZXb/VzL5jZi9kP7cMuq1FM7PQzJ40s29ml/eY2YFsf38tO6V0pZjZZjN70MyeN7PnzOyXRmRffyp7fx80s6+Y2VjV9reZfcnMXjOzg13XrbhvLfVn2bb/2Mze0+vzlTbcuxbi/hBwA/BxM7thsK3aEG3g0+5+A3Ar8IlsO+8FHnH3vcAj2eWq+STwXNflzwCfc/d3AGeAuwfSqo31eeDb7n498C7S7a/0vjazncBvA/vc/UbSU4V/jOrt778EPrjsutX27YeAvdm/e4Av9PpkpQ13uhbidvcmkC/EXSnufsLdn8h+f4P0w76TdFv3Z3fbD9w5mBZuDDPbBfwq8MXssgHvBx7M7lLFbb4M+GXgfgB3b7r761R8X2ciYNzMImACOEHF9re7/wA4vezq1fbtHcBfeepRYLOZXdHL85U53FdaiLv4JcSHiJntBm4CDgA73P1EdtNJYMeAmrVR/hT4PSBf5fpy4HV3b2eXq7i/9wDTwF9k5agvmtkkFd/X7n4c+CzwMmmozwCPU/39Davv277zrczhPlLMbBPwdeB33P1s922ezmetzJxWM/s14DV3f3zQbbnIIuA9wBfc/SZglmUlmKrta4CsznwH6ZfblcAk55cvKq/ofVvmcB+ZhbjNrEYa7F92929kV7+a/zct+/naoNq3Ad4HfNjMDpOW295PWovenP23Haq5v48Bx9z9QHb5QdKwr/K+BrgdeMndp929BXyD9D1Q9f0Nq+/bvvOtzOE+EgtxZ7Xm+4Hn3P1Pum56GLgr+/0u4KGL3baN4u5/4O673H036X79rrv/BvA94CPZ3Sq1zQDufhI4ambXZVfdBjxLhfd15mXgVjObyN7v+XZXen9nVtu3DwO/mc2auRWY6SrfrI27l/Yf8CvAT4CfAv9+0O3ZoG38h6T/Vfsx8FT271dIa9CPAC8A/wvYOui2btD2/xPgm9nv1wKPAS8CfwM0Bt2+DdjedwNT2f7+H8CWUdjXwB8BzwMHgf8GNKq2v4GvkI4ptEj/l3b3avsWMNLZgD8FniadSdTT8+n0AyIiFVTmsoyIiKxC4S4iUkEKdxGRClK4i4hUkMJdRKSCFO4iIhWkcBcRqaD/DzkCyUsNQk+tAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "m_train, m_test = 60, 50  #number of samples\n",
    "\n",
    "n = 100 #number of features\n",
    "\n",
    "theta_true = np.linspace(10, 20, n)\n",
    "\n",
    "mu_x = np.linspace(1, 10, n)\n",
    "approx_rank = 12 ##approx rank because first approx_rank singular values are high. \n",
    "eig_vals = np.concatenate((100 * np.random.rand(approx_rank), 0.1*np.random.rand(n - approx_rank)))\n",
    "\n",
    "#eig_vals is a full rank matrix, but most of its energy is only along \"approx_rank directions\"\n",
    "tmp = np.random.randn(n,n)\n",
    "#eig_vals = 100 * np.random.rand(n)\n",
    "\n",
    "\n",
    "#matlab has a sqrtm function that does matrix square root\n",
    "\n",
    "X_train = np.random.randn(m_train, n) @ np.sqrt(np.diag(eig_vals)) @ tmp.T + mu_x.T ##matrix square root and element-wise square roots are completely different things \n",
    "Y_train = X_train @ theta_true + 1e-3 * np.random.randn(m_train)\n",
    "\n",
    "X_test = np.random.randn(m_test, n) @ np.sqrt(np.diag(eig_vals)) @ tmp.T + mu_x.T\n",
    "Y_test = X_test @ theta_true + 1e-3 * np.random.randn(m_test)\n",
    "\n",
    "plt.plot(eig_vals)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Linear Regression -- all this code is from the second recitation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Final Loss is  137.22443404647188\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAc10lEQVR4nO3df7xVdZ3v8debwy9/gCIcKvkh4JBJaeEcMctGqmti40iXsQT7LY8H2X3YVE4zYc2tO800OXWnpq7ONUryR4Wa2h0ySplS0UYTRDCQkBPkeEgFBH8XCXzuH991OJvN2etszjl7r33OeT8fj/VYe3333mt9vudxOG/Wr+9SRGBmZlbJoKILMDOzxuagMDOzXA4KMzPL5aAwM7NcDgozM8s1uOgCamHMmDExadKkosswM+tTHnzwwR0R0Vze3i+DYtKkSaxataroMszM+hRJj3XW7kNPZmaWy0FhZma5HBRmZpbLQWFmZrkcFGZmlstBYWZmuRwUZmaWy0FR6vrr4ZvfLLoKM7OG4qAodeON8K1vFV2FmVlDcVCUGjkSnnuu6CrMzBqKg6KUg8LM7CAOilIjRjgozMzKOChKjRwJv/89vPxy0ZWYmTUMB0WpkSPT/Pnni63DzKyBOChKOSjMzA7ioCjVHhQ+T2Fmtp+DopSDwszsIA6KUg4KM7ODNPyjUCW9BXgvqdZpEfGmmm3MQWFmdpBC9igkLZa0TdK6svZZkjZKapW0ECAi7omIi4HbgGtrWtiIEWnuoDAz26+oQ0/XALNKGyQ1AVcC5wDTgHmSppV85ELg+zWtynsUZmYHKSQoImIFsLOseQbQGhGbI+KPwA3AbABJE4FnI6K2160eeWSaOyjMzPZrpJPZ44DHS5bbsjaA+cB38r4saYGkVZJWbd++vXsVDBrkYTzMzMo0UlBUFBGfj4j/7OIziyKiJSJampubu7+xkSN9w52ZWYlGCoqtwISS5fFZW315BFkzswM0UlCsBKZKmixpKDAXWFr3KkaNgqefrvtmzcwaVVGXxy4B7gNOkNQmaX5E7AEuAW4HNgA3RcT6uhfX3AzdPcdhZtYPFXLDXUTMq9C+DFhW53IONHYs3H9/oSWYmTWSRjr01BjGjoUdO2DfvqIrMTNrCA6KcmPHwt69sGtX0ZWYmTUEB0W59ktrt20rtg4zswbhoCj3ilek+RNPFFuHmVmDcFCUO/74NG9tLbYOM7MG4aAoN2ECDBsGmzYVXYmZWUNwUJQbNAj+5E/g0UeLrsTMrCE4KDozfTr88pcQUXQlZmaFc1B05owz4KmnfJ7CzAwHRefOPDPNly8vtg4zswbgoOjMCSfAiSfCjTcWXYmZWeEcFJ2RYO5cuOceeOyxoqsxMyuUg6KSD384XQH1jW8UXYmZWaEcFJVMmADveQ9861uws/zx3mZmA4eDIs9ll8GLL8IXvlB0JWZmhXFQ5DnpJJg/H668EtbX/xlKZmaNwEHRlX/8x/R41Pe9D3bvLroaM7O6c1B0ZexYuPpqWLMGLr3Ud2ub2YDjoKjGX/wF/M3fwL/9G3zta0VXY2ZWV4U8M/tQSJoJ/AOwHrghIu4qpJDLL4ctW+Cv/xqGDoVLLimkDDOzeitkj0LSYknbJK0ra58laaOkVkkLs+YAXgCGA231rnW/QYPgu9+F2bPhYx+DL37Rh6HMbEAo6tDTNcCs0gZJTcCVwDnANGCepGnAPRFxDvBp4O/rXOeBhg2DH/wA3vte+Lu/gwsvhJdeKrQkM7NaKyQoImIFUH4X2wygNSI2R8QfgRuA2RGxL3t/FzCs0jolLZC0StKq7du316RuAIYMgeuvT4eibrwRTjsN1q6t3fbMzArWSCezxwGPlyy3AeMkzZH0TeB64IpKX46IRRHREhEtzc3Nta1Ugk9/Gn76U9ixA049NQXHnj213a6ZWQEaKSg6FRG3RsRHIuKCwk5kV/KOd8CvfpXOW1x2GbS0wC9+UXRVZma9qpGCYiswoWR5fNbW2MaMgZtuSucudu5MDz364Adha+OXbmZWjUYKipXAVEmTJQ0F5gJLC66pOhKcfz5s2JD2LJYsSc/d/tSn0qEpM7M+rKjLY5cA9wEnSGqTND8i9gCXALcDG4CbIqJvDbB0xBHwT/8EGzfCBRekm/OmTIHPfQ5qeYLdzKyGFP3wXoCWlpZYtWpV0WXAI4+kkLjlFjjsMLjoojQMyJQpRVdmZnYQSQ9GREt5eyMdeup/pk2Dm29OgTFvHixaBFOnpudc3HWXb9gzsz7BQVEPJ56YBhb87W/TECDLl8Nb35qC5Otfh127iq7QzKwiB0U9HXssfPnL6Yqo73wHjjoKPvEJGDcO3v9+uOMO34thZg3HQVGEww+HD30I7r8fVq+GD3wAbrsNzj47PYL10kvhoYd8aMrMGoKDomjTp8NVV8ETT6ST3qefDldcAaeckg5ZfeYzsGqVQ8PMCuOgaBTDh8OcOXDrrfDkkyk8JkxIh6pOPRWOOw4+/nG4+24fnjKzuvLlsY1u50740Y9SgNxxB/zhD3D00XDWWTBrVpqOPbboKs2sH6h0eayDoi954YU0EOFPfpLmv/tdaj/55BQYZ5+dDl0ddlixdZpZn+Sg6G8i0oCE7cFx773pkNTQoWno85kz0+TgMLMqOSj6u+efT+cv7r473cy3ejXs25eCY8YMOPNMeNObUoiMHl10tWbWgBwUA82zz6Yhz++6qyM49u5N702dCm98Y5pOPx1OOgkGN/zj082sxhwUA90LL8CDD6Z7N+6/H+67D556Kr13+OHpctzp09N0yinprvEhQ4qt2czqqlJQ+L+RA8WRR6bDT2eemZYj4L/+KwXG/ffDypWweDG8+GJ6f+hQeN3rOsJj+vS0PHJkcX0ws0J4j8I67N0Lra3prvDVq9P8oYfg6ac7PjNhArz2tQdO06alIDKzPs17FNa1piY44YQ0zZ2b2iKgrS0Fxrp1sH59mu68E3bv7vjuxIkpNE48EV796nQeZOrUNI7VIN/XadaXOSgsn5T2IiZMgPPO62jfuxc2b+4Ijvbp5z8/MEAOOwyOP/7A8GifXvnKtH4za2gOCuuepqaOP/jveldH+759aQ9k0yZ49NE037QphciPfgQvv9zx2eHDYdKkNE2e3PG6fXnMGAeJWQNwUFjvGjQoHYaaOBHe/vYD39uzBx57LAVHa2t6PseWLWn+wANpuJJShx/eERrHHZf2asaP75jGjfPNhGZ10PBBIelE4OPAGOBnEfF/Cy7Jumvw4HQY6vjjO3//uedSaJRPW7ake0Keeebg74wefWB4lE/jxsGIETXrktlAUEhQSFoMnAtsi4jXlbTPAr4ONAHfjojLI2IDcLGkQcB1gIOivxo5Mo1bdfLJnb//4ovpoU9tbZ1PDzwA27cf/L0jjkjnQ171qvx5c3M6pGZmB+gyKCS9GfhfwHHZ5wVEREzpwXavAa4g/eFv304TcCVwFtAGrJS0NCIekXQe8FHg+h5s0/q6I45IJ8Vf/erKn/nDH9JgiaUB8uST6XkfTz6ZxsdavjzduV5u0CAYO/bgAHnlK1N76XTMMQ4VGzCq2aO4Gvgk8CCwtzc2GhErJE0qa54BtEbEZgBJNwCzgUciYimwVNKPge93tk5JC4AFABMnTuyNMq0vGj4cpkxJU56XXkp3prcHSOm8/fWaNekz+/Yd/H0pnWxvD47m5gODpHS5uTkNDe8T89ZHVRMUz0bET2peCYwDHi9ZbgNOkzQTmAMMA5ZV+nJELAIWQbrhrnZlWr9w+OHpJPnkyfmf27s33XC4fTts29YxlS+vWZPmnZ1HgTQcSnt4lM7HjEnnWcaMOfD16NHp7nizBlBNUNwp6SvArcD+C+QjYnXNqioREXcBd9VjW2YHaWrq2DN47Wu7/vwf/5hCpDRIOguZTZtgx440BlclI0YcHCBdvXa4WA1UExSnZfPS27oDeFsv17IVmFCyPD5rM+s7hg5NV1qNG1fd53fvTnssO3Z0zCu9/vWv0+vnn6+8vvZwyQuTY445cDrySB8Ws1xdBkVEvLUehQArgamSJpMCYi5wYZ22bVaMYcPSo2wP5XG27eHSVbDs2AEbN6Z5XrgMGXJweLRPnQVLe/uIEQ6YAaKaq56OAj4P/FnWdDfwhYjo5LKR6khaAswExkhqAz4fEVdLugS4nXR57OKIWN/dbZj1Wz0Jl507D5w6a3v8cVi7Nr3XPppwZ5qaDj1cjjkmXQbt8b/6lC5Hj5V0C7AOuDZrej/w+oiYU+Paus2jx5r1kt27Ydeu6gKmtD1vD2bQIBg1qvMQGTUqf/Kd+DXVk9Fjj4+IvyxZ/ntJa3qvNDNrWMOGddxLcihefvnggKkULtu3d5x/ee65ruupJlAcMr2qmqD4vaQzIuJe2H8D3u9rW5aZ9WlDhnRcLXYo9u5Nlxjv2lXdtHVrGv5+587qQqZ9T8Yhc0iqCYqPAtdm5yoE7AQ+VMuizGyAampKh6FGjz7073Y3ZHbt6vxO/VIDPGSqueppDfB6SSOz5S5i28ysAI0eMkcfffC80uv2+VFHNcRQMRWDQtL7IuK7ki4tawcgIr5a49rMzOqjpyHz7LPp8FelYHnmmY4g2r49PaulvW1vFyMjjRx5aOHyp3+aRh7oRXl7FEdk887GaPYQGWZmcOBlwocqIt2dXxokpfPO2rZsSc+0f+aZzq8u27gxf+DMbqgYFBHxzezlf0TEL0rfy05om5lZT0jpxsURI9KDuQ7Vnj1pb6Y0SLqzni5UczL7/wCnVNFmZmb1NHhw9w+ZHcpmKr0h6XTgTUBz2XmKkaQ7p83MbADI26MYChyZfab0PMVzwPm1LMrMzBpH3jmKu4G7JV0TEY/VsSYzM2sg1YzM9W1JR7cvSBol6fYa1mRmZg2kmqAYExH7H9sVEbuAQ7wv38zM+qpqgmKfpP0PoZZ0HL6PwsxswKjm8tjPAvdKups01tNbgAU1rcrMzBpGNWM9/VTSKcAbs6ZPRMSO2pZlZmaNopo9CoBhpFFjBwPTJBERK2pXlpmZNYpqHoX6z8AFwHpgX9YcgIPCzGwAqGaP4l3ACRGxu9bFdEbSFNJ5kqMiwjf6mZnVWTVXPW0GhvTmRiUtlrRN0rqy9lmSNkpqlbQQICI2R8T83ty+mZlVr5o9ipeANZJ+Buzfq4iIv+rBdq8BrgCua2+Q1ARcCZwFtAErJS2NiEd6sB0zM+uhaoJiaTb1mohYIWlSWfMMoDUiNgNIugGYDTgozMwKVM3lsdfWoxBgHPB4yXIbcJqk0cAXgemSLouIL3X2ZUkLyO7vmDhxYmcfMTOzbqjmqqctdHIndkRMqUlFB2/naeDiKj63CFgE0NLS4jvHzcx6STWHnlpKXg8H3g1045l/XdoKlD6aaXzWZmZmBeryqqeIeLpk2hoR/wr8eQ1qWQlMlTRZ0lBgLr18bsTMzA5dNYeeSh95Ooi0h1HtHd2V1rkEmAmMkdQGfD4irpZ0CXA76Ql6iyNifU+2Y2ZmPVfNH/x/KXm9B9gCvKcnG42IeRXalwHLerJuMzPrXXnPzP54RHwd+J8RcW8dazIzswaSd47iw9n8G/UoxMzMGlPeoacNkjYBx0p6uKRdQETEybUtzczMGkHFoIiIeZJeSTq5fF79SjIzs0aSezI7Ip4EXl+nWszMrAFVM3qsmZkNYA4KMzPL1WVQSHp3NW1mZtY/VbNHcVmVbWZm1g/l3XB3DvBOYJyk0nspRpLu0DYzswEg76qn3wGrSJfGPljS/jzwyVoWZWZmjSPvPoq1wFpJ34+IlwEkjQImRMSuehVoZmbFquYcxXJJIyUdA6wGviXpazWuy8zMGkQ1QXFURDwHzAGui4jTgLfXtiwzM2sU1QTFYEmvIg0tfluN6zEzswZTTVB8gTTe028iYqWkKcCm2pZlZmaNossHF0XED4AflCxvBv6ylkWZmVnjqObO7PGSfihpWzbdIml8PYozM7PiVXPo6TvAUuDYbPpR1mZmZgNANUHRHBHfiYg92XQN0FzjuvaTNEXS1ZJurtc2zcysQzVB8bSk90lqyqb3AU9Xs3JJi7PDVevK2mdJ2iipVdLCvHVExOaImF/N9szMrPdVExQXkS6NfRJ4Ajifjudpd+UaYFZpg6Qm4ErgHGAaME/SNEknSbqtbBpb5XbMzKxGqrnq6TG6+SjUiFghaVJZ8wygNbt6Ckk3ALMj4kvAud3ZTraeBcACgIkTJ3Z3NWZmVqaaq56ulXR0yfIoSYt7sM1xwOMly21ZW6Xtj5Z0FTBdUsXhzSNiUUS0RERLc3PdTqGYmfV7Xe5RACdHxDPtCxGxS9L0GtZ0gIh4Gri4XtszM7MDVXOOYlA2aiwA2eCA1QRMJVuBCSXL47M2MzNrQNX8wf8X4D5J7Xdnvxv4Yg+2uRKYKmkyKSDmAhf2YH1mZlZDXe5RRMR1pJFjn8qmORFxfTUrl7QEuA84QVKbpPkRsQe4hDR+1AbgpohY390OmJlZbVV1CCkiHgEeOdSVR8S8Cu3LgGWHuj4zM6u/as5RmJnZAOagMDOzXA4KMzPL5aAwM7NcDgozM8vloDAzs1wOCjMzy+WgMDOzXA4KMzPL5aAwM7NcDgozM8vloDAzs1wOCjMzy+WgMDOzXA4KMzPL5aAwM7NcDgozM8vloDAzs1wNHxSSTpR0laSbJX206HrMzAaamgaFpMWStklaV9Y+S9JGSa2SFuatIyI2RMTFwHuAN9eyXjMzO1it9yiuAWaVNkhqAq4EzgGmAfMkTZN0kqTbyqax2XfOA34MLKtxvWZmVmZwLVceESskTSprngG0RsRmAEk3ALMj4kvAuRXWsxRYKunHwPc7+4ykBcACgIkTJ/ZK/WZmVuOgqGAc8HjJchtwWqUPS5oJzAGGkbNHERGLgEUALS0t0RuFmplZMUFxSCLiLuCugsswMxuwirjqaSswoWR5fNZmZmYNqIigWAlMlTRZ0lBgLrC0gDrMzKwKtb48dglwH3CCpDZJ8yNiD3AJcDuwAbgpItbXsg4zM+u+Wl/1NK9C+zJ8qauZWZ/Q8Hdmm5lZsRwUZmaWy0FhZma5HBRmZpbLQWFmZrkcFGZmlstBYWZmuRwUZmaWy0FhZma5HBRmZpbLQWFmZrkcFGZmlstBYWZmuRwUZmaWy0FhZma5HBRmZpbLQWFmZrkcFGZmlqvhg0LSTEn3SLpK0syi6zEzG2hqGhSSFkvaJmldWfssSRsltUpa2MVqAngBGA601apWMzPr3OAar/8a4ArguvYGSU3AlcBZpD/8KyUtBZqAL5V9/yLgnoi4W9IrgK8C761xzWZmVqKmQRERKyRNKmueAbRGxGYASTcAsyPiS8C5OavbBQyrRZ1mZlZZrfcoOjMOeLxkuQ04rdKHJc0BzgaOJu2dVPrcAmABwMSJE3ulUDMzKyYoDklE3ArcWsXnFgGLAFpaWqLWdZmZDRRFXPW0FZhQsjw+azMzswZURFCsBKZKmixpKDAXWFpAHWZmVoVaXx67BLgPOEFSm6T5EbEHuAS4HdgA3BQR62tZh5mZdV+tr3qaV6F9GbCslts2M7Pe0fB3ZpuZWbEcFGZmlstBYWZmuRwUZmaWy0FhZma5HBRmZpbLQWFmZrkafqynevrIR+Dee4uuwsys+1asgNGje3edDooSEyfCtGlFV2Fm1n2Da/BX3UFR4rOfLboCM7PG43MUZmaWy0FhZma5HBRmZpbLQWFmZrkcFGZmlstBYWZmuRwUZmaWy0FhZma5FBFF19DrJG0HHuvm18cAO3qxnL7AfR4Y3Of+r6f9PS4imssb+2VQ9ISkVRHRUnQd9eQ+Dwzuc/9Xq/760JOZmeVyUJiZWS4HxcEWFV1AAdzngcF97v9q0l+fozAzs1zeozAzs1wOCjMzy+WgKCFplqSNklolLSy6nu6StFjSNknrStqOkbRc0qZsPiprl6RvZH1+WNIpJd/5YPb5TZI+WERfqiVpgqQ7JT0iab2kj2ft/bbfkoZLekDS2qzPf5+1T5b0y6xvN0oamrUPy5Zbs/cnlazrsqx9o6Szi+lRdSQ1SXpI0m3Zcr/uL4Ck30r6laQ1klZlbfX73Y4IT+k8TRPwG2AKMBRYC0wruq5u9uXPgFOAdSVtXwYWZq8XAv+cvX4n8BNAwBuBX2btxwCbs/mo7PWoovuW0+dXAadkr0cAjwLT+nO/s9qPzF4PAX6Z9eUmYG7WfhXw0ez1/wCuyl7PBW7MXk/Lft+HAZOzfwdNRfcvp9+XAt8HbsuW+3V/s5p/C4wpa6vb77b3KDrMAFojYnNE/BG4AZhdcE3dEhErgJ1lzbOBa7PX1wLvKmm/LpL7gaMlvQo4G1geETsjYhewHJhV++q7JyKeiIjV2evngQ3AOPpxv7PaX8gWh2RTAG8Dbs7ay/vc/rO4GXi7JGXtN0TE7ojYArSS/j00HEnjgT8Hvp0ti37c3y7U7XfbQdFhHPB4yXJb1tZfvCIinshePwm8Intdqd999ueRHWKYTvofdr/ud3YYZg2wjfQP/zfAMxGxJ/tIaf37+5a9/ywwmr7V538F/hbYly2Ppn/3t10Ad0h6UNKCrK1uv9uDu1u19V0REZL65XXRko4EbgE+ERHPpf9AJv2x3xGxF3iDpKOBHwKvKbikmpF0LrAtIh6UNLPoeursjIjYKmkssFzSr0vfrPXvtvcoOmwFJpQsj8/a+ounst1Psvm2rL1Sv/vcz0PSEFJIfC8ibs2a+32/ASLiGeBO4HTSoYb2/wSW1r+/b9n7RwFP03f6/GbgPEm/JR0afhvwdfpvf/eLiK3ZfBvpPwQzqOPvtoOiw0pganYFxVDSya+lBdfUm5YC7Vc5fBD495L2D2RXSrwReDbbnb0deIekUdnVFO/I2hpSduz5amBDRHy15K1+229JzdmeBJIOA84inZu5Ezg/+1h5n9t/FucDP490lnMpMDe7SmgyMBV4oD69qF5EXBYR4yNiEunf588j4r300/62k3SEpBHtr0m/k+uo5+920WfzG2kiXS3wKOk472eLrqcH/VgCPAG8TDoOOZ90bPZnwCbgP4Bjss8KuDLr86+AlpL1XEQ60dcKfLjofnXR5zNIx3EfBtZk0zv7c7+Bk4GHsj6vAz6XtU8h/eFrBX4ADMvah2fLrdn7U0rW9dnsZ7EROKfovlXR95l0XPXUr/ub9W9tNq1v/9tUz99tD+FhZma5fOjJzMxyOSjMzCyXg8LMzHI5KMzMLJeDwszMcjkozMpI+s9sPknShb287s90ti2zRubLY80qyIaJ+FREnHsI3xkcHeMOdfb+CxFxZG/UZ1Yv3qMwKyOpfUTWy4G3ZM8A+GQ2AN9XJK3Mxvn/SPb5mZLukbQUeCRr+3/ZAG7r2wdxk3Q5cFi2vu+Vbiu7i/YrktZlzx24oGTdd0m6WdKvJX0vuwsdSZcrPX/jYUn/u54/IxtYPCigWWULKdmjyP7gPxsRp0oaBvxC0h3ZZ08BXhdp2GqAiyJiZza0xkpJt0TEQkmXRMQbOtnWHOANwOuBMdl3VmTvTQdeC/wO+AXwZkkbgP8OvCYion0oD7Na8B6FWfXeQRpDZw1pCPPRpHGCAB4oCQmAv5K0FrifNBDbVPKdASyJiL0R8RRwN3BqybrbImIfaWiSSaQhs/8AXC1pDvBSj3tnVoGDwqx6Aj4WEW/IpskR0b5H8eL+D6VzG/8NOD0iXk8aj2l4D7a7u+T1XqD9PMgM0gN5zgV+2oP1m+VyUJhV9jzpsartbgc+mg1njqRXZ6N5ljsK2BURL0l6DelxlO1ebv9+mXuAC7LzIM2kx9lWHNE0e+7GURGxDPgk6ZCVWU34HIVZZQ8De7NDSNeQnn0wCVidnVDeTsfjJ0v9FLg4O4+wkXT4qd0i4GFJqyMNkd3uh6RnSawljYL7txHxZBY0nRkB/Luk4aQ9nUu710WzrvnyWDMzy+VDT2ZmlstBYWZmuRwUZmaWy0FhZma5HBRmZpbLQWFmZrkcFGZmluv/A+TaC2vJhI/tAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def cost_function(theta, X, y):   \n",
    "    \"\"\"\n",
    "    Inputs: \n",
    "        theta: this is the parameter vector of size (n x 1)\n",
    "        x: this is the input training sample of size (m x n)\n",
    "        y: this is the output training sample of size (m x 1)\n",
    "        \n",
    "    Outputs: \n",
    "        cost_function(theta, x, y): the \"loss\" incurred under current value of parameter\n",
    "    \"\"\"\n",
    "    m, n = X.shape\n",
    "    yhat = np.dot(X, theta)\n",
    "    return 1/m * np.linalg.norm( np.squeeze(y) - np.squeeze(yhat), 2) ** 2 #1/m\\| y - theta^T x\\|_2^2\n",
    "\n",
    "def cost_gradient(theta, x, y):   \n",
    "    \"\"\"\n",
    "    Inputs: \n",
    "        theta: this is the parameter vector of size (n x 1)\n",
    "        x: this is the input training sample of size (1 x n)\n",
    "        y: this is the output training sample of size (1 x 1)\n",
    "        \n",
    "    Outputs: \n",
    "        cost_gradient(theta, x, y): the gradient of cost evaluated at parameter theta\n",
    "    \"\"\"\n",
    "    #np squeeze removes all dimensions whose value is 1: so if you have (1, n, 1, 1) ---> (n, )\n",
    "    yhat = np.squeeze(x).dot(np.squeeze(theta)) #using squeeze because sometimes when I use this in code, i end up passing (n,1) (n,), (1, n)\n",
    "    return 2 * (yhat  - y) * x.T  #\\Del_{\\theta} J(\\theta) = (theta^T x - y)x^T\n",
    "\n",
    "\n",
    "def gradient_descent(X, y, theta_init, max_iter=1000, step_size = 1e-2): #hyperparameters for gradient descent are max iter and step sizr\n",
    "    theta_list = [theta_init]\n",
    "    m, n = X.shape\n",
    "    for _ in range(max_iter):\n",
    "        grad = 0\n",
    "        for ii in range(m):\n",
    "            grad += cost_gradient(theta_init, X[ii:ii+1, :], y[ii]) #evaluates gradient for EACH sample and not the entire dataset! \n",
    "            \n",
    "        theta_new  = theta_init - step_size * grad / m #update rule for theta\n",
    "        theta_init = theta_new\n",
    "        #theta -= step_size * grad / m\n",
    "        theta_list.append(theta_new) \n",
    "    return np.squeeze(np.array(theta_list))\n",
    "\n",
    "\n",
    "#function to compute errrors (nice plots)\n",
    "def error_vals(X, y, theta_list):\n",
    "    max_iter = theta_list.shape[0]\n",
    "    err = []\n",
    "    for ii in range(max_iter):\n",
    "        err_tmp = cost_function(theta_list[ii], X, y) #remember that cost_function i am using entire dataset\n",
    "        err.append(err_tmp)\n",
    "    return err\n",
    "\n",
    "\n",
    "max_iter = 5000\n",
    "step_size = 1e-5\n",
    "\n",
    "\n",
    "#error plateus at approximately 1e3\n",
    "\n",
    "theta_list = gradient_descent(X_train, Y_train, np.random.randn(n,1), max_iter, step_size)\n",
    "cost_decay = error_vals(X_train, Y_train, theta_list)\n",
    "\n",
    "print(\"Final Loss is \", cost_decay[-1])\n",
    "plt.plot(cost_decay, 'r', range(max_iter), [cost_function(theta_true, X_train, Y_train)] * max_iter, 'b')\n",
    "plt.yscale(\"log\")\n",
    "plt.xlabel('iterations')\n",
    "plt.ylabel('cost function')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It seems like we can't do better than this. Let us try \"model selection\". Let us try PCA! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEGCAYAAAB2EqL0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3yV9d3/8dcnJ4vsHQIhZDNEZhCQDbXFgVRrXa1bcVTv2nH37vq1teO2y7baWr1x1HJXUdvqLbiLBRwM2coKhLBXBishEDI+vz/OQQOScJKck+uck8/z8bgenHOd9b7gkE+u67tEVTHGGGPOJczpAMYYY4KDFQxjjDFesYJhjDHGK1YwjDHGeMUKhjHGGK+EOx3AH9LS0jQ3N9fpGMYYE1RWrlxZparprT0eUgVDRKYD0wsLC1mxYoXTcYwxJqiIyI62Hg+pS1KqOk9VZyYmJjodxRhjQk5IFQxjjDH+YwXDGGOMV6xgGGOM8YoVDGOMMV6xgmGMMcYrVjCMMcZ4xQrGOew6WMfrH+9zOoYxxjjOCkYb1u89whV//oB7nl3Fpv1HnY5jjDGOsoLRiuXbD3LtrKVEuMIIDxNeXr3H6UjGGOMoKxhnsbC0ghueWkZ6XBT/uPtCJhan88rqvTQ32+qExpjuywrGGV79aC93zF5BQXocL941ht5JPfjisN7sP3qCpeXVTsczxhjHWMFoYc6HO7lvzmqG9UlmzszRpMVFAXDRwEziosLtspQxplsLqYIhItNFZNaRI0fa/VpVZWl5NZOK0/nrrReQEB3xyWPRES4uHtSTN9bt5/jJJl9GNsaYoBFSBaMzs9WKCL/98hBm3VhCj0jXZx6/Ynhvausbmb/xgC+iGmNM0AmpgtFZEa4wIlxn/ysZnZdKVmK0XZYyxnRbVjC8FBYmzBjam0WbK6murXc6jjHGdDkrGO1wxbDeNDUr89budTqKMcZ0OSsY7dCvZzwDsxJ4eY0VDGNM92MFo52uGNabtbsOs7Wy1ukoxhjTpaxgtNPlQ3sRJvCKNX4bY7oZKxjtlJkQzdjCNF5eswdVmyrEGNN9WMHogCuG9WbXweOs3nXY6SjGGNNlrGB0wMTidACWbzvocBJjjOk6VjA6IDUuir6pMazeaWcYxpjuwwpGBw3PSWbVzkPWjmGM6TasYHTQsJwkKmrq2XvkhNNRjDGmS1jB6KBhfZIBWL3zkMNJjDGma1jB6KD+WfFEhYdZO4YxptuwgtFBEa4wBmcn2hmGMabbsILRCcNyklm39yj1jbaokjEm9FnB6IRhfZI42djMxn01Tkcxxhi/s4LRCcP7WsO3Mab7CPiCISL5IvKUiPzD6SxnykyIplditDV8G2O6Bb8WDBF5WkQqRGTdGfuniUipiJSJyHfbeg9VLVfV2/yZszOGeQbwGWNMqPP3GcYzwLSWO0TEBTwKXAwMBK4TkYEicr6IvHrGluHnfJ02LCeJ3YeOU1FjA/iMMaEt3J9vrqrvikjuGbsvAMpUtRxARJ4HZqjqg8BlHf0sEZkJzATIycnp6Nu027CcJADW7DzM58/r2WWfa4wxXc2JNozewK4W93d79p2ViKSKyOPAMBH5XmvPU9VZqlqiqiXp6em+S3sO5/VKJMIlNtW5MSbk+fUMwxdUtRq4y+kcrYmOcDEwK8F6ShljQp4TZxh7gD4t7md79nWaiEwXkVlHjhzxxdt5bVhOMh/tPkJjU3OXfq4xxnQlJwrGcqBIRPJEJBK4FpjrizdW1XmqOjMxMdEXb+e1YTlJ1J1sYvOB2i79XGOM6Ur+7lY7B1gC9BOR3SJym6o2AvcCbwEbgRdVdb0/c/jb8BzPAL5ddlnKGBO6/N1L6rpW9r8OvO7Pz+5K2ck9SIuLZNWOw3xlVF+n4xhjjF8E/Ejv9nCqDUNEGNon2c4wjDEhLaQKhlNtGOBuxyivPMbhupNd/tnGGNMVQqpgOOmTAXw2HsMYE6KsYPjIkOwkwgRW7rDLUsaY0GQFw0dio8IZ0ieJ97ZUOR3FGGP8IqQKhlON3qdMLE5n7e7DHDpm7RjGmNATUgXDyUZvcBcMVXi/zM4yjDGhJ6QKhtMGZyeRFBPBos2VTkcxxhifs4LhQ64wYVxhGos2V6KqTscxxhifsoLhYxOL06msqWfjvhqnoxhjjE+FVMFwutEb3AUDsMtSxpiQE1IFw+lGb4CMhGgGZCXwrhUMY0yICamCESgmFKexYsdBausbnY5ijDE+YwXDDyYWp9PQpCzZWu10FGOM8RkrGH5Q0jeFmEgXizZXOB3FGGN8ptWCISJfbXF77BmP3evPUMEuMjyMCwuse60xJrS0dYbxzRa3/3jGY7f6IUunBUIvqVMmFqex6+BxtlfXOR3FGGN8oq2CIa3cPtv9gBAIvaROmVicAcCiUrssZYwJDW0VDG3l9tnumzPkpMaQlxZr4zGMMSGjrTW9+4vIR7jPJgo8t/Hcz/d7shAwoSiNF1bs4kRDE9ERLqfjGGNMp7RVMAZ0WYoQNbFfOn9dsoMV2w8xrijN6TjGGNMprV6SUtUdLTegFhgOpHnum3MYnZ9KpCvMutcaY0JCW91qXxWRQZ7bWcA63L2j/ldE7u+ifEEtJjKcC/JSWFBq7RjGmODXVqN3nqqu89y+BfiXqk4HRhGg3WoD0ZT+GZRV1LKj+pjTUYwxplPaKhgNLW5PBV4HUNUaoNmfoToqkMZhnDJ1gLt77Tsb7bKUMSa4tVUwdonIfSJyBe62izcBRKQHENEV4dorkMZhnNI3NZbCjDje2XTA6SjGGNMpbRWM24DzgJuBa1T1sGf/aOAvfs4VUqYOyGBZ+UGOnmg495ONMSZAtdVLqkJV71LVGar6dov9C1T1t10TLzR8bkAmjc1qa2QYY4Jaq+MwRGRuWy9U1ct9Hyc0Dc9JJjkmgn9vrOCywb2cjmOMMR3S1sC9McAuYA6wjACdPyoYuMKEyf0yWFBaQVOz4gqzv0pjTPBpqw2jJ/B9YBDwMHARUKWqi1R1UVeECyVTBmRwqK6BVTsPOR3FGGM6pK02jCZVfVNVb8Ld0F0GLLS1MDpmQnE64WHC/I3WW8oYE5zaXHFPRKJE5Ergb8DXgEeAl7siWKhJiI5gVH6KjccwxgSttqYGmQ0swT0G4wFVHamqP1PVPV2Wrp0CceBeS1P7Z9qob2NM0GrrDOOrQBHwdWCxiBz1bDUicrRr4rVPIA7ca8lGfRtjgllbbRhhqhrv2RJabPGqmtCVIUOFjfo2xgSzNtswjO/ZqG9jTLCygtHFbNS3MSZYWcHoYi1HfRtjTDA5Z8EQkV95s89459So7/kbD7D5QI3TcYwxxmvenGFcdJZ9F/s6SHdy2/g8IsPDmP7H93l22Q5U1elIxhhzTm2Nw7hbRD4G+onIRy22bcBHXRcx9JzXK5E3vj6BUfmp/ODlddzz7CqO1FkjuDEmsElrv92KSCKQDDwIfLfFQzWqerALsnVYSUmJrlixwukY59TcrDz5fjm/frOUjPgoHr5uGCNzU5yOZYzppkRkpaqWtPr4uS6HiEgBsFtV60VkEjAYmN1iQaWAEywF45S1uw7zH8+vZkd1HenxUeSmxpCTEkvf1Bj6psbwhfN6Eh3hcjqmMSbE+aJgrAFKgFzc63q/Apynqpf4MKdPBVvBAKg50cBzy3a6pw45WMfO6jr2Hz0BwJXDevO7a4Y6nNAYE+rOVTDaWg/jlGZVbfRMQvhHVf2jiKz2XUQDEB8dwZ0TC07bd6KhiV++sYnZS7Zzz+RCCjPinAlnjDF410uqQUSuA24EXvXsi/BfJHNKdISL+6YUEh3h4pF3tjgdxxjTzXlTMG7BvfreL1R1m4jkAf/r31gdE+iz1XZEalwUN12Yy7yP9rLFxm0YYxx0zoKhqhuA/wJWee5vU9WAHLgX6LPVdtTM8fnERLj4g51lGGMc5M1I7+nAGuBNz/2hIjLX38HMp5JjI7llbB6vf7yP0v12lmGMcYY3l6R+AlwAHAZQ1TVAvh8zmbO4fXwecZHhPPzOZqejGGO6Ka8avVX1zEaBZn+EMa1LionklrG5vP7xfjbsDcj1q4wxIc6bgrFeRK4HXCJSJCJ/BBb7OZc5i9vG5RMfbWcZxhhneFMw7gPOA+qB54AjuJdtNV0sMSaC28bl8db6A6zbEzo9wYwxwcGbgnGpqv5AVUd6th8Cl/s7mDm7W8flkRAdzg/+bx0VNSecjmOM6Ua8KRjf83Kf6QIJ0RH86kuDKd1/lEsefp8lW6udjmSM6SZanRpERC4GLgF6i8gjLR5KABr9Hcy07uLzs8hPj+PuZ1fylSeX8q3P9+PuiQWEhYnT0YwxIaytM4y9wArgBLCyxTYX+IL/o5m29OsZz7x7xzF9SC9+81YptzyznIPHTjodyxgTwryZrTZCVYNqdZ9gnK22o1SV5z7cyQNzN5CVFM38b04kwmVLtRtj2u9cs9V6MzVIUBWL7kZE+Mqovvzmy4PZUV3H8u0BvbaVMSaI2a+iIWLqgEwiXWEs2FThdBRjTIiyghEi4qLCGZWfwjtWMIwxftJWL6l5QKsNHKpqYzECzJT+GTwwbwPbq46RmxbrdBxjTIhp6wzjt8BDwDbgOPCEZ6sFtvo/mmmvKf0zAPi3nWUYY/yg1TMMVV0EICIPndFqPk9EukcXpCDTNzWWgvRYFpRWcOu4PKfjGGNCjDdtGLEi8sl05p4V9+x6R4CaOiCTpeXV1Nbb2EpjjG95UzC+ASwUkYUisghYANzv31imoyb3y6ChSXl/S5XTUYwxIabVS1KnqOqbIlIE9Pfs2qSq9f6N9SkR+SJwKe4pSZ5S1be76rODUUluMvHR4fx70wGmDerpdBxjTAjxZonWGOA/gXtVdS2QIyKXefPmIvK0iFSIyLoz9k8TkVIRKROR77b1Hqr6f6p6B3AXcI03n9udRbjCmFCczoLSSpqb2x7Fb4wx7eHNJam/ACeBMZ77e4Cfe/n+zwDTWu4QERfwKHAxMBC4TkQGisj5IvLqGVtGi5f+0PM6cw5T+2dQWVPPur22ZoYxxnfOeUkKKFDVa0TkOgBVrRMRr6ZFVdV3RST3jN0XAGWqWg4gIs8DM1T1QeAzZy6ez/ol8Iaqrmrts0RkJjATICcnx5t4IWtSvwxE3N1rB2cnOR3HGBMivDnDOCkiPfAM4hORAtyr73VUb2BXi/u7Pftacx/wOeAqEbmrtSep6ixVLVHVkvT09E7EC34psZEM65Nk4zGMMT7lzRnGT4A3gT4i8iwwFrjZj5lOo6qPAI+c84nmNFMHZPKbt0qpqDlBRny003GMMSHAm9lq3wauxF0k5gAlqrqwE5+5B+jT4n62Z5/xocn93M0/CzdVOpzEGBMqvOkl9Q4wSlVfU9VXVbVKRGZ14jOXA0UikicikcC1uBdl6jQRmS4is44cscbeAVnxZCVG22UpY4zPeNOGkQf8l4j8uMW+VhfYaElE5gBLgH4isltEblPVRuBe4C1gI/Ciqq5vZ+6zUtV5qjozMTHRF28X1ESEyf0zeG9LJfWNTU7HMcaEAG8KxmFgKpApIvNExOufxqp6napmqWqEqmar6lOe/a+rarGqFqjqLzqY3ZzD1P4ZHDvZxPJth5yOYowJAd4UDFHVRlW9B/gn8D6QcY7XOMIuSZ1uTEEqkeFhLCi1y1LGmM7zpmA8fuqGqj6Du/E7IKfnsEtSp4uJDGd0fioLrWAYY3yg1YIhIgmem38XkZRTG+71Mb7dJelMp03ul87WymPsrK5zOooxJsi1dYbxnOfPlcAKz58rW9w3QeCT7rWb7SzDGNM5rRYMVb3M82eequZ7/jy15bf2OhNYctNiyU2NYYF1rzXGdFJba3oPb+uFbc3r5BQRmQ5MLywsdDpKQJnUL4M5H+7kREMT0REup+MYY4JUW1ODPNTGYwpM8XGWTlPVecC8kpKSO5zOEkgm98/gmcXbWVpezaR+AdnBzRgTBNpa03tyVwYx/jMqL4XoiDAWllZawTDGdJg33WoRkUEicrWI3Hhq83cw4zvRES4uLEjj35sqULVFlYwxHePNXFI/Bv7o2SYDvwYu93Mu42OT+6Wz82Ad26qOOR3FGBOkvDnDuAr31CD7VfUWYAgQkCPjbKR3605dilpYarPXGmM6xpuCcVxVm4FGz2C+Ck6fnjxg2Ejv1vVJiaEwI86mCTHGdJg3BWOFiCQBT+AetLcK9wy0JshMKk5nWflB6k42Oh3FGBOEvFlA6R5VPayqjwMXATd5Lk2ZIDO5fwYnm5pZXFbtdBRjTBDytpfUYBG5HBgOFIrIlf6NZfyhJDeZ2EiXXZYyxnTIOdf0FpGngcHAeqDZs1uBl/yYy/hBVLiLsYVpLCytRFUREacjGWOCyDkLBjBaVQf6PYkP2NQg5zapXwZvbzhAWUUtRZnxTscxxgQRby5JLRGRoCgY1kvq3Cb1Swdg7tq9DicxxgQbbwrGbNxFo1REPhKRj0XkI38HM/7RK6kHl56fxeOLtrJ+r41XMcZ4z5uC8RRwAzANmA5c5vnTBKmffXEQiT0i+eYLa6lvbHI6jjEmSHhTMCpVda6qblPVHac2vyczfpMSG8mvrzqf0gM1/O7tzU7HMcYECW8avVeLyHPAPKD+1E5VtV5SQWxK/0yuuyCHWe+VM6V/BqPyU52OZIwJcN6cYfTAXSg+j/tS1KnLUibI/fDSAeSkxPCtv6+l5kSD03GMMQGuzYIhIi7gI1W95Yzt1i7KZ/woNiqc3109hL2Hj/OzVzc4HccYE+DaLBiq2gRc10VZOs1mq22/EX1TuHtSAS+u2M3b6/c7HccYE8C8uST1gYj8SUTGi8jwU5vfk3WAjcPomK9PLWZAVgI/fXUDjU3N536BMaZb8qZgDAXOA36Ke53vh4Df+jOU6VqR4WF886Jidh86zmsf73M6jjEmQJ2zl5St7d09TO2fQVFGHI8t3MrlQ3rZPFPGmM/wZonWRBH5nYis8GwPiYhd8wkxYWHCzAn5bNpfw8LNtiqfMeazvLkk9TRQA1zt2Y4Cf/FnKOOMGUN7k5UYzeMLtzodxRgTgLwpGAWq+mNVLfdsDwD5/g5mul5keBi3jctj2baDrN55yOk4xpgA49Wa3iIy7tQdERkLHPdfJOOkay/IISE6nMcX2VmGMeZ03kwNchcw29NuIcBB4GZ/hjLOiYsK58YxuTy6sIytlbUUpMc5HckYEyC8WdN7raoOwb3q3vmqOkxV1/o/WvvZwD3fuHlsLpGuMGYtKnc6ijEmgHjTSypKRK4H7gXuF5EficiP/B+t/Wzgnm+kxUVxdUkfXlq9m/1HTjgdxxgTILxpw3gFmAE0AsdabCaE3TE+n6Zm5ekPtjkdxRgTILxpw8hW1Wl+T2ICSk5qDJcO7sUzi7eTEhvJbePyiHB58/uFMSZUefMTYLGInO/3JCbg/OiygUwoSueXb2zikoffY2l5tdORjDEO8qZgjANW2pre3U96fBRP3lTCkzeWcLyhiWtnLeUbL6yhosbaNYzpjry5JHWx31OYgPa5gZmMLUzjzwvL+J9F5czfcIBLzs/i0sFZjClItUtVxnQToqpOZ/C5kpISXbFihdMxQlJ5ZS1//HcZ/9pwgNr6RpJiIvjCwJ5cMjiL8YVphIXZpIXGBCsRWamqJa0+bgXDdMSJhibe21LF6x/vY/6GA9TUN/Kti4q5b2qR09GMMR10roLhzSUpYz4jOsLFRQMzuWhgJvWNTdz+1xX8bdkO7p5UQLhdojImJNn/bNNpUeEuvjKqLweO1rPIpkY3JmRZwTA+MXVABmlxkbywfJfTUYwxfmIFw/hEhCuML43I5p1NFdbt1pgQZQXD+MzVJX1oalb+uXKP01GMMX4QUgXDZqt1VkF6HBfkpvDiil2EYu87Y7q7kCoYNlut864Z2YdtVcf4cNvBTr3PiYYmKo6esMJjTACxbrXGpy45P4ufzF3PC8t3MSo/1avXVNSc4OH5W9h5sI6Ko/UcqDnB4boGAK4p6cMvv3Q+IjYg0BinWcEwPtUj0sXlQ3vxj5W7+fHl55HYI6LN51fW1HPdrKXsPnScAVkJ9E2N4YK8FDIToth96DjPL99FUWYct4+3ZeTP5mRjM4fqTpKZEO10FNMNWMEwPnftyByeXbaTuWv2cMOY3FafV1Vbz/VPLGXv4RPMvvWCz5yRNDcrR4438N+vb6QoM56Jxel+Tu47q3ceoqr2JJP7pfttIGN9YxM3P72cD7cf5IbRffnGRcXnLNDGdEZItWGYwDCodwIDsxJ4YUXrYzKqa+v5yhPL2HWojqdvHnnWy1dhYcJDVw+hX88E7n1uFeWVtf6M7TOqyr3PreaO2SuY8OsFPLqgjOraep9+RnOz8s0X1rKkvJrJ/dKZvWQ7U367kBeX76K52dp9jH/YXFLGL/66eDs/nrueV+8bx6Dep3dCOHjsJNc/sZRtVcd4+uaRjC1Ma/O9dh+qY8afPiAxJoKX7xkb8L9Fbz5Qw+d//y7XlPRh9+E6PiirJtIVxmVDsrhyWDa9kqJJi48iPiq8Q20zqsoD8zbwzOLtfP+S/sycUMC6PUf48dz1rNxxiCF9kvjOF/rR2KyUVdSytbKWrRW1bK8+RnhYGEkxESTHRJIYE0FyTARXDMtmRN9kP/xNmGBjkw8aRxypa2Dkf8+nKCOO4TnJJMdEkBgTSXJMBE++t42yylqeuqmE8UXeXWb6cNtBrn9iKWML03j65pG4AnhW3FnvbuW/X9/E4u9OoVdSD7YcqGH2kh38c9Vu6k42ffK8qPAw0uKi6J3Ugy+XZDNjaG8iw8990v/Ywq386s1N3DYujx9eOuCToqOqvLx6Dw++sYnKmk/PaJJiIihMjyM3LZZmVQ7XNXCo7iRH6hrYf/QEPSJcLPrOZOKi7Ap1d2cFwzjmkXe28PLqPe4fTscbOPVVi3SF8cRNJe1uk5jz4U6+99LHfG1yAf/5hf5+SOwb1z+xlKraet7+xsTT9tecaGDNrsNU1dZTWVNPVe1JqmrqWb/3KKUHashKjOb28flcO7IPsa388P7Hyt18++9ruXxIL/5wzdCzTidfc6KBhaWVZCZEU5AeS0psZKtnMqt2HuLKPy/mmxcV8x8203C3ZwXDBISmZuXocfdvtvHREaTHR3Xofb714lrmrd3Lgv+cRO+kHj5O2Xm19Y0M++nb3DI2j+9fMsCr16gqizZX8ueFW/lw20GSYiK4cUwu/XvG09SsNKt7q6o5yS/f3MSY/FSevnmkV2cj3pg5ewWLt1bz7ncmkxIb6ZP3NMHJpjc3AcEVJiTHRpLcyR9I3/p8MfPW7uWR+Vv41VWDfZTOdxaXVdHQpExqx9mTiDCpXwaT+mWwcsdBHltYziPvbDnrc8/vnchjXx3us2IB8O0v9GPaH97lzwvK+OFlA332vib0WMEwQaVXUg++MjqH2Ut2cNekAvLSYp2OdJqFmyuJjXRRkpvSodeP6JvCkzelsPfwcY6eaMAlQliYECaCS4ReSdE+76ZbnBnPlcOzmb10B7eOy6NXAJ65mcBg3WpN0LlnUiGRrjD+MH+z01FOo6osKq3kwsK0Tp8B9ErqQf+eCRRlxlOQHkdeWiw5qTF+G9Nx/+eKQOHh+Wc/szEGrGCYIJQeH8XNY3OZu3YvpftrnI7zibKKWvYcPs6kfsEzwPCU7OQYvjq6L39fuYuyiuAY79KdzV27l5+9uqHLP9cKhglKd07IJy4ynIfeLnU6yicWlrpXG5zUL8PhJB3ztckF9IhwBdTfqfms5mbl129u4qn3t3Ho2Mku/WwrGCYoJcVEcseEfN7ecIC1uw47HQeAhZsrKM6MC8jeW95IjYvijgn5vLFuf8D8nZrPer+sit2HjgOwbFt1l362FQwTtG4dl0dKbCS/DYDfiI/VN7J826GgPbs45fbx+aTERvKzVzdwvMUgQxM4nl++k+SYCHpEuFiy1QqGMV6Jiwrn7okFvLelimXlXfsf50yLt1Zzsqm5Xd1pA1FcVDjfv2QAK3ce4ouPfhA083f505wPdzLuV/+mrML59rLKmnreXn+AK4dnU5KbzGIrGMZ474YxfcmIj+L7L3/Mg29sZNa7W/nnyt0sKK1g476jXbYA08LSik51pw0kV43I5plbLuBAzQku/9MHvLlun9ORHHNqdoHdh47zo1fWO76g1z9X7aaxWbnugj5cWJDGlora06aB8Tcbh2GCWnSEi5/OGMRP563nL+9v52RT82mPD8lO5N4pRXxuQIbfFmFSVRb6qDttoJhYnM5r/zGee55dxV1/W8Ud4/P4zrT+RLjCUFUqaurZWlnLzuo68tPjGJaTRISfuvw65XlPsZjUL51xhWn8/LWNvP7xfi4dnOVIHlXlheW7GJmbTGFGPLX17kuGS8urmT6kV5dkCPiCISIDgK8DacA7qvqYw5FMgJk2qCfTBvVEVampb+Rg7Umqj51kw94jzHqvnDtmr6B/z3jum1LEtEE9fT5x4dZKd3faeyYX+PR9ndY7qQcv3jmaX7y2kSfe28Z7W6qIDA+jvPIYtfWNpz03PiqcCwtTmVCczoSidPqkxDiU+rOampXnPtxJfloso/NTvfr3f3H5Lr770sdMLE7n8a+OIMIVxkur9vDz1zYwqV96q3N9+dPS8oNsqzrGvZMLARjUK4G4qHAWbw2RgiEiTwOXARWqOqjF/mnAw4ALeFJVf9nae6jqRuAuEQkDZgNWMMxZiQgJ0REkREeQmxbLiL7JXHdBDnPX7uVPC8r42nOrKEiPZdqgngzMSmRgrwT6psScdQK/9gj27rRtiQp3n8GN6JvMYwu3ktgjgqtGZJOfHktBehzZyT3YsPco726p5N3NVby1/gAAM4b24vdXn31yxK7U2NTMt/6+llfW7AUgKzGaLw7rzZXDelOUGX/W17y4Yhf/9dJHTCxO539uGEF0hAuAn844j6seX8KjC8r4zrSun/zy+eU7iY8O55Lz3Wc44a4wRuWlsLQL2+/8XSafAf6E+wc9ACLiAh4FLgJ2A8tFZC7u4vHgGa+/VVUrRORy4G7gf/2c14SYcFcYVw53Tx3+xrp9PPneNv5nUTmNnkWGYmR0FjEAAA0gSURBVCJdDMhKoDgznvy0WPLSYslLj6VPcozXl5cWllZSlBG83Wm9MWNob2YM7X3Wx/qmxnLx+VmoKlsrj/HC8p088d42+qbE8M3P9+vipJ9qbGrmG57JKr/9+WL6psby0qrdzHq3nMcWbuX83okMzEqgoamZ+qZmGhqbOdHYzHtbKhlXmHZasQAoyU3hyuG9eeK9ck/RjOuyYzl07CRvrNvPtSP70CPy00xjClJ5Z1MF+44cJyvR/98/vxYMVX1XRHLP2H0BUKaq5QAi8jwwQ1UfxH02crb3mQvMFZHXgOfO9hwRmQnMBMjJyfFJfhM6XGHCZYN7cdngXtQ3NrHlQC0b9h1lw1739tb6/RxsMQjKFeaet6lnQjQZCdFkxkeTmRBFYo8I9h45wY7qY2yvrmNH9TEO1zVwx/g8B48uMIgIhRlxfP+SARyua+CRf5dR3DOeywZ3zeWSlhqamrn/hTW89tE+vndxf+6c6L5cOH1ILypr6pm7di+vrNnDos2VRIQLEa4wIl1hRIaH8eUR2fx0xqDTisUp3724P/9af4CfzNvAX28Z6bd2sTO9tHoPJxubuXbk6T/bRntWqlyytZorh2f7PYcTbRi9gZZrd+4GRrX2ZBGZBFwJRAGvt/Y8VZ0FzAL39Oa+CGpCU1S4i0G9Ez+zEuDhupNsqzr2ybajuo6KmhNs3HuUBUcrPln8SMR9fT83NZZLz88iNzWWL43w/3/WYCEi/PyKQWyrOsa3/76W3NTYz/xdt+VEQxNlFbWUVdSy+UANWyrcKwbmpMZw35Sic64O2NDUzH/MWc0b6/bzg0sGcMeE/NMeT4+P4rZxedw2rv1FPiM+mvsvKuZnr27g7Q0H+MJ5Pdv9Hu2lqjz/4U6G9EliYK+E0x4bmJVAYo+IkC4Y7aKqC4GFDscw3UBSTCTDciIZlnP2H0i19Y0cOnaSjIQoosI/+9un+VRUuIvHbxjBjD99wB2zV/DKvWPJiI9u8zX7j5xg1rvlzPlwJ8cb3MU5PEzIS4ulKDOO5dsP8aXHFjO+KI37P1fEiL6nd2Fubla2VR/j129u4q31B/jhpQO4fXz+2T6qU24a05cXl+/ip/M2MLE4/axnIr60auchtlTU8ssrz//MY2Fhwuj8FJZ0UTuGEwVjD9Cnxf1szz5jAlpcVLgtY9oOaXFRzLpxBFc9toQ7/3clc+4YfdYfrrsO1vHYoq38Y8VumlSZMaQXUwdkUpwZR9/U2E/ako7VN/K3pTuY9W45X3psCeMK0/jCoJ6UHahh/d6jbNh39JOzwB9PH8gtY/1zmTDcFcYDM87j2llL+c1bpfw/P68hMufDXcRGulrtCTUmP5W31h9g18E6v/dOc+LbvxwoEpE83IXiWuB6X7yxiEwHphcWFvri7YwxnXRer0R+f80Q7vrbKm58+kOKM+MIDwsjPEwId4Wx78hxXv1oHy4RrirJ5u6JBa3+0IuNCufOiQXcMKbvJ4Xj/bIqYiNdDOyVwNUlfRjYK4HhOUkUZpy9B5SvjM5P5YbRfXnq/W2MLUxlSv9Mv3zO7kN1zF2zly+XZLfalffCwjTA3Y7h74Lh1yVaRWQOMAn3GIoDwI9V9SkRuQT4A+6eUU+r6i98+bm2RKsxgeWp97fxxLvlnGxqpqGpmaZmpbFJiXAJ14zMYeaEfHomtn3J6kzHTzZRWVNPdnIPR7rvnmho4ouPfkBFTT1vfH08mQnty++N7/xjLf+3uu0liVWVkb+Yz/iidH5/zdBOfZ6jS7Sq6nWt7H+dNhqwjTGhpaONzG3pEekiJ9W5AYLRES7+dP0wpv/xA+5/fg1/u32UTweFllfW8s9Ve7hxTN82u2yLCKPyU1m8tQpV9WvPrZAayy8i00Vk1pEjR5yOYozpBgoz4nlgxnksKa/mzwvKfPrev5+/hUhXGPdMOvcl9jH5qRw4Ws+2qmM+zXCmkCoYqjpPVWcmJnrfhc8YYzrjyyOyuXxIL34/fzPLtx/0yXtu3HeUeWv3csvYXNLjo875/AsLPOMx/NxbKqQKhjHGdDUR4RdXDCI7OYavz1nN4brOr4L30NubiY8O584J3s1PlpcWS2ZClN+nO7eCYYwxnRQfHcGfrh9GZW09X3tuFScaOr741Kqdh5i/8QB3TsgnMSbCq9eICGPyU1lWXu3XKditYBhjjA8Mzk7iwSsHs3hrNbf9dXmHVyx86O1SUmMj2z2O5MKCNKpqT7Klwn+LXoVUwbBGb2OMk64akc1vrhrC4q3V3PrMcupONp77RS0s3lrFB2XV3DO5sN1TqE8dkMGzt48ix49jMUKqYFijtzHGaVeNyOZ3Vw9h2TbvioaqcuDoCd7fUsWDr28iKzGar4xq/wSqqXFRjC1M8+tUJTbPgTHG+NgVw7IJE+EbL6zh5r8s5y83jyRMhJ0H69hefYyd1XWUVx1ji2dyxSPHGwD3xJa/v3qo3+en6igrGMYY4wczhvZGRLj/+dVc8Iv5HDujTSM5JoKijHguG5xFUUYcxZnxFPeMJy3u3N1onWIFwxhj/OTyIb2Ii3Lx5rr99EmOoW9aLLmpMfRNifW6B1QgCamCYZMPGmMCzZT+mX6bnLCrWaO3McYYr4RUwTDGGOM/VjCMMcZ4xQqGMcYYr1jBMMYY45WQKhg2NYgxxvhPSBUM6yVljDH+E1IFwxhjjP+IP+dOd4qIVAI7OvjyNKDKh3ECQagdkx1P4Au1Ywq144GzH1NfVU1v7QUhWTA6Q0RWqGqJ0zl8KdSOyY4n8IXaMYXa8UDHjskuSRljjPGKFQxjjDFesYLxWbOcDuAHoXZMdjyBL9SOKdSOBzpwTNaGYYwxxit2hmGMMcYrVjCMMcZ4xQpGCyIyTURKRaRMRL7rdJ72EpGnRaRCRNa12JciIv8SkS2eP5OdzNgeItJHRBaIyAYRWS8iX/fsD+ZjihaRD0VkreeYHvDszxORZZ7v3gsiEul01vYQEZeIrBaRVz33g/14tovIxyKyRkRWePYF8/cuSUT+ISKbRGSjiIzpyPFYwfAQERfwKHAxMBC4TkQGOpuq3Z4Bpp2x77vAO6paBLzjuR8sGoFvqepAYDTwNc+/STAfUz0wRVWHAEOBaSIyGvgV8HtVLQQOAbc5mLEjvg5sbHE/2I8HYLKqDm0xViGYv3cPA2+qan9gCO5/q/Yfj6ra5m74HwO81eL+94DvOZ2rA8eRC6xrcb8UyPLczgJKnc7YiWN7BbgoVI4JiAFWAaNwj7gN9+w/7bsY6BuQ7fmBMwV4FZBgPh5P5u1A2hn7gvJ7ByQC2/B0curM8dgZxqd6A7ta3N/t2RfsMlV1n+f2fiAoFxcWkVxgGLCMID8mz+WbNUAF8C9gK3BYVRs9Twm2794fgO8AzZ77qQT38QAo8LaIrBSRmZ59wfq9ywMqgb94Lhs+KSKxdOB4rGB0I+r+VSLo+lGLSBzwT+B+VT3a8rFgPCZVbVLVobh/M78A6O9wpA4TkcuAClVd6XQWHxunqsNxX6L+mohMaPlgkH3vwoHhwGOqOgw4xhmXn7w9HisYn9oD9GlxP9uzL9gdEJEsAM+fFQ7naRcRicBdLJ5V1Zc8u4P6mE5R1cPAAtyXbJJEJNzzUDB998YCl4vIduB53JelHiZ4jwcAVd3j+bMCeBl3YQ/W791uYLeqLvPc/wfuAtLu47GC8anlQJGnd0ckcC0w1+FMvjAXuMlz+ybc7QBBQUQEeArYqKq/a/FQMB9TuogkeW73wN0msxF34bjK87SgOSZV/Z6qZqtqLu7/M/9W1a8QpMcDICKxIhJ/6jbweWAdQfq9U9X9wC4R6efZNRXYQEeOx+kGmUDagEuAzbivKf/A6TwdyD8H2Ac04P6t4jbc15PfAbYA84EUp3O243jG4T5N/ghY49kuCfJjGgys9hzTOuBHnv35wIdAGfB3IMrprB04tknAq8F+PJ7saz3b+lM/C4L8ezcUWOH53v0fkNyR47GpQYwxxnjFLkkZY4zxihUMY4wxXrGCYYwxxitWMIwxxnjFCoYxxhivWMEwJgCIyE9E5NtO5zCmLVYwjPExcbP/Wybk2JfaGB8QkVzPWiqzcQ/Ie0pEVrRc88LzvO0i8oCIrPKst/CZeaRE5A4RecMzEtyYgBF+7qcYY7xUBNykqktFJEVVD3rWWXlHRAar6kee51Wp6nARuQf4NnD7qTcQkXtxTxfyRVWt7/IjMKYNdoZhjO/sUNWlnttXi8gq3NOAnId7Ua5TTk2iuBL3+iWn3Ih7dtSrrFiYQGQFwxjfOQbu5UlxnzlMVdXBwGtAdIvnnSoGTZx+lv8x7gKS7fekxnSAFQxjfC8Bd/E4IiKZuM8avLEauBOYKyK9/BXOmI6ygmGMj6nqWtw//DcBzwEftOO17+M+O3lNRNL8k9CYjrHZao0xxnjFzjCMMcZ4xQqGMcYYr1jBMMYY4xUrGMYYY7xiBcMYY4xXrGAYY4zxihUMY4wxXvn/34t/21XqoSQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#define a function for normalized MSE \n",
    "\n",
    "def normalized_mse(y_true, y_hat):\n",
    "    return np.linalg.norm(y_true - y_hat, 2) ** 2 / (np.linalg.norm(y_true, 2) **2)\n",
    "\n",
    "\n",
    "err_norm_mse = np.zeros(m_train-1)\n",
    "\n",
    "for r in range(1, m_train):\n",
    "    \n",
    "    ##PCA in general mean subtraction is very important. \n",
    "    #if mean is not subtracted -- then the first principal component (eigenvector corresponding to the largest eigenvalue ) will not give great results\n",
    "    #largest eigenvalue corresponds to the mean in this case! \n",
    "    \n",
    "    mu_hat = np.mean(X_train, axis=0) #empirical mean computed from training data\n",
    "    Z = X_train - mu_hat ##using broadcasting rules\n",
    "\n",
    "\n",
    "    [U_full, S_full, V_full] = np.linalg.svd(Z) #since the \"data\" is along rows, V_full corresponds to the principal components\n",
    "    #and U * S will correspond to \"scores\" -- these are a metric of how close the true data is to the \"compressed data\"\n",
    "    \n",
    "    \n",
    "    V = V_full[:, :r] #top r principal components since SVD orders singular values in decreasing order\n",
    "    B = X_train @ V ##kind of like a compressed version of the data. \n",
    "    \n",
    "    # X_train * V \\approx u_full * S_full * V_full.T * V  \\approx U_full * S_full\n",
    "\n",
    "    B_tilde = np.hstack((np.ones((60,1)), B)) ##proxy/compressed version of X\n",
    "\n",
    "    #linear regression on the compressed data\n",
    "    #could use gradient descent here as well but have to tune hyperparameters \n",
    "    #pseudo inverse is fast enough for small data\n",
    "    \n",
    "    theta_hat_full = np.linalg.pinv(B_tilde) @ Y_train\n",
    "\n",
    "\n",
    "    #extract bias part and the rest\n",
    "    bias_hat = theta_hat_full[0]\n",
    "    theta_hat = theta_hat_full[1:]\n",
    "\n",
    "    #finally predict\n",
    "    #original linear regression -- y_test = X_test * theta_hat, where theta_hat was learned from training data (X_train). \n",
    "    y_hat_test = X_test @ V @ theta_hat + bias_hat ##kind of undoing the compression (verify through SVDs!)\n",
    "    \n",
    "    #another way to think about this is that prediction is also projecting the test data onto subspace V. \n",
    "    err_norm_mse[r-1] = normalized_mse(Y_test, y_hat_test) \n",
    "\n",
    "plt.plot(err_norm_mse)\n",
    "plt.yscale(\"log\")\n",
    "plt.xlabel('rank')\n",
    "plt.ylabel('normalized test MSE')\n",
    "plt.show()\n",
    "\n",
    "\n",
    "#to get a better plot, incease m and also run monte carlo\n",
    "#the plot should look like a U shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD5CAYAAADLL+UrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxU9bnH8c9DAEFAQEBkBxW1uFIptS51R1TULlbFtdWWYmsvVmyL195q29t71bbWBa2lGtEuIMVal1LRohXstUhAVBaRFEEStihhD0vC7/7xTF6ZQALJzCRn5sz3/Xr9XidzMstzkplnfuc5v/M7FkJARETyQ4uoAxARkeajpC8ikkeU9EVE8oiSvohIHlHSFxHJI0r6IiJ5pGXUAdSla9euoX///lGHISKSU+bOnftxCKHbvu6TlUm/f//+FBUVRR2GiEhOMbMV+7tPxpO+mR0G3AF0DCFcllh3OnB14vUGhRBOyfTriojI/jWopm9mhWa2zswW7LF+uJktMbNiMxsHEEJYFkK4Mfl+IYRZIYTRwIvAk5kKXkREGqehB3InAsOTV5hZAfAwcAEwCBhpZoP28zxXAX9sZIwiIpIhDUr6IYSZwPo9Vg8FihM9+53AZODS+p7DzPoCG0MIm1MNVkRE0pPOkM1ewMqk2yVALzPrYmaPAoPN7Pak398IPFHfk5nZKDMrMrOisrKyNMISEZH6ZPxAbgjhE2B0Hevv3M/jJgATAIYMGaKpP0VEmkA6Pf1SoE/S7d6JdSIikqXS6enPAQaa2QA82V+JH6gVEckNK1fCli1QWQm7dvmy+ufkVlkJVVUQAuzeXXsJNUvw9VVVNc+V3KqqardWreDAA6Ft25pl165wzjlNtskNSvpmNgk4E+hqZiXAnSGEx83sZmA6UAAUhhAWNlmkIiKZ9Mtfwm23RR3F3o45BhYs2P/9UtSgpB9CGFnP+mnAtIxGJCLS1IqL4Yc/hPPPh699DVq29NaqFRQU+DK5tWzp61u0ALOaZXWDmiXUPKa6FRTULKtbixa+F1FR4W3bNl82saychkFEpMmEAN/8JrRuDYWF0LNndLFUf6kcdFCzvaSSvojkl6eegldfhV//OtqEHxFNrSwi+aOsDG69FU49FUaNijqaSCjpi0j++O53YfNmmDDBa+p5KD+3WkTyz/Tp8Ic/wO23w6D9TRMWX0r6IhJ/W7fC6NFw1FGe9POYDuSKSPz94hewfDn84x/Qpk3U0URKPX0Rib+5c+G44+CMM6KOJHJK+iISfxs2QJcuUUeRFZT0RST+ysuhc+eoo8gKSvoiEn/l5dCpU9RRZAUlfRGJvw0b1NNPUNIXkXjbtcuHbKqnDyjpi0jclZf7Uj19QElfROJuwwZfKukDSvoiEnfVPX2VdwAlfRGJO/X0a1HSF5F4U02/FiV9EYk3lXdqUdIXkXhTeacWJX0Ribfycp9ZM89n16ympC8i8aYpGGpR0heReNMUDLUo6YtIvGmGzVqU9EUk3lTeqUVJX0TiTeWdWjKe9M3sMDN73Mym7rG+nZkVmdmITL+miEi91NOvpUFJ38wKzWydmS3YY/1wM1tiZsVmNg4ghLAshHBjHU/zA2BK+iGLiDTQ7t3q6e+hoT39icDw5BVmVgA8DFwADAJGmtmguh5sZucBi4B1KUcqItJYmzdDCEr6SVo25E4hhJlm1n+P1UOB4hDCMgAzmwxciif3PZ0JtMO/HCrMbFoIYXeKMYuINIymYNhLOjX9XsDKpNslQC8z62JmjwKDzex2gBDCHSGEW4A/Ar+tK+Gb2ahEzb+orKwsjbBERBI02dpeGtTTb4wQwifA6Hp+N3Efj5sATAAYMmRIyHRcIpKHNO/OXtLp6ZcCfZJu906sExHJDirv7CWdpD8HGGhmA8ysNXAl8HxmwhIRyQCVd/bS0CGbk4A3gaPMrMTMbgwhVAI3A9OBxcCUEMLCpgtVRKSRVN7ZS0NH74ysZ/00YFpGIxIRyZTycmjRAtq3jzqSrKFpGEQkvqrPxm2hVFdNfwkRiS+djbsXJX0RiS/Nu7MXJX0RiS/19PeipC8i8aWe/l6U9EUkvnTVrL0o6YtIfKm8sxclfRGJp4oK2LFD5Z09KOmLSDxpCoY6KemLSDxpCoY6KemLSDxphs06KemLSDypvFMnJX0RiSeVd+qkpC8i8aTyTp2U9EUknpT066SkLyLxtGGDz6PfqlXUkWQVJX0RiSfNu1MnJX0RiSfNu1MnJX0RiacNG9TTr4OSvojEk3r6dVLSF5F4UtKvk5K+iMSTyjt1UtIXkfiprITNm9XTr4OSvojEj6ZgqJeSvojET3XSV3lnLxlP+mZ2mJk9bmZT97VORKTJaIbNejUo6ZtZoZmtM7MFe6wfbmZLzKzYzMYBhBCWhRBuTL5fXetERJqMyjv1amhPfyIwPHmFmRUADwMXAIOAkWY2KKPRiYikQpOt1atBST+EMBNYv8fqoUBxohe/E5gMXJrh+EREGk/lnXqlU9PvBaxMul0C9DKzLmb2KDDYzG4HqGvdnsxslJkVmVlRWVlZGmGJSN5TeadeLTP9hCGET4DR+1tXx+MmABMAhgwZEjIdl4jkkfJyn1K5bduoI8k66fT0S4E+Sbd7J9aJiESregoGs6gjyTrpJP05wEAzG2BmrYErgeczE5aISBo0BUO9GjpkcxLwJnCUmZWY2Y0hhErgZmA6sBiYEkJY2HShiog0kCZbq1eDavohhJH1rJ8GTMtoRCIi6Sovh65do44iK2kaBhGJH5V36qWkLyLxo/JOvZT0RSReQvCevpJ+nZT0RSRetmyBqiqVd+qhpC8i8aIpGPZJSV9E4kVJf5+U9EUkXnQBlX1S0heReFFPf5+U9EUkXpT09ynjs2yKSMPs2gVLlsA778CBB8LgwdCvX91zhFVVwbJlUFzsOa283KsYGzbA9u0wahQcd1zzb0NWSirvbN8Or70GL7wABQUwbhz06hVteFFT0hdJU2UlrF8PLVp4Kyjwtnu3r//kE28ffwzr1sGCBTB/vi937Kj9XJ07w6c/7V8AnTrBokWwcCG8//7e9wWfOTgEKCz0dsUVzbPN2WxrSTkHmnHZDR2Z/gps3Qrt2vmXbGEh3HYbfO970L59hDFuhZtvhjZt4Jpr4JRTmm9CUAsh+6auHzJkSCgqKoo6DMlzH33kibtHDzjggNq/W74cXn7Z24wZNZ3LhujSxZP6iSf68oQTPAnMmwdvv+3L997zJN+3LxxzTE078kifUqZTJ+jY0eNavRq+8hX45z9h7Fi4+25omWfduY0b4dlnYfJkGPHyd7gm/I7jem/g4ovhkkvgzDP97zRuHEyZ4v/T//5vuP56/4IG/1IoL4dt2+rf48pUrBddBG++6f+/igro3x+uusq/AD71qdSf28zmhhCG7PM+SvoitYUA993nPcJqXbt6WaBnTy+xLF3q63v3hvPP9wQOXoapqvIvC4CDD/Yk37Vr7eX+EsquXZ70G9ob3bkTbr0VHn4Yzj7bk1+3bo3b7mwVgn8RzprlCbp165q2cyc8/zxMm+Y/DxgAz7S7lkHr36B1yYd1/p3/7//8y/Ff//L/X4sWvke2ZUvNfT7/ebj/fv9SzqSPP/b3y7vvwh//CMOHw1/+An/4A7zyir9vTj8dXn89tS8dJX2RRqqqgltugfHj4ctfhgsugNLSmrZqlfcSzz8fhg2Do47Krut0TJwIo0dD9+7wn//pewq9enk7+ODmibWy0nvL7dt7Qk3V++/DpEneqr9k69KjB1x+OYwcCUOHgl08wv9R8+bV+5gQvMc/ZQp06OB/m86dfbl1K/z8516Su+EG+NnP/O+ZrlWr4Lzz/NjMM8/AhRfW/v2aNfD0077XeOedqb2Gkn6GrFnjPbR822XON1u3+i728897L/+ee9JLWlEpKoLLLoMVK2qvb9PGE2SPHnDooTWtZ08YNMhbx457P9+OHX7AedEiWLsWNm+u3TZu9B5sdasePFNQ4Im0SxdvBx/sVzCs/uIx81Z9LCS5vfeel7rM4KyzPKGPGOGP37nTY9q503vGAwfWlGgAOO003w149dWU/4YbNsBPfwoPPujHTe64A770Jc8BLVv667Vs6a9fUVG77d7t+eKQQ3ybW7SADz+Ec8/1YzovvODlpqagpJ+iykrfBXzxRf8Hvf++19xuucW/+Tt0iCw0aSJr18LFF8PcufDAA36QLZdVVtbeQ0neU1m71uvba9Z4WSNZ7941xw5KS/0gcnGx7wEla93aPwcdOvgXRbdutUtY7dr5l0H1QexPPvHXqqz0x4dQu+3eXbtV994vv9y/lBrlmGPg6KO9O52mDz7wDsALL6T2+IIC/3tUVPjPL73keyNNRUm/kRYvhv/5H68Prl/vvYozzvCext/+Bm+84W/wb34T/uM/ohn6FYLXax94wD+IrVp5q65xdu/uB6H69vXWr59/MPf8kHXu7MME80UIXnLYtMn/t2vW1LS1a+FPf/Ll5Ml+4C9f7NgBJSU1o4QWLvRRRUuXerJNPoh8zDH+nu/Qwd9rWamy0r+BvvIVmDAhY0/75ptelqmsrGlVVb4n0rZt7WbmX3Lr1tW0bdvg+99v+mG1SvoNtH27J/u77/Z/2he+4L2+YcPgoINq7jd7Nvzyl96BaNECzjnHh9dVD7E77LCaXdctW7ynVFLiZYPzzvPnTkdpqddrX3zRP4B9+/ou7q5dNQf+Vq/2tj+dO/vBowsuSC+mqJWXex3797/3nmV1yaC6bLB9u6/ftGnv3mq1Aw6AI47w4XxN2QuTZjBtmg+NefZZ/yDnGSX9BpgxA266yXs2V1/tSX1/B20+/BAeesgfu2hRzS5rx47eO1q1yhNNsu7dfcTA6NF1l4dWroTnnvPnOvlk/xKpHiYYAjz2mO9m7trlQ83GjNmjjplkxw7/glixwtu2bbWTIcAjj3jd9Cc/8QN+9dWuS0q8N7xli7fNm/1L7IwzvJaaKUVFvivdrZvXQg85xHeLW7Wq+/7z5vlIlUmTfNf55JP9Sze5XBCC17EPOsj/N9XLzp1r17QPOii7DsZKGi6/3Gv5q1Zl8e5I01HSr8euXb4be9998LvfweGHw69/7b3xxtq+3XeHq8dXr11bM1qid29v27f7l8nf/+4JZ8wY+M53vMzwzDPe5syp/bytW/swwJNP9ud/9VU/+PPb33qvNF3btsE3vuHDxi65BJ56quYgXmWl1zDHj6//WFinTn6mY/VQxXQUFvoZpXX1xDt2rKkdH3SQL9ev97/1gQf6uOZvfcvHukueW7/eDwaMHu31zzykpJ+wfr3X5OfM8TZvnifiVq3gBz/wnm66pZeGeOstLyM991zNGGOAIUN8eOCXv+wHwGbP9jHEs2d7vC1b+hCyr389s6NJQvA9lrFjfXzzY4/5WOjf/Mb3PPr08eMXxx/vybZ9e28VFXDppb58/XUf9ZHq6//v//rIiGHD/IuxvLymDlpW5rXRzZu9PFM9WgR8NMd119U92kTy1COPwLe/7R/wTA+wzxFK+nhiOekk74m3bev19898xmu3p5/uPfHm9t57nmD79/dhYP361X/fykovVTTlnuqsWX7ca+1av33uuf7ZGTGi/mGqS5f6CSxmMHNm4/c+du/20VAPPeRltcLCvNwbl0z67Ge9Nzd/ft7W6xqS9AkhZF076aSTQqbMnu1V3nvvDWHXrow9beyUlobw4IMhLF7c8McsWBBCly4h9O0bwvLlDX/c9u0hXH65/1/Gjg2hqqrx8YrUsnChv6Huuy/qSCIFFIX95NccPPWkcQoLvYc/apROrtqXnj39OMPRRzf8Mccc46eOb9zoewf7GjW0Y4eXt8aP9yGwU6bAL37hLRdPgJIs8+ST/gG/+uqoI8l6sU6DFRU+7vrLX1btt6kMHuwnnJx3np/Q06tXzRmYXbr4cZP5873t2uWP6dHDh1jq8ykZUVnpIzIuvNCHfck+xTrpP/us90K/9rWoI4m3k0/2UT6FhTVnX370kR9HqajwA8G33urHUYYO9S+GPC25SlN45RXfzfzqV6OOJCc0S9I3s77Ag8B64IMQwt3N8bpPPOEHS5tqngup8ZnPeBNpdhMn+m7lRRdFHUlOSLmaamaFZrbOzBbssX64mS0xs2IzG5dYfRwwNYRwA9AsY6lWrPCTp776VdWMRWKrvNzHQF99tYZ/NVA66XAiMDx5hZkVAA8DFwCDgJFmNgj4F3Cjmb0KvJTGazbYk0/6cM3rr2+OVxORSDz9tI8SUGmnwVJO+iGEmXi5JtlQoDiEsCyEsBOYDFwKfA24M4RwNtDk+2C7d/se39lne3lHRGJq4kQ/aJSJU8PzRKYLH72AlUm3SxLrXgL+w8weBZbX9UAzG2VmRWZWVFZWllYQM2f6/Dg33JDW04hItqqq8iuNzJ7tvXyNDGiwZjmQG0JYAFy2n/tMACaAn5GbzusVFvo8LV/8YjrPIiJZadUqv9rN6697wr/ppqgjyimZTvqlQJ+k270T65rNpk0wdSpce21+zRcvkhdeesk/3Nu2+YG7666LOqKck+nyzhxgoJkNMLPWwJXA8xl+jX2aMsXHhqu0IxIjlZVw++1+AYhDD/W5uJXwU5LOkM1JwJvAUWZWYmY3hhAqgZuB6cBiYEoIYWFmQm2YJ56AT31KF8MQiY3SUp+74+67fT7wt97yD7mkJOXyTghhZD3rpwHTUo4oDatX+1TEP/uZjuuIxMLLL/sY/IoKv9TbVVdFHVHOi9U0DD16+DGe+q62JCI5oqoK7rrLe3DHHOMXMW7MbIBSr1glffBL7IlIDluzxnv0r73mE2eNH69RGRkUu6QvIjls0SI/WFtW5gfodKZtxinpi0h2mDXLL9jcpg288YZf5k4yTlORiUj0pk71izJ07w5vvqmE34SU9EUkWvffD5df7hez/uc/NWFWE1N5R0Sa3urVfnWr9ev9EmrVbfVqnxr5i1/0IZlt20Ydaewp6YtI0ykuhp//3GfD3LnT57xv2dLHVbdq5bfHjoV77oGCgqijzQtK+iKSeW+/7WfQTp3qyf2GG+C22+Dww6OOLO8p6YtI6nbt8gS/YAG89563BQtg7Vqf6vb734cxY3y+HMkKSvoikprXX4dRo+CDD/x227Z+9uyFF/pB2WuugY4do41R9qKkLyKNU17uPfjHHoMBA+D3v4fPftZ/Vl0+6ynpi0jDhADPPAPf+Q6sW+c1+rvugnbtoo5MGkFJX0T2b9Uq+Na3fHjl4MHw17/qBKocpZOzRKR+IXgZZ9AgmD4d7r3X57NXws9Z6umLSN2WLfMDtTNmwBlnePI/4oioo5I0qacvIrWtWgW/+AUcd5z36n/zG3j1VSX8mFBPXyTfbdgA//iH9+hnzIDFi339RRfBo49C796RhieZpaQvkq927oQf/hDuu8+vVHXggfD5z8ONN8K558Lxx+u6ozGkpC+SjxYv9mvPvv22T5Fw/fVw8sk+F47EmpK+SD4JwUs2Y8f6+PrnnvMLl0je0IFckXxQVQULF3qC/9a3vIzz3ntK+HlIPX2RXBcCbN8OGzfCpk2+3LjR58R55x2YP98TfEUFHHAAPPAA3HwztFCfLx8p6Yvkiq1bvQb/wQfeli71ZXGxJ/26dO7sZ9DedBOceKL38Pv1a964Jaso6Ytkq08+8QuEz5rlbd48qKz037Vq5XPTH3kkDBsG3br5VMYdO9Ys+/eHPn00AkdqaZakb2ZnAj8FFgKTQwj/aI7XFclqlZUwd65PUfzhh1BWBh9/XLNct87v17o1DB0K3/senHKKT4nQt69fgUqkkVJ+15hZITACWBdCODZp/XDgAaAAeCyEcDcQgC1AG6AkrYhFctWmTbBoEcyc6SdDzZoFW7b477p189a1Kxx9tC/794fTToPPfAbatIkycomRdLoKE4HxwFPVK8ysAHgYOA9P7nPM7HlgVgjhdTPrDtwHXJ3G64pkp+oLfa9aBaWlsHIlLFkC77/vy9Wra+77qU/BtdfCWWf5vDaHHBJd3JJXUk76IYSZZtZ/j9VDgeIQwjIAM5sMXBpCWJT4fTlwQKqvKZJ1Kirg9tth0qSackyyzp29537++b486ig/CUqXD5SIZLoo2AtYmXS7BPismX0JOB/ohO8diOS+d9+Fq67y8e9XXOG99549oVcvX/buDV266ECqZJVmORIUQvgz8Od93cfMRgGjAPr27dscYYmkJgR46CG/ZGDnzj7P/LBhUUcl0iCZPjujFOiTdLt3Yt1+hRAmhBCGhBCGdOvWLcNhiWTAjh0+yuaii2DMGDjvPO/tK+FLDsl0T38OMNDMBuDJ/krgqgy/hkj6Kirgo49qWmlpzZms1We1btrk0w5Xt+oToNq0gfHjfToDlW4kx6QzZHMScCbQ1cxKgDtDCI+b2c3AdHzIZmEIYWFGIhVpqMpK2LzZW2mpn7Ga3KrHxO+pffuak5uqW58+0KlTTevY0Xv4Rx7Z/NslkgHpjN4ZWc/6acC0lCMS2Zfdu2H58pppCKrbv/8N5eWe6Hfs2PtxLVr4CU2HHw5f+IJPRdCvn6/r188PvLZq1eybI9LcdEqfZK8QvPQyZw4UFfly7lwvvVRr3x4GDvR5Zbp2hQ4dfF31skcPv8xf//4+2ZhInlPSl+yyebNfj/Vvf4OXXoIVK3x9q1Z+JaeRI+Gkk3y8+8CB0L276uoijaCkL9Hbtg2eeAKeecYnGNu1y3vp55wDt93m886ccIJ66iIZoKQv0Skvh4cf9vndP/4Yjj0WvvtdGD4cTj1Vl+4TaQJK+tL8Vq+GX/3KL9u3ebOPex83zicXE5EmpaQvTa+01Ms21e3dd339FVd4sj/++GjjE8kjSvqSvu3bfRjl8uU+s2RyW7rUR+CAX4j7c5+DH/0IrrnGh0+KSLNS0peG2bbNE/iSJTWX6Fu2zFvpHjNtmPm49z59/KIft97qpZsTTtCFP0Qipk9gvgvBz1AtKvLZIjdt8gt7bN3qy02bPLGvXFn7cT17ek/93HPhsMO8VV+eTyc6iWQtJf18sGOHz/W+di2sWePL4mJP9HPn+iga8B56hw5ehmnfvmZ5xhk+7cBRR9WMjz/wwGi3SURSoqQfZ0uX+tWZZs/e+3ctW/oB1K98BYYM8ROejj1WwyRFYk5JP65+/3u46SZP4nfe6Rf2OPRQP4P10EO9KcGL5B0l/bjZuhVuvhkmToTTT4c//tGv4CQiQuYvoiJReucdL9M8+aQPi3z1VSV8EalFPf1cFgK8/Ta8+CL89a/w1ls+q+SMGXDWWVFHJyJZSEk/F5WVwX/9F7zwAqxa5aNuhg6Fn/wERo8GXW5SROqhpJ9rNm+GCy6ABQvg4othxAi/fcghUUcmIjlAST+X7NgBX/wizJ8Pzz3nE5WJiDSCkn6u2L0brrvO6/VPPqmELyIp0eidXBACjBkDU6bAz3/uyV9EJAVK+rngZz+D8eNh7Fi/kpSISIqU9LPd44/7SJ1rr4V77406GhHJcUr62WzGDB+Cef75nvxb6N8lIulRFslW778Pl13ms1o+/bSmKhaRjFDSz0Yff+zj71u18rNtO3aMOiIRiQkN2cw2O3bAl74EJSXw2mt+YRIRkQxplqRvZl8ALgIOAh4PIbzcHK+bc0KAUaNg1iyfHfNzn4s6IhGJmZTLO2ZWaGbrzGzBHuuHm9kSMys2s3EAIYS/hBC+AYwGrkgv5JgKwee9f+opuOsuGDky6ohEJIbSqelPBIYnrzCzAuBh4AJgEDDSzAYl3eWHid9LsvJyn17hpz+F66/3aZFFRJpAykk/hDATWL/H6qFAcQhhWQhhJzAZuNTcPcDfQgjzUg83ht56CwYPhmnT4Fe/giee8FkzRUSaQKZH7/QCVibdLkms+w5wLnCZmY2u64FmNsrMisysqKysLMNhZaEQ4P774bTT/PYbb8Attyjhi0iTapYDuSGEB4EH93OfCcAEgCFDhoTmiCsSW7fCyy/DhAnw0ktw6aXeu+/cOerIRCQPZDrplwJ9km73TqzLb2vW+AVPnnsO/v53H5bZqZOXc8aMUe9eRJpNppP+HGCgmQ3Ak/2VwFUZfo3sFQJ89JHPd5/cli/33/fv79MqXHqpl3V0lq2INLOUk76ZTQLOBLqaWQlwZwjhcTO7GZgOFACFIYSFGYk0m1RVwYoVPlXCokW12+bNfh8zOPJIv4zh6NF+davjjlOvXkQilXLSDyHUOZA8hDANmJZyRNlm2zYvy7zzDixZAh98AMXFsHNnzX169IBBg+CrX/XliSd6gm/XLrKwRUTqomkY6rN0KTz6qB9kLS/3UswRR3jvfcQIXx59tCd5HYQVkRyhpJ+sstLHyz/yCEyfDi1b+jw4N93kNfiW+nOJSG5TFgOf3Oyxx3zO+pIS6NkTfvxj+PrX/WcRkZjI36S/e7ePk3/0UfjrX33kzbBh8MADcPHFGlkjIrGUn0k/BLjmGpg0Cbp3hx/8AL7xDRgwIOrIRESaVH4m/Z/8xBP+j34Ed9wBrVtHHZGISLPIv6T/9NM+dfF11/lS4+ZFJI/k1+US33rLx9KfeqrPfaOELyJ5Jn+SfkmJT39w6KHw7LNwwAFRRyQi0uzyo7yzdStccokvX3kFunWLOiIRkUjEP+mH4CWdd97xmS6PPTbqiEREIhP/pD9+PEydCvfcAxdeGHU0IiKRindNv6gIbrvN58q57baooxERiVx8k/7GjXDFFX7y1cSJ0CK+myoi0lDxLO+E4PPmrFgBM2dCly5RRyQikhXimfQfecTr+PfeC6ecEnU0IiJZI341j3nz4NZb/aDt2LFRRyMiklXilfQ3bYLLL4dDDoEnn1QdX0RkD/HKirt3+6UKJ0+Grl2jjkZEJOvEq6bfqZPX8kVEpE7x6umLiMg+KemLiOQRJX0RkTyipC8ikkeU9EVE8oiSvohIHlHSFxHJI0r6IiJ5xEIIUcewFzMrA1ak8RRdgY8zFE42iNv2QPy2KW7bA/HbprhtD+y9Tf1CCPu8HmxWJv10mVlRCGFI1HFkSty2B+K3TXHbHojfNsVteyC1bVJ5R0Qkjyjpi4jkkbgm/QlRB5BhcdseiPbk/5UAAAOmSURBVN82xW17IH7bFLftgRS2KZY1fRERqVtce/oiIlKHWCV9MxtuZkvMrNjMxkUdTyrMrNDM1pnZgqR1B5vZK2a2NLHsHGWMjWFmfczsNTNbZGYLzWxMYn0ub1MbM3vLzN5JbNOPE+sHmNnsxPvvaTNrHXWsjWFmBWb2tpm9mLid69uz3MzeM7P5ZlaUWJfL77tOZjbVzN43s8Vm9rlUtic2Sd/MCoCHgQuAQcBIMxsUbVQpmQgM32PdOGBGCGEgMCNxO1dUAmNDCIOAk4FvJ/4vubxNO4CzQwgnACcCw83sZOAe4FchhCOAcuDGCGNMxRhgcdLtXN8egLNCCCcmDWvM5ffdA8BLIYSjgRPw/1XjtyeEEIsGfA6YnnT7duD2qONKcVv6AwuSbi8BeiR+7gEsiTrGNLbtOeC8uGwTcCAwD/gsfpJMy8T6Wu/HbG9A70TSOBt4EbBc3p5EzMuBrnusy8n3HdAR+JDEcdh0tic2PX2gF7Ay6XZJYl0cdA8hrE78vAboHmUwqTKz/sBgYDY5vk2JUsh8YB3wCvBvYEMIoTJxl1x7/90PfB/YnbjdhdzeHoAAvGxmc81sVGJdrr7vBgBlwBOJEtxjZtaOFLYnTkk/LwT/Ss+5IVdm1h54BrglhLAp+Xe5uE0hhKoQwol4D3kocHTEIaXMzEYA60IIc6OOJcNOCyF8Gi/5ftvMPp/8yxx737UEPg38OoQwGNjKHqWchm5PnJJ+KdAn6XbvxLo4WGtmPQASy3URx9MoZtYKT/h/CCH8ObE6p7epWghhA/AaXv7oZGYtE7/KpfffqcAlZrYcmIyXeB4gd7cHgBBCaWK5DngW/3LO1fddCVASQpiduD0V/xJo9PbEKenPAQYmRhy0Bq4Eno84pkx5Hrg+8fP1eF08J5iZAY8Di0MI9yX9Kpe3qZuZdUr83BY/RrEYT/6XJe6WM9sUQrg9hNA7hNAf/9y8GkK4mhzdHgAza2dmHap/BoYBC8jR910IYQ2w0syOSqw6B1hEKtsT9QGKDB/suBD4AK+v3hF1PCluwyRgNbAL/3a/Ea+vzgCWAn8HDo46zkZsz2n4Lue7wPxEuzDHt+l44O3ENi0AfpRYfxjwFlAM/Ak4IOpYU9i2M4EXc317ErG/k2gLq/NBjr/vTgSKEu+7vwCdU9kenZErIpJH4lTeERGR/VDSFxHJI0r6IiJ5RElfRCSPKOmLiOQRJX0RkTyipC8ikkeU9EVE8sj/AzXiXs/A7eTGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Bias Variance Tradeoff\n",
    "\n",
    "biases = np.zeros(m_train-1)\n",
    "variances = np.zeros(m_train-1)\n",
    "for r in range(1, m_train):\n",
    "    \n",
    "    mu_hat = np.mean(X_train, axis=0)\n",
    "    Z = X_train - mu_hat ##using broadcasting rules\n",
    "\n",
    "\n",
    "    [U_full, S_full, V_full] = np.linalg.svd(Z)\n",
    "    V = V_full[:, :r]\n",
    "    B = X_train @ V\n",
    "\n",
    "    B_tilde= np.hstack((np.ones((60,1)), B))\n",
    "\n",
    "    theta_hat_full = np.linalg.pinv(B_tilde) @ Y_train\n",
    "\n",
    "\n",
    "    bias_hat = theta_hat_full[0]\n",
    "    theta_hat = theta_hat_full[1:]\n",
    "\n",
    "    y_hat_test = X_test @ V @ theta_hat + bias_hat\n",
    "    \n",
    "    tmp_bias = 0\n",
    "    tmp_var = 0\n",
    "    \n",
    "    pinv_matrix = X_test[:, :r].T @ X_test[:, :r]\n",
    "    [u_tmp, s_tmp, v_tmp] = np.linalg.svd(pinv_matrix)\n",
    "    \n",
    "    for ii in range(m_test):\n",
    "        x_test = X_test[ii:ii+1, :r]\n",
    "        y_test = Y_test[ii]\n",
    "        tmp_bias += np.dot(np.squeeze(x_test), np.squeeze(theta_true[:r]) - np.squeeze(theta_hat)) ** 2\n",
    "        tmp_var += np.linalg.norm(x_test, 2) ** 2 / s_tmp[-1]  \n",
    "    biases[r-1] = tmp_bias\n",
    "    variances[r-1] = tmp_var\n",
    "\n",
    "    \n",
    "plt.figure\n",
    "plt.plot(biases, 'b', variances, 'r')\n",
    "plt.yscale(\"log\")\n",
    "plt.show()"
   ]
  }
 ],
 "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}