how to invert the matrix with quite small eigenvalues
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 ------------------------------