[R] matrix inversion using solve() and matrices containing large/small values

From: gerardus vanneste <gerard.vanneste_at_gmail.com>
Date: Wed, 05 Mar 2008 14:21:20 +0100


Hello

I've stumbled upon a problem for inversion of a matrix with large values, and I haven't found a solution yet... I wondered if someone could give a hand. (It is about automatic optimisation of a calibration process, which involves the inverse of the information matrix)

code:



> macht=0.8698965
> coeff=1.106836*10^(-8)

> echtecoeff=c(481.46,19919.23,-93.41188,0.5939589,-0.002846272,8.030726e-6
,-1.155094e-8,6.357603e-12)/10000000
> dosis=c(0,29,70,128,201,290,396)

> dfdb <-

array(c(1,1,1,1,1,1,1,dosis,dosis^2,dosis^3,dosis^4,dosis^5,dosis^6,dosis^7),dim=c(7,8))

> dfdbtrans = aperm(dfdb)
> sigerr=sqrt(coeff*dosis^macht)
> sigmadosis = c(1:7)
> for(i in 1:7){
sigmadosis[i]=ifelse(sigerr[i]<2.257786084*10^(-4),2.257786084*10^(-4),sigerr[i])

}
> omega = diag(sigmadosis)
> infomatrix = dfdbtrans%*%omega%*%dfdb


I need the inverse of this information matrix, and

> infomatrix_inv = solve(infomatrix, tol = 10^(-43))

does not deliver adequate results (matrixproduct of infomatrix_inv and infomatrix is not 1). Regular use of solve() delivers the error 'system is computationally singular: reciprocal condition number: 2.949.10^(-41)'

So I went over to an eigendecomposition using eigen(): (so that infomatrix = V D V^(-1) ==> infomatrix^(-1)= V D^(-1) V^(-1) ) in the hope this would deliver better results.)



> infomatrix_eigen = eigen(infomatrix)
> infomatrix_eigen_D = diag(infomatrix_eigen$values)
> infomatrix_eigen_V = infomatrix_eigen$vectors
> infomatrix_eigen_V_inv = solve(infomatrix_eigen_V)

however, the matrix product of these are not the same as the infomatrix itself, only in certain parts:

> infomatrix_eigen_V %*% infomatrix_eigen_D %*% infomatrix_eigen_V_inv
> infomatrix

Therefore, I reckon the inverse of eigendecomposition won't be correct either.

As far as I understand, the problem is due to the very large range of data, and therefore results in numerical problems, but I can't come up with a way to do it otherwise.

Would anyone know how I could solve this problem?

(PS, i'm running under linux suse 10.0, latest R version with MASS libraries (RV package))

F. Crop
UGent -- Medical Physics

        [[alternative HTML version deleted]]



R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Wed 05 Mar 2008 - 13:24:41 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Wed 05 Mar 2008 - 14:00:18 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive