Cython没有提高速度

Question

我正在尝试定义一个包含内部循环以模拟积分的函数。

问题在于速度。在我的计算机上，评估该函数一次可能需要长达30秒的时间。由于我的最终目标是将该函数最小化，因此一些额外的速度会很不错。

因此，我尝试使用Cython来加速，但我肯定犯了严重的错误（可能有许多！）。按照Cython文档的指示，我尝试对变量进行类型定义。但是，在这样做之后，代码的速度与纯Python一样慢。这似乎很奇怪。

以下是我的代码：

import numpy as np 
cimport cython
cimport numpy as np
import minuit

data = np.genfromtxt('q6data.csv', usecols = np.arange(1, 24, 1), delimiter = ',')  

cdef int ns    = 1000                 # Number of simulation draws
cdef int K     = 5                    # Number of observed characteristics, including            constant
cdef int J     = len(data[:,1])       # Number of products, including outside
cdef double tol   = 0.0001            # Inner GMM loop tolerance
nu = np.random.normal(0, 1, (6, ns))  # ns random deviates

@cython.boundscheck(False)
@cython.wraparound(False)


def S(np.ndarray[double, ndim=1] delta, double s1, double s2, double s3, double s4,  double s5, double a):
    """Computes the simulated integrals, one for each good.
    Parameters: delta is an array of length J containing mean product specific utility levels
    Returns: Numpy array with length J."""
    cdef np.ndarray[double, ndim=2] mu_ij = np.dot((data[:,2:7]*np.array([s1, s2, s3, s4, s5])), nu[1:K+1,:])
    cdef np.ndarray[double, ndim=2] mu_y  = a * np.log(np.exp(data[:,21].reshape(J,1) +  data[:,22].reshape(J,1)*nu[0,:].reshape(1, ns)) - data[:,7].reshape(J,1))
    cdef np.ndarray[double, ndim=2] V = delta.reshape(J,1) + mu_ij + mu_y
    cdef np.ndarray[double, ndim=2] exp_vi = np.exp(V)
    cdef np.ndarray[double, ndim=2] P_i = (1.0 / np.sum(exp_vi[np.where(data[:,1] == 71)], 0)) * exp_vi[np.where(data[:,1] == 71)] 
    cdef int yrs = 19
    cdef int yr
    for yr in xrange(yrs):
        P_yr = (1.0 / np.sum(exp_vi[np.where(data[:,1]== (yr + 72))], 0)) *   exp_vi[np.where(data[:,1] == (yr + 72))]
        P_i  = np.concatenate((P_i, P_yr)) 
    cdef np.ndarray[double, ndim=1] S = np.zeros(dtype = "d", shape = J)
    cdef int j
    for j in xrange(ns):
        S += P_i[:,j]
    return (1.0 / ns) * S

def d_infty(np.ndarray[double, ndim=1] x, np.ndarray[double, ndim=1] y):
    """Sup norm."""
    return np.max(np.abs(x - y)) 

def T(np.ndarray[double, ndim=1] delta_exp, double s1, double s2, double s3, double s4,  double s5, double a):
    """The contraction operator.  This function takes the parameters and the exponential
    of the starting value of delta and returns the fixed point.""" 
    cdef int iter = 0
    cdef int maxiter = 200
    cdef int i
    for i in xrange(maxiter): 
        delta1_exp = delta_exp * (data[:, 8] / S(np.log(delta_exp), s1, s2, s3, s4, s5, a))                    
        print i
        if d_infty(delta_exp, delta1_exp) < tol:                                       
            break
        delta_exp = delta1_exp
    return np.log(delta1_exp)


def Q(double s1, double s2, double s3, double s4, double s5, double a):
    """GMM objective function."""  
    cdef np.ndarray[double, ndim=1] delta0_exp = np.exp(data[:,10])                                                     
    cdef np.ndarray[double, ndim=1] delta1 = T(delta0_exp, s1, s2, s3, s4, s5, a)
    delta1[np.where(data[:,10]==0)] = np.zeros(len(np.where(data[:,10]==0)))            
    cdef np.ndarray[double, ndim=1] xi =  delta1 - (np.dot(data[:,2:7],   np.linalg.lstsq(data[:,2:7], delta1)[0]))   
    cdef np.ndarray[double, ndim=2] g_J = xi.reshape(J,1) * data[:,11:21]
    cdef np.ndarray[double, ndim=1] G_J = (1.0 / J) * np.sum(g_J, 0) 
    return np.sqrt(np.dot(G_J, G_J))

我已经对代码进行了分析，看起来是函数S，即积分模拟器，导致了性能问题。无论如何，我本来预期通过对变量进行输入会有一些速度上的提升。但是由于没有产生任何提升，我认为我可能犯了一些基本错误。

有没有人在Cython代码中看到可能导致这种结果的明显错误？

哦，由于我很新手，肯定有很多不好的风格和减慢代码速度的东西。如果你有时间，可以告诉我这些问题的正确方式。

- Randall J

---

data有多大？你能把它分成若干块，其中第一列=={0,1,...18,71}，并在内部循环之外进行吗？ - denis

因此，数据有2237行和21列。这些行可以分成子集，对应于1971年到1990年的年份。对我来说，难点在于找到一种更快的方法来执行该操作，因为我必须仅将给定年份的每个条目除以对应于该年份的条目总和。肯定有比我实现的更好的方法，但我未能找到它！ - Randall J

Answer 1

Cython可以生成一个HTML文件来帮助解决这个问题：

cython -a MODULE.py

这显示了源代码的每一行的颜色从白色到各种黄色的渐变。颜色越深，表示该行上仍在执行更多动态的Python操作。对于每一行包含某些黄色的地方，需要添加更多静态类型声明。

当我这样做时，我喜欢将我遇到麻烦的源代码部分拆分成许多单独的行，每个表达式或运算符一个，以获得最细粒度的视图。

如果没有这个功能，很容易忽视一些变量、函数调用或运算符的静态类型声明。(例如，索引运算符x[y]仍然是完全动态的Python操作，除非你声明它不是)

- Jonathan Hartley

---

谢谢这个提示。看起来非常有用。我已经搜索了如何将其添加到我的setup.py中，但没有运气。有人有关于如何完成这个任务的建议吗？ - Randall J

你不需要将它添加到你的setup.py文件中；你只需要在命令行中运行它。 - Peter Wang

我喜欢为像这样的小随机命令创建一个小的Makefile。然后，“make cython”将执行cythoning（通常使用您的setup.py），然后运行上面的命令。在Windows上，我喜欢安装Cygwin并将bin目录放在我的Path中，这样我就可以从DOS提示符中访问“make”和其他类Unix命令了。 - Jonathan Hartley

谢谢你提供的索引提示 - 它真的是一个救命稻草。有关此内容的相关信息可以在http://docs.cython.org/src/userguide/numpy_tutorial.html#efficient-indexing找到。 - Gerald Senarclens de Grancy

Answer 2

Cython并不会自动提高性能，你需要了解它的内部结构并检查生成的C代码。

特别是如果你想提高循环的性能，你需要避免在循环中调用Python函数，而这种情况在本例中经常发生（所有的 np. 调用都是Python调用，切片和其他操作也一样）。

请参阅此页面了解使用Cython进行性能优化的一般指南（当进行优化时，-a开关确实非常方便），以及此页面了解优化numpy代码的具体方法。

- Luper Rouch

---

我认为 'cimport numpy' 告诉 Cython 使用 C 函数来代替 Python 函数处理 numpy。这意味着 'np' 调用不是 Python 调用，而其他所有内容都是。我是否理解有误？ - Jonathan Hartley

是的，你误解了，cimport 中的 numpy 声明仅在 cdef 语句中使用，并且还提供了访问 numpy 的 C API 的权限，该 API 与常规的 Python API 非常不同：http://docs.scipy.org/doc/numpy/reference/c-api.html - Luper Rouch

Answer 3

你可以通过更多地使用Numpy的功能来加速你的代码。

例如：

cdef np.ndarray[double, ndim=1] S = np.zeros(dtype = "d", shape = J)
cdef int j
for j in xrange(ns):
    S += P_i[:,j]

如果这样写会更快，更易读：

S = P_i.sum(axis=1)

您也重复了一些计算，因此所需时间是必要的两倍。例如

np.where(data[:,1]==(yr + 72))

此计算可以仅计算一次并存储在变量中以便重复使用。

您还进行了大量的重塑和切片操作：从一开始就将变量保持在更简单的格式中可能会有所帮助。如果可能的话，您的代码会更清晰，优化也会更加明显。

- Eric O Lebigot

---

感谢您关于计算 S 的建议，这样更加简洁！此外，为切片数据数组定义变量确实加速了一些操作。非常感谢！ - Randall J

实际上，我不同意你的开场白。通过在C中使用更直接的处理方式并减少对numpy例程的调用，代码将变得更快。此外，根据cython文档，通过切片访问数组（P_i [:，j]）是低效的。http://docs.cython.org/src/tutorial/numpy.html - DaveP

@DaveP：我不确定我是否理解你的意思：我的观点正是如你所说，“切片是低效的”。许多Numpy函数允许您绕过任何需要切片的操作，就像我第一个示例中所示。此外，由于Numpy在C中进行了大量“直接处理”，因此很难编写比其等效的Numpy函数运行更快的代码。 - Eric O. Lebigot

正如我提供的链接所解释的那样，类型P_i[:, j]的索引速度较慢（这是我们希望避免的切片），而编写显式循环并使用类型P_i[k, j]的索引速度较快。据我所知，这种差异是由于第一种情况需要函数调用，而第二种情况可以生成内联代码。当您生成内联代码时，编译器可以更轻松地将整个算法优化为一个，而不是将函数作为瓶颈。 - DaveP

@DaveP：感谢您的解释。因此，我们都同意要避免在i上使用P [:，i]的单个循环。您建议使用完整循环和Cython，而我建议根本不使用循环和Numpy。您的方法可能是最快的，而我的更简单/更短的解决方案除了函数调用（希望与计算时间相比可以忽略不计）之外，可能几乎一样快。 - Eric O. Lebigot

Answer 4

使用性能分析器可以帮助您确定哪一部分速度较慢吗？我喜欢使用标准库的性能分析器来运行程序：

python -O -m cProfile -o profile.out MYAPP.py

然后在“RunSnakeRun”图形界面中查看输出：

runsnake profile.out

您可以从这里安装RunSnakeRun： http://www.vrplumber.com/programming/runsnakerun/ 下面是RunSnakeRun的截图：

- Jonathan Hartley

---

这个问题是关于Cython代码而不是Python代码的，所以你的回答并没有直接帮助到。然而，关于性能分析的观点是正确的，在Cython中有关性能分析的详细信息可以在http://docs.cython.org/src/tutorial/profiling_tutorial.html找到。 - DaveP

+1 对于性能分析：性能分析确实是优化的第一步。还要感谢您提到runsnake：它看起来非常有前途！ - Eric O. Lebigot

Answer 5

采纳这里给出的建议，我花了更多时间对上述代码进行分析。为了希望能够使事情变得更清晰一些，我定义了

我对代码进行了更详细的分析，并且更好地了解了哪些部分是最慢的。此外，我还定义了

X = data[:, 2:7]
m_y = data[:, 21].reshape(J,1)
sigma_y = 1.0
price = data[:, 7].reshape(J, 1)
shares_data = data[:,8]

接下来的这几行代码占用了大部分的总时间。

mu_ij = np.dot((X*np.array([s1, s2, s3, s4, s5])), nu[1:K+1,:])
mu_y  = a * np.log(np.exp(m_y + sigma_y*nu[0,:].reshape(1,ns)) - price)
V = delta.reshape(J,1) + mu_ij + mu_y
exp_vi = np.exp(V)
P_i = (1.0 / np.sum(exp_vi[np.where(data[:,1]==71)], 0)) *  exp_vi[np.where(data[:,1]==71)] 
for yr in xarange(19):
    P_yr = (1.0 / np.sum(exp_vi[np.where(data[:,1]==yr)], 0)) * exp_vi[np.where(data[:,1]==yr)]
P_i  = np.concatenate((P_i, P_yr))

我有一种印象，这是一种过于繁琐的方式来实现我的目标。我希望有人能够提供一些关于如何加速这些代码行的建议。也许我错过了一些Numpy的功能？如果这个问题没有足够明确的规定，以便您提供帮助，我很乐意提供更多关于我的问题背景的细节。谢谢！

- Randall J

---

请不要用更多的问题来“回答”自己的问题。相反，如果相关，请更新您的原始问题，否则提出一个新问题。 - John Machin

抱歉，我之前不清楚更新问题的正确流程。感谢您的澄清。 - Randall J

Answer 6

仅在您于11月28日发表评论后分割数据：

import sys
import time
import numpy as np

def splitdata( data, n, start=1971 ):
    """ split data into n pieces, where col 1 == start .. start + n """
        # not fancy, fast enough for small n
    split = n * [None]
    for j in range(n):
        split[j] = data[ data[:,1] == start + j ]
    return split  # [ arrays: col1 0, col1 1 ... ]

#...........................................................................
N = 2237
ncol = 21
start = 1971
n = 20
seed = 1
exec "\n".join( sys.argv[1:] )  # run this.py N= ...
np.set_printoptions( 2, threshold=100, suppress=True )  # .2f
np.random.seed(seed)

print "N=%d  ncol=%d  n=%d" % (N, ncol, n)
data = np.random.uniform( start, start + n, (N,ncol) )
data[:,1] = data[:,1].round()
t0 = time.time()

split = splitdata( data, n, start )  # n pieces

print "time: %.0f us splitdata" % ((time.time() - t0) * 1e6)
for y, yeardata in enumerate(split):
    print "%d  %d  %g" % (start + y, len(yeardata), yeardata[:,0].sum())

-->

time: 27632 us splitdata  # old mac ppc
1971  69  136638
1972  138  273292
...

- denis

Answer 7

“根本错误”在于您期望Python在长循环中具有良好的性能。它是一种解释性语言，切换实现和ctyping对此毫无作用。有一些用于快速计算的数值Python库，大多数是用C编写的。例如，如果您已经使用numpy进行数组操作，为什么不进一步使用scipy进行高级数学计算呢？这将提高可读性和速度。

- atomizer

---

Cython不是Python的实现，它是一种类似Python的语言，可以编译成C。同时，Python也不是解释型语言。 - Luper Rouch

我有一些使用Scipy的经验，但我不确定在这里如何利用它来提高性能。你有任何具体的建议，在哪里我可以调用它的一些功能吗？ - Randall J

Cython的整个意义在于它专门设计用于加速Python的性能，在许多情况下可能实现类似于C的性能。因此，在这个问题中没有“根本性错误”。 - DaveP

说 Python 不是解释型语言虽然在技术上是正确的，但没有给出上下文是会误导人的。当你运行一个 Python 程序时，“解释器”将把你的源代码转换为字节码，并执行它。从黑盒子的角度来看，它的行为与传统的解释器完全相同，包括性能特征。 - Jonathan Hartley

反思之后，我并不是对性能特征的所有说法都正确。但其中有一些是正确的！ :-) - Jonathan Hartley