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

[(Ted Harding)]

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 ------------------------------