Contents

swp()

Usage:

swp(x, n1 [, n2, ...] [, quiet:T, diag:d, tolerance:tol, keepswept:T]),
  x a REAL matrix, n1, n2, ... positive integers, or vectors of positive
  integers, tol > 0 a REAL scalar, d a REAL vector with positive
  elements

Keywords: matrix algebra, glm

swp(x,Cols) uses a form of the Beaton SWP operator on the REAL matrix x
to compute a REAL result of the same size.  If x is m by n, Cols should
be a vector whose elements are positive integers <= min(m,n).

A single SWP of x on row and column k produces a matrix y of the same
size as x with

     y[i,j] = x[i,j] - x[i,k]x[k,j]/x[k,k], for i and j not equal to k
     y[i,k] = x[i,k]/x[k,k], for i not equal to k
     y[k,j] = -x[k,j]/x[k,k], for j not equal to k
     y[k,k] = 1/x[k,k]

The element x[k,k] is sometimes called a "pivot".

swp(x,Cols) does successive SWPs of x on Cols[1], Cols[2], ... rows and
columns.  If, attempting to SWP row and column M, abs(pivot) is found to
be too small in comparison with the abs(x[M,M]), the message
  WARNING: tolerance failure pivoting column M; not swept
is printed and that SWP is skipped.

The threshold for finding tolerance failure may be modified by keyword
'tolerance'; see below.

If x is n by n and non-singular, swp(x,run(n)) computes the inverse of x
(it can fail for certain non-positive definite but invertible matrices).

swp(x,Cols1, Cols2, ...), where Cols1, Cols2, ..., are vectors of
positive integers is equivalent to swp(x,vector(Cols1, Cols2, ...)).
For example, swp(cp,1,2,3,4) is equivalent to swp(cp,run(4)).

swp() can be very useful in regression and analysis of variance.

swp(x,Cols,...,tolerance:Tol), where Tol is a small REAL scalar does the
same, but Tol is used in the test for what constitutes a small pivot.
Column M will be skipped if abs(pivot) < Tol*abs(x[M,M]), where pivot is
the value of x[M,M] just before a SWP of row and column M is attempted.

The default value of Tol is 10^(-12).

swp(x,Cols,...,quiet:T) does the same, but no warning message is printed
when a too small pivot is found.

swp(x,Cols,...,diag:d [,quiet:T]), where d is a REAL vector, does the
same, but the pivot for row and clumn M is compared with abs(d[M]) to
test whether a row and column will be swept.  The length of d must
always be min(nrows(x), ncols(x)).  This feature allows the sequence,
say,
  Cmd> d <- diag(x); x1 <- swp(x,1); x1 <- swp(x1,2,diag:d)
to be completely equivalent to
  Cmd> x1 <- swp(x,1,2)
Without diag:d, the comparison element for swp(x1,2) would be x1[2,2]
instead of x[2,2].

swp(x,Cols,..., keepswept:T [,tolerance:Tol, ,diag:d, quiet:F]) does the
same, except that the result is a structure with components 'matrix' and
'sweptcols'.  Component 'matrix' is the swept version of x.  Component
'sweptcols' is a vector of integers, with abs(sweptcols[i]) the number
of the i-th row and column swept.  sweptcols[i] < 0 if and only if row
and column abs(sweptcols[i]) failed the tolerance check.  No warning
message is printed unless 'quiet:F' is an argument.

The use of 'keepswept:T' may allow you to compute a generalized inverse
of a singular square matrix.

  Cmd> tmp <- swp(a, run(nrows(a)), keepswept:T)

  Cmd> b <- tmp$matrix; J <- (-tmp$sweptcols)[tmp$sweptcols < 0]

  Cmd> if(!isnull(J)){b[J,] <- b[,J] <- 0;;}

Now, even if a is singular, a %*% b %*% a should be a within rounding
error and b %*% a %*% b should be b within rounding error.

See also solve(), qr().