R help - Mar 2008 - matrix inversion using solve() and matrices containing large/small values

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

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.

Sorry, I meant to send this to the whole list.

On Mar 5, 2008, at 8:46 AM, Charilaos Skiadas wrote:
> The problem doesn't necessarily have to do with the range of data.  
> At first level, it has to do with the simple fact that dfdb has  
> rank 6 at most, (7 at most in general, though in your case since  
> 290=10*29, it is at most 6). Since matrix rank only goes down when  
> multiplying, your infomatrix is an 8x8 matrix, with rank at most 6  
> (7 if you were more lucky with that 290, still not good enough), so  
> it is certainly not invertible, even forgetting the computational  
> issues of computing the inverse.
>
> You would need either power smaller than 7, or a longer dosis  
> vector, I think.
>
> If you manage to make your problem in a case where dfdb is square,  
> then you just have to invert that, which might be easier, seeing  
> how it's a Vandermonde matrix.
>
> Haris Skiadas
> Department of Mathematics and Computer Science
> Hanover College
>
> On Mar 5, 2008, at 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)
>>
>> 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]]
>>
>
>
>
>
Haris Skiadas
Department of Mathematics and Computer Science
Hanover College