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

From: Duncan Murdoch <murdoch_at_stats.uwo.ca>
Date: Wed, 05 Mar 2008 08:43:23 -0500

On 3/5/2008 8:21 AM, gerardus vanneste wrote:
> 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)

If the matrix actually isn't singular, then a rescaling of the parameters should help a lot. I see the diagonal of infomatrix as

 > diag(infomatrix)
[1] 5.930720e-03 3.872339e+02 4.562529e+07 6.281634e+12 9.228140e+17 [6] 1.398687e+23 2.154124e+28 3.345598e+33

so multiplying the parameters by numbers on the order of the square roots of these entries (e.g. 10^c(-1, 1, 4, 6, 9, 12, 14, 17)), and redoing the rest of the calculations on that scale, should work.

Duncan Murdoch
>
> 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.


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:52:42 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 - 15:30: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