[R] Precision in R

David Winsemius dwinsemius at comcast.net
Thu Jan 15 06:04:12 CET 2009


I am seeing different behavior than don Reis on my installation as well:
mtx is the same as his WBtWB

 > mtx <- matrix(c(1916061939, 2281366606,  678696067, 2281366606,  
3098975504, 1092911209, 678696067, 1092911209,  452399849), ncol=3)
 >
 > mtx
            [,1]       [,2]       [,3]
[1,] 1916061939 2281366606  678696067
[2,] 2281366606 3098975504 1092911209
[3,]  678696067 1092911209  452399849
 > eigen(mtx)
$values
[1]  5.208856e+09  2.585816e+08 -4.276959e-01

$vectors
            [,1]       [,2]       [,3]
[1,] -0.5855545 -0.7092633  0.3925195
[2,] -0.7678140  0.3299775 -0.5491599
[3,] -0.2599763  0.6229449  0.7378021

 > rcond(mtx)
[1] 5.33209e-11

Despite a very ill-conditioned matrix, solve still proceeds
 >
 > solve(mtx)
            [,1]       [,2]       [,3]
[1,] -0.3602361  0.5039933 -0.6771204
[2,]  0.5039933 -0.7051189  0.9473348
[3,] -0.6771204  0.9473348 -1.2727543
 >
 > sessionInfo()
R version 2.8.1 Patched (2009-01-07 r47515)
i386-apple-darwin9.6.0

locale:
en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8

snipped package info

-- 
David Winsemius
On Jan 14, 2009, at 11:15 PM, Charles C. Berry wrote:

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




More information about the R-help mailing list