Precision in R
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}}