How to correctly vectorize?

1 次查看（过去 30 天）
Csaba
Csaba 2018-12-20
编辑： Jan 2018-12-21
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 个评论
显示 2更早的评论
Jan
Jan 2018-12-20
@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 2018-12-20
@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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Jan
Jan 2018-12-20
编辑：Jan 2018-12-20
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 个评论
显示 3更早的评论
Csaba
Csaba 2018-12-20
@Stephen Cobeldick
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 2018-12-20
编辑：Jan 2018-12-21
By the way, eps is 2.2e-16.
Maybe this thread clarifies some details: https://www.mathworks.com/matlabcentral/answers/57444-faq-why-is-0-3-0-2-0-1-not-equal-to-zero
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

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Special Values 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by