An embedded and charset-unspecified text was scrubbed... Name: not available Url: https://stat.ethz.ch/pipermail/r-help/attachments/20080305/793a38ef/attachment.pl
matrix inversion using solve() and matrices containing large/small values
3 messages · gerardus vanneste, Duncan Murdoch, Douglas Bates
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.
On Wed, Mar 5, 2008 at 7:43 AM, Duncan Murdoch <murdoch at stats.uwo.ca> wrote:
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...
Someone with experience in numerical linear algebra will immediately focus on the words "inversion of a matrix" in that statement. There is a sort of a meta-result in numerical linear algebra that if the best way you can formulate an answer to your problem involves inversion of a matrix (without specifying a special structure like diagonal or bidiagonal or unit triangular or ...) then you need to think about your problem in more detail. Certainly a formula may be written out in terms of the inverse of an information matrix but it is a bad idea to try to calculate the result that way. If you have an information matrix is it positive definite? If so, you should be working with the Cholesky decomposition of the matrix or perhaps the QR decomposition of a model matrix that generates the information matrix. Think in terms of the factors of a matrix and how you can reexpress the calculation using them.
> 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.