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

[huang min]

Maybe I should state more clear that I define b  to get the orthogonal
matrix bb$vectors.

We also can define diag(b)<-diag(b)+100, which will make the
eigenvalues of b much bigger to make sure the orthogonal matrix is
reliable.

My intention is to invert the covariance matrix to perform some
algorithm which is common in the estimating equations like GEE.

I meet difficulty to invert the covariance matrix. Two possibilities here:

1. The rounding error in defining the covariance matrix make the
eigenvalue to small.

2. The solve function in R can not cope with the matrix with so small
an eigenvalue.

For the first possibility, I think it can not be improved unless we
can  define more precise number than the double precision. So I ask
for the possiblity of coping with the second.

I can not find the default way to invert the matrix with solve().

For the symmetric matrix, I wonder if there are some algorithm which
can naturally make the inverse matrix symmetric and make sure it is
the inverse in the sense that the product is an identity matrix. I
know there are many decompositions which can be used to find the
inverse of a matrix. QR, SVD, Chol, and maybe some iterative method. I
wonder anyone can suggest me which algorithm might be good.

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?

Regards,

Huang


On 5/30/05, Ted Harding <Ted.Harding at nessie.mcc.ac.uk> wrote:
> On 30-May-05 huang min wrote:
> > Dear all,
> >
> > I encounter some covariance matrix with quite small eigenvalues
> > (around 1e-18), which are smaller than the machine precision. The
> > dimension of my matrix is 17. Here I just fake some small matrix for
> > illustration.
> >
> >  a<-diag(c(rep(3,4),1e-18)) # a matrix with small eigenvalues
> >  b<-matrix(1:25,ncol=5) # define b to get an orthogonal matrix
> >  b<-b+t(b)
> 
> NB: b is not *orthogonal*! Each row of b is equal to the preceding
> row plus (row2 - row1) (and similar for the columns), and b has
> rank 2 (though eigen(b) taken literally says that it has 5 non-zero
> eugenvalues; however, 3 of these are snaller that 10^(-14)).
> 
> Perhaps you meant "define b to get a symmetric matrix".
> 
> >  bb<-eigen(b,symmetric=T)
> >  aah<-bb$vectors%*%diag(1/sqrt(diag(a)))
> >  aa<-aah%*%t(aah) # aa should have the same eigenvalues as a and
> > should be #invertable,however,
> >  solve(aa) # can not be solved
> 
> Well, I did get a (non-symmetric) result for solve(aa) ...
> 
> >  solve(aa,tol=1e-19) # can be inverted, however, it is not symmetric
> > and furthermore,
> 
> and an idenitical (to solve(aa)) result for this.
> 
> >  solve(aa,tol=1e-19)%*%aa # deviate much from the identity matrix
> 
> But here I agree with you!
> 
> > I have already define aa to make sure it is symmetric. So the inverse
> > should be symmetric.
> >
> > Does the problem come from the rounding error since the eigenvalue is
> > smaller than the machine precision? In fact, the eigenvalue of aa is
> > negative now, but at least, it is still invertable. How can I get the
> > inverse? Thanks.
> 
> It does indeed, like the eigenvalue result for b above, come from
> the rounding error.
> 
> You should clarify in your mind why you want to ensure that you get
> correct results for matrices like these.
> 
> You are (and in your example deliberately so) treading on the very
> edge of what is capable of being computed, and results are very likely
> to lead to unexpected gross anomalies (such as being unable to
> invert a mathematically invertible matrix, or getting a non-symmetric
> inverse to a symmetric matrix [depending on the algorithm], or
> getting non-zero values for eigenvalues which should be zero, or
> the gross difference from the identity matrix which you expected).
> 
> It is like using a mechanical ditch-digger to peel an apple.
> 
> Exactly what will happen in any particular case can only be
> determined by a very fine-grained analysis of the operation
> of the numerical algorithm, which is beyond your reach in the
> normal usage of R.
> 
> Best wishes,
> Ted.
> 
> 
> --------------------------------------------------------------------
> E-Mail: (Ted Harding) <Ted.Harding at nessie.mcc.ac.uk>
> Fax-to-email: +44 (0)870 094 0861
> Date: 30-May-05                                       Time: 09:43:56
> ------------------------------ XFMail ------------------------------
>