Python 32/64位机器浮点转置矩阵求和不正确？

Question

首先，我不是数学方面的高手，所以大数字的精度问题在我的日常工作中很少涉及。请多包涵。;)

我在用NumPy生成一个从1开始均匀分布的矩阵：

>>> m = numpy.matrix([(1.0 / 1000) for x in xrange(1000)]).T
>>> m
matrix[[ 0.001 ],
       [ 0.001 ],
       ...
       [ 0.001 ]])

在64位的Windows上，使用Python 2.6时，计算总和通常得不到1.0。用math.fsum()来计算这个矩阵的总和是可以的，但如果我把矩阵改成小数字，就不行了。

>>> numpy.sum(m)
1.0000000000000007
>>> math.fsum(m)
1.0
>>> sum(m)
matrix([[ 1.]])
>>> float(sum(m))
1.0000000000000007

在32位的Linux（Ubuntu）上，使用Python 2.6时，计算总和总是能得到1.0。

>>> numpy.sum(m)
1.0
>>> math.fsum(m)
1.0
>>> sum(m)
matrix([[ 1.]])
>>> float(sum(m))
1.0000000000000007

我可以在代码中加一个小的数（epsilon）来判断矩阵的总和是否等于1（比如 -epsilon < sum(m) < +epsilon），但我想先弄明白在Python中造成这种差异的原因，以及有没有更好的方法来正确计算总和。

我的理解是，计算总和时，程序处理数字的机器表示（浮点数）和它们显示的方式不同，而在计算总和时，使用的是内部表示。不过，看我用来计算总和的三种方法，为什么它们在不同平台上会有不同的结果，或者为什么有些结果是相同的，这一点还不太清楚。

计算一个矩阵的总和，最好的方法是什么呢？

如果你想要一个更有趣的矩阵，简单的修改一下就能得到更小的矩阵数字：

>>> m = numpy.matrix([(1.0 / 999) for x in xrange(999)]).T

提前感谢任何帮助！

更新 我想我发现了一些东西。如果我把存储的值改成32位浮点数，结果就和32位Linux的总和一致了。

>>> m = numpy.matrix([(numpy.float32(1.0) / 1000) for x in xrange(1000)]).T
>>> m
matrix[[ 0.001 ],
       [ 0.001 ],
       ...
       [ 0.001 ]])
>>> numpy.sum(m)
1.0

这样会把矩阵中的数字设置为32位浮点数，而不是64位，这样计算的总和就正确了。为什么在32位和64位系统上，0.001这个浮点数的机器表示不相等呢？如果我试图存储很多小数位的非常小的数字，我本来以为它们会不同。

有没有人对此有什么看法？在这种情况下，我应该明确切换到32位浮点数，还是有64位的计算总和的方法？还是我又得加上那个小的数（epsilon）？抱歉如果我听起来很傻，我对大家的看法很感兴趣。谢谢！

Answer 1

首先，如果你使用numpy来存储数据，最好使用numpy提供的方法来处理这些数组或矩阵。也就是说，如果你想相信那些非常厉害的人，他们开发了numpy。

现在，numpy的sum()函数在64位系统下计算的结果，可能无法精确地加到1。这是因为计算机处理浮点数的方式（murgatroid99给你提供了一个链接，还有很多类似的资料）。因此，唯一安全的方法（而且这也能帮助你更好地理解代码中的数学处理，从而更好地解决问题）就是使用一个epsilon值来限制到某个精度。

我为什么觉得这样做有帮助呢？因为计算科学需要处理错误，和实验科学一样，而在这个地方主动处理（也就是：确定这些错误）可以让你在解决代码中的计算错误时迈出第一步。

所以，可能还有其他方法来处理这个问题，但大多数情况下，我会使用一个epsilon值来确定我在特定问题上需要的精度。

Answer 2

我觉得最准确的方法（虽然不是最高效的）是使用decimal模块：

>>> from decimal import Decimal
>>> m = numpy.matrix([(Decimal(1) / 1000) for x in xrange(1000)])
>>> numpy.sum(m)
Decimal('1.000')
>>> numpy.sum(m) == 1.0
True

Answer 3

这是因为你在比较32位浮点数和64位浮点数，正如你已经发现的那样。

如果你在两台机器上都指定使用32位或64位的数据类型，你会看到相同的结果。

Numpy默认的浮点数数据类型（也就是numpy数组的数字类型）和机器的精度是一样的。这就是为什么你在不同的机器上看到不同的结果。

比如说，

32位版本：

m = numpy.ones(1000, dtype=numpy.float32) / 1000
print repr(m.sum())

和64位版本：

m = numpy.ones(1000, dtype=numpy.float64) / 1000
print repr(m.sum())

由于精度不同，它们的结果会不同，但在不同的机器上你会看到相同的结果。（不过，在32位机器上，64位操作会慢很多）

如果你只指定 numpy.float，那么它会根据机器的本地架构变成 float32 或 float64。