[R] qr and Moore-Penrose
Torsten Hothorn
hothorn at amadeus.statistik.uni-dortmund.de
Wed Jun 30 11:12:24 CEST 1999
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 at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list