[R] a numeric problem

[Duncan Murdoch]

On 07/03/2011 8:46 AM, baicaidoufu wrote:
> ### An numeric problem in R  ########
>
> ###I have two matrix one is##########
>
> A<- matrix(c(21.97844, 250.1960, 2752.033, 29675.88, 316318.4, 3349550,
> 35336827,
>                24.89267, 261.4211, 2691.009, 27796.02, 288738.7, 3011839,
> 31498784,
>                21.80384, 232.3765, 2460.495, 25992.77, 274001.6, 2883756,
> 30318645,
>                39.85801, 392.2341, 3971.349, 40814.22, 423126.2, 4409388,
> 46090850,
>                25.55343, 235.5315, 2343.281, 23987.10, 248557.3, 2590821,
> 27089162,
>                16.63100, 195.4592, 2200.264, 24006.66, 257499.5, 2736200,
> 28922851,
>                23.20258, 262.6086, 2824.971, 29973.61, 316369.1, 3331009,
> 35025965
>                ),7,7,byrow=T)
>
> #### the other one is##################
> S<-matrix(c(356.57205,  32.55943, 120.28162,  95.74285,  45.013261,
> 271.557635, 183.9009,
>               32.55943, 299.68103,  81.95644, 246.95280,  96.023270,
> 19.383732, 153.4544,
>              120.28162,  81.95644, 277.09354,  85.42180, 138.751600,
> 138.972234, 140.0991,
>               95.74285, 246.95280,  85.42180, 527.78444, 182.417527,
> 24.946245, 129.1490,
>               45.01326,  96.02327, 138.75160, 182.41753, 328.414655,
> 1.035543,  63.2730,
>              271.55764,  19.38373, 138.97223,  24.94625,   1.035543,
> 322.635669, 158.5317,
>              183.90092, 153.45443, 140.09909, 129.14899,  63.273005,
> 158.531662, 256.1531
>              ),7,7,byrow=T)
>
> ####both of these two matrix are non-singular so their inverse should be
> exists######
> ####while
> KK<- t(A)%*%solve(S)%*%A
> det(KK) ###12553.48 which is not 0, but solve() function refuse to work now.

Look at the eigenvalues of KK:  they range from 10^(-7) to 10^12.  Its 
condition number (according to the rcond definition) is about 10^(-21).
If you could invert it, you would just be seeing noise from rounding error.

So the best advice I can give is:  don't do that.  But if you really 
don't care about the quality of your results, you could use solve(A) %*% 
S %*% t(solve(A)) to avoid the error messages.  If you multiply that by 
KK you'll get something that's a long way from an identity matrix.

Duncan Murdoch

> ##  So I try to use ginv() instead###
> ginv(KK)
>
> ### It is expected that
>
> A%*%ginv(t(A)%*%solve(S)%*%A)%*%t(A)
>
> ### should equal to solve(S)######
>
> solve(S)
>
> ### but it is not!!!
> ### So I am wondering in such situation, do you have any suggestion?
> ### I have tried the argument "tol = sqrt(.Machine$double.eps))" in ginv(),
> but it won't help in
> ### more larger determinant matrix.
>
> --
> View this message in context: http://r.789695.n4.nabble.com/a-numeric-problem-tp3338989p3338989.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> 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.