python3上numpy中float32数字的矩阵乘法不稳定性

Question

当我使用numpy. 简而言之，对于数学上相同的计算，numpy给了我不同的结果。

下面的代码片段说明了这个问题：

import numpy as np
np.set_printoptions(precision=10)
np.random.seed(2)  # for reproducibility

# create matrices:
A = np.random.rand(300, 1)
A = A/np.sum(A)
A = np.repeat(A, 15, 1)
B = np.random.rand(300, 300)

# convert data to float32:
A_star = A.astype("float32")
B_star = B.astype("float32")

# do the matrix multiplication A'BA and take the diagonal:
res_star = np.diag(A_star.transpose().dot(B_star.dot(A_star)))

# print the results:
print(res_star)

在此运行python3.5与numpy1.11.1Windows机器上打印以下数组：

[0.5000683069  0.5000683069  0.5000683069  0.5000683069  0.5000683069
 0.5000683069  0.5000683069  0.5000683069  0.5000683069  0.5000683069
 0.5000683069  0.5000683069  0.5000681877  0.5000681877  0.5000683069]

请注意， 中的值res_star[12:14]与数组的其他元素不同——尽管从数学上讲，人们希望它们重合。

我自己对这些差异的调查并不是特别成功，但我想我缩小了一点范围：

• 使用 dtype float64，结果数组中的值是相同的
• 对于float32值，值的差异通常仅从第 7 个小数位开始（float32精度结束的区域）

然而，对我来说，这些发现并不能真正解释数组中的差异，因为发生了相同的数学运算（如果我要在相同的矩阵上使用 dtype 计算相同的点积，它们确实可以解释值的差异float64）。此外，对于不同的种子和矩阵大小，结果值实际上可能重合 - 因此这种行为不一致。

• 在 UNIX 机器上，上面给出的代码片段返回相同的数字（至少在这个特定配置中）

现在，

1）是什么导致了这种行为和

2）如何在不将 dtype 更改为 的情况下确保数组中的结果相同float64？