在NumPy中，使用向量化代码和标准循环产生不同的结果

Question

我有以下两个函数：

def loop(x):
    a = np.zeros(10)
    for i1 in range(10):
        for i2 in range(10):
            a[i1] += np.sin(x[i2] - x[i1])
    return a

并且

def vectorized(x):    
    b = np.zeros(10)
    for i1 in range(10):
        b += np.sin(np.roll(x, i1) - x)
    return b

然而，当我运行这两个时，发现它们的结果略有不同：

x = np.arange(10)
a, b = loop(x), vectorized(x)
print b - a

我收到：

[  2.22044605e-16   0.00000000e+00   0.00000000e+00   6.66133815e-16
  -2.22044605e-16   2.22044605e-16   0.00000000e+00   2.22044605e-16
   2.22044605e-16   2.22044605e-16]

这个问题很小，但在我的情况下会影响模拟。如果我从函数中删除np.sin，则差异消失。或者，如果使用np.float32代替x，差异也会消失，但这是由float64解决器求解的ode的一部分。有没有办法解决这个差异？

- patatapon

---

很遗憾，当你像你所做的那样重新排列操作时，很难（或者不可能？）保持舍入误差顺序上的差异。如果这实际上对最终解决方案有重大影响，你需要开始考虑是否需要选择一个不太敏感的ODE求解器（或者如果你试图解决的系统表现出混沌行为，从而使某些类型的建模不可能）。 - mgilson

Answer 1

因为您没有按照相同的顺序进行操作。

要获得等效的完全向量化解决方案，请执行c=sin(add.outer(x,-x))).sum(axis=0)。

In [8]: (c==loop(x)).all()
Out[8]: True

您可以获得向量化的全部优势：

In [9]: %timeit loop(x)
1000 loops, best of 3: 750 µs per loop

In [10]: %timeit vectorized(x)
1000 loops, best of 3: 347 µs per loop

In [11]: %timeit sin(x[:,None]-x).sum(axis=0)
10000 loops, best of 3: 46 µs per loop

- B. M.