Sweep() and Inv()

Sweep() may be used for creating a generalized inverse of a square matrix using Gaussian elimination. A generalized inverse of the matrix A, say A+ satisfies the properties, A+AA+ = A+ and AA+A = A. Also, the quantity A+A is an idempotent matrix of zeros and ones, so its rank is its trace. These are handy properties when doing linear models.

The arguments of Sweep() are the rows to be eliminated, from k1 to k2, and the matrix to operate on. If all rows are eliminated, and the matrix has full rank, then the matrix is inverted. The function, Inv(), performs matrix inversion by eliminating all rows and columns. If any pivot is less than 1.0E-8, then the row and column being eliminated is set to zero. This procedure creates A+.

Sweep() also simplifies the calculation of regression coefficients and the residual mean square. Consider the regression of XNxp on YNx1. The steps for calculating the regression coefficients are as follows:

  1. Concatenate X and Y using the statement VMatrix XY = Ch(X,Y); This produces an Nx(p+1) matrix.
  2. Calculate the matrix of crossproducts using VMatrix xypxy = Tran(XY)*XY;
  3. Calculate the coefficients by sweeping out the first p rows of xypxy, VMatrix ixypxy = Sweep(1,p,xypxy);
  4. The regression coefficients are in the right most column: VMatrix betahat = Submat(ixypxy,p,p+1,1,p+1);
  5. The mean square error is calculated from the bottom right corner: double mse = ixypxy.m(p+1,p+1) / ((double)(X.r-X.c));
This results from sweeping out the columns of X'X. The original matrix is partitioned as

xypxy = $\displaystyle \left[\vphantom{
\begin{array}{cc}
X'X & X'Y \\
Y'X & Y'Y
\end{array}
}\right.$$\displaystyle \begin{array}{cc}
X'X & X'Y \\
Y'X & Y'Y
\end{array}$$\displaystyle \left.\vphantom{
\begin{array}{cc}
X'X & X'Y \\
Y'X & Y'Y
\end{array}
}\right]$

The swept matrix is given by

Sweep(1, p, xpyxy) = $\displaystyle \left[\vphantom{
\begin{array}{cc}
(X'X)^{-1} & (X'X)^{-1}X'Y \\
-Y'X(X'X)^{-1} & Y'(I-X(X'X)^{-1}X')Y
\end{array}
}\right.$$\displaystyle \begin{array}{cc}
(X'X)^{-1} & (X'X)^{-1}X'Y \\
-Y'X(X'X)^{-1} & Y'(I-X(X'X)^{-1}X')Y
\end{array}$$\displaystyle \left.\vphantom{
\begin{array}{cc}
(X'X)^{-1} & (X'X)^{-1}X'Y \\
-Y'X(X'X)^{-1} & Y'(I-X(X'X)^{-1}X')Y
\end{array}
}\right]$

The top right matrix is the solution to the least squares normal equations. More information about the sweep operator can be found in Kennedy and Gentle[17], Heiberger[11] and Hocking[12].