Perché x ** 3 è più lento di x * x * x?

Question

In NumPy, x * x * x è un ordine di grandezza più veloce di x ** 3 o persino np.power (x, 3).

x = np.random.rand(1e6) 
%timeit x**3 
100 loops, best of 3: 7.07 ms per loop 

%timeit x*x*x 
10000 loops, best of 3: 163 µs per loop 

%timeit np.power(x, 3) 
100 loops, best of 3: 7.15 ms per loop 

Qualche idea sul perché questo comportamento si verifica? Per quanto ne so, tutti e tre producono lo stesso risultato (verificato con np.allclose).

Answer 1

Come per this answer, è perché l'implementazione di esponenziazione ha un sovraccarico che la moltiplicazione non ha. Tuttavia, la moltiplicazione ingenua diventerà più lenta e più lenta man mano che l'esponente aumenta. Una dimostrazione empirica:

In [3]: x = np.random.rand(1e6) 

In [15]: %timeit x**2 
100 loops, best of 3: 11.9 ms per loop 

In [16]: %timeit x*x 
100 loops, best of 3: 12.7 ms per loop 

In [17]: %timeit x**3 
10 loops, best of 3: 132 ms per loop 

In [18]: %timeit x*x*x 
10 loops, best of 3: 27.2 ms per loop 

In [19]: %timeit x**4 
10 loops, best of 3: 132 ms per loop 

In [20]: %timeit x*x*x*x 
10 loops, best of 3: 42.4 ms per loop 

In [21]: %timeit x**10 
10 loops, best of 3: 132 ms per loop 

In [22]: %timeit x*x*x*x*x*x*x*x*x*x 
10 loops, best of 3: 137 ms per loop 

In [24]: %timeit x**15 
10 loops, best of 3: 132 ms per loop 

In [25]: %timeit x*x*x*x*x*x*x*x*x*x*x*x*x*x*x 
1 loops, best of 3: 212 ms per loop 

Annotare il tempo di elevamento a potenza rimane più o meno costante, fatta eccezione per il caso x**2 che ho il sospetto è speciale prima lettera maiuscola, mentre la moltiplicazione diventa più lento e più lento. Sembra che si potrebbe sfruttare questo per ottenere più velocemente intero elevamento a potenza ... ad esempio:

In [26]: %timeit x**16 
10 loops, best of 3: 132 ms per loop 

In [27]: %timeit x*x*x*x*x*x*x*x*x*x*x*x*x*x*x*x 
1 loops, best of 3: 225 ms per loop 

In [28]: def tosixteenth(x): 
    ....:  x2 = x*x 
    ....:  x4 = x2*x2 
    ....:  x8 = x4*x4 
    ....:  x16 = x8*x8 
    ....:  return x16 
    ....: 

In [29]: %timeit tosixteenth(x) 
10 loops, best of 3: 49.5 ms per loop 

Sembra si potrebbe applicare questa tecnica genericamente suddividendo ogni intero in una somma delle potenze di due, calcolando ogni potenza di due come sopra, e di somma:

In [93]: %paste 
def smartintexp(x, exp): 
    result = np.ones(len(x)) 
    curexp = np.array(x) 
    while True: 
     if exp%2 == 1: 
      result *= curexp 
     exp >>= 1 
     if not exp: break 
     curexp *= curexp 
    return result 
## -- End pasted text -- 

In [94]: x 
Out[94]: 
array([ 0.0163407 , 0.57694587, 0.47336487, ..., 0.70255032, 
     0.62043303, 0.0796748 ]) 

In [99]: x**21 
Out[99]: 
array([ 3.01080670e-38, 9.63466181e-06, 1.51048544e-07, ..., 
     6.02873388e-04, 4.43193256e-05, 8.46721060e-24]) 

In [100]: smartintexp(x, 21) 
Out[100]: 
array([ 3.01080670e-38, 9.63466181e-06, 1.51048544e-07, ..., 
     6.02873388e-04, 4.43193256e-05, 8.46721060e-24]) 

In [101]: %timeit x**21 
10 loops, best of 3: 132 ms per loop 

In [102]: %timeit smartintexp(x, 21) 
10 loops, best of 3: 70.7 ms per loop 

e 'veloce per le piccole potenze anche di due:

In [106]: %timeit x**32 
10 loops, best of 3: 131 ms per loop 

In [107]: %timeit smartintexp(x, 32) 
10 loops, best of 3: 57.4 ms per loop 

Ma diventa più lento come esponente diventa più grande:

In [97]: %timeit x**63 
10 loops, best of 3: 133 ms per loop 

In [98]: %timeit smartintexp(x, 63) 
10 loops, best of 3: 110 ms per loop 

E non più veloce per le grandi peggiore dei casi:

In [115]: %timeit x**511 
10 loops, best of 3: 135 ms per loop 

In [114]: %timeit smartintexp(x, 511) 
10 loops, best of 3: 192 ms per loop 

Answer 2

Questo perché i poteri in python vengono eseguiti come operazioni float (questo vale anche per numpy, perché usa C).

In C, il pow function prevede 3 metodi:

doppia pow (double x, double y)

lungo POWL (long double x, long double y)

galleggiante powf (float x, float y)

Ognuna di queste operazioni è in virgola mobile.

Answer 3

Mi aspetto che sia perché x**y deve gestire il caso generico in cui sia x sia sono float. Matematicamente possiamo scrivere x**y = exp(y*log(x)). Seguendo il tuo esempio trovo

x = np.random.rand(1e6) 
%timeit x**3 
10 loops, best of 3: 178 ms per loop 

%timeit np.exp(3*np.log(x)) 
10 loops, best of 3: 176 ms per loop 

Non ho controllato il codice vero e proprio NumPy ma devo fare qualcosa di simile internamente.

Answer 4

timeit np.multiply(np.multiply(x,x),x) 

volte lo stesso di x*x*x. La mia ipotesi è che np.multiply stia usando un pacchetto di algebra lineare Fortran veloce come BLAS. So da un altro problema che numpy.dot utilizza BLAS per determinati casi.

---

devo prendere quella posteriore. np.dot(x,x) è 3 volte più veloce di np.sum(x*x). Quindi il vantaggio di velocità su np.multiply non è coerente con l'utilizzo di BLAS.

---

Con il mio NumPy (orari varia a seconda della macchina e librerie disponibili)

np.power(x,3.1) 
np.exp(3.1*np.log(x)) 

prendere circa lo stesso tempo, ma

np.power(x,3) 

è 2x più veloce. Non veloce come x*x*x, ma ancora più veloce della potenza generale. Quindi sta sfruttando un po 'il potere dell'intero.

Answer 5

Come nota se si sta calcolando poteri e sono preoccupati per la velocità:

x = np.random.rand(5e7) 

%timeit x*x*x 
1 loops, best of 3: 522 ms per loop 

%timeit np.einsum('i,i,i->i',x,x,x) 
1 loops, best of 3: 288 ms per loop 

Perché einsum è più veloce è ancora una questione aperta di mine. Anche se è simile a einsum in grado di utilizzare SSE2 mentre gli ufunc di Numpy non saranno disponibili fino al 1.8.

Al posto è ancora più veloce:

def calc_power(arr): 
    for x in xrange(arr.shape[0]): 
     arr[x]=arr[x]*arr[x]*arr[x] 
numba_power = autojit(calc_power) 

%timeit numba_power(x) 
10 loops, best of 3: 51.5 ms per loop 

%timeit np.einsum('i,i,i->i',x,x,x,out=x) 
10 loops, best of 3: 111 ms per loop 

%timeit np.power(x,3,out=x) 
1 loops, best of 3: 609 ms per loop 

Answer 6

Secondo il spec:

I due argomenti forma pow (x, y) è equivalente usando l'operatore di potenza : x * * y.

Gli argomenti devono avere tipi numerici. Con i tipi di operandi misti, si applicano le regole di coercizione per gli operatori aritmetici binari.

In altre parole: dal x è un galleggiante, l'esponente viene convertito da un int ad un galleggiante, e viene eseguita l'operazione di alimentazione a virgola mobile generico. Internamente, questo di solito è riscritto come:

x**y = 2**(y*lg(x)) 

2**a e lg a (base 2 logaritmo a) sono singole istruzioni su processori moderni, ma è ancora richiede molto più di un paio di moltiplicazioni.