递归方程的高效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 [s]

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 运行在二次时间，因为这是矩阵向量点积的时间复杂度。

对数图清楚地显示了运行时间的差异（线性表示f，二次表示f2）：

绿线的最右边部分（即，当它不表现为直线时）可以解释，因为 numpy 函数通常具有很高的开销，这对于不是很小但占主导地位的数组来说可以忽略不计运行时间当他们还小的时候。

---

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

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

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 还可以让你在时间语句之前写%timeit 来计时代码，真的很方便。

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