¿Por qué ensum de numpy es más rápido que las funciones integradas de numpy?

Vamos a empezar con tres matrices de dtype=np.double . Los tiempos se realizan en una CPU Intel utilizando un número 1.7.1 comstackdo con icc y vinculado a mkl de mkl . También se usó una CPU de AMD con numpy 1.6.1 comstackda con gcc sin mkl para verificar los tiempos. Tenga en cuenta que los tiempos se escalan casi linealmente con el tamaño del sistema y no se deben a la pequeña sobrecarga en las funciones numpy if declaraciones de esta diferencia se mostrarán en microsegundos y 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 exterior:

 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 deben ser comparaciones de manzanas con manzanas, ya que todo es específicamente de dtype=np.double . Yo esperaría la aceleración 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 el doble de 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 a 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 la DGEMM para la integridad:

 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 proviene del comentario de @sebergs de que np.einsum puede hacer uso de SSE2 , pero los ufuncs de numpy no lo harán hasta el numpy 1.8 (consulte el registro de cambios ). Creo que esta es la respuesta correcta, pero no he podido confirmarlo. Se puede encontrar alguna prueba limitada 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.

En primer lugar, ha habido muchas discusiones pasadas sobre esto en la lista de números. Por ejemplo, consulte: 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, está tratando de ser mejor en cuanto a la alineación de la memoria 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 sobre una altamente optimizada. Sin embargo, sólo estoy especulando.


Sin embargo, algo de lo que estás haciendo no es una comparación “manzanas con manzanas”.

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

Por ejemplo, la sum es más cuidadosa al verificar el tipo de entrada y al usar un acumulador adecuado. 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 equivocado:

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

Pero si usamos un dtype menos limitado, todavía obtendremos el resultado que esperaría:

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

Ahora que se ha publicado el número 1.8, donde, según la documentación, todos los ufuncs deberían usar el SSE2, quería verificar que el comentario de Seberg sobre el SSE2 fuera válido.

Para realizar la prueba, se creó una nueva instalación de python 2.7: los números 1.7 y 1.8 se comstackron con icc usando las opciones estándar en un núcleo AMD Opteron que ejecuta Ubuntu.

Esta es la ejecución de prueba antes y 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 la ESS juega un papel importante en las diferencias de tiempo, se debe tener en cuenta que, al repetir estas pruebas, los tiempos están muy cerca de ~ 0.003s. La diferencia restante debe ser cubierta en las otras respuestas a esta pregunta.

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 

Básicamente, tiene una sobrecarga casi constante cuando llama a np.sum través de np.einsum , por lo que básicamente se ejecutan a la misma velocidad, pero uno tarda un poco más en np.einsum . ¿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 lo que está sucediendo exactamente, pero parece que np.einsum está omitiendo 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 en su mayoría constante, no una carrera más rápida una vez que llegan a ello.