python - array - ¿Por qué es el einsum de numpy más rápido que las funciones integradas de numpy?

numpy vectorize

Question

dtype=np.double con tres matrices de dtype=np.double . Los tiempos se realizan en una CPU Intel usando numpy 1.7.1 compilado con icc y vinculado a intel mkl . También se usó una CPU AMD con numpy 1.6.1 compilado con gcc sin mkl para verificar los tiempos. Tenga en cuenta que las temporizaciones se escalan casi linealmente con el tamaño del sistema y no se deben a la pequeña sobrecarga incurrida en las funciones numpy, if declaraciones se mostrarán en microsegundos, no en milisegundos:

arr_1D=np.arange(500,dtype=np.double)
large_arr_1D=np.arange(100000,dtype=np.double)
arr_2D=np.arange(500**2,dtype=np.double).reshape(500,500)
arr_3D=np.arange(500**3,dtype=np.double).reshape(500,500,500)

Primero veamos la función np.sum :

np.all(np.sum(arr_3D)==np.einsum(''ijk->'',arr_3D))
True
%timeit np.sum(arr_3D)
10 loops, best of 3: 142 ms per loop
%timeit np.einsum(''ijk->'', arr_3D)
10 loops, best of 3: 70.2 ms per loop

Potestades:

np.allclose(arr_3D*arr_3D*arr_3D,np.einsum(''ijk,ijk,ijk->ijk'',arr_3D,arr_3D,arr_3D))
True
%timeit arr_3D*arr_3D*arr_3D
1 loops, best of 3: 1.32 s per loop
%timeit np.einsum(''ijk,ijk,ijk->ijk'', arr_3D, arr_3D, arr_3D)
1 loops, best of 3: 694 ms per loop

Producto externo:

np.all(np.outer(arr_1D,arr_1D)==np.einsum(''i,k->ik'',arr_1D,arr_1D))
True
%timeit np.outer(arr_1D, arr_1D)
1000 loops, best of 3: 411 us per loop
%timeit np.einsum(''i,k->ik'', arr_1D, arr_1D)
1000 loops, best of 3: 245 us per loop

Todo lo anterior es dos veces más rápido con np.einsum . Estas deberían ser comparaciones de manzanas a manzanas, ya que todo es específicamente de dtype=np.double . Esperaría la velocidad en una operación como esta:

np.allclose(np.sum(arr_2D*arr_3D),np.einsum(''ij,oij->'',arr_2D,arr_3D))
True
%timeit np.sum(arr_2D*arr_3D)
1 loops, best of 3: 813 ms per loop
%timeit np.einsum(''ij,oij->'', arr_2D, arr_3D)
10 loops, best of 3: 85.1 ms per loop

Einsum parece ser al menos dos veces más rápido para np.inner , np.outer , np.kron y np.sum independientemente de la selección de axes . La excepción principal es np.dot ya que llama DGEMM desde una biblioteca BLAS. Entonces, ¿por qué np.einsum es más rápido que otras funciones numpy que son equivalentes?

El caso de DGEMM para completar:

np.allclose(np.dot(arr_2D,arr_2D),np.einsum(''ij,jk'',arr_2D,arr_2D))
True
%timeit np.einsum(''ij,jk'',arr_2D,arr_2D)
10 loops, best of 3: 56.1 ms per loop
%timeit np.dot(arr_2D,arr_2D)
100 loops, best of 3: 5.17 ms per loop

La teoría principal es del comentario de @sebergs de que np.einsum puede hacer uso de SSE2 , pero ufuncs de numpy no lo hará hasta numpy 1.8 (ver el registro de cambios ). Creo que esta es la respuesta correcta, pero no he podido confirmarlo. Se pueden encontrar algunas pruebas limitadas cambiando el tipo de matriz de entrada y observando la diferencia de velocidad y el hecho de que no todos observan las mismas tendencias en los tiempos.

Answer 1

Ahora que se lanzó numpy 1.8, donde de acuerdo con los documentos todos los ufuncs deberían usar SSE2, quería verificar que el comentario de Seberg sobre SSE2 fuera válido.

Para realizar la prueba, se creó una nueva instalación de Python 2.7. Numpy 1.7 y 1.8 fueron compilados con icc usando opciones estándar en un núcleo opteron de AMD ejecutando Ubuntu.

Esta es la ejecución de prueba tanto antes como después de la actualización 1.8:

import numpy as np
import timeit
arr_1D=np.arange(5000,dtype=np.double)
arr_2D=np.arange(500**2,dtype=np.double).reshape(500,500)
arr_3D=np.arange(500**3,dtype=np.double).reshape(500,500,500)
print ''Summation test:''
print timeit.timeit(''np.sum(arr_3D)'',
''import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D'',
number=5)/5
print timeit.timeit(''np.einsum("ijk->", arr_3D)'',
''import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D'',
number=5)/5
print ''----------------------/n''
print ''Power test:''
print timeit.timeit(''arr_3D*arr_3D*arr_3D'',
''import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D'',
number=5)/5
print timeit.timeit(''np.einsum("ijk,ijk,ijk->ijk", arr_3D, arr_3D, arr_3D)'',
''import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D'',
number=5)/5
print ''----------------------/n''
print ''Outer test:''
print timeit.timeit(''np.outer(arr_1D, arr_1D)'',
''import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D'',
number=5)/5
print timeit.timeit(''np.einsum("i,k->ik", arr_1D, arr_1D)'',
''import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D'',
number=5)/5
print ''----------------------/n''
print ''Einsum test:''
print timeit.timeit(''np.sum(arr_2D*arr_3D)'',
''import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D'',
number=5)/5
print timeit.timeit(''np.einsum("ij,oij->", arr_2D, arr_3D)'',
''import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D'',
number=5)/5
print ''----------------------/n''

Numpy 1.7.1:

Summation test:
0.172988510132
0.0934836149216
----------------------
Power test:
1.93524689674
0.839519000053
----------------------
Outer test:
0.130380821228
0.121401786804
----------------------
Einsum test:
0.979052495956
0.126066613197

Numpy 1.8:

Summation test:
0.116551589966
0.0920487880707
----------------------
Power test:
1.23683619499
0.815982818604
----------------------
Outer test:
0.131808176041
0.127472200394
----------------------
Einsum test:
0.781750011444
0.129271841049

Creo que esto es bastante concluyente de que SSE juega un papel importante en las diferencias de tiempo, se debe tener en cuenta que al repetir estas pruebas los tiempos son muy poco por ~ 0.003s. La diferencia restante debe ser cubierta en las otras respuestas a esta pregunta.

Answer 2

Creo que estos tiempos explican lo que está pasando:

a = np.arange(1000, dtype=np.double)
%timeit np.einsum(''i->'', a)
100000 loops, best of 3: 3.32 us per loop
%timeit np.sum(a)
100000 loops, best of 3: 6.84 us per loop
a = np.arange(10000, dtype=np.double)
%timeit np.einsum(''i->'', a)
100000 loops, best of 3: 12.6 us per loop
%timeit np.sum(a)
100000 loops, best of 3: 16.5 us per loop
a = np.arange(100000, dtype=np.double)
%timeit np.einsum(''i->'', a)
10000 loops, best of 3: 103 us per loop
%timeit np.sum(a)
10000 loops, best of 3: 109 us per loop

Así que básicamente tienes una sobrecarga casi constante de 3us al llamar a np.sum sobre np.einsum , así que básicamente funcionan tan rápido, pero uno tarda un poco más en ponerse en marcha. ¿Por qué podría ser eso? Mi dinero está en lo siguiente:

a = np.arange(1000, dtype=object)
%timeit np.einsum(''i->'', a)
Traceback (most recent call last):
...
TypeError: invalid data type for einsum
%timeit np.sum(a)
10000 loops, best of 3: 20.3 us per loop

No estoy seguro de qué está sucediendo exactamente, pero parece que np.einsum omite algunas comprobaciones para extraer funciones específicas del tipo para hacer las multiplicaciones y adiciones, y va directamente con * y + para los tipos C estándar.

Los casos multidimensionales no son diferentes:

n = 10; a = np.arange(n**3, dtype=np.double).reshape(n, n, n)
%timeit np.einsum(''ijk->'', a)
100000 loops, best of 3: 3.79 us per loop
%timeit np.sum(a)
100000 loops, best of 3: 7.33 us per loop
n = 100; a = np.arange(n**3, dtype=np.double).reshape(n, n, n)
%timeit np.einsum(''ijk->'', a)
1000 loops, best of 3: 1.2 ms per loop
%timeit np.sum(a)
1000 loops, best of 3: 1.23 ms per loop

Por lo tanto, una sobrecarga constante, no una ejecución más rápida una vez que lleguen a ella.

Answer 3

En primer lugar, ha habido mucha discusión pasada sobre esto en la lista numpy. Por ejemplo, ver: http://numpy-discussion.10968.n7.nabble.com/poor-performance-of-sum-with-sub-machine-word-integer-types-td41.html http://numpy-discussion.10968.n7.nabble.com/odd-performance-of-sum-td3332.html

Algunos de ellos se reducen al hecho de que einsum es nuevo y presumiblemente intenta mejorar la alineación del caché y otros problemas de acceso a la memoria, mientras que muchas de las funciones numpy más antiguas se centran en una implementación fácilmente portátil en lugar de una muy optimizada. Solo estoy especulando, allí, sin embargo.

Sin embargo, algo de lo que estás haciendo no es una comparación "de manzanas a manzanas".

Además de lo que ya dijo @Jamie, sum usa un acumulador más apropiado para arreglos

Por ejemplo, sum tiene más cuidado al verificar el tipo de entrada y usar un acumulador apropiado. Por ejemplo, considere lo siguiente:

In [1]: x = 255 * np.ones(100, dtype=np.uint8)
In [2]: x
Out[2]:
array([255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255], dtype=uint8)

Tenga en cuenta que la sum es correcta:

In [3]: x.sum()
Out[3]: 25500

Mientras que einsum dará el resultado incorrecto:

In [4]: np.einsum(''i->'', x)
Out[4]: 156

Pero si utilizamos un tipo de dtype menos limitado, obtendremos el resultado que esperabas:

In [5]: y = 255 * np.ones(100)
In [6]: np.einsum(''i->'', y)
Out[6]: 25500.0