Householder Triangularization

> library(Matrix)
> 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))
+ }
> Q <- matrix(data = 0, nrow = 16, ncol = 4)
> Q[seq(1, 16, 4), 1] <- 1
> Q[seq(2, 16, 4), 2] <- 1
> Q[seq(3, 16, 4), 3] <- 1
> Q[seq(4, 16, 4), 4] <- 1
> 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]
+     Qt[k:m, k:n] <- diag(length(x)) - 2 * Vk %*% t(Vk)
+     Q[(4 * k - 3):(4 * k), ] <- Qt
+ }

> print(round(A, 8))
          [,1]       [,2]       [,3]       [,4]
[1,] -1.628187 -0.7206857 -0.9536335 -0.7608957
[2,]  0.000000  0.2857674 -0.3555428  0.5260942
[3,]  0.000000  0.0000000  0.7054906  0.1160929
[4,]  0.000000  0.0000000  0.0000000  0.1973830

> Qf <- Q[1:4, ] %*% Q[5:8, ] %*% Q[9:12, ] %*% Q[13:16, ]
> print(round(Qf, 8))
           [,1]        [,2]       [,3]       [,4]
[1,] -0.4427649  0.48076393  0.6752776  0.3417974
[2,] -0.5378825 -0.77430980  0.2856081 -0.1719148
[3,] -0.4673803 -0.04107847 -0.6035221  0.6446932
[4,] -0.5442401  0.40941780 -0.3133516 -0.6618085

> print(round(Qf %*% A, 8))
          [,1]      [,2]       [,3]       [,4]
[1,] 0.7209039 0.4564810 0.72770525 0.73568495
[2,] 0.8757732 0.1663718 0.98973694 0.00113659
[3,] 0.7609823 0.3250954 0.03453544 0.39120334
[4,] 0.8861246 0.5092243 0.15237349 0.46249465
> print(round(Y, 8))
          [,1]      [,2]       [,3]       [,4]
[1,] 0.7209039 0.4564810 0.72770525 0.73568495
[2,] 0.8757732 0.1663718 0.98973694 0.00113659
[3,] 0.7609823 0.3250954 0.03453544 0.39120334
[4,] 0.8861246 0.5092243 0.15237349 0.46249465

> for (k in n:1) {
+     w <- as.vector(V[[3]])
+     w %*% t(w) %*% c(1, 0)
+ }