Re: [R] Matrix inversion

From: <Bill.Venables_at_csiro.au>
Date: Tue, 19 Feb 2008 14:37:35 +1000

Ben Domingue asks:

> I am trying to invert a matrix for the purposes of least squares. I
> have tried a number of things, and the variety of results has me
> confused.

Don't be.

> 1. When I try solve() I get the following:
> >Error in solve.default(t(X) %*% X) : system is computationally
> singular: reciprocal condition number = 3.76391e-20

That seems clear enough.

> 2. When I try qr.solve(), I get:
> >Error in qr.solve(t(X) %*% X) : singular matrix 'a' in solve

That backs up what solve() said.

> 3. I can, however, use lm(y~X) to get coefficients. This confuses me
> since I thought that lm() used qr(). So why did qr.solve() not work
> earlier?

It may use the QR decomposition, but that does not necessarily imply that it uses the R function qr() in the same the way you did.

If you look at the help information for lm() you will see that there is an argument singular.OK with the default value of "TRUE". Essentially this agrument allows you to say whether or not you anticipate the model matrix *could be* not of full column rank. If so, the regression coefficients are not unique and what you will see in the coefficients is one particular version. The residuals and the residual sum of squares, however, are unique.

> 4. I have even tried using ginv(). The process works, but I end up
> with a different set of regression coefficients after I finish the
> process than what I had with lm(). To the best of my knowledge, this
> shouldn't happen.

Then you have learned something and it's a win. If the model matrix is not of full column rank, then the regression coefficients are not unique. This means that different solving methods can give different answers. Indeed you should expect this. lm() uses the QR decomposition; ginv() uses the singular value decomposition. You might like to see what aov() gives. It could be the same as lm() or it could be something different again. All three, though, will give you the same residuals up to computational accuracy.

>
> I've been digging around all day and can't figure this out. Thanks,
>
> Ben Domingue
> PhD Student, School of Education
> University of Colorado at Boulder
>



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. Received on Tue 19 Feb 2008 - 04:41:18 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Tue 19 Feb 2008 - 05:30:15 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive