矢量化代码与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,差异就会消失.或者,如果对x使用np.float32,差异也会消失,但这是使用float64的解算器解决的ode的一部分.有没有办法解决这个差异？

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