Re: ML matrices multiplication accuracy

ellieandrogerxyzzy@xxxxxxxxxxxxxxxxxxxxxx (Roger Stafford) writes:

In article <ef5b370.7@xxxxxxxxxxxxxxxxxxxxxxx>, "Yi Cao"
<y.cao@xxxxxxxxxxxxxxx> wrote:

Which matlab version did you use? On my machine (matlab r2007a, Intel
Core Duo Processor T2500), I got all 9 cases exactily the same:

format long eng
a = hex2num('bfd729edb7ae53dc');
b = hex2num('3fef2524f4fe4a4a');
c = hex2num('bf88abe7503157c0');
d = hex2num('bfddeefa0ebbddf4');

x=[a*c+b*d;
c*a+b*d;
a*c+d*b;
c*a+d*b;
b*d+a*c;
b*d+c*a;
d*b+a*c;
d*b+c*a;
[a,b]*[c;d]]

x =

-450.855637235782e-003
-450.855637235782e-003
-450.855637235782e-003
-450.855637235782e-003
-450.855637235782e-003
-450.855637235782e-003
-450.855637235782e-003
-450.855637235782e-003
-450.855637235782e-003

norm(diff(x))

ans =

0.00000000000000e-003

hth
Yi
---------------
Yi, I use an admittedly very ancient version 4a of matlab on a Macintosh
clone with PowerPC 604e microprocessor chip, which no doubt accounts for
the difference between our respective results. Nevertheless, it surprises
me that, regardless of the age of my matlab and computer, it could have
gotten different results from something so simple as a*c+b*d and
[a,b]*[c;d] if it used only IEEE 754 floating point CPU instructions which
accept just two floating point inputs and then round the result back down
to 53 bits precision each time. It seems to me that the procedure must be
to somehow multiply a by c and b by d, and afterwards add the two
products. What else is there to do? If the products were each rounded
first in accordance with IEEE 754 rounding specs. before the addition,
there is only one possible result. Changes in ordering here with only two
operands in each case can make no difference. I can't help thinking that
somehow a floating point instruction must have been used on my machine
that accepts three inputs, as with the "multiply-add" instructions, so as
to alter the round off situation. I know that there are such instructions
on my computer. I notice that [b,a]*[d;c] gets a different result from
[a,b]*[c;d] on my machine, which strengthens that suspicion.

Roger-

I can't speak to the PPC 604e. But the modern Intel floating point
unit has 80-bit registers and internal computation (extended
precision). So on these machines even the order of operations is not
enough to guarantee consistency, but one also needs to know at what
point values were stored in memory, at which point they become 64-bit
values. So the precise optimization that the compiler performed on a
certain set of code can introduce different rounding effects.

The SSE2 instruction set referred to earlier adds a different floating
point unit, which is capable of processing 2 64-bit floats at once.
This unit has no extended precision. It also (inexplicably) has
better performance for many operations (like int to float
conversions), so that compilers are increasingly using this unit instead
the 387 "standard" Intel floating point instructions, even for scalar
operations.

Maybe your PPC machine has similar extended precision effects?

-Peter
.