Date: Wed, 30 Jun 1999 11:12:24 +0200 (MET DST)
From: Torsten Hothorn <hothorn@amadeus.statistik.uni-dortmund.de>
To: r-help@stat.math.ethz.ch
Subject: [R] qr and Moore-Penrose
In-Reply-To: <Pine.GSO.4.05.9906300738070.5595-100000@auk.stats>
yesterday I had a little shock using qr (or lm). having a matrix
X <- cbind(1,diag(3))
y <- 1:3
the qr.coef returns one NA (because X is singular). So I computed the
Moore-Penrose inverse of X (just from the definition) and I get a correct
result. Whats wrong about qr in this situation?
here is the code:
mpinv <- function(X)
{
# Moore-Penrose inverse
Eps <- 100 * .Machine$double.eps;
# singular value decomposition
s <- svd(X);
d <- s$d;
m <- length(d);
if (!(is.vector(d)))
return(t(s$v%*%(1/d)%*%t(s$u)));
# remove eigenvalues equal zero
d <- d[d > Eps];
notnull <- length(d);
if (notnull == 1)
{
inv <- 1/d;
} else {
inv <- solve(diag(d));
}
# add rows, columns of zeros if needed
if (notnull != m)
{
inv <- cbind(inv, matrix(0, nrow=notnull, ncol=(m - notnull)));
inv <- rbind(inv, matrix(0, nrow=(m-notnull), ncol=m));
}
# compute Moore-Penrose
mp <- s$v%*%inv%*%t(s$u);
# set very small values to zero
mp[abs(mp) < Eps] <- 0;
return(mp);
};
X <- cbind(1, diag(3)); # singular matrix
y <- 1:3
Xp <- qr(X);
b1 <- qr.coef(Xp, y); # contains NA
b2 <- mpinv(X)%*%y # least square fit using Moore-Penrose
X%*%b2 # == y
Torsten
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._