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

From: huang min <minhuangr_at_gmail.com>
Date: Mon 30 May 2005 - 19:39:52 EST

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@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@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:46:49 2005

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