2

I am summing an array using numpy, and it isn't working for float32 type. What am I doing wrong? Once I do a sum directly, and then I use numpy.sum. See the code below

import struct
import numpy as np
import matplotlib.pyplot as plt
import math
from pylab import *
xpt=128
ypt=128
zpt=256
bx1=np.zeros((xpt,ypt,zpt),dtype=float32)
bx2=np.zeros((xpt,ypt,zpt),dtype=float32)
bx3=np.zeros((xpt,ypt,zpt),dtype=float32)

bx1=bx1+1.0
bx2=bx2+1.5
bx3=bx3+2.0

dummy=0.0
for kxi in range (0,xpt) :
  for kyi in range (0,ypt) :
    for kzi in range (0,zpt) :
      dummy=dummy+(bx1[kxi,kyi,kzi]*bx1[kxi,kyi,kzi]+bx2[kxi,kyi,kzi]*bx2[kxi,kyi,kzi]+bx3[kxi,kyi,kzi]*bx3[kxi,kyi,kzi])
print(dummy)

print(np.sum(bx1**2+bx2**2+bx3**2))

Both outputs should match. This gives the output:
30408704.0
3.1323e+07

The direct sum gives the correct result, whereas np.sum is giving something wrong. However, if I use float64, then np.sum gives the correct result. What is the reason behind this?

Thanks.

user3366891
  • 21
  • 2

3 Answers3

4

It is a problem with the accuracy of float32 with a number that big. I haven't run through how 2.25 would be stored, but the minimal example

x = 2.25 * np.ones((128, 128, 256), dtype = float32)
y = 2.25 * np.ones((128, 128, 256), dtype = float64)
x.sum() # 8854642.0
y.sum() # 9437184.0
2.25 * 128 * 128 * 256 # 9437184.0

Shows that you lose accuracy, but gain it back with float64 (python's standard float).

Community
  • 1
  • 1
colcarroll
  • 3,632
  • 17
  • 25
2

Possible cause of rounding error\precision loss when adding up a lot of small numbers.

If you sum up one axis first, and then sum the result of that. You get the correct answer. 

print(np.sum(bx1**2+bx2**2+bx3**2, axis=0).sum())

Floating point's are a devious creatures, never trust them.

M4rtini
  • 13,186
  • 4
  • 35
  • 42
1

The other responses are excellent specific answers to this question. A general answer to all sorts of questions of this type can be found in the article What Every Computer Scientist Should Know About Floating-Point Arithmetic. The original was great when it was first published in 1991 and the edited reprint online is at least as good.

Simon
  • 10,679
  • 1
  • 30
  • 44