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

[Douglas Bates]

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