R help - May 2005 - how to invert the matrix with quite small eigenvalues

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)
 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
 solve(aa,tol=1e-19) # can be inverted, however, it is not symmetric
and furthermore,
 solve(aa,tol=1e-19)%*%aa # deviate much from the identity matrix

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.

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

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

On 30-May-05 huang min wrote:
> Maybe I should state more clear that I define b to get the
> orthogonal matrix bb$vectors.
OK. Certainly bbv<-bb$vectors is close to orthogonal: bbv%*%bbv
differs from the unit matrix only in that the off-diagonal
terms are O(10^(-16)).
> 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.
Which comes back to my general question: why do you need to ensure
that you get correct results for matrices like these? This matrix
is very nearly singular, and in many contexts you would interpret
it as exactly singular (e.g. if I got such a matrix as the
empirical covariance matrix from a sample, or by computation from
the structure of a model such as a design matrix, I would not
normally want to preserve this very small margin of non-singularity:
I would replace it with the lower-dimensional singular version,
e.g. by decomposing it with eigen() or svd() and reconstructing
the singular version after setting very small eigenvalues to 0).

But then of course there would be no inverse, which might be
required by the further computational expressions you want to
evaluate. However, in that case (depending on your application),
you may be able to proceed satisfactorily by using ginv() (see
package MASS) instead of solve(). But if you take that route,
then you have to bear in mind that a generalised inverse is
not a true inverse (if only because the latter does not exist):
the only property required for G=ginv(M) is that

  M%*%G%*%M = M

If such a matrix G will work for the rest of your calculations,
then you are OK, in particular if what you need is a possible
solution to an under-determined system of linear equations.
If not, then you should seriously consider whether your
computational strategy is suitable for the problem you have
in hand.
> [...]
> 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?
These phenomena are yet another instance of the fact that at the
margins of computability the results will be anomalouus. Your
friend is simply using a narrower margin, but it is still not
exact! Not strange at all.

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: 11:41:28
------------------------------ XFMail ------------------------------