teaching/biostatm280-2017spring/slides/11-sweep/sweep.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sweep operator\n",
    "\n",
    "* We learnt Cholesky decomposition and QR decomposition approaches for solving linear regression.\n",
    "\n",
    "* The popular statistical software SAS uses sweep operator for linear regression and matrix inversion.\n",
    "\n",
    "* Assume $\\mathbf{A}$ is symmetric and positive semidefinite.\n",
    "\n",
    "* **Sweep** on the $k$-th diagonal entry $a_{kk} \\ne 0$ yields $\\widehat A$ with entries\n",
    "$$\n",
    "\\begin{eqnarray*}\n",
    "\t\\widehat a_{kk} &=& - \\frac{1}{a_{kk}} \\\\\n",
    "\t\\widehat a_{ik} &=& \\frac{a_{ik}}{a_{kk}} \\\\\n",
    "\t\\widehat a_{kj} &=& \\frac{a_{kj}}{a_{kk}} \\\\\n",
    "\t\\widehat a_{ij} &=& a_{ij} - \\frac{a_{ik} a_{kj}}{a_{kk}}, \\quad i \\ne k, j \\ne k.\n",
    "\\end{eqnarray*}\n",
    "$$\n",
    "$n^2$ flops (taking into account of symmetry).\n",
    "\n",
    "* **Inverse sweep** sends $\\mathbf{A}$ to $\\check A$ with entries\n",
    "$$\n",
    "\\begin{eqnarray*}\n",
    "\t\\check a_{kk} &=& - \\frac{1}{a_{kk}} \\\\\n",
    "\t\\check a_{ik} &=& - \\frac{a_{ik}}{a_{kk}} \\\\\n",
    "\t\\check a_{kj} &=& - \\frac{a_{kj}}{a_{kk}} \\\\\n",
    "\t\\check a_{ij} &=& a_{ij} - \\frac{a_{ik}a_{kj}}{a_{kk}}, \\quad i \\ne k, j \\ne k.\n",
    "\\end{eqnarray*}\n",
    "$$\n",
    "$n^2$ flops (taking into account of symmetry).\n",
    "\n",
    "* $\\check{\\hat{\\mathbf{A}}} = \\mathbf{A}$.\n",
    "\n",
    "* Successively sweeping all diagonal entries of $\\mathbf{A}$ yields $- \\mathbf{A}^{-1}$.\n",
    "\n",
    "* Exercise: invert a $2 \\times 2$ matrix, say \n",
    "$$\n",
    "\\mathbf{A} = \\begin{pmatrix} 4 & 3 \\\\ 3  & 2 \\end{pmatrix},\n",
    "$$\n",
    "on paper using sweep operator.\n",
    "\n",
    "* Block form of sweep: Let the symmetric matrix $\\mathbf{A}$ be partitioned as \n",
    "$$\n",
    "    \\mathbf{A} = \\begin{pmatrix} \\mathbf{A}_{11} & \\mathbf{A}_{12} \\\\ \\mathbf{A}_{21} & \\mathbf{A}_{22} \\end{pmatrix}.\n",
    "$$\n",
    "If possible, sweep on the diagonal entries of $\\mathbf{A}_{11}$ yields  \n",
    "\n",
    "$$\n",
    "\t\\begin{pmatrix} \n",
    "    - \\mathbf{A}_{11}^{-1} & \\mathbf{A}_{11}^{-1} \\mathbf{A}_{12} \\\\\n",
    "\t\\mathbf{A}_{21} \\mathbf{A}_{11}^{-1} & \\mathbf{A}_{22} - \\mathbf{A}_{21} \\mathbf{A}_{11}^{-1} \\mathbf{A}_{12}\n",
    "\t\\end{pmatrix}.\n",
    "$$\n",
    "Order dose **not** matter. The block $\\mathbf{A}_{22} - \\mathbf{A}_{21} \\mathbf{A}_{11}^{-1} \\mathbf{A}_{12}$ is recognized as the **Schur complement** of $\\mathbf{A}_{11}$.\n",
    "\n",
    "* Pd and determinant: $\\mathbf{A}$ is pd if and only if each diagonal entry can be swept in succession and is positive until it is swept. When a diagonal entry of a pd matrix $\\mathbf{A}$ is swept, it becomes negative and remains negative thereafter. Taking the product of diagonal entries just before each is swept yields the determinant of $\\mathbf{A}$. \n",
    "\n",
    "## Applications\n",
    "\n",
    "* Invert a matrix _in place_. $n^3$ flops.  \n",
    "Recall that inversion by Cholesky costs $(1/3)n^3 + (4/3) n^3 = (5/3) n^3$ flops and needs to allocate a matrix of same size.\n",
    "\n",
    "* Stepwise regression. \n",
    "\n",
    "* Conditional multivariate normal density calculation.\n",
    "\n",
    "* MANOVA.\n",
    "\n",
    "* Linear regression (as done in SAS).\n",
    "\n",
    "Sweep on the Gram matrix \n",
    "$$\n",
    "\\begin{pmatrix} \n",
    "    \\mathbf{X}^T \\mathbf{X} & \\mathbf{X}^T \\mathbf{y} \\\\ \n",
    "    \\mathbf{y}^T \\mathbf{X} & \\mathbf{y}^T \\mathbf{y} \n",
    "\\end{pmatrix}\n",
    "$$  \n",
    "yields  \n",
    "\n",
    "$$\n",
    "\\begin{eqnarray*}\n",
    "\\begin{pmatrix}\n",
    "- (\\mathbf{X}^T \\mathbf{X})^{-1} & (\\mathbf{X}^T \\mathbf{X})^{-1} \\mathbf{X}^T \\mathbf{y} \\\\\n",
    "\\mathbf{y}^T \\mathbf{X} (\\mathbf{X}^T \\mathbf{X})^{-1} & \\mathbf{y}^T \\mathbf{y} - \\mathbf{y}^T \\mathbf{X} (\\mathbf{X}^T \\mathbf{X})^{-1} \\mathbf{X}^T \\mathbf{y}\n",
    "\\end{pmatrix} = \n",
    "\\begin{pmatrix}\n",
    "- \\sigma^{-2} \\text{Cov}(\\beta) & \\beta \\\\\n",
    "\\beta^T & \\|\\mathbf{y} - \\hat y\\|_2^2\n",
    "\\end{pmatrix}.\n",
    "\\end{eqnarray*}\n",
    "$$  \n",
    "In total $np^2 + p^3$ flops.\n",
    "\n",
    "## Implementation\n",
    "\n",
    "* [SweepOperator.jl](https://github.com/joshday/SweepOperator.jl) package implements the sweep operator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×4 Array{Float64,2}:\n",
       "  9.589    -4.98208  -5.70052  -1.83933\n",
       " -4.98208   4.94476   3.81663   2.82462\n",
       " -5.70052   3.81663   5.45442   2.46008\n",
       " -1.83933   2.82462   2.46008   4.98343"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using SweepOperator\n",
    "\n",
    "srand(280)\n",
    "\n",
    "X = randn(5, 3) # predictor matrix\n",
    "y = randn(5)    # response vector\n",
    "\n",
    "# form the augmented Gram matrix\n",
    "G = [X y]' * [X y]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×4 Array{Float64,2}:\n",
       " -0.312747  -0.136595  -0.231278  0.379547\n",
       " -4.98208   -0.499383   0.206676  0.650887\n",
       " -5.70052    3.81663   -0.569668  0.39225 \n",
       " -1.83933    2.82462    2.46008   2.87807 "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sweep!(G, 1:3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " 0.379547\n",
       " 0.650887\n",
       " 0.39225 "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# least squares solution by QR\n",
    "X \\ y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×4 Array{Float64,2}:\n",
       "  9.589    -4.98208  -5.70052  -1.83933\n",
       " -4.98208   4.94476   3.81663   2.82462\n",
       " -5.70052   3.81663   5.45442   2.46008\n",
       " -1.83933   2.82462   2.46008   4.98343"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# inverse sweep to restore original matrix\n",
    "sweep!(G, 1:3, true)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Further reading\n",
    "\n",
    "* Section 7.4-7.6 of [Numerical Analysis for Statisticians](http://ucla.worldcat.org/title/numerical-analysis-for-statisticians/oclc/793808354&referer=brief_results) by Kenneth Lange (2010). Probably the best place to read about sweep operator.\n",
    "\n",
    "* The paper [A tutorial on the SWEEP operator](http://www.jstor.org/stable/2683825) by James H. Goodnight (1979). Note this version (anti-symmetry matrix) is slightly different from ours. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.5.2",
   "language": "julia",
   "name": "julia-0.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.5.2"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "68px",
    "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
}