Numerical inconsistency between loop and builtin function
I'm trying to compute the sum of an array of random numbers. But there seems to be an inconcistancy between the results when I do it one element at a time and when I use the built-in function. Furthermore, the error seems to increase when I decrease the data precision.
import torch
columns = 43*22
rows    = 44
torch.manual_seed(0)
array = torch.rand([rows,columns], dtype = torch.float64)
array_sum = 0
for i in range(rows):
    for j in range(columns):
        array_sum += array[i, j]

torch.abs(array_sum - array.sum())

results in:
tensor(3.6380e-10, dtype=torch.float64)

using dtype = torch.float32 results in:
tensor(0.1426)

using dtype = torch.float16 results in (a whooping!):
tensor(18784., dtype=torch.float16)

I find it hard to believe no one has ever asked about it. Yet, I haven't found a similar question in SO.
Can anyone please help me find some explanation or the source of this error?

The first mistake is this:
you should change the summation line to
array_sum += float(array[i, j])

For float64 this causes no problems, for the other values it is a problem, the explenation will follow.
To start with: when doing floating point arithmetic, you should always keep in mind that there are small mistakes due to rounding errors. The most simple way to see this is in a python shell:
&gt;&gt;&gt; .1+.1+.1-.3
5.551115123125783e-17

But how do you take these errors into account?
When summing n positive integers to a total tot, the analysis is fairly simple and it the rule is:
error(tot) &lt; tot * n * machine_epsilon
Where the factor n is usually a gross over-estimation and the machine_epsilon is dependant on the type (representation size) of floating point-number.
And is approximatly:
float64: 2*10^-16
float32: 1*10^-7
float16: 1*10^-3

And one would generally expect as an error approximately within a reasonable factor of tot*machine_epsilon.
And for my tests with float16 we get (with always +-40000 variables summing to a total of +- 20000):
error(float64) = 3*10^-10 ≈ 80* 20000 * 2*10^-16
error(float32) = 1*10^-1  ≈ 50* 20000 * 1*10^-7

which is acceptable.
Then there is another problem with the float 16. There is the machine epsilon = 1e-4 and you can see the problem with
&gt;&gt;&gt; ar = torch.ones([1], dtype=float16)
&gt;&gt;&gt; ar
tensor([2048.], dtype=torch.float16)
&gt;&gt;&gt; ar[0] += .5
&gt;&gt;&gt; ar
tensor([2048.], dtype=torch.float16)

Here the problem is that when the value 2048 is reached, the value is not precise enough to be able to add a value 1 or less. More specifically: with a float16 you can 'represent' the value 2048, and you can represent the value 2050, but nothing in between because it has too little bits for that precision. By keeping the sum in a float64 variable, you overcome this problem. Fixing this we get for float16:
error(float16) = 16  ≈ 8* 20000 * 1*10^-4

Which is large, but acceptable as a value relative to 20000 represented in float16.
If you ask yourself, which of the two methods is 'right' then the answer is none of the two, they are both approximations with the same precision, but a different error.
But as you probably guessed using the sum() method is faster, better and more reliable.
