Sweep operator¶

versioninfo()

Julia Version 1.4.1
Commit 381693d3df* (2020-04-14 17:20 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin18.7.0)
  CPU: Intel(R) Core(TM) i7-6920HQ CPU @ 2.90GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-8.0.1 (ORCJIT, skylake)
Environment:
  JULIA_EDITOR = code
  JULIA_NUM_THREADS = 4

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¶

SweepOperator.jl package (by Josh Day) implements the sweep operator.

using SweepOperator

A = [9. 2 -2; 2 1 0; -2 0 4]

3×3 Array{Float64,2}:
  9.0  2.0  -2.0
  2.0  1.0   0.0
 -2.0  0.0   4.0

B = copy(A)
sweep!(B, 1) # sweep (1, 1) entry

3×3 Array{Float64,2}:
 -0.111111  0.222222  -0.222222
  2.0       0.555556   0.444444
 -2.0       0.0        3.55556

sweep!(B, 2) # sweep (2, 2) entry

3×3 Array{Float64,2}:
 -0.2   0.4  -0.4
  2.0  -1.8   0.8
 -2.0   0.0   3.2

sweep!(B, 3) # sweep (3, 3) entry, we are left with -inv(B)

3×3 Array{Float64,2}:
 -0.25   0.5  -0.125
  2.0   -2.0   0.25
 -2.0    0.0  -0.3125

# check correctness
inv(A)

3×3 Array{Float64,2}:
  0.25   -0.5    0.125
 -0.5     2.0   -0.25
  0.125  -0.25   0.3125

using LinearAlgebra

# sweep! function only changes the upper triangular part
UpperTriangular(inv(A)) ≈ UpperTriangular(- B)

true

# inverse sweep to bring negative inverse back to original matrix
sweep!(B, 1:3, true)

3×3 Array{Float64,2}:
  9.0  2.0  -2.0
  2.0  1.0  -1.11022e-16
 -2.0  0.0   4.0