singular_value_decomposition.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.datasets import load_digits\n",
    "from sklearn.linear_model import LinearRegression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Singular Value Decomposition\n",
    "\n",
    "A more **general form of eigendecomposition** that generalizes decomposition from a square matrix to any $m \\times n$ matrix via an extension of **polar decomposition**. \n",
    "\n",
    "If $A$ is an $m \\times n$ matrix with rank $r,$ then there exists a factorization of $A$ as\n",
    "\n",
    "\\begin{equation}\n",
    "A=U \\Sigma V^{T}\n",
    "\\end{equation}\n",
    "\n",
    "where $U(m \\times r)$ and $V(n \\times r)$ are orthogonal matrices, and $\\sum(r \\times r)$ is a matrix with only non-zero diagonal elements.\n",
    "\n",
    "$\\Sigma=\\left(\\sigma_{i j}\\right)$, where $\\sigma_{i i}=\\sigma_{i}$, for $i=1, \\ldots, r$ are the singular values of $A$ such that $\\sigma_{1} \\geq \\sigma_{2} \\geq \\ldots \\geq \\sigma_{r} \\geq 0$, and all non-diagonal entries of $\\sum$ are zero. \n",
    "\n",
    "The singular vector corresponding to the largest singular value represents **the 1-D subspace that retains the most information**, then the second largest and so on.."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Calculation of singular vectors and singular values**:\n",
    "\n",
    "$C=A A^{T}$\n",
    "$=\\left(U \\Sigma V^{T}\\right)\\left(U \\Sigma V^{T}\\right)^{T}$\n",
    "$=U \\Sigma V^{T} V \\Sigma^{T} U^{T}$\n",
    "$=U \\Sigma^{2} U^{T}$\n",
    "\n",
    "Since $C$ is symmetric which means it is diagonalizable. Also, we know that $U$ and $V$ are orthogonal matrices. \n",
    "\n",
    "Hence, $C$ can be represented as:\n",
    "\n",
    "$C= U \\Sigma^{2} U^{-1}$\n",
    "\n",
    "- This implies that eigenvectors of $C$ are left-singular vectors $U$ of $A$. \n",
    "- Likewise, eigenvectors of $A^{T}A$ are right-singular vectors $V$ of $A$. \n",
    "- And, square root of eigenvalues of $C$ are singular values of $A$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5, 5)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X=np.random.normal(0,1,25).reshape(5,-1)\n",
    "\n",
    "#symmetric matrix\n",
    "X = (X + np.transpose(X)) / 2\n",
    "X.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((5, 5), (5, 5), (5, 5))"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U = np.linalg.svd(X)[0]\n",
    "\n",
    "S = np.diag(np.linalg.svd(X)[1])\n",
    "\n",
    "VT = np.linalg.svd(X)[2]\n",
    "\n",
    "U.shape, S.shape, VT.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[2.45202372, 0.        , 0.        , 0.        , 0.        ],\n",
       "       [0.        , 1.44188024, 0.        , 0.        , 0.        ],\n",
       "       [0.        , 0.        , 1.17519512, 0.        , 0.        ],\n",
       "       [0.        , 0.        , 0.        , 0.66631124, 0.        ],\n",
       "       [0.        , 0.        , 0.        , 0.        , 0.33384111]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "C = X@np.transpose(X)\n",
    "\n",
    "eigvalue = np.linalg.eig(C)[0]\n",
    "eigvector = np.linalg.eig(C)[1]\n",
    "\n",
    "eigvalue_vector_map = dict(zip(eigvalue, eigvector))\n",
    "\n",
    "#sorted in eigenvalues\n",
    "eigvalue_vector_map = {k: eigvalue_vector_map[k] for k in np.sort(eigvalue)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[2.45202372, 0.        , 0.        , 0.        , 0.        ],\n",
       "       [0.        , 1.44188024, 0.        , 0.        , 0.        ],\n",
       "       [0.        , 0.        , 1.17519512, 0.        , 0.        ],\n",
       "       [0.        , 0.        , 0.        , 0.66631124, 0.        ],\n",
       "       [0.        , 0.        , 0.        , 0.        , 0.33384111]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.diag(np.sort(np.sqrt(eigvalue))[::-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.82391094, -0.31558874,  0.39671631, -0.25250596, -0.02077056],\n",
       "       [ 0.27556414, -0.82249527, -0.06294942, -0.01066225,  0.49344668],\n",
       "       [-0.07676445,  0.20590547, -0.48247752, -0.79018774,  0.30745559],\n",
       "       [ 0.26967162,  0.33441897,  0.72994658, -0.19183933,  0.49579919],\n",
       "       [-0.40819065,  0.26395486, -0.27027203,  0.52433584,  0.64477374]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.82391094,  0.31558874, -0.39671631,  0.25250596, -0.02077056],\n",
       "       [ 0.27556414,  0.82249527,  0.06294942,  0.01066225,  0.49344668],\n",
       "       [-0.07676445, -0.20590547,  0.48247752,  0.79018774,  0.30745559],\n",
       "       [ 0.26967162, -0.33441897, -0.72994658,  0.19183933,  0.49579919],\n",
       "       [-0.40819065, -0.26395486,  0.27027203, -0.52433584,  0.64477374]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.array([eigvalue_vector_map[k] for k in np.sort(eigvalue)[::-1]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Notice that the magnitudes are same but direction of eigenvectors can be different**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.82391094, -0.31558874, -0.39671631, -0.25250596,  0.02077056],\n",
       "       [-0.27556414, -0.82249527,  0.06294942, -0.01066225, -0.49344668],\n",
       "       [ 0.07676445,  0.20590547,  0.48247752, -0.79018774, -0.30745559],\n",
       "       [-0.26967162,  0.33441897, -0.72994658, -0.19183933, -0.49579919],\n",
       "       [ 0.40819065,  0.26395486,  0.27027203,  0.52433584, -0.64477374]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.transpose(VT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.82391094,  0.31558874, -0.39671631,  0.25250596, -0.02077056],\n",
       "       [ 0.27556414,  0.82249527,  0.06294942,  0.01066225,  0.49344668],\n",
       "       [-0.07676445, -0.20590547,  0.48247752,  0.79018774,  0.30745559],\n",
       "       [ 0.26967162, -0.33441897, -0.72994658,  0.19183933,  0.49579919],\n",
       "       [-0.40819065, -0.26395486,  0.27027203, -0.52433584,  0.64477374]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C = np.transpose(X)@X\n",
    "\n",
    "eigvalue = np.linalg.eig(C)[0]\n",
    "eigvector = np.linalg.eig(C)[1]\n",
    "\n",
    "eigvalue_vector_map = dict(zip(eigvalue, eigvector))\n",
    "\n",
    "#sorted in eigenvalues\n",
    "eigvalue_vector_map = {k: eigvalue_vector_map[k] for k in np.sort(eigvalue)}\n",
    "\n",
    "np.array([eigvalue_vector_map[k] for k in np.sort(eigvalue)[::-1]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Again, magnitudes are same but direction of eigenvectors can be different**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Low Rank Approximation of SVD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((5, 5), (5, 5), (5, 5))"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U = np.linalg.svd(X)[0]\n",
    "\n",
    "S = np.diag(np.linalg.svd(X)[1])\n",
    "\n",
    "VT = np.linalg.svd(X)[2]\n",
    "\n",
    "U.shape, S.shape, VT.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_approximated_data(n_components):\n",
    "    k = n_components\n",
    "    \n",
    "    return U[:,:k] @ S[:k, :k] @ VT[:k]\n",
    "\n",
    "def rmse(X,Y):\n",
    "    return np.sqrt(np.mean(np.square(X - Y)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rmse</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>n_components</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.633339</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.400775</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.278317</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.149053</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.066768</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  rmse\n",
       "n_components          \n",
       "0             0.633339\n",
       "1             0.400775\n",
       "2             0.278317\n",
       "3             0.149053\n",
       "4             0.066768"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "RMSE=[]\n",
    "\n",
    "for n_components in range(min(X.shape)):\n",
    "    Y = find_approximated_data(n_components)\n",
    "    \n",
    "    RMSE.append(rmse(X,Y))\n",
    "    \n",
    "df = pd.DataFrame(RMSE, index=range(n_components+1), columns = ['rmse'])\n",
    "df.index.name='n_components'\n",
    "\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Elbow method to find optimal number of components**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAhk0lEQVR4nO3deXRU9f3/8ec7kUURJi6pKAFRk0ADsviNiFjQLxBAg4CCAgrIUtAvUEAWf6yKLFWL4gYuWEmwWlkKIkECJqhVqK0EBCESQkRQWRQXwmKR7fP7g0gjjTCBSe7M5PU4J+dwl8x9nTnw4p47c9/XnHOIiEjoi/A6gIiIBIYKXUQkTKjQRUTChApdRCRMqNBFRMLEOV4d+OKLL3Y1a9b06vAiIiFp9erV3zrnoova5lmh16xZk6ysLK8OLyISksxs269t0yUXEZEwoUIXEQkTKnQRkTChQhcRCRMqdBGRMKFCFxEJEyp0EZEwEXKFnpeXx6hRozh27JjXUUREgkrIFfrChQt59NFHGTp0KJrlLiLyH57dKXqmhg0bxvbt23nqqaeIiopi/PjxXkcSEQkKIVfoZsYTTzxBfn4+Dz/8MD6fj/vvv9/rWCIingu5QgeIiIhgxowZ7N27l6FDh+Lz+ejdu7fXsUREPBWShQ5wzjnn8Nprr7Fv3z769u1LlSpV6NSpk9exREQ8E3IfihZWoUIFFixYQOPGjbnrrrtYtmyZ15FERDwT0oUOUKlSJd566y3q1KnDbbfdxooVK7yOJCLiiZAvdICoqCiWLVtG9erVSU5O5uOPP/Y6kohIqQuLQgf4zW9+Q0ZGBj6fj9atW7Np0yavI4mIlCq/Ct3M2pjZJjPLM7ORv7LPnWb2qZllm9lfAxvTPzVq1CAzMxMzIykpiS+++MKLGCIinjhtoZtZJDAduBlIALqaWcJJ+8QBo4AbnHN1gCGBj+qf+Ph43n77bfbu3UvLli35+uuvvYoiIlKq/DlDbwTkOee2OOcOAbOB9ift0xeY7pz7AcA5901gYxZP/fr1WbJkCdu3b6d169b88MMPXsYRESkV/hR6NeDLQstfFawrLB6IN7OVZvZPM2tT1AuZWT8zyzKzrN27d59ZYj81adKEhQsXsnHjRpKTkzlw4ECJHk9ExGuB+lD0HCAOuAnoCrxkZlEn7+Scm+GcS3TOJUZHRwfo0L8uKSmJ119/nX/961/cdttt/PTTTyV+TBERr/hT6NuB6oWWYwrWFfYVsMg5d9g59zmQy/GC99ztt9/Oyy+/TEZGBl27duXIkSNeRxIRKRH+FPoqIM7MrjCz8kAXYNFJ+yzk+Nk5ZnYxxy/BbAlczLPTs2dPnnrqKd544w1+//vfa5a6iISl085ycc4dMbOBwDIgEpjpnMs2swlAlnNuUcG2Vmb2KXAUGOGc+64kgxfX4MGDyc/P56GHHqJKlSo8/fTTmJnXsUREAsav4VzOuSXAkpPWPVjozw4YWvATtMaNG8eePXt48sknueCCC3j44Ye9jiQiEjAhO23xTBSepT5hwgR8Ph9Dhwb1/0EiIn4rU4UOx0v951nqw4YNw+fz0adPH69jiYictTJX6ACRkZG8+uqr7Nu3j379+lGlShXuuOMOr2OJiJyVsBnOVVwVKlRg/vz5XH/99dx9990sXbrU60giImelzBY6HJ+lvnjxYurUqcPtt9+uWeoiEtLKdKHDf2ap16hRg+TkZNasWeN1JBGRM1LmCx3+M0s9KiqK1q1bk5OT43UkEZFiU6EXqF69OpmZmURGRpKUlMS2bdu8jiQiUiwq9ELi4uJ4++232b9/v2api0jIUaGfpF69eixZsoQdO3bQqlUrzVIXkZChQi/C9ddfz8KFC8nJySE5OZn9+/d7HUlE5LRU6L8iKSmJ2bNnn5ilfvDgQa8jiYickgr9FG677TZmzpxJZmamZqmLSNBToZ/GPffcw9NPP83ChQvp06ePZqmLSNAqk7NcimvQoEHk5+fz4IMP4vP5NEtdRIKSCt1PY8eOZc+ePUydOpWoqCgmTJjgdSQRkV9QofvJzHj88cfJz89n4sSJ+Hw+hg0b5nUsEZETVOjFYGa8+OKL7N27l+HDh+Pz+fj973/vdSwREUCFXmxFzVK/8847vY4lIqJvuZyJ8uXLM3/+fG644Qa6detGenq615FERFToZ+q8885j8eLF1K1bl44dO/LBBx94HUlEyjgV+lnw+XwsW7aMyy+/nLZt22qWuoh4SoV+lqKjo8nIyOCCCy7QLHUR8ZQKPQBiYmLIyMggMjKSli1bsnXrVq8jiUgZpEIPkJ9nqR84cICkpCR27drldSQRKWNU6AFUr1490tPT2blzJ61ateL777/3OpKIlCF+FbqZtTGzTWaWZ2Yji9je08x2m9nagp8ye7dN48aNWbhwIZs2beKWW27RLHURKTWnLXQziwSmAzcDCUBXM0soYtc5zrkGBT9/DnDOkNKyZUtmz55NVlYWHTp00Cx1ESkV/pyhNwLynHNbnHOHgNlA+5KNFfp+nqW+fPlyunTpolnqIlLi/Cn0asCXhZa/Klh3so5m9omZ/c3Mqhf1QmbWz8yyzCxr9+7dZxA3tPTo0YNnnnmGN998k969e2uWuoiUqEB9KJoG1HTO1QMygFlF7eScm+GcS3TOJUZHRwfo0MHtD3/4AxMnTuQvf/kLgwcPxjnndSQRCVP+DOfaDhQ+444pWHeCc+67Qot/Bv509tHCx5gxY9izZw9PPPEEUVFRTJw40etIIhKG/Cn0VUCcmV3B8SLvAtxVeAczu9Q5t7NgsR2wMaApQ5yZMWXKFPLz85k0aRI+n4/hw4d7HUtEwsxpC905d8TMBgLLgEhgpnMu28wmAFnOuUXAIDNrBxwBvgd6lmDmkGRmvPDCC+zdu5cRI0bg8/no27ev17FEJIyYV9d0ExMTXVZWlifH9tKhQ4fo0KEDS5cu5fXXX6dz585eRxKREGJmq51ziUVt052ipax8+fL87W9/43e/+x3dunVjyZIlXkcSkTChQvfAeeedR1paGvXq1aNjx468//77XkcSkTCgQveIz+dj6dKl1KxZk7Zt27J69WqvI4lIiFOhe+jnWeoXXnghbdq0YeNGfTlIRM6cCt1jMTExZGZmEhkZSVJSkmapi8gZU6EHgdjYWDIyMvjxxx9p2bIlO3fuPP0viYicRIUeJK6++mqWLFnCrl27NEtdRM6ICj2ING7cmDfffJPc3Fxuvvlm9u3b53UkEQkhKvQg06JFC+bMmcPq1as1S11EikWFHoQ6dOhASkoK77zzDl26dOHw4cNeRxKREKBCD1Ldu3dn2rRpmqUuIn7zZ9qieGTAgAHk5+czZswYfD4fzz77LGbmdSwRCVIq9CA3atQo9uzZw5QpU4iKimLSpEleRxKRIKVCD3JmxmOPPcaePXuYPHkyPp+PESNGeB1LRIKQCj0EmBnPP/88e/fu5YEHHiAqKkqz1EXkv6jQQ0RkZCSvvPIK+/bt495776VKlSqapS4iv6BvuYSQ8uXLM2/ePJo2bapZ6iLyX1ToIebnWer169fXLHUR+QUVegiqUqUKS5cu5YorrtAsdRE5QYUeoi6++GIyMjK46KKLaN26NZ9++qnXkUTEYyr0EFatWjUyMjIoV64cSUlJfP75515HEhEPqdBDXGxsLG+//Tb//ve/NUtdpIxToYeBq6++mvT0dL7++muSkpL47rvvvI4kIh5QoYeJ6667jkWLFpGXl8ctt9yiWeoiZZAKPYw0b978xCz19u3ba5a6SBnjV6GbWRsz22RmeWY28hT7dTQzZ2aJgYsoxdG+fXtmzZrFe++9R+fOnTVLXaQMOW2hm1kkMB24GUgAuppZQhH7VQYGA/8KdEgpnrvvvptp06axaNEievXqpVnqImWEP7NcGgF5zrktAGY2G2gPnPzF54nAY4BGAQaB/v37k5+fz+jRo6lcuTLPPvss55yj0T0i4cyfSy7VgC8LLX9VsO4EM7sGqO6ceyuA2eQsjRw5kgceeIAXXniBevXqsWTJEpxzXscSkRJy1h+KmlkEMBUY5se+/cwsy8yydu/efbaHltMwMx599FHeeOMNjhw5QnJyMq1bt+aTTz7xOpqIlAB/Cn07UL3QckzBup9VBuoC75nZVqAxsKioD0adczOcc4nOucTo6OgzTy1+MzM6dOjAhg0beOqpp8jKyqJhw4b07duXXbt2eR1PRALIn0JfBcSZ2RVmVh7oAiz6eaNzLt85d7FzrqZzribwT6Cdcy6rRBLLGSlfvjyDBw8mLy+PwYMHM2vWLGJjY5k0aRI//vij1/FEJABOW+jOuSPAQGAZsBGY65zLNrMJZtaupANKYF144YVMnTqV7OxsWrVqxbhx46hVqxavvvqqvg0jEuLMqw/JEhMTXVaWTuK99v777zN06FBWr15NYmIiU6dOpWnTpl7HEpFfYWarnXNF3uujO0XLuGbNmvHRRx/xyiuvsHPnTpo1a0bHjh3Jy8vzOpqIFJMKXYiIiKB79+7k5uYyYcIEli1bRkJCAkOHDuWHH37wOp6I+EmFLiecd955jBs3js2bN9OjRw+eeuopYmNjeeaZZzRCQCQEqNDlv1x66aX8+c9/5uOPP6Zhw4YMHjyYunXr8uabb+rGJJEgpkKXX1W/fn0yMjJYvHgxERERdOjQgebNm7NmzRqvo4lIEVTockpmRnJyMp988gnTpk1j/fr1JCYm0rNnT7Zv3376FxCRUqNCF7+UK1eOAQMGkJeXx/Dhw3n99deJj49n/PjxHDhwwOt4IoIKXYopKiqKP/3pT2zcuJHk5GQefvhh4uLiSElJ4ejRo17HEynTVOhyRq688krmzp3LihUrqF69Or179yYxMZF33nnH62giZZYKXc7KDTfcwIcffshf//pXvv/+e1q0aEH79u3ZtGmT19FEyhwVupy1iIgIunbtSk5ODo888gjvvvsudevWZdCgQXz33XdexxMpM1ToEjDnnnsuI0eOZPPmzfTp04fp06cTGxvL1KlT+emnn7yOJxL2VOgScJdccgkvvPAC69ato3HjxgwbNow6deowf/583ZgkUoJU6FJi6tatS3p6Ounp6VSsWJFOnTrRrFkzVq1a5XU0kbCkQpcS16ZNG9auXcuLL75Ibm4ujRo1onv37nz55Zen/2UR8ZsKXUrFOeecQ79+/di8eTOjRo1i3rx5xMfHM3bsWPbt2+d1PJGwoEKXUlWlShX++Mc/smnTJm6//XYmT55MXFwcL730km5MEjlLKnTxxOWXX85rr73GP//5T6666ir69etHw4YNycjI8DqaSMhSoYunrrvuOlasWMG8efPYv38/rVq1Ijk5mU8//dTraCIhR4UunjMzOnXqxMaNG5kyZQorV66kXr169O/fn2+++cbreCIhQ4UuQaNChQoMHz6cvLw8/u///o8ZM2YQFxfHY489xsGDB72OJxL0VOgSdC6++GKeffZZNmzYQLNmzRg5ciS1a9dmzpw5ujFJ5BRU6BK0ateuTVpaGpmZmURFRdGlSxeaNGnChx9+6HU0kaCkQpeg16JFC1avXs3LL7/M1q1badKkCV26dGHr1q1eRxMJKip0CQmRkZH07t2bzZs38+CDD7Jo0SJq167NyJEjyc/P9zqeSFBQoUtIOf/883n44YfJzc2lc+fOPPbYY8TGxvL8889z5MgRr+OJeMqvQjezNma2yczyzGxkEdvvM7P1ZrbWzFaYWULgo4r8R0xMDLNmzWLVqlUkJCTQv39/6tWrx5IlS/TBqZRZpy10M4sEpgM3AwlA1yIK+6/Ouaudcw2APwFTAx1UpCiJiYm89957vPHGGxw+fJjk5GRat27N+vXrvY4mUur8OUNvBOQ557Y45w4Bs4H2hXdwzu0ttFgJ0CmSlBozo0OHDmRnZ/Pkk0+SlZVFgwYN6NevH7t27fI6nkip8afQqwGF55x+VbDuF8xsgJl9xvEz9EGBiSfiv/LlyzNkyBDy8vIYNGgQKSkpxMXFMXnyZP797397HU+kxAXsQ1Hn3HTn3FXA/wPGFrWPmfUzsywzy9q9e3egDi3yCxdeeCFPPvkk2dnZtGzZkrFjx1KrVi1ee+01jh075nU8kRLjT6FvB6oXWo4pWPdrZgMditrgnJvhnEt0ziVGR0f7HVLkTMTHx/PGG2/w3nvvER0dTbdu3WjcuDErVqzwOppIifCn0FcBcWZ2hZmVB7oAiwrvYGZxhRaTgc2Biyhydm688UZWrVrFrFmz2LFjB02bNqVTp0589tlnXkcTCajTFrpz7ggwEFgGbATmOueyzWyCmbUr2G2gmWWb2VpgKHBPSQUWORMRERH06NGD3NxcJkyYQHp6Or/97W8ZPnw4P/zwg9fxRALCvPrObmJiosvKyvLk2CI7duxg3LhxpKSkcMEFFzB+/Hjuu+8+ypUr53U0kVMys9XOucSitulOUSmTLrvsMl5++WXWrFlDgwYNGDRoEHXr1mXRokW6MUlClgpdyrQGDRqQmZlJWloaZkb79u1p0aKFbkySkKRClzLPzGjbti3r169n2rRprFu3joYNGzJ48GD27NnjdTwRv6nQRQqUK1eOAQMGkJubS79+/Xj22WeJj49n5syZ+v66hAQVushJLrroIp577jmysrKIi4ujT58+XH/99axatcrraCKnpEIX+RXXXHMNK1as4JVXXmHbtm1cd9119O3bF93lLMFKhS5yCmZG9+7dyc3N5f777yc1NZX4+HimT5+u+esSdFToIn6oUqUKTzzxBOvWreN//ud/GDhwIImJiRojIEFFhS5SDAkJCWRkZDBv3jy+//57mjZtSrdu3dixY4fX0URU6CLFZWZ06tSJjRs3MnbsWObNm0etWrV4/PHHOXTokNfxpAxToYucoUqVKjFx4kSys7O56aabGDFiBPXr1ycjI8PraFJGqdBFzlJsbCxpaWmkpaVx+PBhWrVqRceOHdm2bZvX0aSMUaGLBEjbtm3ZsGEDkydPJj09ndq1azNhwgQ9LUlKjQpdJIAqVqzI6NGjycnJ4dZbb+Whhx6iTp06GvolpUKFLlICatSowdy5c1m+fDnnnnsu7du355ZbbiE3N9fraBLGVOgiJah58+asXbuWqVOnsnLlSurWrcuoUaPYv3+/19EkDKnQRUpYuXLluP/++8nNzeWuu+7i0UcfpXbt2syePVuXYSSgVOgipaRq1aqkpqaycuVKfvOb39C1a1f+93//lw0bNngdTcKECl2klDVp0oRVq1bxwgsvsH79eho0aMCQIUM0e13OmgpdxAORkZHce++95Obm0rdvX5555hlq1apFSkqKZq/LGVOhi3jooosu4vnnnycrK4urrrqK3r1706RJE/QAdTkTKnSRIPDz7PVZs2axdetWGjVqRL9+/fj222+9jiYhRIUuEiQiIiLo0aMHmzZt4v7772fmzJnEx8fz3HPPcfToUa/jSQhQoYsEGZ/Pd2L2esOGDRkwYIBmr4tfVOgiQapOnTpkZmYyd+5cvv32W5o2bUr37t3ZuXOn19EkSKnQRYKYmXHHHXeQk5PDmDFjmDt3LvHx8Zq9LkXyq9DNrI2ZbTKzPDMbWcT2oWb2qZl9YmbLzezywEcVKbsqVarEpEmTyM7O5sYbb9TsdSnSaQvdzCKB6cDNQALQ1cwSTtrtYyDROVcP+Bvwp0AHFZHjs9cXL178i9nrnTp10ux1Afw7Q28E5DnntjjnDgGzgfaFd3DOveuc+7Fg8Z9ATGBjikhhP89enzRpEkuWLOG3v/0tkyZN4uDBg15HEw/5U+jVgC8LLX9VsO7X9AHSi9pgZv3MLMvMsnbv3u1/ShH5LxUrVmTMmDHk5OTQtm1bxo0bR506dUhLS9PQrzIqoB+Kmlk3IBGYUtR259wM51yicy4xOjo6kIcWKbN+nr2emZlJhQoVaNeuHW3btmXz5s1eR5NS5k+hbweqF1qOKVj3C2bWEhgDtHPO/RSYeCLirxYtWrBu3TqmTp3KBx98QN26dRk9ejQHDhzwOpqUEn8KfRUQZ2ZXmFl5oAuwqPAOZtYQeJHjZf5N4GOKiD8Kz17v0qULjzzyCLVr12bOnDm6DFMGnLbQnXNHgIHAMmAjMNc5l21mE8ysXcFuU4DzgXlmttbMFv3Ky4lIKahatSqzZs1i5cqVREdH06VLF5o3b67Z62HOvPpfOzEx0WminEjJO3r0KC+99BJjxowhPz+fP/zhD4wfPx6fz+d1NDkDZrbaOZdY1DbdKSoS5iIjI7nvvvtOzF5/+umniY+PJzU1VbPXw4wKXaSMKDx7/corr6RXr17ccMMNrF692utoEiAqdJEy5pprrmHlypWkpqayZcsWrr32Wu69917NXg8DKnSRMigiIoJ77rmH3NxchgwZwssvv0x8fDzPP/+8Zq+HMBW6SBnm8/mYOnUq69ato0GDBvTv35/ExERWrlzpdTQ5Ayp0EaFOnTosX76cOXPm8O233/K73/2OHj16aPZ6iFGhiwhwfPb6nXfeSU5ODqNHj2bOnDnUqlWLqVOncvjwYa/jiR9U6CLyC5UqVWLy5MlkZ2fTtGlThg0bRv369cnMzPQ6mpyGCl1EihQbG8tbb71FWloaP/30E0lJSdxxxx188cUXXkeTX6FCF5FTatu2LdnZ2UycOJG33nqL2rVra/Z6kFKhi8hpVaxYkbFjx5KTk0NycvKJ2euLFy/2OpoUokIXEb/VqFGDefPmnZi9fuutt5KcnMzf//53TXMMAip0ESm2n2evP/HEE6xYsYKbbrqJ2NhYJkyYoOebekiFLiJnpFy5cgwdOpQdO3bwyiuvULNmTR566CGuuOIKWrZsyauvvsqPP/54+heSgFGhi8hZqVSpEt27d2f58uV8/vnnjB8/ni1bttC9e3eqVq1K3759+cc//qFLMqVA89BFJOCOHTvG+++/T2pqKvPmzePHH38kPj6enj170qNHD6pVO9Vz5uVUNA9dREpVREQEN910E6mpqezatYuZM2dyySWXMHr0aGrUqEGbNm2YM2eOvvoYYCp0ESlRlStXplevXrz//vts3ryZ0aNH8+mnn9KlSxcuvfRS+vfvz6pVq3RJJgB0yUVESt3Ro0d59913SUlJYcGCBRw8eJA6derQs2dPunXrRtWqVb2OGLR0yUVEgkpkZCQtW7bktddeY9euXbz44otUrlyZESNGEBMTw6233sqCBQs4dOiQ11FDigpdRDzl8/no168fH374IRs3bmT48OGsXr2ajh07ctlllzF48GDWrl3rdcyQoEsuIhJ0jhw5QkZGBikpKbz55pscOnSI+vXr06tXL+666y6io6O9jugZXXIRkZByzjnncPPNNzN37lx27tzJtGnTKFeuHEOGDKFatWrcfvvtpKWlaU77SXSGLiIhY/369aSmpvLqq6/yzTffcMkll9CtWzd69uxJ3bp1vY5XKk51hq5CF5GQc/jwYdLT00lNTSUtLY0jR46QmJhIz5496dq1KxdeeKHXEUuMLrmISFgpV64c7dq1Y8GCBezYsYMnn3ySQ4cOMXDgQC699FI6d+5Meno6R48e9TpqqfKr0M2sjZltMrM8MxtZxPZmZrbGzI6YWafAxxQRKVp0dDRDhgxh3bp1rFmzhvvuu4/ly5dzyy23UKNGDUaOHMmmTZu8jlkqTlvoZhYJTAduBhKArmaWcNJuXwA9gb8GOqCIiL8aNmzI008/zfbt25k/fz7XXHMNjz/+OLVr16ZJkybMmDGD/Px8r2OWGH/O0BsBec65Lc65Q8BsoH3hHZxzW51znwDHSiCjiEixVKhQ4cQ3Yb766iumTJlCfn4+9957L1WrVuXuu+8mMzOTY8fCq7L8KfRqwJeFlr8qWFdsZtbPzLLMLGv37t1n8hIiIsVStWpVhg8fzoYNG/joo4/o1asXS5YsISkpiZo1azJu3Dg+++wzr2MGRKl+KOqcm+GcS3TOJZblGwNEpPSZGddeey3PPfccO3fuZPbs2SQkJDB58mRiY2Np1qwZKSkp7N+/3+uoZ8yfQt8OVC+0HFOwTkQkJFWsWJHOnTuzdOlSvvjiC/74xz/y9ddf07t3b6pWrUrPnj35+9//HnKXZPwp9FVAnJldYWblgS7AopKNJSJSOmJiYhg1ahQ5OTmsXLmSrl27smDBgl88J3Xr1q1ex/TLaQvdOXcEGAgsAzYCc51z2WY2wczaAZjZtWb2FXAH8KKZZZdkaBGRQDMzmjRpwksvvcSuXbv4y1/+wpVXXnniOaktWrQI+uek6k5REZFT2LZtG7NmzSI1NZXPP/+cypUr07lzZ3r27EmTJk0ws1LNo1v/RUTO0rFjx/jggw9OPCf1wIEDJ56T2r17d2JiYkolh279FxE5SxEREdx4442kpKSwa9cuUlJSqFq1KqNHj+byyy+nTZs2zJ4929PnpKrQRUSK6fzzzz/xTZi8vDzGjBnDxo0b6dq164nnpH700Uel/pxUXXIREQmAY8eOnXhO6vz58zl48CAJCQknLskE6jmpuuQiIlLCIiIiTnwTZteuXcyYMQOfz8cDDzxw4jmp8+fPL9HnpKrQRUQCzOfz0bdvX/7xj3+Qk5PDiBEjWLNmDZ06deKyyy7j9ddfL5HjqtBFREpQrVq1eOSRR/jiiy9IT0+nZcuW1KhRo0SOpWvoIiIhRNfQRUTKABW6iEiYUKGLiIQJFbqISJhQoYuIhAkVuohImFChi4iECRW6iEiY8OzGIjPbDWw7w1+/GPg2gHHCnd6v4tH7VXx6z4rnbN6vy51z0UVt8KzQz4aZZf3anVLy3/R+FY/er+LTe1Y8JfV+6ZKLiEiYUKGLiISJUC30GV4HCDF6v4pH71fx6T0rnhJ5v0LyGrqIiPy3UD1DFxGRk6jQRUTCRMgVupm1MbNNZpZnZiO9zhPMzGymmX1jZhu8zhIKzKy6mb1rZp+aWbaZDfY6UzAzs4pm9pGZrSt4vx72OlMoMLNIM/vYzBYH+rVDqtDNLBKYDtwMJABdzSzB21RBLRVo43WIEHIEGOacSwAaAwP09+uUfgKaO+fqAw2ANmbW2NtIIWEwsLEkXjikCh1oBOQ557Y45w4Bs4H2HmcKWs6594Hvvc4RKpxzO51zawr+vI/j/+iqeZsqeLnj9hcsliv40bcsTsHMYoBk4M8l8fqhVujVgC8LLX+F/sFJCTCzmkBD4F8eRwlqBZcP1gLfABnOOb1fp/YU8ABwrCRePNQKXaTEmdn5wHxgiHNur9d5gplz7qhzrgEQAzQys7oeRwpaZtYW+MY5t7qkjhFqhb4dqF5oOaZgnUhAmFk5jpf5a865BV7nCRXOuT3Au+gzm1O5AWhnZls5frm4uZm9GsgDhFqhrwLizOwKMysPdAEWeZxJwoSZGfAysNE5N9XrPMHOzKLNLKrgz+cCSUCOp6GCmHNulHMuxjlXk+Pd9Y5zrlsgjxFShe6cOwIMBJZx/AOruc65bG9TBS8zex34EKhlZl+ZWR+vMwW5G4DuHD9zWlvwc4vXoYLYpcC7ZvYJx0+2MpxzAf8qnvhPt/6LiISJkDpDFxGRX6dCFxEJEyp0EZEwoUIXEQkTKnQRkTChQhcRCRMqdBGRMPH/AWXpToWAMEluAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "components = range(n_components+1)\n",
    "\n",
    "plt.plot(components, RMSE, c='k')\n",
    "plt.xticks(components)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rmse</th>\n",
       "      <th>rmse_shifted</th>\n",
       "      <th>gradient</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>n_components</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.400775</td>\n",
       "      <td>0.633339</td>\n",
       "      <td>-0.232563</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.278317</td>\n",
       "      <td>0.400775</td>\n",
       "      <td>-0.122459</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.149053</td>\n",
       "      <td>0.278317</td>\n",
       "      <td>-0.129264</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.066768</td>\n",
       "      <td>0.149053</td>\n",
       "      <td>-0.082285</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  rmse  rmse_shifted  gradient\n",
       "n_components                                  \n",
       "1             0.400775      0.633339 -0.232563\n",
       "2             0.278317      0.400775 -0.122459\n",
       "3             0.149053      0.278317 -0.129264\n",
       "4             0.066768      0.149053 -0.082285"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df['rmse_shifted']=df['rmse'].shift(1)\n",
    "df['gradient']=df['rmse']-df['rmse_shifted']\n",
    "\n",
    "df.dropna(inplace=True)\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>rmse</th>\n",
       "      <th>rmse_shifted</th>\n",
       "      <th>gradient</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>n_components</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.400775</td>\n",
       "      <td>0.633339</td>\n",
       "      <td>-0.232563</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  rmse  rmse_shifted  gradient\n",
       "n_components                                  \n",
       "1             0.400775      0.633339 -0.232563"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df[abs(df['gradient'])==max(abs(df['gradient']))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}