{ "cells": [ { "cell_type": "markdown", "metadata": { "toc": "true" }, "source": [ "# Table of Contents\n", "

1  Sweep operator
1.1  Definition
1.2  Applications
1.2.1  Linear regression (as done in SAS)
1.2.2  Invert a matrix in place
1.2.3  Conditional multivariate normal density calculation
1.2.4  Stepwise regression
1.2.5  MANOVA
1.3  Implementation
1.4  Further reading
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Sweep operator\n", "\n", "## Definition\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 " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\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}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* 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", "### 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", "### Invert a matrix _in place_\n", "\n", "$n^3$ flops. 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", "### Conditional multivariate normal density calculation\n", "\n", "### Stepwise regression\n", "\n", "### MANOVA\n", "\n", "\n", "## Implementation\n", "\n", "* [SweepOperator.jl](https://github.com/joshday/SweepOperator.jl) package (by Josh Day) implements the sweep operator." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "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": [ "# Pkg.add(\"SweepOperator\")\n", "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": {}, "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": {}, "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": {}, "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. " ] }, { "cell_type": "code", "execution_count": 5, "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": "68px", "width": "252px" }, "navigate_menu": true, "number_sections": true, "sideBar": true, "threshold": 4, "toc_cell": true, "toc_section_display": "block", "toc_window_display": true, "widenNotebook": false } }, "nbformat": 4, "nbformat_minor": 2 }