python - 递归方程的高效python方法

Question

我正在尝试优化一段代码中的循环。我认为以更 NumPy 的方式编写它会使它更快，但现在更慢了!方程将长度为 n 的 numpy.array vec 作为输入:

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 [秒]

5.66514706612 [s]

最后我决定将整个函数插入我的循环中，获得 3 倍的性能提升。我真的在寻找 100 倍的性能提升，但不知道还能做什么。这是我的最终功能:

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

当线性算法被二次算法取代时会发生这种情况:无论执行多快，更好的算法总是获胜(对于足够大的问题)。

很明显，f 以线性时间运行，而 f2 以二次时间运行，因为这是矩阵向量点积的时间复杂度。

log-log 图清楚地显示了运行时间的差异(线性指 f，二次方指 f2):

可以解释绿线最右边的部分(即，当它不是一条直线时)，因为 numpy 函数通常有很高的开销，对于不是很小但支配运行时间的数组来说可以忽略不计当它们很小的时候。

---

在已经使用快速算法的 Python 中加速代码的“标准”方法是获取已编译代码并编写扩展。 Cython允许您通过使用一些类型注释来注释 Python 源代码来做到这一点，并且它理解 numpy 数组。

通过告诉 Cython vec 是一个 double 组，aux 是一个 double 组，i 是一个整数，它能够生成一个 C对我来说快 400 倍的扩展。

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

如果您碰巧使用 IPython ，您可以只运行 %load_ext cythonmagic，然后将该函数复制到以 %%cython 行为前缀的单元格中进行尝试。 Cython documentation 中解释了构建和编译它的其他方法。 .顺便说一句，IPython 还允许您通过在语句 to time 之前编写 %timeit 来计时代码，这真的很方便。

一个完全不同的选择是使用 PyPy ，一个带有 JIT 的 Python 2.7 实现，并具有一些基本的 numpy 支持。它可以通过将 import numpy 替换为 import numpy 来运行这个小片段，但它可能无法运行您的整个程序。它比 Cython 慢一点，但它不需要编译器，也不需要注释代码。