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.
On Aug 16, 2011, at 10:15 AM, Frank Paetzold wrote:> i am trying to use matinv from the Design packageWhich 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 1It 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=9The 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'XSee 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
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.