¿Es posible vectorizar el cálculo recursivo de una matriz NumPy donde cada elemento depende de la anterior?

T(i) = Tm(i) + (T(i-1)-Tm(i))**(-tau(i)) 

Tm y tau son vectores NumPy de la misma longitud que se han calculado previamente, y el deseo es crear un nuevo vector T La i se incluye solo para indicar el índice del elemento para lo que se desea.

¿Es necesario un bucle for para este caso?

Podrías pensar que esto funcionaría:

 import numpy as np n = len(Tm) t = np.empty(n) t[0] = 0 # or whatever the initial condition is t[1:] = Tm[1:] + (t[0:n-1] - Tm[1:])**(-tau[1:]) 

pero no lo hace: en realidad no se puede hacer la recursión de forma entrecortada (ya que numpy calcula el RHS completo y luego lo asigna a la LHS).

Por lo tanto, a menos que pueda encontrar una versión no recursiva de esta fórmula, tendrá un bucle explícito:

 tt = np.empty(n) tt[0] = 0. for i in range(1,n): tt[i] = Tm[i] + (tt[i-1] - Tm[i])**(-tau[i]) 

En algunos casos, es posible tener este tipo de recursión, es decir, si hay un NumPy ufunc para la fórmula de recursión, por ejemplo

 c = numpy.arange(10.) numpy.add(c[:-1], c[1:], c[1:]) 

Esto calcula las sums acumulativas sobre c en su lugar, utilizando el parámetro de salida de numpy.add . No puede ser escrito como

 c[1:] = c[:-1] + c[1:] 

porque ahora el resultado de la adición es un temporal que se copia en c[1:] finalizado el cálculo.

Lo más natural que puedes intentar ahora es definir tu propio ufunc:

 def f(T, Tm, tau): return Tm + (T - Tm)**(-tau) uf = numpy.frompyfunc(f, 3, 1) 

Pero por razones que están más allá de mí, lo siguiente no funciona:

 uf(T[:-1], Tm[1:], tau[1:], T[1:]) 

Al parecer, el resultado no se escribe directamente en T[1:] , sino que se almacena de forma temporal y se copia una vez finalizada la operación. Incluso si funcionara, no esperaría ninguna aceleración de esto en comparación con un bucle ordinario, ya que necesita llamar a una función de Python para cada entrada.

Si realmente quiere evitar el bucle de Python, probablemente deba elegir Cython o ctypes.

Realicé algunos puntos de referencia y en 2018 utilizando Numba es la primera opción que las personas deberían intentar acelerar las funciones recursivas en Numpy (propuesta de Aronstef). Numba ya está preinstalado en el paquete Anaconda y tiene uno de los tiempos más rápidos (aproximadamente 20 veces más rápido que cualquier Python). En 2018, Python admite las anotaciones de @numba sin pasos adicionales (al menos las versiones 3.6 y 3.7). Aquí hay dos puntos de referencia: uno realizado en 2018-10-20 y el otro en 2016-05-18.

Y, como lo mencionó Jaffe, en 2018 todavía no es posible vectorizar las funciones recursivas. Revisé la vectorización por Aronstef y NO funciona.

Puntos de referencia ordenados por tiempo de ejecución:

 ----------------------------------- |Variant |2018-10 |2016-05 | ----------------------------------- |Pure C | na | 2.75 ms| |C extension | na | 6.22 ms| |Cython float32 | 1.01 ms| na | |Cython float64 | 1.05 ms| 6.26 ms| |Fortran f2py | na | 6.78 ms| |Numba float32 | 2.81 ms| na | |(Aronstef) | | | |Numba float64 | 5.28 ms| na | |Append to list |48.2 ms|91.0 ms| |Using a.item() |58.3 ms|74.4 ms| |np.fromiter() |60.0 ms|78.1 ms| |Loop over Numpy|71.9 ms|87.9 ms| |(Jaffe) | | | |Loop over Numpy|74.4 ms| na | |(Aronstef) | | | ----------------------------------- 

El código correspondiente se proporciona al final de la respuesta.

No verifiqué Pure C en 2018, pero supongo que sigue siendo el más rápido según el índice de referencia anterior.

Tampoco verifiqué la extensión C en 2018 y creo que tiene casi el mismo tiempo que Cython en base al punto de referencia anterior.

Fortran fue muy difícil de depurar y comstackr, así que no revisé la versión f2py en 2018. Y de todos modos era peor que Cython.

Tengo la siguiente configuración en 2018:

 Processor: Intel i7-7500U 2.7GHz Versions: Python: 3.7.0 Numba: 0.39.0 Cython: 0.28.5 Numpy: 1.15.1 

El código Numba recomendado que usa float32 (de Aronstef):

 @numba.jit("float32[:](float32[:], float32[:])", nopython=False, nogil=True) def calc_py_jit32(Tm_, tau_): tt = np.empty(len(Tm_),dtype="float32") tt[0] = Tm_[0] for i in range(1, len(Tm_)): tt[i] = Tm_[i] - (tt[i-1] + Tm_[i])**(-tau_[i]) return tt[1:] 

Todos los demás códigos:

Creación de datos (como comentario de Aronstef + Mike T):

 np.random.seed(0) n = 100000 Tm = np.cumsum(np.random.uniform(0.1, 1, size=n).astype('float64')) tau = np.random.uniform(-1, 0, size=n).astype('float64') ar = np.column_stack([Tm,tau]) Tm32 = Tm.astype('float32') tau32 = tau.astype('float32') Tm_l = list(Tm) tau_l = list(tau) 

El código en 2016 fue ligeramente diferente, ya que usé la función abs () para evitar nans y no la variante de Mike T. En 2018, la función es exactamente la misma que la que escribió OP (Póster original).

Cython float32 usando Jupyter %% magic. La función se puede utilizar directamente en Python . Cython necesita un comstackdor de C ++ en el que se compiló Python. La instalación de la versión correcta del comstackdor de Visual C ++ (para Windows) podría ser problemática:

 %%cython import cython import numpy as np cimport numpy as np from numpy cimport ndarray cdef extern from "math.h": np.float32_t exp(np.float32_t m) @cython.boundscheck(False) @cython.wraparound(False) @cython.infer_types(True) @cython.initializedcheck(False) def cy_loop32(np.float32_t[:] Tm,np.float32_t[:] tau,int alen): cdef np.float32_t[:] T=np.empty(alen, dtype=np.float32) cdef int i T[0]=0.0 for i in range(1,alen): T[i] = Tm[i] + (T[i-1] - Tm[i])**(-tau[i]) return T 

Cython float64 usando Jupyter %% magic. La función se puede utilizar directamente en Python :

 %%cython cdef extern from "math.h": double exp(double m) import cython import numpy as np cimport numpy as np from numpy cimport ndarray @cython.boundscheck(False) @cython.wraparound(False) @cython.infer_types(True) @cython.initializedcheck(False) def cy_loop(double[:] Tm,double[:] tau,int alen): cdef double[:] T=np.empty(alen) cdef int i T[0]=0.0 for i in range(1,alen): T[i] = Tm[i] + (T[i-1] - Tm[i])**(-tau[i]) return T 

Numba float64:

 @numba.jit("float64[:](float64[:], float64[:])", nopython=False, nogil=True) def calc_py_jit(Tm_, tau_): tt = np.empty(len(Tm_),dtype="float64") tt[0] = Tm_[0] for i in range(1, len(Tm_)): tt[i] = Tm_[i] - (tt[i-1] + Tm_[i])**(-tau_[i]) return tt[1:] 

Anexar a la lista . La solución no comstackda más rápida:

 def rec_py_loop(Tm,tau,alen): T = [Tm[0]] for i in range(1,alen): T.append(Tm[i] - (T[i-1] + Tm[i])**(-tau[i])) return np.array(T) 

Utilizando a.item ():

 def rec_numpy_loop_item(Tm_,tau_): n_ = len(Tm_) tt=np.empty(n_) Ti=tt.item Tis=tt.itemset Tmi=Tm_.item taui=tau_.item Tis(0,Tm_[0]) for i in range(1,n_): Tis(i,Tmi(i) - (Ti(i-1) + Tmi(i))**(-taui(i))) return tt[1:] 

np.fromiter ():

 def it(Tm,tau): T=Tm[0] i=0 while True: yield T i+=1 T=Tm[i] - (T + Tm[i])**(-tau[i]) def rec_numpy_iter(Tm,tau,alen): return np.fromiter(it(Tm,tau), np.float64, alen)[1:] 

Loop over Numpy (basado en la idea de Jaffe):

 def rec_numpy_loop(Tm,tau,alen): tt=np.empty(alen) tt[0]=Tm[0] for i in range(1,alen): tt[i] = Tm[i] - (tt[i-1] + Tm[i])**(-tau[i]) return tt[1:] 

Bucle sobre Numpy (código de Aronstef). En mi computadora, float64 es el tipo predeterminado para np.empty .

 def calc_py(Tm_, tau_): tt = np.empty(len(Tm_),dtype="float64") tt[0] = Tm_[0] for i in range(1, len(Tm_)): tt[i] = (Tm_[i] - (tt[i-1] + Tm_[i])**(-tau_[i])) return tt[1:] 

Pure C sin usar Python en absoluto. Versión del año 2016 (con función fabs ()):

 #include  #include  #include  #include  #include  double randn() { double u = rand(); if (u > 0.5) { return sqrt(-1.57079632679*log(1.0 - pow(2.0 * u - 1, 2))); } else { return -sqrt(-1.57079632679*log(1.0 - pow(1 - 2.0 * u,2))); } } void rec_pure_c(double *Tm, double *tau, int alen, double *T) { for (int i = 1; i < alen; i++) { T[i] = Tm[i] + pow(fabs(T[i - 1] - Tm[i]), (-tau[i])); } } int main() { int N = 100000; double *Tm= calloc(N, sizeof *Tm); double *tau = calloc(N, sizeof *tau); double *T = calloc(N, sizeof *T); double time = 0; double sumtime = 0; for (int i = 0; i < N; i++) { Tm[i] = randn(); tau[i] = randn(); } LARGE_INTEGER StartingTime, EndingTime, ElapsedMicroseconds; LARGE_INTEGER Frequency; for (int j = 0; j < 1000; j++) { for (int i = 0; i < 3; i++) { QueryPerformanceFrequency(&Frequency); QueryPerformanceCounter(&StartingTime); rec_pure_c(Tm, tau, N, T); QueryPerformanceCounter(&EndingTime); ElapsedMicroseconds.QuadPart = EndingTime.QuadPart - StartingTime.QuadPart; ElapsedMicroseconds.QuadPart *= 1000000; ElapsedMicroseconds.QuadPart /= Frequency.QuadPart; if (i == 0) time = (double)ElapsedMicroseconds.QuadPart / 1000; else { if (time > (double)ElapsedMicroseconds.QuadPart / 1000) time = (double)ElapsedMicroseconds.QuadPart / 1000; } } sumtime += time; } printf("1000 loops,best of 3: %.3f ms per loop\n",sumtime/1000); free(Tm); free(tau); free(T); } 

Fortran f2py. La función se puede utilizar desde Python . Versión del año 2016 (con función abs ()):

 subroutine rec_fortran(tm,tau,alen,result) integer*8, intent(in) :: alen real*8, dimension(alen), intent(in) :: tm real*8, dimension(alen), intent(in) :: tau real*8, dimension(alen) :: res real*8, dimension(alen), intent(out) :: result res(1)=0 do i=2,alen res(i) = tm(i) + (abs(res(i-1) - tm(i)))**(-tau(i)) end do result=res end subroutine rec_fortran 

Actualización: 21-10-2018 He corregido mi respuesta en base a los comentarios.

Es posible vectorizar operaciones en vectores siempre que el cálculo no sea recursivo. Debido a que una operación recursiva depende del valor calculado anterior, no es posible procesar en paralelo la operación. Esto por lo tanto no funciona:

 def calc_vect(Tm_, tau_): return Tm_[1:] - (Tm_[:-1] + Tm_[1:]) ** (-tau_[1:]) 

Dado que (el procesamiento en serie / un bucle) es necesario, el mejor rendimiento se obtiene al acercarse lo más posible al código de máquina optimizado, por lo que Numba y Cython son las mejores respuestas aquí.

Un enfoque Numba se puede lograr de la siguiente manera:

 init_string = """ from math import pow import numpy as np from numba import jit, float32 np.random.seed(0) n = 100000 Tm = np.cumsum(np.random.uniform(0.1, 1, size=n).astype('float32')) tau = np.random.uniform(-1, 0, size=n).astype('float32') def calc_python(Tm_, tau_): tt = np.empty(len(Tm_)) tt[0] = Tm_[0] for i in range(1, len(Tm_)): tt[i] = Tm_[i] - pow(tt[i-1] + Tm_[i], -tau_[i]) return tt @jit(float32[:](float32[:], float32[:]), nopython=False, nogil=True) def calc_numba(Tm_, tau_): tt = np.empty(len(Tm_)) tt[0] = Tm_[0] for i in range(1, len(Tm_)): tt[i] = Tm_[i] - pow(tt[i-1] + Tm_[i], -tau_[i]) return tt """ import timeit py_time = timeit.timeit('calc_python(Tm, tau)', init_string, number=100) numba_time = timeit.timeit('calc_numba(Tm, tau)', init_string, number=100) print("Python Solution: {}".format(py_time)) print("Numba Soltution: {}".format(numba_time)) 

Comparación de Timeit de las funciones de Python y Numba:

 Python Solution: 54.58057559299999 Numba Soltution: 1.1389029540000024 

Para aprovechar la respuesta de NPE, estoy de acuerdo en que tiene que haber un bucle en alguna parte ¿Quizás su objective es evitar la sobrecarga asociada con un bucle de Python for? En ese caso, numpy.fromiter vence un bucle for, pero solo un poco:

Usando la relación de recursión muy simple,

 x[i+1] = x[i] + 0.1 

yo obtengo

 #FOR LOOP def loopit(n): x = [0.0] for i in range(n-1): x.append(x[-1] + 0.1) return np.array(x) #FROMITER #define an iterator (a better way probably exists -- I'm a novice) def it(): x = 0.0 while True: yield x x += 0.1 #use the iterator with np.fromiter def fi_it(n): return np.fromiter(it(), np.float, n) %timeit -n 100 loopit(100000) #100 loops, best of 3: 31.7 ms per loop %timeit -n 100 fi_it(100000) #100 loops, best of 3: 18.6 ms per loop 

Curiosamente, la asignación previa de una matriz numpy produce una pérdida sustancial en el rendimiento. Esto es un misterio para mí, aunque supongo que debe haber más sobrecarga asociada con el acceso a un elemento de matriz que con la adición a una lista.

 def loopit(n): x = np.zeros(n) for i in range(n-1): x[i+1] = x[i] + 0.1 return x %timeit -n 100 loopit(100000) #100 loops, best of 3: 50.1 ms per loop