Re: [R] how to invert the matrix with quite small eigenvalues

From: Thomas Lumley <tlumley_at_u.washington.edu>
Date: Tue 31 May 2005 - 02:24:43 EST

On Mon, 30 May 2005, huang min wrote:>
> My intention is to invert the covariance matrix to perform some
> algorithm which is common in the estimating equations like GEE.

In that case there is no benefit in being able to invert very extreme covariance matrices. The asymptotic approximations to the distribution of regression coefficients will break down really badly with such extreme working covariance matrices.

I think in a case like this you should either 1) report an error and stop
2) shrink the covariance matrix towards a diagonal one, eg increase the diagonal entries until the condition number becomes reasonable. 3) Use a one-step estimator from the independence working model (which is asymptotically equivalent to the full solution and better behaved).

Remember that in GEE the matrix V^{-1} is just a set of weights, chosen to get good efficiency. Your matrix solve(a) is not a good set of weights.

I think in an earlier thread on this topic Brian Ripley recommended using the singular value decomposition if you really have to compute something like D^TV^{-1}D. In your example this still isn't good enough.

>
> Another strange point is that my friend use the LU decomposition in
> Fortran to solve the inverse matrix of aa for me. He used double
> precision of course, otherwise no inverse matrix in Fortran too. He
> show that the product of the two matrix is almost identity(The biggest
> off-digonal element is about 1e-8). I copy his inverse matrix(with 31
> digits!) to R and read in aa I sent to him(16 digits). The product is
> also not an identity matrix. It is fairly strange! Any suggestion?
>

It isn't that strange. The system is computationally singular, meaning that you should expect to get different answers with apparently similar computations on different systems.

Also remember that what you care about for GEE is the result of solve(a,y-mu), rather than solve(a,a). Getting one of these right is no guarantee of getting the other one right.

If you really had to work with this matrix in double precision you would need to track very carefully the error bounds on all your computations, which would be very difficult. Fortunately this is almost never necessary in statistics, and I don't think it's necessary in your case.

A good habit to get into when you aren't tracking error bounds carefully is to think of the last couple of digits of the result of any calculation as random.

         -thomas



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Tue May 31 02:29:35 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:32:16 EST