[R] Matrix inversion

[Bill.Venables at csiro.au]

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
>