chapter-03/3.4.4-least-angle-regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3.4.4 Least Angle Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Least angle regression (LAR) uses a similar strategy to Forwarf stepwise regression, but only enters \"as much\" of a predictor as it deserves. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Algorithm 3.2**\n",
    "\n",
    "1. Standardize the predictors to have mean zero and unit norm. Start with the residual $\\mathbf{r} = \\mathbf{y} - \\mathbf{\\overline{y}}$ and $\\beta_1,...,\\beta_p = 0$\n",
    "\n",
    "2. Find the predictor $\\mathbf{x}_j$ most correlated with $\\mathbf{r}$.\n",
    "\n",
    "3. Move $\\beta_j$ from 0 towards its least-squares coefficient $\\langle \\mathbf{x}_j, \\mathbf{r} \\rangle$, until some other competitor $\\mathbf{x}_k$ has as much correlation with the current residual as does $\\mathbf{x}_j$.\n",
    "\n",
    "4. Move $\\beta_j$ and $\\beta_k$ in the direction defined by their joint least squares coefficient of the current residual on $\\langle \\mathbf{x}_j, \\mathbf{x}_k \\rangle$, until some other competitor $\\mathbf{x}_l$ has as much correlation with the current residual.\n",
    "\n",
    "5. Continue in this way until all $p$ predictors have been entered. After min(N - 1, p) steps, we arrive at the full least-squares solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose at the beginning of the kth step:\n",
    "\n",
    "- $\\mathcal{A}_k$ is the active set of variables\n",
    "\n",
    "- $\\beta_{\\mathcal{A}_k}$ be the coefficients \n",
    "\n",
    "- $\\mathbf{r}_k=\\mathbf{y} - \\mathbf{X}_{\\mathcal{A}_k}\\beta_{\\mathcal{A}_k}$ is the current residual,\n",
    "\n",
    "then the direction for this step is (3.55):\n",
    "\n",
    "$$\\delta_k = (\\mathbf{X}_{\\mathcal{A}_k}^T\\mathbf{X}_{\\mathcal{A}_k})^{-1}\\mathbf{X}_{\\mathcal{A}_k}^T\\mathbf{r}_k$$\n",
    "\n",
    "The coefficient profile then evolves as $\\beta_{\\mathcal{A}_k}(\\alpha)=\\beta_{\\mathcal{A}_k} + \\alpha \\cdot \\delta_k$ and the fit vector evolves as $\\hat{f}_k(\\alpha)=\\hat{f}_k + \\alpha \\cdot \\mathbf{u}_k$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "\n",
    "df = pd.read_csv('../data/prostate/prostate.data', delimiter='\\t', index_col=0)\n",
    "mask_train = df.pop('train')\n",
    "df_y = df.pop('lpsa')\n",
    "\n",
    "train_x = df[mask_train == 'T']\n",
    "train_y = df_y[mask_train == 'T']\n",
    "\n",
    "train_x_centered = train_x - train_x.mean(axis = 0)\n",
    "train_x_centered /= np.linalg.norm(train_x_centered, axis=0)\n",
    "train_y_centered = train_y - train_y.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lars(X, y):\n",
    "    n, p = X.shape\n",
    "    mu = np.zeros_like(y)\n",
    "    beta = np.zeros(p)\n",
    "\n",
    "    for _ in range(p):\n",
    "        c = X.T @ (y - mu) \n",
    "        c_abs = np.abs(c)\n",
    "        c_max = c_abs.max() \n",
    "\n",
    "        active = np.isclose(c_abs, c_max)\n",
    "        signs = np.where(c[active] > 0, 1, -1)\n",
    "      \n",
    "        X_active = signs * X[:, active]\n",
    "\n",
    "        G = X_active.T @ X_active\n",
    "        Ginv = np.linalg.inv(G)\n",
    "\n",
    "        A = Ginv.sum() ** (-0.5)\n",
    "\n",
    "        w = A * Ginv.sum(axis = 1)\n",
    "        u = X_active @ w\n",
    "\n",
    "        gamma = c_max / A\n",
    "\n",
    "        if not np.all(active):\n",
    "            a = X.T @ u\n",
    "            complement = np.invert(active)\n",
    "            cc = c[complement]\n",
    "            ac = a[complement]\n",
    "            candidates = np.concatenate([(c_max - cc) / (A - ac),\n",
    "                                         (c_max + cc) / (A + ac)])\n",
    "            gamma = candidates[candidates >= 0].min()\n",
    "        mu += gamma * u\n",
    "        beta[active] += gamma * signs\n",
    "    return mu, beta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Beta:  [10.06592433  6.58925225 -1.95834047  4.00665484  5.62572779 -1.65081245\n",
      " -0.20495795  3.9589639 ]\n",
      "train error:  0.43919976805833433\n"
     ]
    }
   ],
   "source": [
    "y_fit, beta = lars(train_x_centered.as_matrix(), train_y_centered.as_matrix())\n",
    "train_error = np.mean((y_fit - train_y_centered) ** 2)\n",
    "print ('Beta: ', beta)\n",
    "print ('train error: ', train_error)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Algorithm 3.2a**\n",
    "\n",
    "4a. If a non-zero coefficient hits zero, drop its variable from the active set of variables and recompute the current joint least squares direction.\n",
    "\n",
    "The LAR(lasso) algorithm is extremely efficient, requiring the same order of computation as that of a single least squares fit using the p predictors.\n",
    "\n",
    "**Heuristic argument why LAR and Lasso are similar**\n",
    "\n",
    "Suppose $\\mathcal{A}$ is the active set of variables at some stage. We can express as (3.56):\n",
    "$$\\mathbf{x}_j^T(\\mathbf{y}-\\mathbf{X}\\beta)=\\lambda \\cdot s_j, j \\in \\mathcal{A}$$\n",
    "\n",
    "also $|\\mathbf{x}_j^T(\\mathbf{y}-\\mathbf{X}\\beta)| \\le \\lambda, j \\notin \\mathcal{A}$. Now consider the lasso criterian (3.57):\n",
    "\n",
    "$$R(\\beta)=\\frac{1}{2}||\\mathbf{y}-\\mathbf{X}\\beta||_2^2 + \\lambda||\\beta||_1$$\n",
    "\n",
    "Let $\\mathcal{B}$ be the active set of variables in the solution for a given value of $\\lambda$, and $R(\\beta)$ is differentiable, and the stationarity conditions give (3.58):\n",
    "\n",
    "$$\\mathbf{x}_j^T(\\mathbf{y}-\\mathbf{X}\\beta)=\\lambda \\cdot sign(\\beta_j), j \\in \\mathcal{B}$$\n",
    "\n",
    "Comparing (3.56) and (3.58), we see that they are identical only if the sign of $\\beta{j}$ matches the sign of the inner product. That is why the LAR algorithm and lasso starts to differ when an active coefficient passes through zero; The stationary conditions for the non-active variable require that (3.59):\n",
    "\n",
    "$$|\\mathbf{x}_j^T(\\mathbf{y}-\\mathbf{X}\\beta)|\\le \\lambda, j \\notin \\mathcal{B}$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Degrees-of-Freedom Formula for LAR and Lasso\n",
    "\n",
    "We define the degrees of freedom of the fitted vector $\\hat{y}$ as:\n",
    "\n",
    "$$\n",
    "df(\\hat{y})=\\frac{1}{\\sigma^2}\\sum_{i=1}^N Cov(\\hat{y}_i,y_i)\n",
    "$$\n",
    "\n",
    "This makes intuitive sense: the harder that we fit to the data, the larger this covariance and hence $df(\\hat{\\mathbf{y}})$."
   ]
  }
 ],
 "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
}