From: Charilaos Skiadas <cskiadas_at_gmail.com>

Date: Wed, 05 Mar 2008 08:48:51 -0500

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:58:04 GMT

Date: Wed, 05 Mar 2008 08:48:51 -0500

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

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:58:04 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.
*