Last active
December 18, 2015 03:19
-
-
Save grayclhn/5717763 to your computer and use it in GitHub Desktop.
Code for LU vs. QR decomposition in regression. Idempotence of projection matrices using the LU decomposition fails faster.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # 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