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.

Answer 2

Given the size of each number, the rounding error in simple operations on this number is of the order 10^-14

You have 5*10^6 numbers, hence if you are really unlucky the rounding error can become anything upto 5*10^-8.

Your observed error is of size 10^10, which is well in the expected range.

Note that the difference is not caused by the complex transpose, but by the product of the sum vs the matrix product.

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

B1 = conj(B); 
B2 = B';

isequal(B1(:),B2(:)) % This returns true

A*transpose(conj(B)) - A*B' % Hence this returns zero

sum(A.*transpose(B')) - A*B' % But this returns something like 1e-10

A similar effect occurs for non complex A and B:

N=1e6;
A = 1:N; 
B=1:N;

(N * (N + 1) * (2*N + 1))/6 % This will give exactly the right answer
A*B' 
fliplr(A)*fliplr(B)'

Note that the two lowest answers only vary a few hundred from eachother, whilst they are actually over 2000 from the correct answer. If this is a problem consider using the symbolic toolbox. That allows you to calculate with arbitrary precision.