GramSchmidt - Householder Transformation

> gs <- function(A) {
+     m <- dim(A)[1]
+     n <- dim(A)[2]
+     Q <- matrix(data = NA, nrow = m, ncol = n)
+     Q[, 1] <- A[, 1]
+     Q[, 1] <- Q[, 1]/sqrt(sum(Q[, 1]^2))
+     for (j in 2:n) {
+         vj <- A[, j]
+         for (i in 1:(j - 1)) {
+             rij <- t(Q[, i]) %*% vj
+             vj <- vj - rij * Q[, i]
+         }
+         rjj <- sqrt(sum(vj^2))
+         qj <- vj/rjj
+         Q[, j] <- qj
+     }
+     return(Q)
+ }

> gs.mod <- function(A) {
+     m <- dim(A)[1]
+     n <- dim(A)[2]
+     Q <- matrix(data = NA, nrow = m, ncol = n)
+     for (j in 1:n) {
+         Q[, j] <- A[, j]
+     }
+     j <- 1
+     for (j in 1:n) {
+         rjj <- sqrt(sum(Q[, j]^2))
+         Q[, j] <- Q[, j]/rjj
+         if (j < n) {
+             for (i in (j + 1):n) {
+                 rij <- t(Q[, j]) %*% Q[, i]
+                 Q[, i] <- Q[, i] - rij * Q[, j]
+             }
+         }
+     }
+     return(Q)
+ }

> A <- matrix(c(1, 2, 3, 0, 1, 0.5, 0, 0.1, 1), nrow = 3, ncol = 3,
+     T)
> print(A)
     [,1] [,2] [,3]
[1,]    1  2.0  3.0
[2,]    0  1.0  0.5
[3,]    0  0.1  1.0
> print(gs(A))
     [,1]       [,2]        [,3]
[1,]    1 0.00000000  0.00000000
[2,]    0 0.99503719 -0.09950372
[3,]    0 0.09950372  0.99503719
> print(gs.mod(A))
     [,1]       [,2]        [,3]
[1,]    1 0.00000000  0.00000000
[2,]    0 0.99503719 -0.09950372
[3,]    0 0.09950372  0.99503719

> HHT <- function(A) {
+     A.bkup <- A
+     m <- dim(A)[1]
+     n <- dim(A)[2]
+     norm2 <- function(x) {
+         sqrt(sum(x^2))
+     }
+     Q <- matrix(data = 0, nrow = (n * m), ncol = m)
+     k <- 1
+     for (k in 1:n) {
+         Qt <- diag(m)
+         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)
+         Qt[k:m, k:n] <- diag(length(x)) - 2 * Vk %*% t(Vk)
+         A[k:m, k:n] <- A[k:m, k:n] - 2 * Vk %*% t(Vk) %*% A[k:m,
+             k:n]
+         Q[(4 * k - 3):(4 * k), ] <- Qt
+     }
+     Q.f <- diag(m)
+     for (i in 1:m) {
+         Q.f <- Q[(4 * i - 3):(4 * i), ] %*% Q.f
+     }
+     result <- list(A = A.bkup, A.t = A, Q = Q.f)
+ }

> set.seed(12345)
> A <- matrix(runif(16), nrow = 4, ncol = 4)
> print(HHT(A))
$A
          [,1]      [,2]       [,3]        [,4]
[1,] 0.7209039 0.4564810 0.72770525 0.735684952
[2,] 0.8757732 0.1663718 0.98973694 0.001136587
[3,] 0.7609823 0.3250954 0.03453544 0.391203335
[4,] 0.8861246 0.5092243 0.15237349 0.462494654

$A.t
              [,1]       [,2]          [,3]       [,4]
[1,] -1.628187e+00 -0.7206857 -9.536335e-01 -0.7608957
[2,]  0.000000e+00  0.2857674 -3.555428e-01  0.5260942
[3,] -1.110223e-16  0.0000000  7.054906e-01  0.1160929
[4,]  0.000000e+00  0.0000000 -1.110223e-16  0.1973830

$Q
           [,1]       [,2]        [,3]       [,4]
[1,] -0.4427649 -0.5378825 -0.46738026 -0.5442401
[2,]  0.4807639 -0.7743098 -0.04107847  0.4094178
[3,]  0.6752776  0.2856081 -0.60352213 -0.3133516
[4,]  0.3417974 -0.1719148  0.64469317 -0.6618085