Velocidad Numpy vs Cython


Tengo un código de análisis que hace algunas operaciones numéricas pesadas usando numpy. Solo por curiosidad, traté de compilarlo con cython con pequeños cambios y luego lo reescribí usando bucles para la parte numpy.

Para mi sorpresa, el código basado en bucles era mucho más rápido (8x). No puedo publicar el código completo, pero armé un cálculo sin relación muy simple que muestra un comportamiento similar (aunque la diferencia de tiempo no es tan grande):

Versión 1 (sin cython)

import numpy as np

def _process(array):

    rows = array.shape[0]
    cols = array.shape[1]

    out = np.zeros((rows, cols))

    for row in range(0, rows):
        out[row, :] = np.sum(array - array[row, :], axis=0)

    return out

def main():
    data = np.load('data.npy')
    out = _process(data)
    np.save('vianumpy.npy', out)

Versión 2 (construyendo un módulo con cython)

import cython
cimport cython

import numpy as np
cimport numpy as np

DTYPE = np.float64
ctypedef np.float64_t DTYPE_t

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
cdef _process(np.ndarray[DTYPE_t, ndim=2] array):

    cdef unsigned int rows = array.shape[0]
    cdef unsigned int cols = array.shape[1]
    cdef unsigned int row
    cdef np.ndarray[DTYPE_t, ndim=2] out = np.zeros((rows, cols))

    for row in range(0, rows):
        out[row, :] = np.sum(array - array[row, :], axis=0)

    return out

def main():
    cdef np.ndarray[DTYPE_t, ndim=2] data
    cdef np.ndarray[DTYPE_t, ndim=2] out
    data = np.load('data.npy')
    out = _process(data)
    np.save('viacynpy.npy', out)

Versión 3 (construyendo un módulo con cython)

import cython
cimport cython

import numpy as np
cimport numpy as np

DTYPE = np.float64
ctypedef np.float64_t DTYPE_t

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
cdef _process(np.ndarray[DTYPE_t, ndim=2] array):

    cdef unsigned int rows = array.shape[0]
    cdef unsigned int cols = array.shape[1]
    cdef unsigned int row
    cdef np.ndarray[DTYPE_t, ndim=2] out = np.zeros((rows, cols))

    for row in range(0, rows):
        for col in range(0, cols):
            for row2 in range(0, rows):
                out[row, col] += array[row2, col] - array[row, col]

    return out

def main():
    cdef np.ndarray[DTYPE_t, ndim=2] data
    cdef np.ndarray[DTYPE_t, ndim=2] out
    data = np.load('data.npy')
    out = _process(data)
    np.save('vialoop.npy', out)

Con una matriz de 10000x10 guardada en datos.npy, los tiempos son:

$ python -m timeit -c "from version1 import main;main()"
10 loops, best of 3: 4.56 sec per loop

$ python -m timeit -c "from version2 import main;main()"
10 loops, best of 3: 4.57 sec per loop

$ python -m timeit -c "from version3 import main;main()"
10 loops, best of 3: 2.96 sec per loop

¿ Es esto esperado o hay una optimización que me falta? El hecho de que las versiones 1 y 2 den el mismo resultado es de alguna manera esperado, pero ¿por qué la versión 3 es más rápida?

Ps.- Este NO es el cálculo que necesito hacer, solo un simple ejemplo que muestra lo mismo cosa.

Author: Hernan, 2011-10-18

5 answers

Como se mencionó en las otras respuestas, la versión 2 es esencialmente la misma que la versión 1, ya que cython no puede profundizar en el operador de acceso a la matriz para optimizarlo. Hay 2 razones para esto

  • En primer lugar, hay una cierta cantidad de sobrecarga en cada llamada a una función numpy, en comparación con el código C optimizado. Sin embargo, esta sobrecarga será menos significativa si cada operación trata con arreglos grandes

  • En segundo lugar, existe la creación de matriz. Esto es más claro si considera una operación más compleja como out[row, :] = A[row, :] + B[row, :]*C[row, :]. En este caso, se debe crear una matriz completa B*C en memoria, y luego agregarla a A. Esto significa que el caché de la CPU está siendo golpeado, ya que los datos se leen y escriben en la memoria en lugar de mantenerse en la CPU y usarse de inmediato. Es importante destacar que este problema empeora si se trata de arreglos grandes.

Particularmente porque usted afirma que su código real es más complejo que su ejemplo, y muestra una aceleración mucho mayor, sospecho que la segunda razón es probable que sea el factor principal en su caso.

Como un aparte, si sus cálculos son lo suficientemente simples, puede superar este efecto utilizando numexpr, aunque por supuesto cython es útil en muchas más situaciones por lo que puede ser el mejor enfoque para usted.

 34
Author: DaveP,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/ajaxhispano.com/template/agent.layouts/content.php on line 61
2018-04-03 06:46:47

Con una ligera modificación, la versión 3 se vuelve el doble de rápida:

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
def process2(np.ndarray[DTYPE_t, ndim=2] array):

    cdef unsigned int rows = array.shape[0]
    cdef unsigned int cols = array.shape[1]
    cdef unsigned int row, col, row2
    cdef np.ndarray[DTYPE_t, ndim=2] out = np.empty((rows, cols))

    for row in range(rows):
        for row2 in range(rows):
            for col in range(cols):
                out[row, col] += array[row2, col] - array[row, col]

    return out

El cuello de botella en su cálculo es el acceso a la memoria. Su matriz de entrada está ordenada en C, lo que significa que moverse a lo largo del último eje hace el salto más pequeño en la memoria. Por lo tanto, su bucle interno debe estar a lo largo del eje 1, no del eje 0. Hacer este cambio reduce el tiempo de ejecución a la mitad.

Si necesita usar esta función en matrices de entrada pequeñas, puede reducir la sobrecarga usando np.empty en lugar de np.ones. Reducir la sobrecarga utiliza PyArray_EMPTY desde la API de numpy C.

Si utiliza esta función en matrices de entrada muy grandes (2**31), entonces los enteros utilizados para la indexación (y en la función range) se desbordarán. Para un uso seguro:

cdef Py_ssize_t rows = array.shape[0]
cdef Py_ssize_t cols = array.shape[1]
cdef Py_ssize_t row, col, row2

En lugar de

cdef unsigned int rows = array.shape[0]
cdef unsigned int cols = array.shape[1]
cdef unsigned int row, col, row2

Momento:

In [2]: a = np.random.rand(10000, 10)
In [3]: timeit process(a)
1 loops, best of 3: 3.53 s per loop
In [4]: timeit process2(a)
1 loops, best of 3: 1.84 s per loop

Donde process es su versión 3.

 39
Author: kwgoodman,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/ajaxhispano.com/template/agent.layouts/content.php on line 61
2012-05-07 22:41:39

Recomendaría usar la bandera-a para que cython genere el archivo html que muestra lo que se está traduciendo a c vs puro llamando a la API de python:

Http://docs.cython.org/src/quickstart/cythonize.html

La versión 2 da casi el mismo resultado que la Versión 1, porque todo el trabajo pesado lo está haciendo la API de Python (a través de numpy) y cython no está haciendo nada por ti. De hecho en mi máquina, numpy está construido contra MKL, así que cuando compilo el cython el código c generado usando gcc, la versión 3 es en realidad un poco más lenta que las otras dos.

Cython brilla cuando se está haciendo una manipulación de matriz que numpy no puede hacer de una manera 'vectorizada', o cuando se está haciendo algo intensivo en memoria que le permite evitar la creación de una matriz temporal grande. He conseguido 115x aceleraciones usando cython vs numpy para algunos de mi propio código:

Https://github.com/synapticarbors/pylangevin-integrator

Parte de eso era llamar el directorio randomkit al nivel del código c en lugar de llamarlo a través de numpy.random, pero la mayor parte de eso fue cython traduciendo los bucles for computacionalmente intensivos a c puro sin llamadas a python.

 7
Author: JoshAdel,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/ajaxhispano.com/template/agent.layouts/content.php on line 61
2011-10-18 00:21:55

La diferencia puede deberse a que las versiones 1 y 2 hacen una llamada a nivel de Python a np.sum() para cada fila, mientras que la versión 3 probablemente compila un bucle C apretado y puro.

Estudiar la diferencia entre la versión 2 y la fuente C generada por Cython 3 debería ser esclarecedor.

 3
Author: Pi Delport,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/ajaxhispano.com/template/agent.layouts/content.php on line 61
2011-10-17 22:17:04

Supongo que la sobrecarga principal que está guardando son los arrays temporales creados. Se crea una gran matriz array - array[row, :], luego se reduce a una matriz más pequeña usando sum. Pero construir esa gran matriz temporal no será gratis, especialmente si necesita asignar memoria.

 1
Author: wisty,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/ajaxhispano.com/template/agent.layouts/content.php on line 61
2011-10-17 23:14:58