sessionInfo()

library(Matrix)

set.seed(2020) # seed

n <- 5
p <- 3
(X <- matrix(rnorm(n * p), nrow=n)) # predictor matrix
(y <- rnorm(n))    # response vector

# find the (minimum L2 norm) least squares solution
lm(y ~ X - 1)

cgs <- function(X) {  # not in-place
    n <- dim(X)[1]
    p <- dim(X)[2]    
    Q <- Matrix(0, nrow=n, ncol=p, sparse=FALSE)
    R <- Matrix(0, nrow=p, ncol=p, sparse=FALSE)
    for (k in seq_len(p)) {
        Q[, k] <- X[, k]
        for (j in seq_len(k-1)) {
            R[j, k] = sum(Q[, j] * X[, k])  # dot product
            Q[, k] <- Q[, k] - R[j, k] * Q[, j]
        }
        R[k, k] <- Matrix::norm(Q[, k, drop=FALSE], "F")
        Q[, k] <- Q[, k] / R[k, k]
    }
    list(Q=Q, R=R)
}

e <- .Machine$double.eps
(A <- t(Matrix(c(1, 1, 1, e, 0, 0, 0, e, 0, 0, 0, e), nrow=3)))

res <- cgs(A)
res$Q

with(res, t(Q) %*% Q)

mgs <- function(X) { #  in-place
    n <- dim(X)[1]
    p <- dim(X)[2]
    R <- Matrix(0, nrow=p, ncol=p)
    for (k in seq_len(p)) {
        for (j in seq_len(k-1)) {
            R[j, k] <- sum(X[, j] * X[, k])  # dot product
            X[, k] <- X[, k] - R[j, k] * X[, j]
        }
        R[k, k] <- Matrix::norm(X[, k, drop=FALSE], "F")
        X[, k] <- X[, k] / R[k, k]
    }
    list(Q=X, R=R)
}

mres <- mgs(A)
mres$Q

with(mres, t(Q) %*% Q)

(B <- t(Matrix(c(0.7, 1/sqrt(2), 0.7 + e, 1/sqrt(2)), nrow=2)))

mres2 <- mgs(B)
mres2$Q

with(mres2, t(Q) %*% Q)

X  # to recall
y

coef(lm(y ~ X - 1)) # least squares solution by QR

# same as
qr.solve(X, y)  # or solve(qr(X), y)

(qrobj <- qr(X, LAPACK=TRUE))

qr.Q(qrobj)  # extract Q

qr.R(qrobj) # extract R

qr.qty(qrobj, y)  # this uses the "compact form"

t(qr.Q(qrobj)) %*% y  # waste of resource

solve(qr.R(qrobj), qr.qty(qrobj, y)[1:p])  # solves R * beta = Q' * y

Xchol <- Matrix::chol(Matrix::symmpart(t(X) %*% X))  # solution using Cholesky of X'X
solve(Xchol, solve(t(Xchol), t(X) %*% y))

qrA <- qr(A, LAPACK=TRUE)
qr.Q(qrA)

t(qr.Q(qrA)) %*% qr.Q(qrA)   # orthogonality preserved

qrB <- qr(B, LAPACK=TRUE)
qr.Q(qrB)

t(qr.Q(qrB)) %*% qr.Q(qrB)  # orthogonality preserved