risultato diverso con il codice di vettorializzare al ciclo standard NumPy

Question

Ho le seguenti due funzioni:

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 

e

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

Tuttavia, quando corro sia, trovo che i loro risultati differiscono leggermente:

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

ottengo:

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

che è molto piccolo, ma nel mio caso, influisce sulla simulazione. Se rimuovo np.sin dalle funzioni, la differenza scompare. In alternativa la differenza scompare anche se si usa np.float32 per x, ma questo fa parte di un'ode che viene risolta da un risolutore che usa float64. C'è un modo per risolvere questa differenza?

Answer 1

È perché non si effettua l'operazione nello stesso ordine.

Per la soluzione equivalente totalmente vettoriale, fare c=sin(add.outer(x,-x))).sum(axis=0).

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

E si vince la piena avantage di vettorizzazione:

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