thr3ads.net - R help - [R] generalized inverse using matinv (Design) [Aug 2011]

If this information is useful, please help other people find it:
Share via:

Frank Paetzold

2011-Aug-16 14:15 UTC

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

David Winsemius

2011-Aug-16 15:13 UTC

head link

[R] generalized inverse using matinv (Design)

On Aug 16, 2011, at 10:15 AM, Frank Paetzold wrote:
> i am trying to use matinv from the Design package
Which has been replaced by the rms package (about 2 years ago).
> 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
It would have been courteous to supply paste()-able code. If you use  
this code and then paste in the above matrix, then others would not  
need to do it themselves:

XtX = matrix(scan(), byrow=TRUE, ncol=16)
>
> this matrix has rank=9
The rankMatrix::Matrix function agrees. And in essence, so does the  
eigen function, since the last 7 eigenvalues are effectively zero:

 > eigen(XtX)
$values
  [1]  1.600000e+01  4.000000e+00  4.000000e+00  4.000000e+00   
4.000000e+00
  [6]  1.000000e+00  1.000000e+00  1.000000e+00  1.000000e+00   
1.065814e-14
[11]  4.440892e-15  4.173627e-16  2.220446e-16  1.717376e-16  
-1.193166e-16
[16] -9.593521e-16

>
> however, matinv(X'X) falsely returns rank=4,
Well, the current version of matinv returns a different value (= 6) ,  
although it is not 9.
> 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
See the FAQ regarding issues related to "==" and numerical accuracy.
>
> is not satisfied.
>
> if i use qr (from the base package) the rank is
> correctly determined as 9.
Since matinv is not described as being designed to deal with rank- 
deficient matrices, why have you chosen it over a function designed  
for the sort of problem you are dealing with. You remember the old  
joke about the guy who goes to the doctor and complains: "Doc, it  
hurts a lot when I do ....".

I guess the more modern Zen-oriented question might be: "And how has  
that method been working out for you?"

??"generalized inverse"

I see two "generalized inverse"  functions in the packages I have  
installed: ginv::MASS and pinv::pracma

require(MASS)
 > all.equal(XtX %*% ginv(XtX) %*% XtX ,  XtX)
[1] TRUE


-- 
David Winsemius, MD
West Hartford, CT

Doran, Harold

2011-Aug-16 15:33 UTC

head link

[R] generalized inverse using matinv (Design)

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.

Maybe Matching Threads

Search for more seemingly similar threads

R help - Aug 2011 - generalized inverse using matinv (Design)

[R] generalized inverse using matinv (Design)

[R] generalized inverse using matinv (Design)

[R] generalized inverse using matinv (Design)

Maybe Matching Threads