Algo Comparison for Least Squares

> m <- 100
> n <- 15
> t <- seq(0, (m - 1))/(m - 1)
> A <- matrix(data = 0, nrow = m, ncol = n)
> for (j in 1:n) {
+     A[, j] <- t^(j - 1)
+ }
> b <- exp(sin(4 * t))
> b <- b/2006.78745308021

> 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:m] <- 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]
+         start.ind <- k * m - m + 1
+         end.ind <- k * m
+         Q[start.ind:end.ind, ] <- Qt
+     }
+     Q.f <- diag(m)
+     for (i in 1:n) {
+         start.ind <- i * m - m + 1
+         end.ind <- i * m
+         Q.f <- Q.f %*% Q[start.ind:end.ind, ]
+     }
+     result <- list(A = A.bkup, A.t = A, Q = Q.f)
+ }
> HHResult <- HHT(A)
> b.new <- t(HHResult$Q) %*% b
> R <- HHResult$A.t
> x <- backsolve(R, b.new)
> print(x[15], digits = 15)
[1] 1.00000012504902
> print(paste("relative error - hht", x[15] - 1))
[1] "relative error - hht 1.25049022692281e-07"

> gs.mod <- function(A) {
+     m <- dim(A)[1]
+     n <- dim(A)[2]
+     Q <- matrix(data = 0, nrow = m, ncol = n)
+     R <- matrix(data = 0, nrow = m, ncol = n)
+     V <- A
+     j <- 1
+     for (i in 1:n) {
+         rii <- sqrt(sum(V[, i]^2))
+         R[i, i] <- rii
+         Q[, i] <- V[, i]/rii
+         if (i < n) {
+             for (j in (i + 1):n) {
+                 rij <- t(Q[, i]) %*% V[, j]
+                 V[, j] <- V[, j] - rij * Q[, i]
+                 R[i, j] <- rij
+             }
+         }
+     }
+     result <- list(Q = Q, R = R)
+     return(result)
+ }
> Q <- gs.mod(A)$Q
> R <- gs.mod(A)$R
> b.new <- t(Q) %*% b
> x <- backsolve(R, b.new)
> print(x[15], digits = 15)
[1] 0.969854877113548
> print(paste("relative error - hht", x[15] - 1))
[1] "relative error - hht -0.0301451228864515"

> U <- svd(A)$u
> V <- svd(A)$v
> d <- diag(svd(A)$d)
> b.new <- t(U) %*% b
> w <- backsolve(d, b.new)
> x <- V %*% w
> print(x[15], digits = 15)
[1] 1.00000013199462
> print(paste("relative error - hht", x[15] - 1))
[1] "relative error - hht 1.31994622343257e-07"