chapter-03/3.4.1-ridge-regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3.4.1 Ridge Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ridge regression shrinks the regression coefficients by imposing a penalty on their size. The ridge coefficients minimize a penalized RSS (3.41):\n",
    "$$\n",
    "\\hat{\\beta}^{ridge} = \\underset{\\beta}{argmin} \\left\\{\n",
    "\\sum_{i=1}^N \\left( y_i-\\beta_0-\\sum_{j=1}^p{x_{ij}\\beta_j} \\right)^2\n",
    "+\\gamma\\sum_{j=1}^p{\\beta_j^2}\n",
    "\\right\\}\n",
    "$$\n",
    "Here $\\gamma \\ge 0$ is a complexity parameter: $\\gamma \\rightarrow \\infty$, the coefficients are shrunk toward zero (and each other). This idea also used in neural networks (known as *weight decay*)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An equivalent way to write the ridge problem is (3.42):\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\hat{\\beta}^{ridge} = \\underset{\\beta}{argmin} \\left\\{\n",
    "\\sum_{i=1}^N \\left( y_i-\\beta_0-\\sum_{j=1}^p{x_{ij}\\beta_j} \\right)^2\n",
    "\\right\\},\\\\\n",
    "\\text{subject to } \\sum_{j=1}^p {\\beta_j^2 \\le t}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "There is a one-to-one correspondence between $\\gamma \\text{ and } t$. A large positive coefficient on one variable can be canceled by a similarly large negative coefficient on its correlated cousin. By imposing constraints as in (3.42), this problem is alleviated."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ridge solutions are not equivariant under scaling of the inputs, so one normally standardizes the inputs before solving (3.41). \n",
    "\n",
    "Notice that the intercepts has been left out of the penalty term: because, adding a constant to each target $y_i$ would not simply result in a shift of the prediction by the same amount. It can be shown that the solution to (3.41) can be separated into two parts, after reparametrization using *centered* inputs, i.e replacing $x_{ij} \\text{ by } x_{ij}-\\overline{x}_j$:\n",
    "1. we estimate $\\beta_0 \\text{ as } \\overline{y}$,\n",
    "2. The remaining coefficients get estimated by a ridge regression.\n",
    "\n",
    "*Proof*:\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\hat{\\beta}^{ridge} = \\underset{\\beta}{argmin} \\left\\{\n",
    "\\sum_{i=1}^N \\left( y_i- \\beta_0-\\sum_{j=1}^p\\overline{x}_j\\beta_j -\\sum_{j=1}^p{(x_{ij} - \\overline{x}_j)\\beta_j} \\right)^2\n",
    "+\\gamma\\sum_{j=1}^p{\\beta_j^2}\n",
    "\\right\\}\\\\\n",
    "\\hat{\\beta}^{ridge} = \\underset{\\beta}{argmin} \\left\\{\n",
    "\\sum_{i=1}^N \\left( y_i- \\beta_0^C  -\\sum_{j=1}^p{(x_{ij} - \\overline{x}_j)\\beta_j} \\right)^2\n",
    "+\\gamma\\sum_{j=1}^p{\\beta_j^2}\n",
    "\\right\\}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "where $\\beta_0^C = \\beta_0-\\sum_{j=1}^p\\overline{x}_j\\beta_j$. Now the stationary point w.r.t $\\beta_0^C$ will be:\n",
    "$$\n",
    "\\begin{align}\n",
    "\\beta_0^C &= \\frac{1}{N}(\\sum_{i=1}^N y_i- \\sum_{i=1}^N\\sum_{j=1}^p{x_{ij}\\beta_j +\n",
    "\\sum_{i=1}^N\\sum_{j=1}^p\\overline{x}_j\\beta_j})\n",
    "\\\\\n",
    "&=\\frac{1}{N}(\\sum_{i=1}^N y_i- \\sum_{j=1}^p\\beta_j\\sum_{i=1}^N x_{ij} + \n",
    "N\\sum_{j=1}^p\\overline{x}_j\\beta_j)\n",
    "\\\\\n",
    "&=\\frac{1}{N}(\\sum_{i=1}^N y_i- \\sum_{j=1}^p\\beta_j\\overline{x}_jN +\n",
    "N\\sum_{j=1}^p\\overline{x}_j\\beta_j)\n",
    "\\\\\n",
    "&=\\overline{y}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "This completes proof."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We assume that the centering has been done, so that the matrix **X** has p columns (3.43):\n",
    "\n",
    "$$\n",
    "RSS(\\gamma)=(\\mathbf{y} - \\mathbf{X}\\beta)^T(\\mathbf{y} - \\mathbf{X}\\beta) + \\gamma\\beta^T\\beta\n",
    "$$\n",
    "\n",
    "the ridge regression solutions are easily seen to be:\n",
    "$$\n",
    "\\hat{\\beta}^{ridge}=(\\mathbf{X}^T\\mathbf{X}+\\gamma\\mathbf{I})^{-1}\\mathbf{X}^T\\mathbf{y}\n",
    "$$\n",
    "\n",
    "where **I** is the $p\\times p$ identity matrix. Notice that the solution is again a linear function of **y**.\n",
    "\n",
    "*Proof*:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\cfrac{\\partial{RSS(\\gamma)}}{\\partial{\\beta}}\n",
    "&=\\cfrac{\\partial{((\\mathbf{y} - \\mathbf{X}\\beta)^T(\\mathbf{y} - \\mathbf{X}\\beta) + \\gamma\\beta^T\\beta)}}{\\partial{\\beta}}\\\\\n",
    "&=\\cfrac{\\partial{\\left(\n",
    "  \\mathbf{y}^T\\mathbf{y} - 2\\mathbf{y}^T\\mathbf{X}\\beta \n",
    "  + (\\mathbf{X}\\beta)^T\\mathbf{X}\\beta \n",
    "  + \\gamma\\beta^T\\beta \n",
    "\\right)}}{\\partial{\\beta}}\\\\\n",
    "&=-2\\mathbf{y}^T\\mathbf{X}+2\\beta^T\\mathbf{X}^T\\mathbf{X}+2\\gamma\\beta^T\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "We set the first derivative to zero:\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\cfrac{\\partial{RSS(\\gamma)}}{\\partial{\\beta}} = 0\\\\\n",
    "-\\mathbf{y}^T\\mathbf{X}+\\beta^T\\mathbf{X}^T\\mathbf{X}+\\gamma\\beta^T=0\\\\\n",
    "\\beta^T(\\mathbf{X}^T\\mathbf{X}+\\gamma\\mathbf{I})=\\mathbf{y}^T\\mathbf{X}\\\\\n",
    "\\beta^T=\\mathbf{y}^T\\mathbf{X}(\\mathbf{X}^T\\mathbf{X}+\\gamma\\mathbf{I})^{-1}\\\\\n",
    "\\beta=(\\mathbf{X}^T\\mathbf{X}+\\gamma\\mathbf{I})^{-1}\\mathbf{X}^T\\mathbf{y}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Figure 3.8 shows the ridge coefficient estimates for the prostate cancer example, plotted as functions of $df(\\gamma)$, the *effective degrees of freedom* implied by the penalty $\\gamma$. In the case of orthonormal inputs, the ridge estimates are: $\\hat{\\beta}^{ridge}=\\hat{\\beta}/(1+\\gamma)$\n",
    "\n",
    "**TODO:** Implement Figure 3.8.\n",
    "\n",
    "Ridge expression can be also derived as the mean or mode of a posterior distribution, with a prior distribution. Suppose, $y_i \\sim N(\\beta_0+x_i^T\\beta, \\sigma^2)$ and the $\\beta_j \\sim N(0, \\mathcal{T}^2)$, independently of one another. Then the log-posterior density of $\\beta$, is equal to the expression in (3.41), with $\\gamma=\\sigma^2/\\mathcal{T}^2$.\n",
    "\n",
    "**TODO**: proof."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The *singular value decomposition* (SVD) of the centered input matrix **X** gives us some insight into the nature of ridge regression. The SVD of the $N \\times p$ matrix **X** (3.45):\n",
    "$$\n",
    "\\mathbf{X}=\\mathbf{UD}\\mathbf{V}^T\n",
    "$$\n",
    "- **U** - $N \\times p$ orthogonal matrix, the columns of **U** span the column space of **X**.\n",
    "- **V** - $N \\times p$ orthogonal matrix, the columns of **V** span the row space of **X**.\n",
    "- **D** - $p \\times p$ diagonal matrix, with diagonal entries $d_1 \\ge d_2 \\ge ... \\ge d_p \\ge 0$, called singular values of **X**. If one or more values $d_j = 0$, **X** is singular.\n",
    "\n",
    "Using SVD we can write the least squares fitted vector as (3.46):\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathbf{X}\\hat{\\beta}^{ls} &= \\mathbf{X}(\\mathbf{X}^T\\mathbf{X})^{-1}\\mathbf{X}^T\\mathbf{y}\\\\\n",
    "&= \\mathbf{UD}\\mathbf{V}^T(\\mathbf{VD}\\mathbf{U}^T\\mathbf{UD}\\mathbf{V}^T)^{-1}\\mathbf{VD}\\mathbf{U}^T\\mathbf{y}\\\\\n",
    "&= \\mathbf{UD}\\mathbf{V}^T(\\mathbf{VD}\\mathbf{D}\\mathbf{V}^T)^{-1}\\mathbf{VD}\\mathbf{U}^T\\mathbf{y}\\\\\n",
    "&= \\mathbf{UD}\\mathbf{V}^T[\\mathbf{V}^T]^{-1}\\mathbf{D}^{-1}\\mathbf{D}^{-1}\\mathbf{V}^{-1}\\mathbf{VD}\\mathbf{U}^T\\mathbf{y}\\\\\n",
    "&= \\mathbf{UD}\\mathbf{D}^{-1}\\mathbf{D}^{-1}\\mathbf{D}\\mathbf{U}^T\\mathbf{y}\\\\\n",
    "&= \\mathbf{U}\\mathbf{U}^T\\mathbf{y}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the ridge solutions are (3.47):\n",
    "    \n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathbf{X}\\hat{\\beta}^{ridge}&=\\mathbf{X}(\\mathbf{X}^T\\mathbf{X}+\\gamma\\mathbf{I})^{-1}\\mathbf{X}^T\\mathbf{y}\\\\\n",
    "&=\\mathbf{UD}\\mathbf{V}^T\n",
    "(\\mathbf{V}\\mathbf{D}^2\\mathbf{V}^T+\\gamma\\mathbf{V}\\mathbf{V}^T)^{-1}\n",
    "\\mathbf{VD}\\mathbf{U}^T\\mathbf{y}\\\\\n",
    "&=\\mathbf{UD}\\mathbf{V}^T\n",
    "(\\mathbf{V}(\\mathbf{D}^2+\\gamma\\mathbf{I})\\mathbf{V}^T)^{-1}\n",
    "\\mathbf{VD}\\mathbf{U}^T\\mathbf{y}\\\\\n",
    "&=\\mathbf{UD}\\mathbf{V}^T\n",
    "\\mathbf{V}(\\mathbf{D}^2+\\gamma\\mathbf{I})^{-1}\\mathbf{V}^T\n",
    "\\mathbf{VD}\\mathbf{U}^T\\mathbf{y}\\\\\n",
    "&=\\mathbf{UD}(\\mathbf{D}^2+\\gamma\\mathbf{I})^{-1}\\mathbf{D}\\mathbf{U}^T\\mathbf{y}\\\\\n",
    "&=\\sum_{j=1}^p\\mathbf{u}_j\\cfrac{d_j^2}{d_j^2+\\gamma}\\mathbf{u}_j^T\\mathbf{y}\n",
    "\\end{align}\n",
    "$$\n",
    "- Computes the coordinates of **y** w.r.t the orthonormal basis **U**. \n",
    "- Then shrinks these coordinates of **y** by the factors  $d_j^2/(d_j^2+\\gamma)$. \n",
    "- A greater amount of shrinkage is applied to the coordinates of basis vectors with smaller $d_j^2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Principal Components**\n",
    "\n",
    "What does a small value of $d_j^2$ mean? The sample covariance matrix is $\\mathbf{S}=\\mathbf{X}^T\\mathbf{X}/N$ (3.46):\n",
    "$$\n",
    "\\mathbf{X}^T\\mathbf{X} = \\mathbf{VD}^2\\mathbf{V}^T\n",
    "$$\n",
    "\n",
    "which is *eigen decomposition* of $\\mathbf{X}^T\\mathbf{X}$, aka the *principal components* directions of **X**.\n",
    "\n",
    "The first principal component (PC) direction $v_1$ has the largest variance, that is:\n",
    "$$\n",
    "\\begin{align}\n",
    "Var(\\mathbf{z}_1) &= Var(\\mathbf{X}v_1)\\\\\n",
    "&= \\frac{1}{N} v_1^T\\mathbf{VD}^2\\mathbf{V}^Tv_1\\\\\n",
    "&= \\frac{1}{N} \n",
    "   \\begin{pmatrix} 1&0 & \\dots & 0 \\end{pmatrix} \n",
    "   \\mathbf{D}^2\n",
    "   \\begin{pmatrix} 1 \\\\ 0 \\\\ \\vdots \\\\ 0 \\end{pmatrix}\\\\\n",
    "&=\\frac{d_1^2}{N}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "and in fact $\\mathbf{z}_1 = \\mathbf{X}v_1 = \\mathbf{UD}\\mathbf{V}^{T}v_1= \\mathbf{u}_1d_1$ and it is called the first PC. Ridge regression shrinks  PC directions having small variance the most."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Effective degrees of freedom**\n",
    "\n",
    "In Figure 3.7 we have plotted the estimated prediction error versus the quantity:\n",
    "$$\n",
    "\\begin{align}\n",
    "df(\\gamma) &= tr[\\mathbf{X}(\\mathbf{X}^T\\mathbf{X}+\\gamma\\mathbf{I})^{-1}\\mathbf{X}^T]\\\\\n",
    "&= tr(\\mathbf{H_{\\gamma}})\\\\\n",
    "&= tr[\\mathbf{UD}(\\mathbf{D}^2+\\gamma\\mathbf{I})^{-1}\\mathbf{D}\\mathbf{U}^T]\\\\\n",
    "&= \\sum_{j=1}^p \\cfrac{d_j^2}{d_j^2+\\gamma}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "This monotone decreasing function of $\\gamma$ is the *effective degrees of freedom* of the ridge regression fit. Note that: \n",
    "- $\\gamma = 0 \\text{  then  } df(\\gamma)=p$\n",
    "- $\\gamma \\rightarrow \\infty \\text{  then  } df(\\gamma) \\rightarrow 0$\n",
    "- Although all p coefficients in a ridge fit will be non-zero, they are fit in a restricted fashion controlled by $\\gamma$."
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}