Efektywna metoda pythonowa dla równań rekurencyjnych

Question

Próbuję zoptymalizować pętlę, którą mam w kawałku mojego kodu. Pomyślałem, że pisanie tego w bardziej odręczny sposób spowoduje, że będzie to szybsze, ale teraz jest wolniejsze! równania wprowadzany jest vec numpy.array o długości n:

from numpy import * 

def f(vec): 
    n=len(vec) 
    aux=0 
    for i in range(n): 
     aux = aux + (1- aux)*vec[i] 
    return aux 

def f2(vec): 
    n=len(vec) 
    G=tril(array([-vec]*n),-1)+1    #numpy way! 
    aux=dot(G.prod(1),vec) 
    return aux 


if __name__ == '__main__': 
    import timeit 
    print(timeit.timeit("f(ones(225)+4)", setup="from __main__ import f\nfrom numpy import ones",number=1000)) 
    print(timeit.timeit("f2(ones(225)+4)", setup="from __main__ import f2\nfrom numpy import ones,tril,dot",number=1000)) 

0,429496049881 [s]

5,66514706612 [s]

końcu postanowił inserte całą funkcję mi pętli uzyskanie 3-krotnego zwiększenia wydajności. Naprawdę szukam 100-krotnego zwiększenia wydajności, ale nie wiem, co jeszcze można zrobić. To moja ostateczna funkcja:

def CALC_PROB_LOC2(int nSectors, int nZones,double[:] beta, double[:] thetaLoc,np.ndarray[double, ndim=2] h, np.ndarray[double, ndim=2] p, np.ndarray[np.float64_t, ndim=3] U_nij, np.ndarray[double, ndim=2] A_ni): 
    cdef np.ndarray[double, ndim=3] Pr_nij =np.zeros((nSectors,nZones,nZones),dtype="d") 
    cdef np.ndarray[double, ndim=2] U_ni =np.zeros((nSectors,nZones),dtype="d") 
    #cdef np.ndarray[np.float64_t, ndim=1] A_ni_pos 
    cdef Py_ssize_t n,i,opt 
    cdef int aux_bool,options 
    cdef np.ndarray[np.float64_t, ndim=1] prob,attractor,optionCosts 
    cdef np.ndarray[np.float64_t, ndim=1] eq23,utilities 
    cdef double disu 
    cdef double eq22 
    cdef double aux17 
    for n in range(nSectors): 
     aux_bool=1 
     if n in [0,2,9,10,11,12,13,14,18,19,20]: 
      for i in xrange(nZones): 
       U_ni[n,i]=p[n,i]+h[n,i] 
       Pr_nij[n,i,i]=1 
      aux_bool=0 
     if aux_bool==1: 
      if beta[n]<=0: 
       for i in xrange(nZones): 
        U_ni[n,i]=U_nij[n,i,i] 
      else: 
       A_ni_pos=A_ni[n,:]>0 
       options=len(A_ni[n,:][A_ni_pos]) 
       attractor=A_ni[n,:][A_ni_pos] 
       if options>0: 
        for i in xrange(nZones): 
         optionCosts=U_nij[n,i,A_ni_pos] 
         disu=0 
         eq22=0 
         aux17=0 
         prob=np.ones(options)/options #default value 
         if beta[n]==0: 
          Pr_nij[n,i,A_ni_pos],U_ni[n,i]= prob,0 
         if options==1: 
          Pr_nij[n,i,A_ni_pos],U_ni[n,i]= prob,optionCosts 
         else: 
          if thetaLoc[n]<=0: 
           cmin=1 
          else: 
           cmin=(optionCosts**thetaLoc[n]).min() 
           if cmin==0: 
            cmin=100 
          utilities=optionCosts/cmin 
          eq23=-beta[n]*utilities 
          eq23=np.exp(eq23) 
          aux17=np.dot(attractor,eq23) 
          if aux17==0: 
           Pr_nij[n,i,A_ni_pos],U_ni[n,i]= 0*prob,0 
          else: 
           for opt in range(options): 
            eq22=eq22+(1-eq22)*eq23[opt] 
           prob=attractor*eq23/aux17 
           disu=cmin*(-np.log(eq22)/beta[n]) 
           Pr_nij[n,i,A_ni_pos],U_ni[n,i]= prob,disu 


    return Pr_nij,U_ni 

Answer 1

To co się dzieje, gdy liniowy algorytmy zostanie zastąpiony przez jednego kwadratowego: Bez względu na to, jak szybko to jest wykonywane, tym lepszy algorytm zawsze wygrywa (dla problemu wystarczająco duży).

Jest całkiem jasne, że f działa w czasie liniowym, a f2 działa w trybie kwadratu, ponieważ jest to złożoność czasowa iloczynu macierzy-wektora.

Działka logarytmiczny wyraźnie pokazuje różnicę w czasie jazdy (liniowy dotyczy f, quadractic do f2):

Prawo część najbardziej zielonej linii (czyli gdy nie zachowuje się jak linia prosta) można wyjaśnić, ponieważ funkcje numpy zwykle mają wysokie koszty ogólne, które są nieistotne dla tablic, które nie są małe, lecz dominują w czasie działania, gdy są małe.

---

„Standardowe” sposobem na przyspieszenie kodu w Pythonie, który jest już za pomocą szybki algorytm jest dotarcie do skompilowanego kodu i napisać rozszerzenie. Cython pozwala to zrobić, dodając adnotacje do kodu źródłowego Pythona za pomocą kilku adnotacji typu i rozumie numpy arrays.

Mówiąc Cython że vec jest tablicą podwaja, aux jest podwójne i i liczbą całkowitą, to jest w stanie generować rozszerzenie C, który jest 400x szybciej dla mnie.

def f(double[:] vec): 
    n = len(vec) 
    cdef double aux = 0 
    cdef int i 
    for i in range(n): 
     aux = aux + (1- aux)*vec[i] 
    return aux 

Jeśli zdarzy się przy użyciu IPython ty, możesz po prostu uruchomić %load_ext cythonmagic i skopiuj tę funkcję do komórki poprzedzony linii %%cython go wypróbować. Inne metody budowania i kompilowania tego są wyjaśnione w Cython documentation. Nawiasem mówiąc, IPython pozwala również na timeit kodu, pisząc %timeit przed stwierdzeniem czasu, jest to naprawdę przydatne.

Zupełnie inną opcją jest użycie PyPy, implementacji Pythona 2.7 dostarczanej z JIT i posiadającej pewną podstawową obsługę numpy. Może uruchomić ten mały fragment zastępując import numpypy dla import numpy, ale możliwe, że nie będzie w stanie uruchomić całego programu. Jest odrobinę wolniejszy niż Cython, ale nie wymaga kompilatora ani adnotacji kodu.