# housed.r         Demonstration of Householder for QR
 #                     step by step
 #                  response vector
  y <-c( 4.5, 5.5, 6.5, 8, 10, 12)
 #                  design matrix -- intercept, linear, quadratic
  X <- matrix(1,6,3)          # get shape and intercept, fix later
  X[,2] <- 1:6
  X[,3] <- (1:6)*(1:6)
 #                  put X and y together for easy printing
  Xy <- cbind(X,y)  # put them together
  Xy                # write it out
  I6 <- diag(c(1,1,1,1,1,1))      # will need this
 #
 #                  first iteration
 #
  u1 <- Xy[,1]      # construct first Householder tfm
  s2 <- crossprod(u1,u1)
  s <- sqrt(s2)*sign(u1[1])   # choose sign to avoid cancellation
  d <- 1/(s2 + s*u1[1])
  u1[1] <- u1[1] + s 
  U1 <- I6 - outer(d*u1,u1)
  Xy <- U1 %*% Xy   # first step
  Xy                # write it out
 # coefficients and error ss
  Xy[1,4]/Xy[1,1]   # should be ybar
  crossprod(Xy[2:6,4],Xy[2:6,4])  # should be css
 # modify for later use
  u1 <- u1/s
 #
 #                  second iteration
 #
  u2 <- Xy[2:6,2]      # construct second Householder tfm
  s2 <- crossprod(u2,u2)
  s <- sqrt(s2)*sign(u2[1])   # choose sign to avoid cancellation
  d <- 1/(s2 + s*u2[1])
  u2[1] <- u2[1] + s
  U2 <- I6             # start with this and modify 
  U2[2:6,2:6] <- I6[2:6,2:6] - outer(d*u2,u2)
  Xy <- U2 %*% Xy   # first step
  Xy                # write it out
 # coefficients and error ss
  backsolve(Xy[1:2,1:2],Xy[1:2,4])   # should be two coefs
  crossprod(Xy[3:6,4],Xy[3:6,4])  # should be css
 # modify for later use
  u2 <- u2/s
 #
 #                  third iteration
 #
  u3 <- Xy[3:6,3]      # construct third Householder tfm
  s2 <- crossprod(u3,u3)
  s <- sqrt(s2)*sign(u3[1])   # choose sign to avoid cancellation
  d <- 1/(s2 + s*u3[1])
  u3[1] <- u3[1] + s
  U3 <- I6             # start with this and modify 
  U3[3:6,3:6] <- I6[3:6,3:6] - outer(d*u3,u3)
  Xy <- U3 %*% Xy   # first step
  Xy                # write it out
 # coefficients and error ss
  backsolve(Xy[1:3,1:3],Xy[1:3,4])   # should be three coefs
  crossprod(Xy[4:6,4],Xy[4:6,4])  # should be css
 # modify for later use
  u3 <- u3/s
 #
 #
 #*** to understand qr output, first multiply U3*U2*U1 = Q'
 #
  U3 %*% U2 %*% U1
 # and print these vectors
  u1
  u2 
  u3
 #
 #**** now do all of this using qr function
  qrstuff <- qr(X)
  qrstuff
 # get Q -- not usually needed -- compare to Q' above
 #
  qr.Q( qrstuff, complete=TRUE )
 #
 # use qr.qty  to multiply Q'y
  Qty <- qr.qty(qrstuff,y)
  Qty               # write it out
 #
 # use qr.qy to multiply by Q and get fitted values 
  yhat <- qr.qy(qrstuff,c(Qty[1:3],rep(0,3)) )
  yhat              # write out fitted values
 #
 # reverse zeros to get residuals
  ehat <- qr.qy(qrstuff,c(rep(0,3),Qty[4:6]) )
  ehat              # write out residuals 
 #                  done
  rm(list=ls())     # clean up
 q()