[R] Matrix Multiplication, Floating-Point, etc.

Daniel Nordlund res90sx5 at verizon.net
Tue Jul 31 04:05:58 CEST 2007


> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch [mailto:r-help-bounces at stat.math.ethz.ch]
> On Behalf Of Talbot Katz
> Sent: Monday, July 30, 2007 10:55 AM
> To: cberry at tajo.ucsd.edu
> Cc: HDoran at air.org; r-help at stat.math.ethz.ch
> Subject: Re: [R] Matrix Multiplication, Floating-Point, etc.
> 
> Thank you for responding!
> 
> I realize that floating point operations are often inexact, and indeed, the
> difference between the two answers is within the all.equal tolerance, as
> mentioned in FAQ 7.31 (cited by Charles):
> 
> >(as.numeric(ev1%*%ev2))==(sum(ev1*ev2))
> [1] FALSE
> >all.equal((as.numeric(ev1%*%ev2)),(sum(ev1*ev2)))
> [1] TRUE
> >
> 
> I suppose that's good enough for numerical computation.  But I was still
> surprised to see that matrix multiplication (ev1%*%ev2) doesn't give the
> exact right answer, whereas sum(ev1*ev2) does give the exact answer.  I
> would've expected them to perform the same two multiplications and one
> addition.  But I guess that's not the case.
> 
> However, I did find that if I multiplied the two vectors by 10, making the
> entries integers (although the class was still "numeric" rather than
> "integer"), both computations gave equal answers of 0:
> 
> >xf1<-10*ev1
> >xf2<-10*ev2
> >(as.numeric(xf1%*%xf2))==(sum(xf1*xf2))
> [1] TRUE
> >
> 
> Perhaps the moral of the story is that one should exercise caution and keep
> track of significant digits.
> 
> --  TMK  --
> 212-460-5430	home
> 917-656-5351	cell
> 
There may other issues involved here besides R version, floating point precision, and OS version.  On my WinXP system running R-2.5.1 binary from CRAN, I get what you expected:

> ev2<-c(0.8,-0.6)
> ev1<-c(0.6,0.8)
> ev1%*%ev2
     [,1]
[1,]    0
>

There could be differences in OS release, service packs installed, cpu, etc.  But the moral you draw is probably a reasonable one.  

Dan

Daniel Nordlund
Bothell, WA



More information about the R-help mailing list