How are vector - vector products calculated?
1 view (last 30 days)
Show older comments
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.
0 Comments
Accepted Answer
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)
See Also
Categories
Find more on Logical in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!