[R] Precision in R

[dos Reis, Marlon]

Thanks Nath,

I get the following results:
system.time(tcp1 <- tcrossprod(WB))
   user  system elapsed 
      0       0       0 
> system.time(tcp2 <- crossprod(t(WB)))
   user  system elapsed 
      0       0       0 
> all.equal(tcp1, tcp2,check.attributes=FALSE)
[1] TRUE 

But, when I try the inverse: 

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
>


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) : 
  s


> tcp1
           [,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




-----Original Message-----
From: Nathan S. Watson-Haigh [mailto:nathan.watson-haigh at csiro.au] 
Sent: Thursday, 15 January 2009 3:25 p.m.
To: dos Reis, Marlon
Cc: r-help at r-project.org
Subject: Re: [R] Precision in R

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

I've recently been working with cross products etc.

You should try the following function:
tcp1 <- tcrossprod(WB)
or
tcp2 <- crossprod(t(WB))

Both should be the same (check for equality accounting for some floating
point imprecision):
all.equal(tcp1, tcp2, check.attributes=FALSE)

You may wish to time how long it takes to do each, since in my recent
experience doing
crossprod(t(WB)) would be faster:
system.time(tcp1 <- tcrossprod(WB))
system.time(tcp2 <- crossprod(t(WB)))
all.equal(tcp1, tcp2, check.attributes=FALSE)

HTH,
Nath


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.


- --
- --------------------------------------------------------
Dr. Nathan S. Watson-Haigh
OCE Post Doctoral Fellow
CSIRO Livestock Industries
Queensland Bioscience Precinct
St Lucia, QLD 4067
Australia

Tel: +61 (0)7 3214 2922
Fax: +61 (0)7 3214 2900
Web: http://www.csiro.au/people/Nathan.Watson-Haigh.html
- --------------------------------------------------------

-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.9 (MingW32)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org

iEYEARECAAYFAklunowACgkQ9gTv6QYzVL7/AwCfcvoeS7QXxa/LCOQQhBrE+JHc
/qoAn24mXmd6fvhKdfiOavzbV1esBycu
=1WL+
-----END PGP SIGNATURE-----
=======================================================================
Attention: The information contained in this message and...{{dropped:12}}