np.sum() and sum don't agree each other

Question

A numpy array contains a list of positive and negative numbers, when do np.sum() and native sum, they don't agree with each other and return extremely small number, even if the original type is float64. Surprisingly, once turn them into float16, they seem like working.

Just wondering what causes this issue, for me , float64 is more precise than float16, there shouldn't be an issue like this.

data = np.array([-0.90569317, -1.46995447, -1.74996384, -1.8867866 , -1.87511954,
   -1.57071542, -1.28009964, -0.93220921, -0.51855901, -0.20460912,
   -0.44537475, -0.2321858 ,  0.22495071, -0.07839277, -0.19347238,
   -0.37325112, -0.19082078, -0.14362223,  0.17350959,  0.29018016,
    0.67731432,  0.40419913,  0.16290318, -0.30590038, -0.97622584,
   -1.19896056, -1.76534314, -1.90852975, -1.84276997, -1.59776178,
   -1.40631598, -1.03456112, -0.05293738,  0.15388772, -0.36582663,
   -0.30218814,  0.71072453,  0.13691746, -0.16483506, -0.15157704,
   -0.14892544,  0.20161659,  0.27904342,  0.99126422,  1.30786572,
    0.55003734,  0.78496944,  0.1353265 , -0.47825463, -0.91682991,
   -1.62692942, -1.77435859, -1.76799475, -1.68844663, -1.37343609,
   -1.03562176, -0.11074234,  0.27798278, -0.12877325, -0.03172455,
    0.77754495,  0.47685307, -0.012633  , -0.01952717, -0.02376973,
    0.06108159,  0.37980437,  1.32324502,  1.50355408,  1.24634851,
    1.48021997,  1.12914762, -0.04869481, -0.74023309, -1.32464658,
   -1.56541221, -1.5574574 , -1.54632066, -1.43866555, -0.95448269,
   -0.12612164,  0.28593759,  0.07911249,  0.00857983,  0.78762104,
    0.55693151,  0.09024923, -0.03596711,  0.41692683,  0.22017782,
    0.62110032,  1.18111905,  1.67378705,  1.31369925,  1.48817478,
    1.61280016,  2.02963229,  0.41586618, -0.54454473, -1.28116028,
   -1.41055855, -1.45881774, -1.40896758, -0.74341502, -0.42097999,
    0.02714106, -0.17915372,  0.56912889,  1.50726633,  1.13922371,
    0.59670557,  0.62216096,  0.67360208,  0.83428927,  0.94353535,
    2.18554659,  2.3032778 ,  1.39748993,  0.95149016,  0.63223706,
    1.51522114,  1.68757539,  1.113238  ,  0.13691746, -0.50477067,
   -0.72591443, -1.1411556 , -1.19789992, -0.90675382, -0.41514646,
   -0.16218346,  0.06479383,  0.69057234,  0.40366881,  0.15123612,
    0.08494602,  0.23555712,  0.5770837 ,  0.79345457,  1.31475989,
    1.56188937,  1.00929512,  1.11376832,  0.9679301 ,  2.25183669,
    2.21577488,  1.60802727,  0.38775918, -0.2067304 , -0.54242345,
   -1.15547426, -1.19153607, -1.1284279 , -0.9751652 , -0.55037826,
   -0.33824995,  0.21169269, -0.00308722, -0.2629444 ,  0.03297459,
   -0.10543913, -0.12028811,  0.30237754,  0.67148079,  1.29195609,
    0.64920732,  0.48427756,  0.14752387])

For_loop accumulate, the last output become a extremely small value

c = 0

for i in data:
    c += i
    print(c)

-0.90569317
-2.37564764
-4.12561148
-6.01239808
-7.88751762
-9.45823304
-10.73833268
-11.670541889999999
-12.1891009
-12.39371002
-12.83908477
-13.07127057
-12.84631986
-12.92471263
-13.11818501
-13.49143613
-13.68225691
-13.82587914
-13.65236955
-13.36218939
-12.684875069999999
-12.280675939999998
-12.117772759999998
-12.423673139999998
-13.399898979999998
-14.598859539999998
-16.36420268
-18.272732429999998
-20.115502399999997
-21.713264179999996
-23.119580159999995
-24.154141279999994
-24.207078659999993
-24.053190939999993
-24.419017569999994
-24.721205709999992
-24.010481179999992
-23.873563719999993
-24.038398779999994
-24.189975819999994
-24.338901259999993
-24.137284669999993
-23.858241249999992
-22.86697702999999
-21.55911130999999
-21.009073969999992
-20.22410452999999
-20.08877802999999
-20.567032659999988
-21.48386256999999
-23.11079198999999
-24.885150579999987
-26.653145329999987
-28.341591959999988
-29.715028049999987
-30.750649809999988
-30.86139214999999
-30.58340936999999
-30.71218261999999
-30.74390716999999
-29.96636221999999
-29.48950914999999
-29.50214214999999
-29.52166931999999
-29.54543904999999
-29.48435745999999
-29.104553089999992
-27.78130806999999
-26.27775398999999
-25.03140547999999
-23.55118550999999
-22.422037889999988
-22.47073269999999
-23.21096578999999
-24.53561236999999
-26.10102457999999
-27.65848197999999
-29.20480263999999
-30.64346818999999
-31.59795087999999
-31.72407251999999
-31.43813492999999
-31.35902243999999
-31.35044260999999
-30.56282156999999
-30.00589005999999
-29.91564082999999
-29.95160793999999
-29.53468110999999
-29.31450328999999
-28.69340296999999
-27.51228391999999
-25.83849686999999
-24.52479761999999
-23.03662283999999
-21.42382267999999
-19.39419038999999
-18.978324209999993
-19.52286893999999
-20.80402921999999
-22.21458776999999
-23.67340550999999
-25.08237308999999
-25.82578810999999
-26.24676809999999
-26.21962703999999
-26.398780759999987
-25.829651869999985
-24.322385539999985
-23.183161829999985
-22.586456259999984
-21.964295299999986
-21.290693219999987
-20.456403949999988
-19.512868599999987
-17.327322009999985
-15.024044209999985
-13.626554279999985
-12.675064119999984
-12.042827059999984
-10.527605919999985
-8.840030529999986
-7.726792529999986
-7.5898750699999855
-8.094645739999985
-8.820560169999984
-9.961715769999984
-11.159615689999983
-12.066369509999983
-12.481515969999984
-12.643699429999984
-12.578905599999985
-11.888333259999985
-11.484664449999986
-11.333428329999986
-11.248482309999986
-11.012925189999986
-10.435841489999987
-9.642386919999987
-8.327627029999988
-6.765737659999988
-5.756442539999988
-4.6426742199999875
-3.6747441199999873
-1.422907429999987
0.792867450000013
2.400894720000013
2.7886539000000132
2.581923500000013
2.0395000500000133
0.8840257900000132
-0.3075102799999867
-1.4359381799999866
-2.411103379999987
-2.961481639999987
-3.299731589999987
-3.088038899999987
-3.091126119999987
-3.3540705199999867
-3.321095929999987
-3.426535059999987
-3.5468231699999873
-3.244445629999987
-2.572964839999987
-1.2810087499999872
-0.6318014299999871
-0.14752386999998712
1.2878587085651816e-14

math.fsum(), also return a unsense extremely small value

math.fsum(data)
-9.367506770274758e-17

sum() and np.sum() for float16, give usuall result

print(sum(data.astype(np.float16)))
print(np.sum(data, dtype = np.float16))
0.0011463165283203125
0.001146

sum() and np.sum() for float64, again give extremely small value

print(sum(data.astype(np.float64)))
print(np.sum(data, dtype = np.float64))
1.2878587085651816e-14
0.0





    

Answer 1

Firstly the correct result is 0. I know that because I've used decimals to compute the sum. So I know that, at least if we consider those numbers as their exact representation as shown in decimal (which they are probably not: they obviously come from some other computations, so are just roundings)

import decimal
s='''copyand paste of the numbers from your post'''
sum(decimal.Decimal(x) for x in s.split(','))

Returns 0. Not a very small number, but exactly 0.

So, obviously float64 works better that float16. I am not sure why you claim that the result are less strange with float16.

Sure, one could claim that from the origin, all we know that each of those numbers are within a [-5e-9, 5e-9] interval from what they are. So, there is a margin for interpretation after all. But even in worst case (all those digits are 0.4999...e-8 over the real value for example), that is at most 8.4e-7 around 0. And any way, when in doubt, if given decimals numbers, there is no reason to assume that they are not what they are.

So, long story short, the sum is 0. So data.sum() is right when saying exactly 0. fsum(data) is right when giving an answer well below the margin of error with those decimals representation.

When data.astype(np.float16).sum() (or same with sum or fsum; they give the same result, see below), are wrong when claiming that result is 0.001146. Which, even for a float16 is not the same thing as 0. But that is because the error accumulate for each addition in float16.

Secondly, there is a rounding behavior in the display only. Python floats and numpy floats do not have the same print function. Numpy rounds more when displaying.

Just try

print(np.sum(data, dtype = np.float16)-0.001146)

And you'll see the rest of your digits displaying. Proof that they were already there.

Or just

np.sum(data, dtype = np.float16)*1.0

There are two reasons why this works:

• The missing information is already there. If there were really some missing precision (0.0000003165283...) in np.sum(data, dtype = np.float16) then there would be no reason to magically find it back
• multiplying by 1.0 (which is a python's float, not a numpy one. And is converted to a numpy.float64 to perform the operation), force the result back into a numpy.float64. And numpy.float64 show more decimals when printed.

As for why numpy truncate the decimal representation when printing, there is another easy way to see it

np.float16(0.001146)*1.0

Displays yet again the same 0.0011463165283203125

In other words, in float16, those are the same things. They have the same exact binary representation.

Which we can also check another way

print(np.float16(0.001146).tobytes())
print(np.float16(0.0011463165283203125).tobytes())

Numpy chooses to print only the decimals that could change the binary representation. Among all the numbers that share the same binary representation (the discrete representation of floats means that an infinite number of reals have the same IEEE float representation), numpy chooses to display the one with the simplest decimal representation. Which makes sense.

But, well, one must be aware that the print function are not the same. And that data.sum(), even tho it is a scalar, is still a numpy.float64 (or numpy.float16, if you astyped it as such), not a python number.

As for what error to expect in float16 for 168 additions of numbers between 1 and 2 (some below 0.5, some other over 2), you can also check, for example, with your first number

d16=data.astype(np.float16)
d16[0].tobytes() ## b'?\xbb'
# Since my machine is little endian (and so is yours, I don't take lot of risk claiming that)
# the least significant byte is the `b'?'` part. If we change it by only 1, 
# we change the float representation by the smallest possible interval
# and one after '?' is chr(ord('?')+1) = @
d16[0] - np.frombuffer(b'@\xbb', dtype=np.float16)[0]

Result is 0.0004883

Note that you can get that more simply by

np.spacing(d16[0])

I just wanted to show "greasy hands", where this come from.

So, at worst, when summing that many numbers, there would be no surprise in accumulating 168 times that as an error. That is 0.08

So, 0.001146 is not even that much (of course, to have 0.08, it would require all the 168 numerical error to be maximal and in the same direction. In reality, errors goes both way and compensate. So, claiming that the error could be 0.08, would be like claiming that drawing 168 coins could give 168 tails...)

So, tl;dr

1. Real answer is 0. Your result only shows that float64 is, unsurprisingly more accurate than float16
2. All results in float16, with .sum() sum(...) of fsum(...) are identical. Only decimal printing changes.
3. Numerical error, for both fsum on float64, and all methods on float16, are way withing what to be expected when summing 168 numbers.

Incidentally, the only method giving the exact 0 result is numpy with float64... the one you were accusing :D

Answer 2

Use math.fsum() instead of np.sum()

np.sum() is for fast calculation, but math.fsum() is slower, but more precise.

math.fsum() returns an accurate floating point sum of values in the iterable. Avoids loss of precision by tracking multiple intermediate partial sums.

The algorithm’s accuracy depends on IEEE-754 arithmetic guarantees and the typical case where the rounding mode is half-even.