[R] Precision in R
David Winsemius
dwinsemius at comcast.net
Thu Jan 15 14:46:01 CET 2009
I'm not getting that difference:
> .Machine$double.eps
[1] 2.220446e-16
# so I don't think that explains the different behavior.
> WB<-as.matrix(read.table(file.choose()))
> tcp1 <- tcrossprod(WB)
> tcp2 <- crossprod(t(WB))
> tcp1
[,1] [,2] [,3]
[1,] 1916061939 2281366606 678696067
[2,] 2281366606 3098975504 1092911209
[3,] 678696067 1092911209 452399849
>>> WBtWB
>> [,1] [,2] [,3]
>> [1,] 1916061939 2281366606 678696067
>> [2,] 2281366606 3098975504 1092911209
>> [3,] 678696067 1092911209 452399849
> tcp2
[,1] [,2] [,3]
[1,] 1916061939 2281366606 678696067
[2,] 2281366606 3098975504 1092911209
[3,] 678696067 1092911209 452399849
> solve(tcp1)
Error in solve.default(tcp1) :
system is computationally singular: reciprocal condition number =
9.60696e-17
> solve(tcp2)
Error in solve.default(tcp2) :
system is computationally singular: reciprocal condition number =
9.60696e-17
That is somewhat interesting since yesterday my machine solved my
input version of your:
>>> WBtWB
>> [,1] [,2] [,3]
>> [1,] 1916061939 2281366606 678696067
>> [2,] 2281366606 3098975504 1092911209
>> [3,] 678696067 1092911209 452399849
I assume that despite those matrices being displayed as the same they
are represented differently in the machine.
Berry wrote:
> 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.
That would seem to be a possibility. You are using an out-of-date
version which may limit people's interest in investigating the problem.
--
David Winsemius
On Jan 15, 2009, at 12:47 AM, dos Reis, Marlon wrote:
> Hi,
> I attached the files I'm using, it may help.
> I'm using Windows XP
>> sessionInfo()
> R version 2.6.0 (2007-10-03)
> i386-pc-mingw32
> locale: LC_COLLATE=English_New Zealand.1252;LC_CTYPE=English_New
> Zealand.1252;LC_MONETARY=English_New
> Zealand.1252;LC_NUMERIC=C;LC_TIME=English_New Zealand.1252
> attached base packages:
> [1] stats graphics grDevices utils datasets methods base
>
> Try for example:
> WB<-as.matrix(read.table("WB.dat"))
> tcp1 <- tcrossprod(WB)
> tcp2 <- crossprod(t(WB))
>
> solve(tcp1)
> [,1] [,2] [,3]
> [1,] -41692.80 58330.89 -78368.17
> [2,] 58330.89 -81608.66 109642.09
> [3,] -78368.17 109642.09 -147305.32
>
> solve(tcp2)
> Error in solve.default(tcp2) :
> system is computationally singular: reciprocal condition number =
> 2.17737e-17
>
> Marlon.
>
> -----Original Message-----
> From: Charles C. Berry [mailto:cberry at tajo.ucsd.edu]
> Sent: Thursday, 15 January 2009 5:16 p.m.
> To: dos Reis, Marlon
> Cc: r-help at r-project.org
> Subject: Re: [R] Precision in R
>
>
> 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
>
>
>
> =
> ======================================================================
> Attention: The information contained in this message and/or
> attachments
> from AgResearch Limited is intended only for the persons or entities
> to which it is addressed and may contain confidential and/or
> privileged
> material. Any review, retransmission, dissemination or other use of,
> or
> taking of any action in reliance upon, this information by persons or
> entities other than the intended recipients is prohibited by
> AgResearch
> Limited. If you have received this message in error, please notify the
> sender immediately.
> =
> ======================================================================
> <WB.dat>
More information about the R-help
mailing list