Table of Contents¶

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 Gram matrix $$ \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{Cov}(\beta) & \beta \\ \beta^T & \|\mathbf{y} - \hat 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.

Further reading¶

• Section 7.4-7.6 of Numerical Analysis for Statisticians by Kenneth Lange (2010). Probably the best place to read about sweep operator.

• The paper A tutorial on the SWEEP operator by James H. Goodnight (1979). Note this version (anti-symmetry matrix) is slightly different from ours.