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

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