How to correctly vectorize?

2018 12 20
1 个回答
回答已采纳
更新时间：2018 12 21
3 次查看（30 天）

0 个投票

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更早的评论 隐藏 2更早的评论

Luna
Luna 2018-12-20
You are adding values in each iteration why do you expect to get all zeros?
c(j,k)=c(j,k)+fj*f(k)
Could you please explain verbally what is your code supposed to do with your f vector? And what is your desired output?
Stephen23
Stephen23 2018-12-20
编辑：Stephen23 2018-12-20
"Why it is not zero all the times?"
Because you are summing floating point numbers.
"Which is correct, vectorized or for loop?"
Why do you think one of them is not "correct" ?
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

3 个投票

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更早的评论 隐藏 3更早的评论

Csaba
Csaba 2018-12-20
编辑：Csaba 2018-12-20
Dear Jan,
Your explanation is not fully correct. What I am doing here is i.e. for c(1,1), that in one case I am addig ten times 0.01 which is 0.1, second time I am multiplying 10*0.01. There is no magnitude of difference between the numbers when I add them, actually they are exactly the same. Try this code:
10*0.01-0.01-0.01-0.01-0.01-0.01-0.01-0.01-0.01-0.01-0.01
and the result is:
ans =
1.0408e-17
Isn't it funny???
Help says:
"Use double-precision to store values greater than approximately 3.4 x 10e38 or less than approximately -3.4 x 10e38. For numbers that lie between these two limits, you can use either double- or single-precision, but single requires less memory."
and
"The range for double is:
-1.79769e+308 to -2.22507e-308 and
2.22507e-308 to 1.79769e+308"
Csaba
Csaba
Csaba 2018-12-20
编辑：Csaba 2018-12-20
@Jan:
I red the help for "eps", and have to correct my previous comment. Yes, for double the difference between two floating point numbers is the same order of magnitude as I mentioned previously as the difference (around 1e-17).
But I still have 1e-13 or even 1e-11 differences between adding and multiplying (in the case of M=10, but more pronounced at M=100).
I would like to know which method is better (meaning closer to the theoretically correct value) when I do vectorization.
Stephen23
Stephen23 2018-12-20
编辑：Stephen23 2018-12-20
"...in one case I am addig ten times 0.01 which is 0.1, second time I am multiplying 10*0.01"
Forget about algebra, in floating point maths those operations are NOT equivalent. Ten times 0.01 is certainly NOT 0.1, because neither of those values can be represented exactly using floating point binary numbers:
>> 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
Floating point operations are not commutative or associative and are not the same as symbolic algebra, so I don't see any reason why your two operations should be the same.
"I would like to know which method is better (meaning closer to the theoretically correct value) when I do vectorization."
Repeated addition (or any operation) can accumulate a lot of error. Best to avoid.
PS: you might like to try this:
https://www.mathworks.com/matlabcentral/fileexchange/22239-num2strexact-exact-version-of-num2str
This is well worth reading as well:
https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.htmlhttp://faculty.tarleton.edu/agapie/documents/cs_343_arch/papers/1991_Goldberg_FloatingPoint.pdf
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 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Logical 的更多信息

标签

提问:

Csaba
Csaba
2018-12-20

编辑：

Jan
Jan
2018-12-21

Community Treasure Hunt

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

Start Hunting!

Translated by