Matrix-multiply giving two different answers

Question

Here's some simple code that shows what I've been seeing:

A = randn(1,5e6)+1i*randn(1,5e6);
B = randn(1,5e6)+1i*randn(1,5e6);

sum(A.*conj(B)) - A*B'
sum(A.*conj(B)) - mtimes(A,B')
A*B' - mtimes(A,B')

Now, the three methods shown on the bottom are supposed to do the same thing, so the answers should be zero, right? Wrong! The differences are small, though not small enough that I would consider them negligible. In addition, the error increases as the length of A and B increases.

Does anyone know what the actual difference between these methods is? I understand that there are probably shortcuts written into the code, but I would like to quantify that if possible. Does Matlab post the differences anywhere? I've looked around, but have not found anything.

Answer 1

It probably has something to do with the order in which the operations are performed. For example,

sum(A.*conj(B)) - fliplr(A)*fliplr(B)'

gives a result different than

sum(A.*conj(B)) - A*B'

Or, more strikingly,

A*B' - fliplr(A)*fliplr(B)'

gives a nonzero result, of the same order as your tests.

So my bet is that depending on the method (sum or *) Matlab internally does the operations in a different order, and that may well account for the different roundoff errors you observed.