How are vector - vector products calculated?

When I run the following code:

n = 10000;

for i = 1:n
    A(i,1) = rand/100;
    B(i,1) = rand/100;
end

D1 = A'*B;
D2 = dot(A,B);

fprintf('%16.15f\n', D1)
fprintf('%16.15f\n', D2)

The least significant digit is sometimes off by one or two[1], which to me indicates that the underlying mechanism is different[2]. Normally such a difference is absolutely negligible, but I was comparing two iterative methods that should give the same answer, but after a number of iterations the differences in the least significant digit propagated to larger decimals, and the final results were different due to this very issue.

I am aware of the non-commutative nature of floating point arithmetic, but this means that the * operator does not compute the vector product in sequential order, whereas dot() does. Can anybody shed some light on why vector-vector multiplication and dot() produce slightly "different" answers?

[1] IEEE 754 specifies 15–17 significant decimal digits precision, so this sort of fine-grained comparison should still be valid.

[2] To quote the online documentation: C = dot(A,B) returns the scalar product of the vectors A and B. A and B must be vectors of the same length. When A and B are both column vectors, dot(A,B) is the same as A'*B.

Answer 1

[kjetil87]

It is different, if you type

type dot

You will see that the actual computation done in dot after all the dimension checks are:

D2=sum(conj(A).*B);

The sum .* and * functions are built in functions so it is hard to give a better answer then; yes the algorithms are different (see mtimes , sum and times)

But, it seems that * (or mtimes) is multithreaded. So it probably splits the vector into chunks, multiplies and sums them, then sum them all together in the end. You can confirm this by typing

test1=A'*B;
LastN=maxNumCompThreads(1);
test2=A'*B;

isequal(test1,test2)

ans =

       0

maxNumCompThreads does not persist through sessions, but to avoid having to restart matlab to enable multiple threads use

maxNumCompThreads(LastN);

[lyreco]

Beautiful answer. Thanks