Re: odd behaviour with 2x2 matrix floating-point multiplication

[s.jo]

Daniel Sebald wrote

On 12/15/2014 09:53 AM, Jordi Gutiérrez Hermoso wrote: > On Sun, 2014-12-14 at 23:05 -0800, s.jo wrote: >> I found some errors in my Cygwin system: > > Hardly an error, > > http://wiki.octave.org/FAQ#Why_is_this_floating_point_computation_wrong.3F > > if you are writing code that depends on checking exact zeros instead > of almost zeros, your code is buggy. Yes, your code may have a bug. > You cannot depend on exact equality except in very narrow cases. > >> Can other users check this example to see if all cygwin release has same >> problem? > > > The fact that you're seeing errors of the order of 1e-18, which is > smaller than eps(1), makes me think that your BLAS compilation is > probably storing some floats in registers, giving you a different > precision. This is probably correct. The _expression_ a*a + b*b is mathematically the same as the 1,1 element of the matrix multiplication, but that doesn't necessarily mean the machine code is compiled in exactly the same way. For example, in the matrix multiplication case it might be that there is some cumulative sum (after all, generally the matrix operation works on N pairs). I would guess though that it is possible to make the underlying code the same in the parser/interpreter somehow with a little work. The cost would likely be a small loss of efficiency given we're looking at a very fundamental feature. Is it worth it? Not sure, but it's probably a more logical thing to attempt than trying to declare equality when floats are next to one another. Dan

After Jordi's comment, I updated my cygwin system from 1.7.32
to 1.7.33-2 (ver  2014-11-13).

Now I have a perfect computation on zero entry, but a*a+b*b term still shows 
mismatch with a probability of 0.24 or so: If I try 10000 pairs of a and b,
about 2300 cases shows mismatch at the a*a+b*b.
I checked that the mismatch occurs at only last bit of hex form, say +/-eps(a*a+b*b),
as Dan pointed out.

I understand that N pairs term floating-point operation of matrix could be different
from the exact math _expression_.
But, in this case, only two terms are involved in machine code that BLAS may generate. 
I am not sure of how BLAS routine can shuffle the operation order or do anything
to show different floating-point results.
I suspect that BLAS routines scale the numbers unnecessarily, even trivial cases.

Is this a problem related to a release version of BLAS and/or cygwin?

-- jo

ps.

The follow script is to check the probability of a*a+b*b mismatch:

%% 
num_trial=10000;
count_err=[0;0];
for ii=1:num_trial
	a = randn;
	b = randn;
	C = [a b; -b a]*[a; b];
	D = [a*a+b*b;-b*a+a*b];
	err=(C ~= D);
	count_err = count_err + err;
	if count_err(1)==100
		num2hex(C(1))
		num2hex(D(1))
		num2hex(D(1)-eps(D(1)))
		num2hex(D(1)+eps(D(1)))
		break
	end
end
ii
count_err
ii/count_err

%% result at 100th mismatch 
ans = 40181cd9a6c0611d
ans = 40181cd9a6c0611e
ans = 40181cd9a6c0611d
ans = 40181cd9a6c0611f
a =    -1.96601051483154
b =     1.47070618400453
ans = 40181cd9a6c0611e
ii =                  415
count_err =

                   100
                     0
ans = 

     0.240963855421687
                     0

---

View this message in context: Re: odd behaviour with 2x2 matrix floating-point multiplication
Sent from the Octave - Maintainers mailing list archive at Nabble.com.