{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Summary of linear regression\n",
    "\n",
    "Methods for solving linear regression $\\widehat \\beta = (\\mathbf{X}^T \\mathbf{X})^{-1} \\mathbf{X}^T \\mathbf{y}$:\n",
    "\n",
    "| Method            | Flops            | Remarks                 | Software | Stability   |\n",
    "|-------------------|------------------|-------------------------|----------|-------------|\n",
    "| Sweep             | $np^2 + p^3$     | $(X^TX)^{-1}$ available | SAS      | less stable |\n",
    "| Cholesky          | $np^2 + p^3/3$   |                         |          | less stable |\n",
    "| QR by Householder | $2np^2 - (2/3)p^3$ |                         | R        | stable      |\n",
    "| QR by MGS         | $2np^2$          | $Q_1$ available         |          | stable |  \n",
    "| QR by SVD         | $4n^2p + 8np^2 + 9p^3$         |  $X = UDV^T$        |          | most stable |  \n",
    "\n",
    "Remarks:\n",
    "\n",
    "0. When $n \\gg p$, sweep and Cholesky are twice faster than QR and need less space.  \n",
    "0. Sweep and Cholesky are based on the **Gram matrix** $\\mathbf{X}^T \\mathbf{X}$, which can be dynamically updated with incoming data. They can handle huge $n$, moderate $p$ data sets that cannot fit into memory.  \n",
    "0. QR methods are more stable and produce numerically more accurate solution.  \n",
    "0. Although sweep is slower than Cholesky, it yields standard errors and so on.  \n",
    "0. MGS appears slower than Householder, but it yields $\\mathbf{Q}_1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "> **There is simply no such thing as a universal 'gold standard' when it comes to algorithms.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "linreg_svd (generic function with 1 method)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using SweepOperator, BenchmarkTools\n",
    "\n",
    "linreg_cholesky(y::Vector, X::Matrix) = cholfact!(X'X) \\ (X'y)\n",
    "\n",
    "linreg_qr(y::Vector, X::Matrix) = X \\ y\n",
    "\n",
    "function linreg_sweep(y::Vector, X::Matrix)\n",
    "    p = size(X, 2)\n",
    "    tableau = [X y]' * [X y]\n",
    "    sweep!(tableau, 1:p)\n",
    "    return tableau[1:p, end]\n",
    "end\n",
    "\n",
    "function linreg_svd(y::Vector, X::Matrix)\n",
    "    xsvd = svdfact(X)\n",
    "    return xsvd[:V] * ((xsvd[:U]'y) ./ xsvd[:S])\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "linreg_cholesky(y, X) = [0.390365, 0.262759, 0.149047]\n",
      "linreg_qr(y, X) = [0.390365, 0.262759, 0.149047]\n",
      "linreg_sweep(y, X) = [0.390365, 0.262759, 0.149047]\n",
      "linreg_svd(y, X) = [0.390365, 0.262759, 0.149047]\n"
     ]
    }
   ],
   "source": [
    "srand(280) # seed\n",
    "\n",
    "n, p = 10, 3\n",
    "X = randn(n, p)\n",
    "y = randn(n)\n",
    "\n",
    "# check these methods give same answer\n",
    "@show linreg_cholesky(y, X)\n",
    "@show linreg_qr(y, X)\n",
    "@show linreg_sweep(y, X)\n",
    "@show linreg_svd(y, X);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  708.34 KiB\n",
       "  allocs estimate:  10\n",
       "  --------------\n",
       "  minimum time:     2.236 ms (0.00% GC)\n",
       "  median time:      2.578 ms (0.00% GC)\n",
       "  mean time:        2.672 ms (1.37% GC)\n",
       "  maximum time:     5.330 ms (24.89% GC)\n",
       "  --------------\n",
       "  samples:          1865\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n, p = 1000, 300\n",
    "X = randn(n, p)\n",
    "y = randn(n)\n",
    "\n",
    "@benchmark linreg_cholesky(y, X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  6.03 MiB\n",
       "  allocs estimate:  922\n",
       "  --------------\n",
       "  minimum time:     9.249 ms (0.00% GC)\n",
       "  median time:      11.248 ms (0.00% GC)\n",
       "  mean time:        11.382 ms (2.68% GC)\n",
       "  maximum time:     16.164 ms (9.21% GC)\n",
       "  --------------\n",
       "  samples:          439\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@benchmark linreg_sweep(y, X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  4.05 MiB\n",
       "  allocs estimate:  2470\n",
       "  --------------\n",
       "  minimum time:     12.299 ms (0.00% GC)\n",
       "  median time:      15.969 ms (0.00% GC)\n",
       "  mean time:        15.946 ms (1.73% GC)\n",
       "  maximum time:     19.512 ms (8.74% GC)\n",
       "  --------------\n",
       "  samples:          314\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@benchmark linreg_qr(y, X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  8.74 MiB\n",
       "  allocs estimate:  46\n",
       "  --------------\n",
       "  minimum time:     54.615 ms (0.00% GC)\n",
       "  median time:      75.000 ms (0.00% GC)\n",
       "  mean time:        72.860 ms (0.64% GC)\n",
       "  maximum time:     83.025 ms (0.00% GC)\n",
       "  --------------\n",
       "  samples:          69\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@benchmark linreg_svd(y, X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Julia Version 0.6.2\n",
      "Commit d386e40c17 (2017-12-13 18:08 UTC)\n",
      "Platform Info:\n",
      "  OS: macOS (x86_64-apple-darwin14.5.0)\n",
      "  CPU: Intel(R) Core(TM) i7-6920HQ CPU @ 2.90GHz\n",
      "  WORD_SIZE: 64\n",
      "  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)\n",
      "  LAPACK: libopenblas64_\n",
      "  LIBM: libopenlibm\n",
      "  LLVM: libLLVM-3.9.1 (ORCJIT, skylake)\n"
     ]
    }
   ],
   "source": [
    "versioninfo()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.2",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.2"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "31px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false,
   "widenNotebook": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}