Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
Code for LU vs. QR decomposition in regression. Idempotence of projection matrices using the LU decomposition fails faster.
# Explicitly calculate the inverse by LU decomposition: LAPACK DGESV,
# then calculate the projection matrix X(X'X)^(-1)X'
projection.LU1 <- function(x)
x %*% solve(crossprod(x)) %*% t(x)
# Use DGESV but without explicitly calculating the inverse
projection.LU2 <- function(x)
crossprod(t(x), solve(crossprod(x), t(x)))
# Calculate the projection matrix using the QR decomposition
projection.QR <- function(x) {
QR <- qr(x)
tcrossprod(qr.Q(QR)[, QR$pivot[seq_len(QR$rank)], drop = FALSE])
}
X <- outer(seq(0, 1, 0.02), 0:10, "^")
all.equal(projection.LU1(X), projection.LU1(X) %*% projection.LU1(X)) # False
all.equal(projection.LU2(X), projection.LU2(X) %*% projection.LU2(X)) # False but less inaccurate
all.equal(projection.QR(X), projection.QR(X) %*% projection.QR(X)) # True
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment