フィルターのクリア

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 件のコメント
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 件のコメント
Csaba
Csaba 2018 年 12 月 20 日
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.
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 件)

Community Treasure Hunt

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

Start Hunting!

Translated by