# 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

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.

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