[R] qr and Moore-Penrose

[Torsten Hothorn]

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._