How to correctly vectorize?

5 views (last 30 days)
Csaba
Csaba on 20 Dec 2018
Edited: Jan on 21 Dec 2018
Hi,
I have this code (Matlab 2018b):
m=5;
f=[.1,-3,.05,70.56,110.32456];
c=zeros(m);
M=10;
for i=1:M
for j=1:m
fj=f(j);
for k=j:m
c(j,k)=c(j,k)+fj*f(k);
c(k,j)=c(j,k);
end
end
end
cc=M*(f'*f);
c-cc
If M=1 (or =2) the result is all zeros. If M=10, the result is not all zeros, but some. If M=100, the result is not zeros at all. I have plenty of this type of code and want to accelerate with vectorization, but I am confused about the results.
What is the correct vectorization of these kind of for loops? Why it is not zero all the times? I migt imagine that the result is around the minimum number of representation but here the difference is 1e-12 - 1e-17. It seems to me way too high.
So what should I do? Which is correct, vectorized or for loop? With for loops it works correctly.
Csaba
  4 Comments
Jan
Jan on 20 Dec 2018
@Luna: I did not understand it directly also. Csaba does not expect the elements of c and cc to be 0, but the difference c-cc.
Csaba
Csaba on 20 Dec 2018
@Luna: The "result" means - in my opinion - the end of the script, so c-cc. It should be zero.
@Stephen Cobeldick: I am summarizing floating point numbers in both cases. In principle the difference should be zero. The difference is however much bigger than the difference which comes from the floating point number representation. That is my problem.

Sign in to comment.

Accepted Answer

Jan
Jan on 20 Dec 2018
Edited: Jan on 20 Dec 2018
The differences between the loop and the linear algebra implementation have the expected range. The matrix multiplication uses highly optimized BLAS routines. The dot products contain a sum and summing is numerical instable in general, see e.g. https://www.mathworks.com/matlabcentral/fileexchange/26800-xsum
If you display the relative errors, you see that the deviations are in the magnitude of eps:
(c - cc) ./ c
This is the typical dimension of errors, e.g. caused by calculating the values one time with floating point commands and the other time with SSE/AVX. Both results are correct.
What do you consider as correct result of:
1 + 1e-17 - 1
  5 Comments
Csaba
Csaba on 20 Dec 2018
Hi,
Just to make clear. I am aware that floating point operations are not algebra. What I was surprised about was that I had too big error (1e-11 vs. 1e-17 for eps).
>> fprintf('%.40f\n',0.01)
0.0100000000000000002081668171172168513294
>> fprintf('%.40f\n',10*0.01)
0.1000000000000000055511151231257827021182
>> fprintf('%.40f\n',0.01+0.01+0.01+0.01+0.01+0.01+0.01+0.01+0.01+0.01)
0.0999999999999999916733273153113259468228
If the eps~1e-17, printing out 40 decimals has no meaning. I wonder however, where the hell Matlab gets the remaining 23(or 22) characters??
"Repeated addition (or any operation) can accumulate a lot of error. Best to avoid."
I agree. But sometimes you cannot avoid it. The example was a very simple example just to show what I mean. In my program there is not just a "repeated" addition, I just simplified the problem. You can imagine a matrix multiplication as well with for loops and vectorized, the question remains the same. How Matlab treats hundreds of additions in matrix multiplication vs. addition in for loops. Which is more accurate?
Jan
Jan on 20 Dec 2018
Edited: Jan on 21 Dec 2018
By the way, eps is 2.2e-16.
0.3 - 0.2 - 0.1
This is not 0.0 also.
There are several methods for creating a sum with reduced rounding effects: https://www.mathworks.com/matlabcentral/fileexchange/26800-xsum

Sign in to comment.

More Answers (0)

Categories

Find more on Creating and Concatenating Matrices 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!