有没有办法提高 cython 代码的速度答案

【问题标题】：Is there a way of improving speed of cython code有没有办法提高 cython 代码的速度
【发布时间】：2019-04-11 15:03:17
【问题描述】：

我正在寻找加快 python numpy 代码的速度：

def fun_np(m,data):
a, b, c = data[:,0], data[:,1], data[:,2] 

M = len(data[:,0]) 
n = round((m+1)*(m+2)*(m+3)/6) 
u =np.zeros((M,n))

C = 0
for i in range(0,m+1):
    for j in range(0,i+1):
        for k in range(0,j+1):
            if ((i-j)!=0):
                u[:,C] = (j-k)*(a)**(i-j)*(b)**(j-k-1)*(c)**k

        C=C+1  
return u

对应的cython代码如下：

%%cython 
import numpy as np
cimport numpy as np
from cython import wraparound, boundscheck, nonecheck

@boundscheck(False)
@wraparound(False)
@nonecheck(False)

cpdef fun_cyt(int m,np.ndarray[np.float64_t, ndim=2] data):

cdef:
    np.ndarray[np.float64_t, ndim=1] a = data[:,0]
    np.ndarray[np.float64_t, ndim=1] b = data[:,1]
    np.ndarray[np.float64_t, ndim=1] c = data[:,2]
    int M, n
    Py_ssize_t i, j, k, s
M = len(data[:,0]) 
n = round((m+1)*(m+2)*(m+3)/6)   
cdef np.ndarray[np.float64_t, ndim=2]  u = np.zeros((M,n), dtype=np.float64)

cdef int C = 0
for i in range(m+1): #range(0,m+1):
    for j in range(i+1):
        for k in range(j+1):
            for s in range(M):
                if (i-j)!=0:
                    u[s,C] = (j-k)*(a[s])**(i-j)*(b[s])**(j-k-1)*(c[s])**k

            C=C+1
return u

这里是时间

z = np.random.randn(6000, 3); m=20;

%timeit fun_np(m,z);

结果：每个循环 1.97 秒 ± 11.2 毫秒（平均值 ± 标准偏差，7 次运行，每次 1 个循环）

%timeit fun_cyt(m,z);

结果：每个循环 1.91 秒 ± 12.7 毫秒（平均值 ± 标准偏差，7 次运行，每个循环 1 个）

如您所见，numpy 和 cython 代码之间没有显着的速度。如果可能的话，如果您能帮助优化 cython 代码，我将不胜感激。

cython 代码的注释 html html

【问题讨论】：

如果您正在寻找代码审查，请在 codereview.stackexchange.com 上试试运气，因为您的问题不在 StackOverflow 上。
@HIlle 不 - 没有理由在这里偏离主题。我想说它属于the help“特定编程问题”中列出的第一个要点。事实上，它也可能在代码审查中成为主题，但这并不意味着它在这里是离题的。
@ForBonder 你看过带注释的 html 输出（来自运行 cython -a yourfile.pyx）吗？如果您遗漏了什么，它通常可以为您提供线索
@DavidW 谢谢，我查看了带注释的 html，但没有提前。我已经添加了带注释的 html 的链接。
看起来不错。这是 Cython 非常擅长的问题，而且看起来你做了正确的事情，所以你没有看到差异有点令人惊讶。

标签： python performance numpy optimization cython

【解决方案1】：

正如 cmets 中已经提到的，您可以尝试使用 numba。我会建议进一步并行化循环：

from numba import prange, jit

@jit(nopython=True, parallel=True)
def fun_numba(m,data):
    a, b, c = data[:,0], data[:,1], data[:,2] 

    M = len(data[:,0]) 
    n = round((m+1)*(m+2)*(m+3)/6) 
    u = np.zeros((M,n))

    C = 0
    for i in range(0,m+1):
        for j in range(0,i+1):
            for k in prange(0,j+1):
                if ((i-j)!=0):
                    u[:,C] = (j-k)*(a)**(i-j)*(b)**(j-k-1)*(c)**k

            C=C+1  
    return u

在我的机器上给我：

In [11]: %timeit fun_np(m,z)                                                                         
642 ms ± 4.13 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [12]: %timeit fun_numba(m,z)                                                                      
101 ms ± 7.15 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

【讨论】：

【解决方案2】：

非常有趣的例子！大多数操作都在 6000 个元素向量上。在大向量功率、乘法和加法方面，Cython 真的不能比 numpy 快。通过在 Cython 中实现这一点，您可能会像 numpy 一样快，甚至可以通过消除一些 numpy 开销来获得 10% 到 20% 的速度。

但是，还有其他方法可以加快计算速度。向量操作是对数据向量的三列的操作，并且您写入输出向量的列。默认情况下，numpy 数组具有行优先排序，即在内存中，行在内存中是连续的。对于这里完成的操作，这很糟糕。延伸阅读：https://en.wikipedia.org/wiki/Row-_and_column-major_order.

这两个函数基本相同，如果输出向量的创建发生在函数之外，它们将是相同的。

请注意以下几点：我将 u[:,C] = ... 替换为 u[:,C] +=，因为否则结果仅由 k=j 定义，因此始终为 0。我不知道这些计算的意义是什么，但可能不是这样。

import numpy as np
def fun_np(m,data):
    a, b, c = data[:,0], data[:,1], data[:,2] 

    M = len(data[:,0]) 
    n = round((m+1)*(m+2)*(m+3)/6) 
    u = np.zeros((M,n))

    C = 0
    for i in range(0,m+1):
        for j in range(0,i+1):
            for k in range(0,j+1):
                if ((i-j)!=0):
                    u[:,C] += (j-k)*(a)**(i-j)*(b)**(j-k-1)*(c)**k

            C=C+1  
    return u

def fun_npF(m,data):
    a, b, c = data[:,0], data[:,1], data[:,2] 

    M = len(data[:,0]) 
    n = round((m+1)*(m+2)*(m+3)/6) 
    u = np.zeros((M,n),order='F')

    C = 0
    for i in range(0,m+1):
        for j in range(0,i+1):
            for k in range(0,j+1):
                if ((i-j)!=0):
                    u[:,C] += (j-k)*(a)**(i-j)*(b)**(j-k-1)*(c)**k

            C=C+1  
    return u

z = np.random.randn(6000, 3); m=20;
print("Numpy Row-major")
%timeit fun_np(m,z);

# Fortran order, because vector operations on columns
print("Numpy Column-major")
zF = np.asarray(z.copy(),order='F')
%timeit fun_npT(m,zF);

# Check if output the same
diff = (max(np.ravel(abs(fun_np(m,z)-fun_npF(m,zF)))))
max_rm = (max(np.ravel(abs(fun_np(m,z)))))
max_cm = (max(np.ravel(abs(fun_npF(m,zF)))))
print("Dffference: %f, Max value Row-major: %f, Max value Column-major: %f"%(diff, max_rm, max_cm))

这给了我

Numpy Row-major
1.64 s ± 12.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
Numpy Column-major
16 ms ± 345 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Dffference: 0.000000, Max value Row-major: 196526643123.792450, Max value Column-major: 196526643123.792450

在考虑将“if”放在哪里并将其与 Cython 结合使用时，您可以获得更多收益，但我猜也只有大约 10% 到 20%。

【讨论】：