From: Ted Harding <Ted.Harding_at_nessie.mcc.ac.uk>

Date: Mon 30 May 2005 - 18:51:12 EST

E-Mail: (Ted Harding) <Ted.Harding@nessie.mcc.ac.uk> Fax-to-email: +44 (0)870 094 0861

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Mon May 30 19:02:48 2005

Date: Mon 30 May 2005 - 18:51:12 EST

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@nessie.mcc.ac.uk> Fax-to-email: +44 (0)870 094 0861

Date: 30-May-05 Time: 09:43:56 ------------------------------ XFMail ------------------------------ ______________________________________________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 Mon May 30 19:02:48 2005

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