Different result with vectorized code to standard loop in numpy

Question

I have the following two functions:

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

and

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

However, when I run both, I find that their results slightly differ:

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

I get:

[  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]

which is very small, but in my case, affects the simulation. If I remove the np.sin from the functions, the difference disappears. Alternatively the difference also goes away if use np.float32 for x, but this is part of an ode which is being solved by by a solver that uses float64. Is there a way to resolve this difference?

Answer 1

It's because you don't make the operation in the same order.

For the equivalent totally vectored solution, do c=sin(add.outer(x,-x))).sum(axis=0).

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

And you win the full avantage of vectorisation :

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