[R] generalized inverse using matinv (Design)

[Doran, Harold]

Frank,

It is very rare that one needs to ever invert a matrix. This is particularly true if you are trying to solve a linear system of equations. Rather than offering advice to you on how to compute the inverse, can you indicate what you're trying to accomplish in the end? Maybe we can offer better solutions using decompositions. 

> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
> Behalf Of Frank Paetzold
> Sent: Tuesday, August 16, 2011 10:16 AM
> To: r-help at r-project.org
> Subject: [R] generalized inverse using matinv (Design)
> 
> i am trying to use matinv from the Design package
> to compute the generalized inverse of the normal equations
> of a 3x3 design via the sweep operator.
> 
> That is, for the linear model
> y = µ + x1 + x2 + x1*x2
> where x1, x2 are 3-level factors and dummy coding is being used
> 
> the matrix to be inverted is
> 
> X'X =
> 
> 9 3 3 3 3 3 3 1 1 1 1 1 1 1 1 1
> 3 3 0 0 1 1 1 1 0 0 1 0 0 1 0 0
> 3 0 3 0 1 1 1 0 1 0 0 1 0 0 1 0
> 3 0 0 3 1 1 1 0 0 1 0 0 1 0 0 1
> 3 1 1 1 3 0 0 1 1 1 0 0 0 0 0 0
> 3 1 1 1 0 3 0 0 0 0 1 1 1 0 0 0
> 3 1 1 1 0 0 3 0 0 0 0 0 0 1 1 1
> 1 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0
> 1 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0
> 1 0 0 1 1 0 0 0 0 1 0 0 0 0 0 0
> 1 1 0 0 0 1 0 0 0 0 1 0 0 0 0 0
> 1 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0
> 1 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0
> 1 1 0 0 0 0 1 0 0 0 0 0 0 1 0 0
> 1 0 1 0 0 0 1 0 0 0 0 0 0 0 1 0
> 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 1
> 
> this matrix has rank=9
> 
> however, matinv(X'X) falsely returns rank=4,
> no matter what the tolerance threshold eps is set to.
> 
> also the defining property of the
> generalized inverse
>              _
> X'X %*% (X'X)  %*% X'X = X'X
> 
> is not satisfied.
> 
> if i use qr (from the base package) the rank is
> correctly determined as 9.
> 
> any ideas?
> 
> thank you
> 
> --
> View this message in context: http://r.789695.n4.nabble.com/generalized-
> inverse-using-matinv-Design-tp3747337p3747337.html
> Sent from the R help mailing list archive at Nabble.com.
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.