# sweep.s                         J F Monahan, June 2001
 # Regression via sweep on Brownlee problem
 Xy <- matrix(scan('browncut.dat'),4,22)
 Xy <- matrix(c(rep(1,22),t(Xy)),22,5)
 "Augmented design matrix with responses"
 Xy
 "Sums of squares & crossproducts matrix -- start"
 A <- t(Xy) %*% Xy
 A
# now here's the sweep operator
  p <- 5
  N <- 22
  list <- seq(1,p)         # create full list
# loop through row/cols 
  for( k in c(1:4,1:4)) {
  print("sweep row/column")
  print(k) 
# sweep operator
  d <- A[k,k]
  A[k,] <- A[k,]/d
  lk <- list[-k]           # drop k-th row/col
  A[lk,lk] <- A[lk,lk] - outer( A[lk,k], A[k,lk] )
  A[lk,k] <- - A[lk,k]/d
  A[k,k] <- 1/d
  print(A) 
                 }
 rm(list=ls())