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? Marlon. -----Original Message----- From: David Winsemius [mailto:dwinsem...@comcast.net] Sent: Thursday, 15 January 2009 6:04 p.m. To: Charles C. Berry Cc: dos Reis, Marlon; r-help@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}} >> >> ______________________________________________ >> R-help@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:cbe...@tajo.ucsd.edu UC San Diego > http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego > 92093-0901 > > ______________________________________________ > R-help@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. ======================================================================= Attention: The information contained in this message and...{{dropped:12}} ______________________________________________ R-help@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.
