[R] Precision in R

[Charles C. Berry]

Marlon,

Are you using a current version of R? sessionInfo()?

It would help if you had something we could _fully_ reproduce.

Taking the _printed_ values you have below (WBtWB) and adding or 
subtracting what you have printed as the difference of the two visually 
equal matrices ( say Delta ) , I am able to run

 	solve( dat3 )
 	solve( WBtWB + Delta )
 	solve( WBtWB - Delta )
 	solve( WBtWB + 2*Delta )
 	solve( WBtWB - 2*Delta )

and get the results to agree to 3 significant digits. And perturbing 
things even more I still get solve() to return a value:

> for ( i in 1:1000 ) solve(WBtWB - tcrossprod(rnorm(3)))
> for ( i in 1:1000 ) solve(WBtWB + tcrossprod(rnorm(3)))


And I cannot get condition numbers anything like what you report:

> range(replicate( 10000, 1/kappa(dat3-tcrossprod(matrix(rnorm(9),3)))))
[1] 5.917764e-11 3.350445e-09
>


So I am very curious that you got the results that you print below.

I suppose that it is possible that the difference between what you report 
and what I see lies in the numerical libraries (LINPACK/LAPACK) that R 
calls upon.

This was done on a windows XP PC. Here is my sessionInfo()

> sessionInfo()
R version 2.8.1 Patched (2008-12-22 r47296)
i386-pc-mingw32

locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United 
States.1252;LC_MONETARY=English_United 
States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base
>

HTH,

Chuck

On Thu, 15 Jan 2009, dos Reis, Marlon wrote:

> Dear All,
> I'm preparing a simple algorithm for matrix multiplication for a
> specific purpose, but I'm getting some unexpected results.
> If anyone could give a clue, I would really appreciate.
> Basically what I want to do is a simple matrix multiplication:
> (WB) %*% t(WB).
> The WB is in the disk so I compared to approaches:
> -	Load 'WB' using 'read.table' (put it in WB.tmp) and then to the
> simple matrix multiplication
> WB.tmp%*%t(WB.tmp)
>
> -	Scan each row of WB and do the cross products 'sum(WB.i*WB.i)'
> and 'sum(WB.i*WB.j)', which proper arrangement leads to WBtWB.
>
> Comparing these two matrices, I get the very similar values, however
> when I tried their inverse, WBtWB leads to a singular system.
> I've tried different tests and my conclusion is that  my precision
> problem is related to cross products 'sum(WB.i*WB.i)' and
> 'sum(WB.i*WB.j)'.
> Does it makes sense?
> Thanks,
> Marlon
>
>
>> WB.tmp%*%t(WB.tmp)
>           WB.i       WB.i       WB.i
> WB.i 1916061939 2281366606  678696067
> WB.i 2281366606 3098975504 1092911209
> WB.i  678696067 1092911209  452399849
>
>> WBtWB
>           [,1]       [,2]       [,3]
> [1,] 1916061939 2281366606  678696067
> [2,] 2281366606 3098975504 1092911209
> [3,]  678696067 1092911209  452399849
>
>
>> WBtWB-WB.tmp%*%t(WB.tmp)
>             WB.i          WB.i          WB.i
> WB.i 2.861023e-06  4.768372e-07  4.768372e-07
> WB.i 4.768372e-07  3.814697e-06 -2.622604e-06
> WB.i 4.768372e-07 -2.622604e-06  5.960464e-08
>
>> solve(WB.tmp%*%t(WB.tmp))
>          WB.i      WB.i       WB.i
> WB.i -41692.80  58330.89  -78368.17
> WB.i  58330.89 -81608.66  109642.09
> WB.i -78368.17 109642.09 -147305.32
>
>> solve(WBtWB)
> Error in solve.default(WBtWB) :
>  system is computationally singular: reciprocal condition number =
> 2.17737e-17
>
>
>
>
>     WB.tmp<-NULL
>     WBtWB<-matrix(NA,n,n)
>      for (i in 1:n)
>      {
>       setwd(Home.dir)
>       WB.i<-scan("WB.dat", skip = (i-1), nlines = 1)
>       WB.tmp<-rbind(WB.tmp,WB.i)
>       WBtWB[i,i]<-sum(WB.i*WB.i)
>       if (i<n)
>        {
>          for (j in (i+1):n)
>           {
>              setwd(Home.dir)
>              WB.j<-scan("WB.dat", skip = (j-1), nlines = 1)
>              WBtWB[i,j]<-sum(WB.i*WB.j)
>              WBtWB[j,i]<-sum(WB.i*WB.j)
>           }
>         }
>       }
>
>
> =======================================================================
> Attention: The information contained in this message and...{{dropped:15}}
>
> ______________________________________________
> 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.
>

Charles C. Berry                            (858) 534-2098
                                             Dept of Family/Preventive Medicine
E mailto:cberry at tajo.ucsd.edu	            UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901