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:
- Concatenate X and Y using the statement
VMatrix XY = Ch(X,Y); This produces an Nx(p+1) matrix.
- Calculate the matrix of crossproducts using
VMatrix xypxy = Tran(XY)*XY;
- Calculate the coefficients by sweeping out the first p
rows of xypxy, VMatrix ixypxy = Sweep(1,p,xypxy);
- The regression coefficients are in the right most
column: VMatrix betahat = Submat(ixypxy,p,p+1,1,p+1);
- 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
The swept matrix is given by
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[#!KG:sc!#],
Heiberger[#!HE:ad!#] and Hocking[#!HK:am!#].