How to correctly vectorize?

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

[Luna]

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]

"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]

@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]

@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.

Answer 1

[Jan]

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

[Csaba]

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]

@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]

"...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]

@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]

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