How are vector - vector products calculated?

1 view (last 30 days)
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.

Accepted Answer

kjetil87
kjetil87 on 24 Aug 2013
Edited: kjetil87 on 24 Aug 2013
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);

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!