In [1]:
versioninfo()
Julia Version 1.1.0
Commit 80516ca202 (2019-01-21 21:24 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin14.5.0)
  CPU: Intel(R) Core(TM) i7-6920HQ CPU @ 2.90GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, skylake)
Environment:
  JULIA_EDITOR = code

Sweep operator

Definition

  • We learnt Cholesky decomposition and QR decomposition approaches for solving linear regression.

  • The popular statistical software SAS uses sweep operator for linear regression and matrix inversion.

  • Assume $\mathbf{A}$ is symmetric and positive semidefinite.

  • Sweep on the $k$-th diagonal entry $a_{kk} \ne 0$ yields $\widehat A$ with entries $$ \begin{eqnarray*} \widehat a_{kk} &=& - \frac{1}{a_{kk}} \\ \widehat a_{ik} &=& \frac{a_{ik}}{a_{kk}} \\ \widehat a_{kj} &=& \frac{a_{kj}}{a_{kk}} \\ \widehat a_{ij} &=& a_{ij} - \frac{a_{ik} a_{kj}}{a_{kk}}, \quad i \ne k, j \ne k. \end{eqnarray*} $$ $n^2$ flops (taking into account of symmetry).

  • Inverse sweep sends $\mathbf{A}$ to $\check A$ with entries $$ \begin{eqnarray*} \check a_{kk} &=& - \frac{1}{a_{kk}} \\ \check a_{ik} &=& - \frac{a_{ik}}{a_{kk}} \\ \check a_{kj} &=& - \frac{a_{kj}}{a_{kk}} \\ \check a_{ij} &=& a_{ij} - \frac{a_{ik}a_{kj}}{a_{kk}}, \quad i \ne k, j \ne k. \end{eqnarray*} $$ $n^2$ flops (taking into account of symmetry).

  • $\check{\hat{\mathbf{A}}} = \mathbf{A}$.

  • Successively sweeping all diagonal entries of $\mathbf{A}$ yields $- \mathbf{A}^{-1}$.

  • Exercise: invert a $2 \times 2$ matrix, say $$ \mathbf{A} = \begin{pmatrix} 4 & 3 \\ 3 & 2 \end{pmatrix}, $$ on paper using sweep operator.

  • Block form of sweep: Let the symmetric matrix $\mathbf{A}$ be partitioned as $$ \mathbf{A} = \begin{pmatrix} \mathbf{A}_{11} & \mathbf{A}_{12} \\ \mathbf{A}_{21} & \mathbf{A}_{22} \end{pmatrix}. $$ If possible, sweep on the diagonal entries of $\mathbf{A}_{11}$ yields

$$ \begin{pmatrix} - \mathbf{A}_{11}^{-1} & \mathbf{A}_{11}^{-1} \mathbf{A}_{12} \\ \mathbf{A}_{21} \mathbf{A}_{11}^{-1} & \mathbf{A}_{22} - \mathbf{A}_{21} \mathbf{A}_{11}^{-1} \mathbf{A}_{12} \end{pmatrix}. $$

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}$.

  • 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}$.

Applications

Linear regression (as done in SAS)

Sweep on the (augmented) Gram matrix $$ \begin{pmatrix} \mathbf{X}, \mathbf{y} \end{pmatrix}^T \begin{pmatrix} \mathbf{X}, \mathbf{y} \end{pmatrix} = \begin{pmatrix} \mathbf{X}^T \mathbf{X} & \mathbf{X}^T \mathbf{y} \\ \mathbf{y}^T \mathbf{X} & \mathbf{y}^T \mathbf{y} \end{pmatrix} $$
yields

$$ \begin{eqnarray*} \begin{pmatrix} - (\mathbf{X}^T \mathbf{X})^{-1} & (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{y} \\ \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} \end{pmatrix} = \begin{pmatrix} - \sigma^{-2} \text{Var}(\widehat{\beta}) & \widehat{\beta} \\ \widehat{\beta}^T & \|\mathbf{y} - \widehat{\mathbf{y}}\|_2^2 \end{pmatrix}. \end{eqnarray*} $$


In total $np^2 + p^3$ flops.

Invert a matrix in place

$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.

Conditional multivariate normal density calculation

Stepwise regression

MANOVA

Implementation

In [2]:
using SweepOperator

A = Float64.([9 2 -2;
    2 1 0;
    -2 0 4])
Out[2]:
3×3 Array{Float64,2}:
  9.0  2.0  -2.0
  2.0  1.0   0.0
 -2.0  0.0   4.0
In [3]:
B = copy(A)
sweep!(B, 1) # sweep (1, 1) entry
Out[3]:
3×3 Array{Float64,2}:
 -0.111111  0.222222  -0.222222
  2.0       0.555556   0.444444
 -2.0       0.0        3.55556 
In [4]:
sweep!(B, 2) # sweep (2, 2) entry
Out[4]:
3×3 Array{Float64,2}:
 -0.2   0.4  -0.4
  2.0  -1.8   0.8
 -2.0   0.0   3.2
In [5]:
sweep!(B, 3) # sweep (3, 3) entry, we are left with -inv(B)
Out[5]:
3×3 Array{Float64,2}:
 -0.25   0.5  -0.125 
  2.0   -2.0   0.25  
 -2.0    0.0  -0.3125
In [6]:
# check correctness
inv(A)
Out[6]:
3×3 Array{Float64,2}:
  0.25   -0.5    0.125 
 -0.5     2.0   -0.25  
  0.125  -0.25   0.3125
In [7]:
# inverse sweep to bring negative inverse back to original matrix
sweep!(B, 1:3, true)
Out[7]:
3×3 Array{Float64,2}:
  9.0  2.0  -2.0        
  2.0  1.0  -1.11022e-16
 -2.0  0.0   4.0        

Further reading