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}}
-----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-----
Yes, computing WB.%*%t(WB) may be the problem, by either method. if the goal is to compute the inverse of WB%*%t(WB), you should consider computing the singular value or QR decomposition for the matrix WB. If WB = Q%*%R, where Q is orthogonal, then WB %*% t(WB) =R %*%t(R), and the inverse of WB%*%t(WB) is inv.t(R)%*% inv.R. computing (WB) %*% t(WB) squares the condition number of the matrix. This is similar to the loss of precision that occurs when you compute the variance via mean(X^2)-mean(X)^2. albyn Quoting "dos Reis, Marlon" <Marlon.dosReis at agresearch.co.nz>:> 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. > >
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