Precision in R

I'm using:
"solve(a, b, tol, LINPACK = FALSE, ...)"
Therefore,tol ==.Machine$double.eps

.Machine$double.eps

[1] 2.220446e-16

It explains why 'solve' works for 'tcp1' but not for 'tcp2':

eigen(tcp1)$values
[1]  5.208856e+09  2.585816e+08 -3.660671e-06

-3.660671e-06/5.208856e+09
[1] -7.027783e-16

eigen(tcp2)$values
[1]  5.208856e+09  2.585816e+08 -6.416393e-08
-6.416393e-08/5.208856e+09
[1] -1.231824e-17

My question would be why 'tcrossprod' and 'crossprod' leads to such
difference?

tcp1-tcp2

Marlon.

-----Original Message-----
From: David Winsemius [mailto:dwinsemius at comcast.net]
Sent: Thursday, 15 January 2009 6:04 p.m.
To: Charles C. Berry
Cc: dos Reis, Marlon; r-help at r-project.org
Subject: Re: [R] Precision in R

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}}