An interesting Matlab bug...

[Jordan Rosenthal]

Hi all,

The following, I think, speaks for itself:

%-----------
>> clear;
>> Y = complex(rand(5),rand(5));
>> Y(:,:,1)* Y(:,:,1)' == Y*Y'
ans =
0 0 0 0 0
1 0 0 0 0
0 0 0 0 1
0 0 0 0 1
0 0 1 0 0
%-----------

1) The differences are on the order of machine precision.

2) After a bit of investigation it seems the problem also occurs if one
tries Y * conj(Y.') == Y * Y'.

3) If you assign Z = Y(:,:,1) or Z = Y + 0 and then do

Z*Y' == Y*Y'

then it also fails. Because Z has to be a physically new matrix, the
left side of the equal sign could potentially be calling a different
low-level function than the right side of the equals sign, for which
Matlab could be calling a more optimized function since it knows that Y
and Y' are sharing the same data. This seems likely since BLAS has a
function that is optimized for calculating Y*Y'.

If my assumption is correct, I'm not sure what Matlab should do in this
case. Fixing the problem would either mean taking out the optimization
or making the parser smarter, which may be difficult.

3) This is machine dependent. The above bug occurs on:
3a) a Linux machine with an Intel processor running Matlab R14.3 (but
works correctly in R2006a) and the Fedora Core 4 OS.
3b) a Linux machine with an Opteron 64-bit processor running Matlab
R14.3 and R2006a and the Fedora Core 3 OS.

By the way, the "fix" is the following

%----
>> clear
>> Y = complex(rand(5),rand(5))
>> Z = Y(:,:,1);
>> isequal(Z*Z', Y*Y')
ans =
1
%----

Then, the two will be equal because (if my assumptions were right) both
operations are using the same optimization.

Jordan

[Scott Seidman]

Jordan Rosenthal <j...@ll.mit.edu> wrote in news:VpXjg.1662$O5.540
@llslave.llan.ll.mit.edu:

> Then, the two will be equal because (if my assumptions were right) both
> operations are using the same optimization.
>
> Jordan
>

This is a shot in the dark, but have you turning the accelerator off?
--
Scott
Reverse name to reply

[ver...@hotmail.com]

What's more troubling is that
abs(Y(:,:,1)* Y(:,:,1)' - Y*Y') < eps(Y*Y')
also produced non-uniform answers. On the other hand
abs(Y(:,:,1)* Y(:,:,1)' - Y*Y') < 2*eps(Y*Y')
did produce only ones.
My feeling is that the error is too large and is a real bug that
Matworks needs to fix.

[Bobby Cheng]

A few release ago, MATLAB has added an optimization that A*A' returns
Hermitian result as expected theoretically.

Ideally, X*Y' will also return a Hermitian result if X == Y. Unfortunately,
this is no fast way to verify this property in general. I hope MATLAB will
get there eventually.

For your case, it will be perfect for MATLAB to extract Y(:,:,1) once and
return a Hermitian result.

I think of this as an enhancement rather than a bug.

---Bob.

<ver...@hotmail.com> wrote in message
news:1150307334.2...@u72g2000cwu.googlegroups.com...