{ "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": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Julia Version 1.1.0\n", "Commit 80516ca202 (2019-01-21 21:24 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", " LIBM: libopenlibm\n", " LLVM: libLLVM-6.0.1 (ORCJIT, skylake)\n", "Environment:\n", " JULIA_EDITOR = code\n" ] } ], "source": [ "versioninfo()" ] }, { "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", "\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 (augmented) Gram matrix \n", "$$\n", "\\begin{pmatrix} \\mathbf{X}, \\mathbf{y} \\end{pmatrix}^T \\begin{pmatrix} \\mathbf{X}, \\mathbf{y} \\end{pmatrix} = \\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{Var}(\\widehat{\\beta}) & \\widehat{\\beta} \\\\\n", "\\widehat{\\beta}^T & \\|\\mathbf{y} - \\widehat{\\mathbf{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": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 9.0 2.0 -2.0\n", " 2.0 1.0 0.0\n", " -2.0 0.0 4.0" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using SweepOperator\n", "\n", "A = Float64.([9 2 -2;\n", " 2 1 0;\n", " -2 0 4])" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " -0.111111 0.222222 -0.222222\n", " 2.0 0.555556 0.444444\n", " -2.0 0.0 3.55556 " ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "B = copy(A)\n", "sweep!(B, 1) # sweep (1, 1) entry" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " -0.2 0.4 -0.4\n", " 2.0 -1.8 0.8\n", " -2.0 0.0 3.2" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sweep!(B, 2) # sweep (2, 2) entry" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " -0.25 0.5 -0.125 \n", " 2.0 -2.0 0.25 \n", " -2.0 0.0 -0.3125" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sweep!(B, 3) # sweep (3, 3) entry, we are left with -inv(B)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 0.25 -0.5 0.125 \n", " -0.5 2.0 -0.25 \n", " 0.125 -0.25 0.3125" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# check correctness\n", "inv(A)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " 9.0 2.0 -2.0 \n", " 2.0 1.0 -1.11022e-16\n", " -2.0 0.0 4.0 " ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# inverse sweep to bring negative inverse back to original matrix\n", "sweep!(B, 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": { "@webio": { "lastCommId": null, "lastKernelId": null }, "kernelspec": { "display_name": "Julia 1.1.0", "language": "julia", "name": "julia-1.1" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.1.0" }, "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, "skip_h1_title": true, "threshold": 4, "toc_cell": true, "toc_section_display": "block", "toc_window_display": true, "widenNotebook": false } }, "nbformat": 4, "nbformat_minor": 2 }