numpy sum 没有给出 float32 类型的正确答案答案

【问题标题】：numpy sum is not giving right answer for float32 typenumpy sum 没有给出 float32 类型的正确答案
【发布时间】：2014-02-28 23:44:15
【问题描述】：

我正在使用 numpy 对数组求和，但它不适用于 float32 类型。我究竟做错了什么？一旦我直接进行求和，然后我使用 numpy.sum。请看下面的代码

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))

两个输出应该匹配。这给出了输出：
30408704.0
3.1323e+07

直接求和给出了正确的结果，而 np.sum 给出了错误的结果。但是，如果我使用 float64，那么 np.sum 会给出正确的结果。这背后的原因是什么？

谢谢。

【问题讨论】：

标签： python numpy

【解决方案1】：

这么大的数字，float32的精度有问题。我还没有介绍 2.25 是如何存储的，但是最小的例子

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

表明您失去了准确性，但使用 float64（python 的 standard float）恢复了准确性。

【讨论】：

【解决方案2】：

将大量小数相加时可能导致舍入误差\精度损失。

如果你先对一个轴求和，然后对结果求和。你得到正确的答案。

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

浮点数是一种狡猾的生物，永远不要相信它们。

【讨论】：

【解决方案3】：

其他回答是对这个问题的很好的具体回答。此类问题的一般答案可以在文章What Every Computer Scientist Should Know About Floating-Point Arithmetic 中找到。原版在 1991 年首次出版时很棒，在线编辑的再版至少也一样好。

【讨论】：