[R] Precision in R
Albyn Jones
jones at reed.edu
Thu Jan 15 04:01:50 CET 2009
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.
>
>
More information about the R-help
mailing list