Vectorizar un bucle de python sobre una matriz numpy

Necesito acelerar el procesamiento de este bucle ya que es muy lento. Pero no sé cómo vectorizarlo, ya que el resultado de un valor depende del resultado de un valor anterior. ¿Alguna sugerencia?

import numpy as np sig = np.random.randn(44100) alpha = .9887 beta = .999 out = np.zeros_like(sig) for n in range(1, len(sig)): if np.abs(sig[n]) >= out[n-1]: out[n] = alpha * out[n-1] + (1 - alpha) * np.abs( sig[n] ) else: out[n] = beta * out[n-1] 

Bajo potencial de vectorización en un código de “bucle dependiente de avance”

La mayoría de su paralelismo de “vectorización” está fuera del juego, una vez que se analiza la dependencia. (El comstackdor JIT tampoco puede vectorizar “contra” dicha barrera de dependencia)

puede pre-calcular algunos valores reutilizados de una manera vectorizada, pero no hay una manera de syntax de Python directa (sin una solución de comstackción JIT externa) para organizar el cálculo del bucle de dependencia de desplazamiento hacia adelante en su CPU vector-registro alineado co- computación paralela:

 from zmq import Stopwatch # ok to use pyzmq 2.11 for [usec] .Stopwatch() aStopWATCH = Stopwatch() # a performance measurement .Stopwatch() instance sig = np.abs(sig) # self-destructive calc/assign avoids memalloc-OPs aConst = ( 1 - alpha ) # avoids many repetitive SUB(s) in the loop for thisPtr in range( 1, len( sig ) ): # FORWARD-SHIFTING-DEPENDENCE LOOP: prevPtr = thisPtr - 1 # prevPtr->"previous" TimeSlice in out[] ( re-used 2 x len(sig) times ) if sig[thisPtr] < out[prevPtr]: # 1st re-use out[thisPtr] = out[prevPtr] * beta # 2nd else: out[thisPtr] = out[prevPtr] * alpha + ( aConst * sig[thisPtr] ) # 2nd 

Se puede ver un buen ejemplo de aceleración vectorizada en los casos en los que la estrategia de cálculo se puede paralelizar / transmitir a lo largo de la estructura 1D, 2D o incluso 3D de la matriz numpy nativa. Para una aceleración de aproximadamente 100x, vea un procesamiento acelerado de matriz RGBA-2D en código Vectorizado para un procesamiento de imagen PNG (una tubería de sombreado OpenGL)

El rendimiento aumentó todavía alrededor de 3x

Incluso esta simple revisión del código de python ha aumentado la velocidad más de aproximadamente 2.8 veces (en este momento, es decir, sin realizar una instalación para permitir el uso de un comstackdor de optimización de JIT ad hoc):

 >>> def aForwardShiftingDependenceLOOP(): # proposed code-revision ... aStopWATCH.start() # 

.start … for thisPtr in range( 1, len( sig ) ): … # |vvvvvvv|————# FORWARD-SHIFTING-LOOP DEPENDENCE … prevPtr = thisPtr – 1 #|vvvvvvv|–STEP-SHIFTING avoids Numpy syntax … if ( sig[ thisPtr] < out[prevPtr] ): ... out[ thisPtr] = out[prevPtr] * beta ... else: ... out[ thisPtr] = out[prevPtr] * alpha + ( aConst * sig[thisPtr] ) ... usec = aStopWATCH.stop() #

.stop … print usec, ” [usec]” >>> aForwardShiftingDependenceLOOP() 57593 [usec] 57879 [usec] 58085 [usec] >>> def anOriginalForLOOP(): … aStopWATCH.start() … for n in range( 1, len( sig ) ): … if ( np.abs( sig[n] ) >= out[n-1] ): … out[n] = out[n-1] * alpha + ( 1 – alpha ) * np.abs( sig[n] ) … else: … out[n] = out[n-1] * beta … usec = aStopWATCH.stop() … print usec, ” [usec]” >>> anOriginalForLOOP() 164907 [usec] 165674 [usec] 165154 [usec]

El comstackdor justo a tiempo de Numba debería lidiar con la sobrecarga de indexación a la que se enfrenta al comstackr la función en código nativo durante la primera ejecución:

 In [1]: %cpaste Pasting code; enter '--' alone on the line to stop or use Ctrl-D. :import numpy as np : :sig = np.random.randn(44100) :alpha = .9887 :beta = .999 : :def nonvectorized(sig): : out = np.zeros_like(sig) : : for n in range(1, len(sig)): : if np.abs(sig[n]) >= out[n-1]: : out[n] = alpha * out[n-1] + (1 - alpha) * np.abs( sig[n] ) : else: : out[n] = beta * out[n-1] : return out :-- In [2]: nonvectorized(sig) Out[2]: array([ 0. , 0.01862503, 0.04124917, ..., 1.2979579 , 1.304247 , 1.30294275]) In [3]: %timeit nonvectorized(sig) 10 loops, best of 3: 80.2 ms per loop In [4]: from numba import jit In [5]: vectorized = jit(nonvectorized) In [6]: np.allclose(vectorized(sig), nonvectorized(sig)) Out[6]: True In [7]: %timeit vectorized(sig) 1000 loops, best of 3: 249 µs per loop 

EDITAR: como se sugiere en un comentario, agregando puntos de referencia jit. jit(nonvectorized) está creando una envoltura ligera y, por lo tanto, es una operación barata.

 In [8]: %timeit jit(nonvectorized) 10000 loops, best of 3: 45.3 µs per loop 

La función en sí misma se comstack durante la primera ejecución (por lo tanto, justo a tiempo ), lo que lleva un tiempo, pero probablemente no tanto:

 In [9]: %timeit jit(nonvectorized)(sig) 10 loops, best of 3: 169 ms per loop