Skip to content
{{ message }}

Instantly share code, notes, and snippets.

# grayclhn/idempotence_test.R

Last active Dec 18, 2015
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
to join this conversation on GitHub. Already have an account? Sign in to comment