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

[Prof Brian Ripley]

On Mon, 30 Jul 2007, Moshe Olshansky wrote:

> After multiplication by 10 you get 6*8 = 48 - the
> result is an exact machine number so there is no
> roundoff, while 0.6*0.8 = 0.48, where neither of the 3
> numbers (0.6, 0.8, 0.48) is an exact machine mumber.
> However, (-0.6)*0.8 should be equal EXACTLY to
> -(0.6*0.8), and in fact you get that sum(ev1*ev2) is
> exactly 0.
> What is strange is that you are not getting this
> result from ev1 %*% ev2. This means that either %^%
> uses some non-straightforward algorithm or it somehow
> sets the rounding control to something different from
> "round to nearest". In the later case (-0.6) does not
> necessarily equal to -(0.6) and the rounding after
> multiplication is not necessarily symetric.

Mr Olshansky seems unaware of the effects of extended-precision 
intermediate arithmetic on ix86 CPUs.

sum() does use a higher-precision accumulator (where available, including 
on Windows), but ev1*ev2 is done in R and so stored to basic precision. 
OTOH, %*% (sic) calls the BLAS routine dgemm and hence may accumulate in 
80-bit floating-point registers.  What result you get will depend on what 
compiler, compiler flags and BLAS is in use, but with the default 
reference BLAS it is very likely that some of the intermediate results are 
stored in FP registers to extended precision.

It is a simple experiment to confirm this: recompile the BLAS with 
-fforce-store and you do get 0 (at least on my Windows build system).

Let's see less speculation and more homework in future.


>
> Regards,
>
> Moshe.
>
> --- Talbot Katz <topkatz at msn.com> wrote:
>
>> 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
>>
>>
>>
>>> From: "Charles C. Berry" <cberry at tajo.ucsd.edu>
>>> To: Talbot Katz <topkatz at msn.com>
>>> CC: r-help at stat.math.ethz.ch
>>> Subject: Re: [R] Matrix Multiplication,
>> Floating-Point, etc.
>>> Date: Mon, 30 Jul 2007 09:27:42 -0700
>>>
>>>
>>>
>>> 7.31 Why doesn't R think these numbers are equal?
>>>
>>> On Fri, 27 Jul 2007, Talbot Katz wrote:
>>>
>>>> Hi.
>>>>
>>>> I recently tried the following in R 2.5.1 on
>> Windows XP:
>>>>
>>>>> ev2<-c(0.8,-0.6)
>>>>> ev1<-c(0.6,0.8)
>>>>> ev1%*%ev2
>>>>              [,1]
>>>> [1,] -2.664427e-17
>>>>> sum(ev1*ev2)
>>>> [1] 0
>>>>>
>>>>
>>>> (I got the same result with R 2.4.1 on a different
>> Windows XP machine.)
>>>>
>>>> I expect this issue is very familiar and probably
>> has been discussed in
>>>> this
>>>> forum before.  Can someone please point me to some
>> documentation or
>>>> discussion about this?  Is there some standard way
>> to get the "correct"
>>>> answer from %*%?
>>>>
>>>> Thanks!
>>>>
>>>> --  TMK  --
>>>> 212-460-5430	home
>>>> 917-656-5351	cell

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595