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.
