La forma más eficiente de asignar la función sobre una matriz numpy

¿Cuál es la forma más eficiente de asignar una función a través de una matriz numpy? La forma en que lo he estado haciendo en mi proyecto actual es la siguiente:

import numpy as np x = np.array([1, 2, 3, 4, 5]) # Obtain array of square of each element in x squarer = lambda t: t ** 2 squares = np.array([squarer(xi) for xi in x]) 

Sin embargo, esto parece que es probablemente muy ineficiente, ya que estoy usando una comprensión de lista para construir la nueva matriz como una lista de Python antes de volver a convertirla en una matriz de números.

¿Podemos hacerlo mejor?

He probado todos los métodos sugeridos más np.array(map(f, x)) con perfplot (un pequeño proyecto mío).

Mensaje # 1: Si puede usar las funciones nativas de numpy, hágalo.

Si la función que está intentando vectorizar ya está vectorizada (como el ejemplo x**2 en la publicación original), el uso es mucho más rápido que cualquier otra cosa (observe la escala de registro):

introduzca la descripción de la imagen aquí

Si realmente necesita vectorización, realmente no importa mucho qué variante utilice.

introduzca la descripción de la imagen aquí


Código para reproducir las plots:

 import numpy as np import perfplot import math def f(x): # return math.sqrt(x) return np.sqrt(x) vf = np.vectorize(f) def array_for(x): return np.array([f(xi) for xi in x]) def array_map(x): return np.array(list(map(f, x))) def fromiter(x): return np.fromiter((f(xi) for xi in x), x.dtype) def vectorize(x): return np.vectorize(f)(x) def vectorize_without_init(x): return vf(x) perfplot.show( setup=lambda n: np.random.rand(n), n_range=[2**k for k in range(20)], kernels=[ f, array_for, array_map, fromiter, vectorize, vectorize_without_init ], logx=True, logy=True, xlabel='len(x)', ) 

¿Qué hay de usar numpy.vectorize .

 >>> import numpy as np >>> x = np.array([1, 2, 3, 4, 5]) >>> squarer = lambda t: t ** 2 >>> vfunc = np.vectorize(squarer) >>> vfunc(x) array([ 1, 4, 9, 16, 25]) 

https://docs.scipy.org/doc/numpy/reference/generated/numpy.vectorize.html

TL; DR

Como lo señaló @ user2357112 , un método “directo” para aplicar la función es siempre la forma más rápida y sencilla de asignar una función a través de matrices Numpy:

 import numpy as np x = np.array([1, 2, 3, 4, 5]) f = lambda x: x ** 2 squares = f(x) 

En general, evite np.vectorize , ya que no tiene un buen desempeño y tiene (o tuvo) varios problemas . Si está manejando otros tipos de datos, es posible que desee investigar los otros métodos que se muestran a continuación.

Comparacion de metodos

Aquí hay algunas pruebas simples para comparar tres métodos para asignar una función, este ejemplo se usa con Python 3.6 y NumPy 1.15.4. Primero, las funciones de configuración para la prueba:

 import timeit import numpy as np f = lambda x: x ** 2 vf = np.vectorize(f) def test_array(x, n): t = timeit.timeit( 'np.array([f(xi) for xi in x])', 'from __main__ import np, x, f', number=n) print('array: {0:.3f}'.format(t)) def test_fromiter(x, n): t = timeit.timeit( 'np.fromiter((f(xi) for xi in x), x.dtype, count=len(x))', 'from __main__ import np, x, f', number=n) print('fromiter: {0:.3f}'.format(t)) def test_direct(x, n): t = timeit.timeit( 'f(x)', 'from __main__ import x, f', number=n) print('direct: {0:.3f}'.format(t)) def test_vectorized(x, n): t = timeit.timeit( 'vf(x)', 'from __main__ import x, vf', number=n) print('vectorized: {0:.3f}'.format(t)) 

Pruebas con cinco elementos (ordenados de más rápido a más lento):

 x = np.array([1, 2, 3, 4, 5]) n = 100000 test_direct(x, n) # 0.265 test_fromiter(x, n) # 0.479 test_array(x, n) # 0.865 test_vectorized(x, n) # 2.906 

Con 100s de elementos:

 x = np.arange(100) n = 10000 test_direct(x, n) # 0.030 test_array(x, n) # 0.501 test_vectorized(x, n) # 0.670 test_fromiter(x, n) # 0.883 

Y con miles de elementos de matriz o más:

 x = np.arange(1000) n = 1000 test_direct(x, n) # 0.007 test_fromiter(x, n) # 0.479 test_array(x, n) # 0.516 test_vectorized(x, n) # 0.945 

Las diferentes versiones de Python / NumPy y la optimización del comstackdor tendrán resultados diferentes, así que haga una prueba similar para su entorno.

 squares = squarer(x) 

Las operaciones aritméticas en las matrices se aplican automáticamente de manera elemental, con bucles de nivel C eficientes que evitan toda la sobrecarga del intérprete que se aplicaría a un bucle o comprensión de nivel de Python.

La mayoría de las funciones que desearía aplicar a una matriz NumPy elementwise simplemente funcionarán, aunque algunas pueden necesitar cambios. Por ejemplo, if no funciona elementwise. Desea convertirlos para usar construcciones como numpy.where :

 def using_if(x): if x < 5: return x else: return x**2 

se convierte en

 def using_where(x): return numpy.where(x < 5, x, x**2) 

Desde que se respondió a esta pregunta, sucedieron muchas cosas: hay numexpr , numba y cython alrededor. El objective de esta respuesta es tener en cuenta estas posibilidades.

Pero primero digamos lo obvio: no importa cómo asigne una función Python a una matriz numpy, sigue siendo una función Python, lo que significa para cada evaluación:

  • El elemento numpy-array se debe convertir en un objeto Python (por ejemplo, un objeto Float ).
  • todos los cálculos se realizan con objetos Python, lo que significa tener la sobrecarga de intérprete, envío dynamic y objetos inmutables.

Por lo tanto, la maquinaria que se utiliza para recorrer realmente la matriz no juega un papel importante debido a la sobrecarga mencionada anteriormente: se mantiene mucho más lenta que con la vectorización de Numpy.

Echemos un vistazo al siguiente ejemplo:

 # numpy-functionality def f(x): return x+2*x*x+4*x*x*x # python-function as ufunc import numpy as np vf=np.vectorize(f) vf.__name__="vf" 

np.vectorize se selecciona como un representante de la clase de funciones de python puros de enfoques. Al usar perfplot (ver código en el apéndice de esta respuesta) obtenemos los siguientes tiempos de ejecución:

introduzca la descripción de la imagen aquí

Podemos ver, que el enfoque numérico es 10x-100x más rápido que la versión de python pura. La disminución del rendimiento para tamaños de matriz más grandes es probablemente porque los datos ya no se ajustan a la memoria caché.

Uno oye a menudo, que el rendimiento numpy es tan bueno como se obtiene, porque es C puro bajo el capó. Sin embargo, hay mucho espacio para mejorar!

La versión numpy vectorizada usa mucha memoria adicional y accesos a memoria. Numexp-library trata de agrupar los arreglos numpy y así obtener una mejor utilización del caché:

 # less cache misses than numpy-functionality import numexpr as ne def ne_f(x): return ne.evaluate("x+2*x*x+4*x*x*x") 

Conduce a la siguiente comparación:

introduzca la descripción de la imagen aquí

No puedo explicar todo en la gráfica anterior: al principio podemos ver una mayor sobrecarga para la biblioteca numexpr, pero como se utiliza mejor la memoria caché, ¡es 10 veces más rápido para matrices más grandes!


Otro enfoque es jit-comstackr la función y así obtener un UFunc puro-C real. Este es el enfoque de Numba:

 # runtime generated C-function as ufunc import numba as nb @nb.vectorize(target="cpu") def nb_vf(x): return x+2*x*x+4*x*x*x 

Es 10 veces más rápido que el método de numpy original:

introduzca la descripción de la imagen aquí


Sin embargo, la tarea es vergonzosamente paralelizable, por lo que también podríamos usar el prange para calcular el bucle en paralelo:

 @nb.njit(parallel=True) def nb_par_jitf(x): y=np.empty(x.shape) for i in nb.prange(len(x)): y[i]=x[i]+2*x[i]*x[i]+4*x[i]*x[i]*x[i] return y 

Como se esperaba, la función paralela es más lenta para entradas más pequeñas, pero más rápida (casi factor 2) para tamaños más grandes:

introduzca la descripción de la imagen aquí


Mientras que numba se especializa en optimizar las operaciones con numerosos arrays, Cython es una herramienta más general. Es más complicado extraer el mismo rendimiento que con numba; a menudo se reduce a llvm (numba) frente al comstackdor local (gcc / MSVC):

 %%cython -c=/openmp -a import numpy as np import cython #single core: @cython.boundscheck(False) @cython.wraparound(False) def cy_f(double[::1] x): y_out=np.empty(len(x)) cdef Py_ssize_t i cdef double[::1] y=y_out for i in range(len(x)): y[i] = x[i]+2*x[i]*x[i]+4*x[i]*x[i]*x[i] return y_out #parallel: from cython.parallel import prange @cython.boundscheck(False) @cython.wraparound(False) def cy_par_f(double[::1] x): y_out=np.empty(len(x)) cdef double[::1] y=y_out cdef Py_ssize_t i cdef Py_ssize_t n = len(x) for i in prange(n, nogil=True): y[i] = x[i]+2*x[i]*x[i]+4*x[i]*x[i]*x[i] return y_out 

Cython da como resultado funciones algo más lentas:

introduzca la descripción de la imagen aquí


Conclusión

Obviamente, probar solo una función no prueba nada. También se debe tener en cuenta que, para el ejemplo de la función elegida, el ancho de banda de la memoria fue el cuello de la botella para tamaños más grandes que 10 ^ 5 elementos, por lo que tuvimos el mismo rendimiento para numba, numexpr y cython en esta región.

Sin embargo, a partir de esta investigación y de mi experiencia hasta ahora, afirmaría que el numba parece ser la herramienta más fácil con el mejor rendimiento.


Trazado de tiempos de ejecución con perfplot -paquete :

 import perfplot perfplot.show( setup=lambda n: np.random.rand(n), n_range=[2**k for k in range(0,24)], kernels=[ f, vf, ne_f, nb_vf, nb_par_jitf, cy_f, cy_par_f, ], logx=True, logy=True, xlabel='len(x)' ) 

Creo que en la versión más nueva (yo uso 1.13) de numpy, simplemente puede llamar a la función pasando la matriz numpy a la función que escribió para el tipo escalar, aplicará automáticamente la llamada de función a cada elemento sobre la matriz numpy y le devolverá otra matriz numpy

 >>> import numpy as np >>> squarer = lambda t: t ** 2 >>> x = np.array([1, 2, 3, 4, 5]) >>> squarer(x) array([ 1, 4, 9, 16, 25]) 

Como se mencionó en esta publicación , solo use expresiones generadoras como:

 numpy.fromiter(((x) for x in ),,) 

Tal vez usar vectorizar es mejor

 def square(x): return x**2 vfunc=vectorize(square) vfunc([1,2,3,4,5]) output:array([ 1, 4, 9, 16, 25])