Skip to content

Instantly share code, notes, and snippets.

Last active December 18, 2015 03:19
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