demo/ridge_regression_demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ridge Regression\n",
    "\n",
    "*Linear Regression using least square loss with $L_2$ regularization penalty*\n",
    "\n",
    "---\n",
    "* [Theory and Derivation](../theory/ridge_regression.ipynb)\n",
    "* [Implementation in Python](../pymlalgo/regression/ridge_regression.py)\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.core.interactiveshell import InteractiveShell\n",
    "InteractiveShell.ast_node_interactivity = \"all\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Tell python to look for modules in the current project directory\n",
    "\n",
    "import os\n",
    "import sys\n",
    "\n",
    "module_path = os.path.abspath(os.path.join('..'))\n",
    "\n",
    "if module_path not in sys.path:\n",
    "    sys.path.append(module_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The demo uses the [Diamonds dataset](https://www.kaggle.com/shivam2503/diamonds/).\n",
    "\n",
    "Price will be used as the response variable other columns would be used as predictors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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>carat</th>\n",
       "      <th>cut</th>\n",
       "      <th>color</th>\n",
       "      <th>clarity</th>\n",
       "      <th>depth</th>\n",
       "      <th>table</th>\n",
       "      <th>price</th>\n",
       "      <th>x</th>\n",
       "      <th>y</th>\n",
       "      <th>z</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.23</td>\n",
       "      <td>Ideal</td>\n",
       "      <td>E</td>\n",
       "      <td>SI2</td>\n",
       "      <td>61.5</td>\n",
       "      <td>55.0</td>\n",
       "      <td>326</td>\n",
       "      <td>3.95</td>\n",
       "      <td>3.98</td>\n",
       "      <td>2.43</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.21</td>\n",
       "      <td>Premium</td>\n",
       "      <td>E</td>\n",
       "      <td>SI1</td>\n",
       "      <td>59.8</td>\n",
       "      <td>61.0</td>\n",
       "      <td>326</td>\n",
       "      <td>3.89</td>\n",
       "      <td>3.84</td>\n",
       "      <td>2.31</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.23</td>\n",
       "      <td>Good</td>\n",
       "      <td>E</td>\n",
       "      <td>VS1</td>\n",
       "      <td>56.9</td>\n",
       "      <td>65.0</td>\n",
       "      <td>327</td>\n",
       "      <td>4.05</td>\n",
       "      <td>4.07</td>\n",
       "      <td>2.31</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.29</td>\n",
       "      <td>Premium</td>\n",
       "      <td>I</td>\n",
       "      <td>VS2</td>\n",
       "      <td>62.4</td>\n",
       "      <td>58.0</td>\n",
       "      <td>334</td>\n",
       "      <td>4.20</td>\n",
       "      <td>4.23</td>\n",
       "      <td>2.63</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.31</td>\n",
       "      <td>Good</td>\n",
       "      <td>J</td>\n",
       "      <td>SI2</td>\n",
       "      <td>63.3</td>\n",
       "      <td>58.0</td>\n",
       "      <td>335</td>\n",
       "      <td>4.34</td>\n",
       "      <td>4.35</td>\n",
       "      <td>2.75</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   carat      cut color clarity  depth  table  price     x     y     z\n",
       "0   0.23    Ideal     E     SI2   61.5   55.0    326  3.95  3.98  2.43\n",
       "1   0.21  Premium     E     SI1   59.8   61.0    326  3.89  3.84  2.31\n",
       "2   0.23     Good     E     VS1   56.9   65.0    327  4.05  4.07  2.31\n",
       "3   0.29  Premium     I     VS2   62.4   58.0    334  4.20  4.23  2.63\n",
       "4   0.31     Good     J     SI2   63.3   58.0    335  4.34  4.35  2.75"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diamonds = pd.read_csv('https://raw.githubusercontent.com/tidyverse/ggplot2/master/data-raw/diamonds.csv')\n",
    "diamonds.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Convert categorical values to category codes**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "diamonds['cut'] = diamonds['cut'].astype('category').cat.codes\n",
    "diamonds['color'] = diamonds['color'].astype('category').cat.codes\n",
    "diamonds['clarity'] = diamonds['clarity'].astype('category').cat.codes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Split into predictors and response**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "predictors = diamonds.drop(columns='price')\n",
    "response = diamonds['price']\n",
    "\n",
    "# Recall that the algorithm needs predictors in shape, (n, d)\n",
    "x = np.asarray(predictors)\n",
    "y = np.asarray(response).reshape(-1, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Train, test split**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((44135, 9), (10788, 9), (44135, 1), (10788, 1))"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from pymlalgo.util.model_selection import train_test_split\n",
    "\n",
    "x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2)\n",
    "x_train.shape, x_test.shape, y_train.shape, y_test.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Normalize the Data**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymlalgo.util.standardization import Standardizer\n",
    "\n",
    "x_standardizer = Standardizer(x_train)\n",
    "x_train = x_standardizer.standardize(x_train)\n",
    "x_test = x_standardizer.standardize(x_test)\n",
    "\n",
    "y_standardizer = Standardizer(y_train)\n",
    "y_train = y_standardizer.standardize(y_train)\n",
    "y_test = y_standardizer.standardize(y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Train and test using `RigdeRegression`**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymlalgo.regression.ridge_regression import RidgeRegression\n",
    "\n",
    "lambd=1e-5\n",
    "\n",
    "ridge_model = RidgeRegression(x_train, y_train, lambd=lambd, min_grad=1e-4, max_iter=1500)\n",
    "ridge_model.train()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Cost history vs number of iterations**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEWCAYAAACXGLsWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XuYXXV97/H3Z+4hCYTLiIEEAopo9KGgk1QtIkcoBh9PaE8ViBekVemRh2OV4wXEBy32OVo59WArKtRaRG4CilKNjRTB2pZLhlsgYDSGSIYAGcSQhNxmMt/zx/pts7Kz9157yKzsyeTzetjP7HX/zgp7f+b3++29liICMzOzRtpaXYCZmY1/DgszMyvksDAzs0IOCzMzK+SwMDOzQg4LMzMr5LCwCU/SmyQta3UdE5XP797BYWFjRtK7JPVL2iDpKUk/lnT8Lu5zpaSTGyw/UdJAjfl3SvoAQET8PCKObuJYn5V0za7Uu7tJmiUpJHWk6ask/U3JxwxJL69MN3t+bc/msLAxIel84DLg/wAHA4cBXwVOa2Vdu1PlDXtPNhF+BytJRPjhxy49gP2ADcA7G6zTTRYmq9PjMqA7LTsI+CGwFngO+DnZHzLfBkaATWn/n6ix3xOBgRrz7wQ+UGsd4JPAk8B6YBlwEjAP2AoMpWM9lNY9BLg11bUc+GBuP58FbgauAdYBnwY2Agfm1nkdMAh0VtV3SPq9DsjNOw54FugEXg78DHg+zftOnfM6CwigAzgn1b81/Q7/kjvWd1MdjwMfbvA7fACYC9yV/j2eAr4CdKX1/z0d74V0jDNqnN9XpfO/FlgKzM8tuwq4HPhROv/3AC9r9f/DfjTxOm91AX7s+Y/0RjsMdDRY5xLgbuAlQC/wX8Dn0rLPA19Pb5KdwJsApWUrgZMb7HeHN6rc/DupERbA0cAq4JA0PavyZpXeOK+p2s/PyFpIPcCx6Q33pNz6Q8CfkIXbJGAh8KHc9v8P+Ic6tf+UHcPnUuDr6fn1wEVpvz3A8XX2MSu9eXek6auAv8ktbwPuAy4GuoAjgRXAWxv8Dq8DXk8WQLOAx4CP5PYZwMtr/Rukf7/lwKfS8d6SQuHoXH3PkQVSB3AtcEOr/x/2o/jhbigbCwcCz0bEcIN13g1cEhFrImIQ+GvgvWnZEDAdODwihiLrAx/NRcsOkbQ2/wDqjZVsI2vlzJbUGRErI+LXtVaUNDPt55MRsTkiHgS+kasb4K6I+H5EjETEJuBbwHvS9u3AArIWUi3XpeVIEnBmmgfZOTmcLNQ2R8R/NHMiapgD9EbEJRGxNSJWAP+YjlXzd4iI+yLi7ogYjoiVwBXAm5s83uuBKcAX0vF+StZqXJBb53sRcW/6/+VashC2cc5hYWPht8BBBf3dhwC/yU3/Js2D7C/q5cBPJK2QdMEoj786IqblH0DNN9eIWA58hOwv6jWSbpB0SK11U33PRcT6qroPzU2vqtrmB2RBdCTwx8DzEXFvnf3fDLwhHf8Esr/Yf56WfQIQcK+kpZL+os4+ihxOVZiS/dV/cL3fQdIrJP1Q0tOS1pGNQx3U5PEOAVZFxEhuXvU5ezr3fCNZuNg457CwsXAXsJmsK6Oe1WRvXBWHpXlExPqI+N8RcSTw34HzJZ2U1hvzyyJHxHURcXyqJ4C/rXOs1cABkqZW1f1kfndV+94M3EjWknov9VsVRMRa4CfA6cC7gOsrLaqIeDoiPhgRhwB/CXw1/wmkRr9e1fQq4PGqMJ0aEW9rsM3XgF8AR0XEvmThoiaODdk5mykp/95Sfc5sD+SwsF0WEc+T9YlfLulPJO0jqVPSqZK+mFa7Hvi0pF5JB6X1rwGQ9HZJL09dMevIuoq2pe2eIetnHxOSjpb0FkndZAG3qepYsypvdBGximxs5fOSeiQdA7yfrOukkauBs4H5pN+xgeuAs4A/Y3sXFJLeKWlGmvwd2Rv6tp0330n1+boXWCfpk5ImSWqX9BpJcxrsYyrZv8MGSa8EPlRwjLx7yAa/P5H+HziR7A+AG5qo3cYxh4WNiYj4EnA+2SeCBsn+oj0P+H5a5W+AfmAJ8DBwf5oHcBTwb2SfrrkL+GpE3JmWfZ4sZNZK+tgYlNoNfIHsE0ZPkw24fyotuyn9/K2k+9PzBWSDvKuBW4DPRMRtjQ4QEf9J9imu+1OffyO3kv3+z0TEQ7n5c4B7JG1I6/xVRDxe+NvBP5F1g62V9P2I2Eb2Zn0s2SehniUbd9mvwT4+RtbSWU82vvGdquWfBb6VjnF6fkFEbCULyVPTsb4KnBURv2iidhvHKp84MbMxJOmnwHUR8Y1W12I2FhwWZmMsdfHcBsysGhw322O5G8psDEn6FlmX2kccFDaRuGVhZmaF3LIwM7NCE+aiYQcddFDMmjWr1WWYme1R7rvvvmcjordovQkTFrNmzaK/v7/VZZiZ7VEk/aZ4LXdDmZlZExwWZmZWyGFhZmaFHBZmZlbIYWFmZoUcFmZmVshhYWZmhfb6sFi/eYgv3fZLHly1ttWlmJmNW3t9WAxvC/7+9l/xwBO/a3UpZmbj1l4fFpO7sy+xv7BluMWVmJmNX3t9WHR1tNHV3saGLc3csdLMbO9UalhImidpmaTlki6osfx8SY9KWiLpdkmHVy3fV9KTkr5SZp37dLezcatbFmZm9ZQWFpLagcvJ7sU7G1ggaXbVag8AfRFxDHAz8MWq5Z8DflZWjRX7dLazcatbFmZm9ZTZspgLLI+IFekm7jcAp+VXiIg7ImJjmrwbmFFZJul1wMHAT0qssXIsfA8oM7P6ygyLQ4FVuemBNK+e9wM/BpDUBvwd8PFGB5B0jqR+Sf2Dg4O7VGzgtDAzq6fMsFCNeTXfkSW9B+gDLk2zzgUWRsSqWuv/fmcRV0ZEX0T09fYW3rujfqGqV5mZmUG5Nz8aAGbmpmcAq6tXknQycBHw5ojYkma/AXiTpHOBKUCXpA0RsdMg+ViQnBVmZo2UGRaLgaMkHQE8CZwJvCu/gqTjgCuAeRGxpjI/It6dW+dsskHwUoICQIjwoIWZWV2ldUNFxDBwHrAIeAy4MSKWSrpE0vy02qVkLYebJD0o6day6mnELQszs8ZKvQd3RCwEFlbNuzj3/OQm9nEVcNVY15Yn8KehzMwa2Ou/wQ3po7OtLsLMbBxzWFBpWTguzMzqcVgAeMzCzKwhhwXpCyFOCzOzuhwWVMYsnBZmZvU4LPCnoczMijgsSN+zcFiYmdXlsDAzs0IOC9LlPjxmYWZWl8MCd0OZmRVxWCTOCjOz+hwW+E55ZmZFHBZU7tLktDAzq8dhgccszMyKOCzw/SzMzIo4LPCd8szMijgscMvCzKyIwwJfG8rMrIjDAsB3yjMza8hhge+UZ2ZWxGFBNmZhZmb1lRoWkuZJWiZpuaQLaiw/X9KjkpZIul3S4Wn+sZLukrQ0LTuj1DrxmIWZWSOlhYWkduBy4FRgNrBA0uyq1R4A+iLiGOBm4Itp/kbgrIh4NTAPuEzStBJr9VVnzcwaKLNlMRdYHhErImIrcANwWn6FiLgjIjamybuBGWn+LyPiV+n5amAN0FtWoW5ZmJk1VmZYHAqsyk0PpHn1vB/4cfVMSXOBLuDXY1rdDsdwWJiZNdJR4r5rDRvXfEuW9B6gD3hz1fzpwLeB90XESI3tzgHOATjssMN2oVB3Q5mZNVJmy2IAmJmbngGsrl5J0snARcD8iNiSm78v8CPg0xFxd60DRMSVEdEXEX29vbvQS+WWhZlZQ2WGxWLgKElHSOoCzgRuza8g6TjgCrKgWJOb3wXcAlwdETeVWGN2PHy5DzOzRkoLi4gYBs4DFgGPATdGxFJJl0ian1a7FJgC3CTpQUmVMDkdOAE4O81/UNKxZdXq71mYmTVW5pgFEbEQWFg17+Lc85PrbHcNcE2Zte180N16NDOzPYq/wY0HuM3Mijgs8EdnzcyKOCzw/SzMzIo4LPCd8szMijgscMvCzKyIwyJxw8LMrD6HBZWrzpqZWT0OC9JFrNy0MDOry2GBxyzMzIo4LPD9LMzMijgs8J3yzMyKOCxwy8LMrIjDAl/uw8ysiMMCAH901sysEYcFlZaF48LMrB6HBbVvFm5mZts5LPCYhZlZEYcFvvmRmVkRhwVuWZiZFXFY4Mt9mJkVcViQdUOZmVl9DovEH501M6uv1LCQNE/SMknLJV1QY/n5kh6VtETS7ZIOzy17n6Rfpcf7yqwTd0OZmTVUWlhIagcuB04FZgMLJM2uWu0BoC8ijgFuBr6Ytj0A+Azwh8Bc4DOS9i+tVnBamJk1UGbLYi6wPCJWRMRW4AbgtPwKEXFHRGxMk3cDM9LztwK3RcRzEfE74DZgXlmF+k55ZmaNlRkWhwKrctMDaV497wd+PJptJZ0jqV9S/+Dg4IsuNLvqrOPCzKyeMsOi1keMar4jS3oP0AdcOpptI+LKiOiLiL7e3t4XX6jHLMzMGiozLAaAmbnpGcDq6pUknQxcBMyPiC2j2Xas+H4WZmaNlRkWi4GjJB0hqQs4E7g1v4Kk44AryIJiTW7RIuAUSfunge1T0rxS+E55ZmaNdZS144gYlnQe2Zt8O/DNiFgq6RKgPyJuJet2mgLcJAngiYiYHxHPSfocWeAAXBIRz5VVq1sWZmaNlRYWABGxEFhYNe/i3POTG2z7TeCb5VWX42tDmZk15G9w48t9mJkVcVjgO+WZmRVxWJDGLFpdhJnZOOawwPezMDMr4rDAd8ozMyvisMAtCzOzIg4LfLkPM7MiDgsA5JaFmVkDDguyloXbFmZm9Tks8OU+zMyKOCzwmIWZWRGHhZmZFXJYkL5n4X4oM7O6HBa4G8rMrIjDAg9wm5kVcViQ7pTntDAzq8thkTgqzMzqc1iQvpTntDAzq6upsJD07Wbm7amyq86amVk9zbYsXp2fkNQOvG7sy2kN3ynPzKyxhmEh6UJJ64FjJK1Lj/XAGuAHu6XC3cC9UGZmjTUMi4j4fERMBS6NiH3TY2pEHBgRFxbtXNI8ScskLZd0QY3lJ0i6X9KwpHdULfuipKWSHpP091J2ub8y+H4WZmaNNdsN9UNJkwEkvUfSlyQd3miD1FV1OXAqMBtYIGl21WpPAGcD11Vt+0bgj4BjgNcAc4A3N1nrqEm+U56ZWSPNhsXXgI2S/gD4BPAb4OqCbeYCyyNiRURsBW4ATsuvEBErI2IJMFK1bQA9QBfQDXQCzzRZ66j5S3lmZo01GxbDkY0AnwZ8OSK+DEwt2OZQYFVueiDNKxQRdwF3AE+lx6KIeKzJWkfPl/swM2uo2bBYL+lC4L3Aj1IXU2fBNrXGGJp6T5b0cuBVwAyygHmLpBNqrHeOpH5J/YODg83suk6hTgszs0aaDYszgC3AX0TE02Rv4JcWbDMAzMxNzwBWN3m8PwXujogNEbEB+DHw+uqVIuLKiOiLiL7e3t4md72z7EKCTgszs3qaCosUENcC+0l6O7A5IorGLBYDR0k6QlIXcCZwa5N1PQG8WVKHpE6ywe3SuqE8ZmFm1liz3+A+HbgXeCdwOnBP9Uddq0XEMHAesIjsjf7GiFgq6RJJ89N+50gaSPu9QtLStPnNwK+Bh4GHgIci4l9G/ds1yZcoNzNrrKPJ9S4C5kTEGgBJvcC/kb2p1xURC4GFVfMuzj1fTNY9Vb3dNuAvm6xtl/nmR2ZmjTU7ZtFWCYrkt6PYdtxzy8LMrLFmWxb/KmkRcH2aPoOqFsOezGMWZmaNNQyL9BHWgyPi45L+B3A82XvrXWQD3hNDeVcSMTObEIq6ki4D1gNExPci4vyI+ChZq+KysoszM7PxoSgsZqXLcewgIvqBWaVU1ALrNw8BsGHLcIsrMTMbn4rCoqfBskljWUgr/fN/rgTg+nueaG0hZmbjVFFYLJb0weqZkt4P3FdOSa3T0e6xCzOzWoo+DfUR4BZJ72Z7OPSRXQ32T8ssrBU62ifMp4HNzMZUw7CIiGeAN0r6b2T3lQD4UUT8tPTKWqCzzS0LM7NamvqeRUTcQXbJ8Amt0y0LM7Oa/O4InPTKlwAwubu9xZWYmY1PDgvg/FNeAWS3VzUzs505LID2NFaxbcTX/DAzq8VhAbTLYWFm1ojDgu0tixFfTdDMrCaHBdvDYnibw8LMrBaHBdBW6YZyy8LMrCaHBbluKI9ZmJnV5LAAOirdUA4LM7OaHBZAmwe4zcwacliwvWWxdXikxZWYmY1PDgtg355OJnW28+TaTa0uxcxsXCo1LCTNk7RM0nJJF9RYfoKk+yUNS3pH1bLDJP1E0mOSHpU0q6w629rE9Gk9rFm3paxDmJnt0UoLC0ntwOXAqcBsYIGk2VWrPQGcDVxXYxdXA5dGxKuAucCasmoFmNrd4duqmpnV0dQlyl+kucDyiFgBIOkG4DTg0coKEbEyLdthsCCFSkdE3JbW21BinQBM6XFYmJnVU2Y31KHAqtz0QJrXjFcAayV9T9IDki5NLZUdSDpHUr+k/sHBwV0qdnJXBxs2OyzMzGopMyxqXe+72c+mdgBvAj4GzAGOJOuu2nFnEVdGRF9E9PX29r7YOgHo7mxny/C2XdqHmdlEVWZYDAAzc9MzgNWj2PaBiFgREcPA94HXjnF9O+hsF0O+NpSZWU1lhsVi4ChJR0jqAs4Ebh3FtvtLqjQX3kJurKMMnW1tDI/4exZmZrWUFhapRXAesAh4DLgxIpZKukTSfABJcyQNAO8ErpC0NG27jawL6nZJD5N1af1jWbUCdLTLV501M6ujzE9DERELgYVV8y7OPV9M1j1Va9vbgGPKrC+vs72NoW1uWZiZ1eJvcCed7fKFBM3M6nBYJB3tbe6GMjOrw2GRdLaJoZERwleeNTPbicMi6WhvIwK2uSvKzGwnDouksz07FR63MDPbmcMimdSZnYpNW/0tbjOzag6LZEpPJwDrfX0oM7OdOCySqT3ZV07WbR5qcSVmZuOPwyKZ2p2FhS9Tbma2M4dF0tOVXQF905DHLMzMqjkskp6OLCy2OCzMzHbisEh60qehNg/5+lBmZtUcFsmk1A212S0LM7OdOCySSjeUw8LMbGcOi6SnszLA7W4oM7NqDouku6MyZuGWhZlZNYdF0tYmujra2DzssDAzq+awyJnU2c4Wd0OZme3EYZHT09nmCwmamdXgsMjp6Wx3N5SZWQ0Oi5yejnYPcJuZ1VBqWEiaJ2mZpOWSLqix/ARJ90salvSOGsv3lfSkpK+UWWdFT2ebv8FtZlZDaWEhqR24HDgVmA0skDS7arUngLOB6+rs5nPAz8qqsVp3p1sWZma1lNmymAssj4gVEbEVuAE4Lb9CRKyMiCXATn/OS3odcDDwkxJr3MEkh4WZWU1lhsWhwKrc9ECaV0hSG/B3wMdLqKsud0OZmdVWZlioxrxocttzgYURsarRSpLOkdQvqX9wcHDUBVbzp6HMzGrrKHHfA8DM3PQMYHWT274BeJOkc4EpQJekDRGxwyB5RFwJXAnQ19fXbBDV5U9DmZnVVmZYLAaOknQE8CRwJvCuZjaMiHdXnks6G+irDooyuBvKzKy20rqhImIYOA9YBDwG3BgRSyVdImk+gKQ5kgaAdwJXSFpaVj3N6Ols921VzcxqKLNlQUQsBBZWzbs493wxWfdUo31cBVxVQnk76elsZ+vwCCMjQVtbrSEXM7O9k7/BnVO5p8WWYXdFmZnlOSxyJqX7cG/cOtziSszMxheHRc7k7qxX7oUtHrcwM8tzWORMqYSFWxZmZjtwWORsb1k4LMzM8hwWOZO7swHuDQ4LM7MdOCxyPGZhZlabwyJncpe7oczManFY5FQGuN0NZWa2I4dFjge4zcxqc1jkdHW00dXexgZ/dNbMbAcOiyqTu9vdsjAzq+KwqDK5u4ON/jSUmdkOHBZVJnd1eIDbzKyKw6LK5O52X+7DzKyKw6LK5O4ONrgbysxsBw6LKlO6OzzAbWZWxWFRZUp3B+s3D7W6DDOzccVhUWXaPp08v8lhYWaW57CoMm2fLjYPjbB5yOMWZmYVDosq++/TBcDajW5dmJlVOCyqTNunE4C1m7a2uBIzs/Gj1LCQNE/SMknLJV1QY/kJku6XNCzpHbn5x0q6S9JSSUsknVFmnXmVsPjdC25ZmJlVlBYWktqBy4FTgdnAAkmzq1Z7AjgbuK5q/kbgrIh4NTAPuEzStLJqzeud0g3A4IYtu+NwZmZ7hI4S9z0XWB4RKwAk3QCcBjxaWSEiVqZlI/kNI+KXueerJa0BeoG1JdYLwEv36wHgqbWbyj6Umdkeo8xuqEOBVbnpgTRvVCTNBbqAX9dYdo6kfkn9g4ODL7rQvKk9nUzp7uCp5zePyf7MzCaCMsNCNebFqHYgTQe+Dfx5RIxUL4+IKyOiLyL6ent7X2SZO5u+Xw8Dv3PLwsysosywGABm5qZnAKub3VjSvsCPgE9HxN1jXFtDrz5kX5YMrCViVNlmZjZhlRkWi4GjJB0hqQs4E7i1mQ3T+rcAV0fETSXWWNNrD9+fNeu3sNpdUWZmQIlhERHDwHnAIuAx4MaIWCrpEknzASTNkTQAvBO4QtLStPnpwAnA2ZIeTI9jy6q12pxZBwDwz//xOEPbdur9MjPb62iidLX09fVFf3//mO3vo995kFseeJKX9U7mmg/8IdP3mzRm+zYzGy8k3RcRfUXr+RvcdXzp9D/gH8/q45l1W/jg1f1s2uprRZnZ3sthUYck/nj2wXz5zGNZunodn/zuEiKCNes2c9evf8svnl7Hsxu2sHV4xAPhZjbhlfmlvAnhpFcdzMdOOZpLFy3jJ48+zeahnccw2gSTOtvZp7uDSZ3tbBsJhraNsG0kkABEm0CCNom2bCZK84TSz2JSM2s1t69mVmruaM3XZWZj71XT9+UfFhxX6jEcFk0498SXcdCULh55ch2HH7gPR790Ks9vGuK5F7aybtMQm4dG2Lh1G5uGhtm0dRvtbW10tou2tuwNNCKIgAgYiSDIfpL9ly1voo5mGzDN7at4rabbS25YmbXUzP3LH1N1WDRBEmfMOYwz5rS6EjOz1vCYhZmZFXJYmJlZIYeFmZkVcliYmVkhh4WZmRVyWJiZWSGHhZmZFXJYmJlZoQlz1VlJg8BvdmEXBwHPjlE5ZRjv9YFrHAvjvT5wjWNhPNV3eEQU3mp0woTFrpLU38xleltlvNcHrnEsjPf6wDWOhfFeXy3uhjIzs0IOCzMzK+Sw2O7KVhdQYLzXB65xLIz3+sA1joXxXt9OPGZhZmaF3LIwM7NCDgszMyu014eFpHmSlklaLumCFtYxU9Idkh6TtFTSX6X5B0i6TdKv0s/903xJ+vtU9xJJr91NdbZLekDSD9P0EZLuSfV9R1JXmt+dppen5bN2U33TJN0s6RfpXL5hHJ7Dj6Z/40ckXS+pp9XnUdI3Ja2R9Ehu3qjPm6T3pfV/Jel9Jdd3afp3XiLpFknTcssuTPUtk/TW3PzSXu+1aswt+5ikkHRQmt7t53CXZbf83DsfQDvwa+BIoAt4CJjdolqmA69Nz6cCvwRmA18ELkjzLwD+Nj1/G/Bjsttkvx64ZzfVeT5wHfDDNH0jcGZ6/nXgQ+n5ucDX0/Mzge/spvq+BXwgPe8Cpo2ncwgcCjwOTMqdv7NbfR6BE4DXAo/k5o3qvAEHACvSz/3T8/1LrO8UoCM9/9tcfbPTa7kbOCK9xtvLfr3XqjHNnwksIvvS8EGtOoe7/Pu1uoCW/vLwBmBRbvpC4MJW15Vq+QHwx8AyYHqaNx1Ylp5fASzIrf/79UqsaQZwO/AW4Ifpf/Rncy/Y35/P9OJ4Q3rekdZTyfXtm96IVTV/PJ3DQ4FV6c2gI53Ht46H8wjMqnozHtV5AxYAV+Tm77DeWNdXtexPgWvT8x1ex5VzuDte77VqBG4G/gBYyfawaMk53JXH3t4NVXnhVgykeS2VuhqOA+4BDo6IpwDSz5ek1VpR+2XAJ4CRNH0gsDYihmvU8Pv60vLn0/plOhIYBP45dZV9Q9JkxtE5jIgngf8LPAE8RXZe7mN8nceK0Z63Vr6e/oLsL3Ua1LHb65M0H3gyIh6qWjRuamzW3h4WqjGvpZ8lljQF+C7wkYhY12jVGvNKq13S24E1EXFfkzW04tx2kHUDfC0ijgNeIOs+qWe315j6/U8j6x45BJgMnNqgjnH3/yj1a2pJrZIuAoaBayuz6tSxu18z+wAXARfXWlynlvH47w04LAbI+hMrZgCrW1QLkjrJguLaiPhemv2MpOlp+XRgTZq/u2v/I2C+pJXADWRdUZcB0yR11Kjh9/Wl5fsBz5VYX+WYAxFxT5q+mSw8xss5BDgZeDwiBiNiCPge8EbG13msGO152+3nMw0Avx14d6R+m3FU38vI/ih4KL1uZgD3S3rpOKqxaXt7WCwGjkqfROkiG0C8tRWFSBLwT8BjEfGl3KJbgconIt5HNpZRmX9W+lTF64HnK10GZYiICyNiRkTMIjtPP42IdwN3AO+oU1+l7nek9Uv9CykingZWSTo6zToJeJRxcg6TJ4DXS9on/ZtXahw35zFntOdtEXCKpP1TC+qUNK8UkuYBnwTmR8TGqrrPTJ8kOwI4CriX3fx6j4iHI+IlETErvW4GyD7E8jTj5ByOSqsHTVr9IPtUwi/JPiVxUQvrOJ6subkEeDA93kbWP3078Kv084C0voDLU90PA327sdYT2f5pqCPJXojLgZuA7jS/J00vT8uP3E21HQv0p/P4fbJPlIyrcwj8NfAL4BHg22Sf2mnpeQSuJxtDGSJ7U3v/izlvZGMHy9Pjz0uubzlZ/37l9fL13PoXpfqWAafm5pf2eq9VY9XylWwf4N7t53BXH77ch5mZFdrbu6HMzKwJDgszMyvksDAzs0IOCzMzK+SwMDOzQg4Lm3AkbUg/Z0l61xjv+1NV0/81hvv+WLqK6iOSHpJ01ljtO+1/mqRzx3KftvdwWNhENgsYVVhIai9YZYewiIg3jrKmesf9n2QXjpwbEa8hu4JprUs/7IppZFexNRs1h4VNZF8A3iTpQWX3kGhP90BYnO4h8JcAkk5Udi+R68i+IIWk70u6T9l9J85J874ATEr7uzbNq7RilPb9iKSHJZ2R2/ed2n6E1knwAAACW0lEQVSPjWvTN7erfQo4N9L1wCLi+Yj4VtrHSenCiA8ru2dCd5q/Utvvj9An6c70/LNpvTslrZD04dz5eFmq/9KxPtk2wbX6W4F++DHWD2BD+nki6Zvmafoc4NPpeTfZN72PSOu9AByRW7fybeVJZN+0PjC/7xrH+jPgNrJ7JhxMdlmP6Wnfz5Nd46cNuAs4vmofU4Hf1fldesi+pfyKNH012UUmYcdvBPcBd6bnnwX+K/2OBwG/BTppcIlvP/woerhlYXuTU8iux/Mg2eXfDyS7bhDAvRHxeG7dD0t6CLib7MJuR9HY8cD1EbEtIp4BfgbMye17ICJGyC5LMatqW1H/yqJHk1148Jdp+ltkXVRFfhQRWyLiWbILAB7cxDZmdXUUr2I2YQj4XxGxw4XZJJ1I1rLIT59MdtOhjal7p6eJfdezJfd8G1Wvu4hYJ+kFSUdGxIpR7HeY7V3J1fU1PKbZaLllYRPZerIunopFwIfSpeCR9AplN0eqth9Zt9BGSa8ku+1lxVBl+yr/DpyRxkV6yf76v3cUtX4euFzSvqm2fdNYyS+AWZJentZ7L1mrBbJuqNel53/WxDGqz4dZ0xwWNpEtAYbTx1A/CnyD7HLg90t6hOyWlbX+4v5XoEPSEuBzZF1RFVcCSyoD3Dm3pOM9BPwU+ERkl6Ju1tfILlO+ONX2M2BjRGwG/hy4SdLDZHcp/Hra5q+BL0v6OVnroaGI+C3wn2kQ3gPcNiq+6qyZmRVyy8LMzAo5LMzMrJDDwszMCjkszMyskMPCzMwKOSzMzKyQw8LMzAr9f8/RPZRomSVxAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(ridge_model.cost_history[1:])\n",
    "plt.xlabel('Iteration Count')\n",
    "plt.ylabel('Cost')\n",
    "plt.title('Cost History vs Iteration');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "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>train_score</th>\n",
       "      <th>test_score</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.885868</td>\n",
       "      <td>0.881184</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   train_score  test_score\n",
       "0     0.885868    0.881184"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame({'train_score' : [ridge_model.r_squared(x_train, y_train)], \n",
    "              'test_score': [ridge_model.r_squared(x_test, y_test)]})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Comparison with `sklearn`**  \n",
    "\n",
    "The cost function in `sklearn` is $$ F(\\beta) = ||Y - X^T\\beta||_2^2 + \\alpha ||\\beta||_2^2$$\n",
    "To make it equivalent to the cost function in current implementation, $$\\lambda = \\frac{\\alpha}{n} \\implies \\alpha = n\\lambda$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "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>coef</th>\n",
       "      <th>sk_coef</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.313943</td>\n",
       "      <td>1.314920</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.019259</td>\n",
       "      <td>0.019255</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-0.115235</td>\n",
       "      <td>-0.115236</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.125249</td>\n",
       "      <td>0.125221</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-0.055442</td>\n",
       "      <td>-0.055505</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>-0.052765</td>\n",
       "      <td>-0.052761</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>-0.332819</td>\n",
       "      <td>-0.334140</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0.007871</td>\n",
       "      <td>0.008149</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>-0.004492</td>\n",
       "      <td>-0.004401</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       coef   sk_coef\n",
       "0  1.313943  1.314920\n",
       "1  0.019259  0.019255\n",
       "2 -0.115235 -0.115236\n",
       "3  0.125249  0.125221\n",
       "4 -0.055442 -0.055505\n",
       "5 -0.052765 -0.052761\n",
       "6 -0.332819 -0.334140\n",
       "7  0.007871  0.008149\n",
       "8 -0.004492 -0.004401"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.linear_model import Ridge\n",
    "\n",
    "alpha = x_train.shape[1] * lambd\n",
    "ridge_sk = Ridge(alpha = alpha)\n",
    "ridge_sk = ridge_sk.fit(x_train, y_train)\n",
    "\n",
    "pd.DataFrame({'coef': ridge_model.w.flatten(), 'sk_coef' : ridge_sk.coef_.flatten()})"
   ]
  },
  {
   "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>train_score</th>\n",
       "      <th>train_score_sk</th>\n",
       "      <th>test_score</th>\n",
       "      <th>test_score_sk</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.885868</td>\n",
       "      <td>0.885868</td>\n",
       "      <td>0.881184</td>\n",
       "      <td>0.881181</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   train_score  train_score_sk  test_score  test_score_sk\n",
       "0     0.885868        0.885868    0.881184       0.881181"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame({'train_score' : [ridge_model.r_squared(x_train, y_train)], \n",
    "              'train_score_sk' : [ridge_sk.score(x_train, y_train)],\n",
    "              'test_score': [ridge_model.r_squared(x_test, y_test)],\n",
    "              'test_score_sk' : [ridge_sk.score(x_test, y_test)]})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "$\\dagger$ For this dataset, $n >> d$ and hence, $\\lambda$ can be set to $0$"
   ]
  }
 ],
 "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}