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