Tentativo di vectorizzare il calcolo iterativo con numpy

Question

Sto cercando di rendere più efficiente una parte di codice utilizzando il modulo vettoriale in numpy. Lascia che ti mostri un esempio in modo da sapere cosa intendo.

Dato il codice seguente:

a = np.zeros([4,4]) 
a[0] = [1., 2., 3., 4.] 
for i in range(len(a)-1): 
    a[i+1] = 2*a[i] 
print a 

Produce

[[ 1. 2. 3. 4.] 
[ 2. 4. 6. 8.] 
[ 4. 8. 12. 16.] 
[ 8. 16. 24. 32.]] 

Quando ora provo vectorize il codice come questo:

a = np.zeros([4,4]) 
a[0] = [1., 2., 3., 4.] 
a[1:] = 2*a[0:-1] 
print a 

ottengo solo la prima iterazione corretta :

[[ 1. 2. 3. 4.] 
[ 2. 4. 6. 8.] 
[ 0. 0. 0. 0.] 
[ 0. 0. 0. 0.]] 

E 'possibile scrivere il codice sopra in modo efficiente in una forma vettoriale (dove l'iterazione successiva accede sempre alla precedente iterazione) o devo mantenere il ciclo for?

Answer 1

Una ricorrenza lineare come questo può essere calcolato usando scipy.signal.lfilter:

In [19]: from scipy.signal import lfilter 

In [20]: num = np.array([1.0]) 

In [21]: alpha = 2.0 

In [22]: den = np.array([1.0, -alpha]) 

In [23]: a = np.zeros((4,4)) 

In [24]: a[0,:] = [1,2,3,4] 

In [25]: lfilter(num, den, a, axis=0) 
Out[25]: 
array([[ 1., 2., 3., 4.], 
     [ 2., 4., 6., 8.], 
     [ 4., 8., 12., 16.], 
     [ 8., 16., 24., 32.]]) 

Vedere quanto segue per maggiori dettagli: python recursive vectorization with timeseries, Recursive definitions in Pandas

---

Si noti che l'utilizzo di lfilter ha davvero senso solo se si risolve un problema non omogeneo lem come x[i+1] = alpha*x[i] + u[i], dove u è un determinato array di input. Per la semplice ricorrenza a[i+1] = alpha*a[i], è possibile utilizzare la soluzione esatta a[i] = a[0]*alpha**i. La soluzione per più valori iniziali può essere vettorizzata utilizzando la trasmissione.Ad esempio,

In [271]: alpha = 2.0 

In [272]: a0 = np.array([1, 2, 3, 4]) 

In [273]: n = 5 

In [274]: a0 * (alpha**np.arange(n).reshape(-1, 1)) 
Out[274]: 
array([[ 1., 2., 3., 4.], 
     [ 2., 4., 6., 8.], 
     [ 4., 8., 12., 16.], 
     [ 8., 16., 24., 32.], 
     [ 16., 32., 48., 64.]]) 

Answer 2

I calcoli vettoriali di Numpy agiscono sul vettore, non come una sequenza di passaggi, quindi è necessario vettorizzare l'intera espressione. Per esempio:

np.multiply(np.arange(1,5), 2**np.arange(0,4)[np.newaxis].T) 

di affrontare la questione "finale", sì, devi tenere il ciclo for se si vuole fare un calcolo sequenziale. Potresti renderlo più efficiente con map o [... for ...] ma l'ottimizzazione in questo modo richiede un sacco di tentativi ed errori. La bellezza del modo di pensare in termini vettoriali e l'uso di Numpy da implementare è che ottieni un risultato efficiente senza tutte le prove e gli errori.

Le funzioni cumsum e cumprod possono fare qualcosa di simile a quello che stai chiedendo. Invece di 2**np.arange(...), è possibile ottenere la stessa cosa da

np.multiply(np.arange(1,5), np.cumprod([1,2,2,2,])[np.newaxis].T) 

Answer 3

guardando il tuo matrice risultato

result = [[ 1, 2, 3, 4,] 
      [ 2, 4, 6, 8,] 
      [ 4, 8, 12, 16,] 
      [ 8, 16, 24, 32,]] 

può essere decostruita in prodotto (elem-saggio) di due matrici come

a = [[1, 2, 3, 4], 
    [1, 2, 3, 4], 
    [1, 2, 3, 4], 
    [1, 2, 3, 4]] 

b = [[1, 1, 1, 1], 
    [2, 2, 2, 2], 
    [4, 4, 4, 4] 
    [8, 8, 8, 8]] 

result = a * b 

è possibile calcolare questo tipo di operazione utilizzando la funzione meshgrid

aa, bb = np.meshgrid(np.array([1.0, 2.0, 3.0, 4.0]), 
        np.array([1.0, 2.0, 4.0, 8.0])) 
result = aa * bb