HH Function

> set.seed(12345)
> A <- matrix(runif(16), nrow = 4, ncol = 4)
> Y <- A
> m <- dim(A)[1]
> n <- dim(A)[2]
> k <- 1
> norm2 <- function(x) {
+     sqrt(sum(x^2))
+ }
> V <- matrix(0, nrow = 4, ncol = 4)
> for (k in 1:n) {
+     Qt <- diag(4)
+     x <- A[k:m, k]
+     e1 <- c(1, rep(0, (length(x) - 1)))
+     Vk <- sign(x[1]) * norm2(x) * e1 + x
+     Vk <- Vk/norm2(Vk)
+     A[k:m, k:n] <- A[k:m, k:n] - 2 * Vk %*% t(Vk) %*% A[k:m,
+         k:n]
+     V[k:m, k] <- Vk
+ }