# How to correctly vectorize?

1 回表示 (過去 30 日間)
Csaba 2018 年 12 月 20 日

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 件のコメント3 件の古いコメントを表示3 件の古いコメントを非表示
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 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 件のコメント4 件の古いコメントを表示4 件の古いコメントを非表示
Jan 2018 年 12 月 20 日

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

サインインしてコメントする。

### カテゴリ

Help Center および File ExchangeLoops and Conditional Statements についてさらに検索

### Community Treasure Hunt

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

Start Hunting!

Translated by