Re: 2x2 orthogonal matrix multiplication error

[Daniel J Sebald]

On 12/14/2014 11:04 PM, s.jo wrote:
> Hello all,
>
> In octave(3.8.2), 2x2 orthogonal matrix M shows numerical error after
> multiplied by own transpose.
>
> Say,
>     M*M' is not equal to eye(2)!
>
> There are errors in off-diagonal entries.
> But, if we try it with element-wise multiplication, it works well.
>
> Say,
> Mt = M';
> [Mt(1, :) * M(:, 1), Mt(1, :) * M(:, 2); Mt(2, :) * M(:, 1), Mt(2, :) * M(:,
> 2)] is equal to eye(2)!!
>
> I check this also with matlab, but it works well.
> I suspect that octave has some unreliable routine with
> matrix-multiplication.
> Or, is this only problem with my cygwin machine?

This sounds like something that came up a couple months back, displaying 
zeros, or equality of two numbers.  The issue is likely that in matrix 
multiplication case there is addition.  In the result you used above, 
assuming orthogonality, there is no addition.  It's addition which is 
the dodgy numerical operation at the register level.  Floating point 
multiplication is rather accurate and a shift of the exponent.

[snip]

> Can somebody help this how it figure out?

From the recent discussion about the related issue, I had found in the 
(I think it was standard) C library some functions that will identify 
two floating point numbers that are next to one another.  With that, I 
think there is the possibility of modifying the equality operation to 
declare floating point numbers next to one another as equal.

That would "solve" some of these "surprising" results where the user is 
expecting the result of different operations to be equal because the 
associated mathematics in theory says they should be equal.  On the 
other hand, it could create equally surprising results like a violation 
of the transitive property, i.e., say 'a' is a floating point number 
next to 'b', and 'b' is a floating point number next to 'c', then:

a = b -> true
b = c -> true
a = c -> false

i.e., it could be that 'a' and 'c' are not floating point numbers next 
to one another.

In some sense, I don't know how important it is, aside from 
"compatibility" with Matlab.  If one is using Octave on practical 
problems, there's going to be a random noise component so "equality" 
might mean some test of norm distance, i.e., |x - y| < alpha, or something.

Could you do more analysis to identify where the equality is failing on 
Cygwin?  Subtract individual elements of the two matrices and see how 
far they deviate from zero?

Dan



> Many thanks,
>
> jo