Cython通过数组广播加速循环

Question

摘要：

你们太棒了。。。我的真正代码工作起来了。我采纳了JoshAdel的建议，即：

1( 已将所有ndarray更改为键入的内存视图2( 手动展开所有numpy数组计算3( 为索引使用静态定义的无符号int4( 禁用边界检查和环绕

同时，也非常感谢韦德拉克的洞察力！

原始帖子：

我知道python做这些代码的速度非常慢：

import numpy as np
def func0():
    x = 0.
    for i in range(1000):
        x += 1.
    return

如果我把它改成Cython，它会更快：

import numpy as np
cimport numpy as np
def func0():
    cdef double x = 0.
    for i in range(1000):
        x += 1.
    return

结果是：

# Python
10000 loops, best of 3: 58.9 µs per loop
# Cython
10000000 loops, best of 3: 66.8 ns per loop

然而，现在我有了这种代码，它不是单个数字的循环，而是数组的循环。(是的……我正在解决一个PDE，所以会发生这种情况(。

我知道下面的例子很愚蠢，但只是试着了解计算类型：

纯python：

def func1():
    array1 = np.random.rand(50000, 4)
    array2 = np.random.rand(50000)
    for i in range(1000):
        array1[:, 0] += array2
        array1[:, 1] += array2
        array1[:, 2] += array2
        array1[:, 3] += array2
    return

Cython:

def func1():
    cdef np.ndarray[np.double_t, ndim=2] array1 = np.random.rand(50000, 4)
    cdef np.ndarray[np.double_t, ndim=1] array2 = np.random.rand(50000)
    for i in range(1000):
        array1[:, 0] += array2
        array1[:, 1] += array2
        array1[:, 2] += array2
        array1[:, 3] += array2
    return

而且几乎没有任何改善。同时，我知道Python不擅长处理这些巨大的循环，因为开销很大。

# Python
1 loops, best of 3: 299 ms per loop
# Cython
1 loops, best of 3: 300 ms per loop

关于我应该如何改进这类代码，有什么建议吗？

Answer 1

在这两个其他实现中，我已经尝试过

• 使用cython编译器指令删除numpy通常必须进行的一些检查
• 使用类型化的内存视图，这样我就可以指定内存布局(有时它们通常比旧的缓冲区接口更快(
• 展开循环，这样我们就不会使用numpy的切片机制：

否则，您只是通过cython使用numpy，而numpy已经在相当快的c代码中实现了这一点。

方法：

import numpy as np
cimport numpy as np
cimport cython
def func1():
    cdef np.ndarray[np.double_t, ndim=2] array1 = np.random.rand(50000, 4)
    cdef np.ndarray[np.double_t, ndim=1] array2 = np.random.rand(50000)
    cdef unsigned int i
    for i in range(1000):
        array1[:, 0] += array2
        array1[:, 1] += array2
        array1[:, 2] += array2
        array1[:, 3] += array2
    return
@cython.boundscheck(False)
@cython.wraparound(False)
def func2():
    cdef np.ndarray[np.double_t, ndim=2] array1 = np.random.rand(50000, 4)
    cdef np.ndarray[np.double_t, ndim=1] array2 = np.random.rand(50000)
    cdef unsigned int i, k
    for i in range(1000):
        for k in xrange(50000):
            array1[k, 0] += array2[k]
            array1[k, 1] += array2[k]
            array1[k, 2] += array2[k]
            array1[k, 3] += array2[k]
    return

@cython.boundscheck(False)
@cython.wraparound(False)
def func3():
    cdef np.ndarray[np.double_t, ndim=2] array1 = np.random.rand(50000, 4)
    cdef np.ndarray[np.double_t, ndim=1] array2 = np.random.rand(50000)
    cdef np.double_t[::1] a2 = array2
    cdef np.double_t[:,::1] a1 = array1
    cdef unsigned int i, k
    for i in range(1000):
        for k in xrange(50000):
            a1[k, 0] += a2[k]
            a1[k, 1] += a2[k]
            a1[k, 2] += a2[k]
            a1[k, 3] += a2[k]
    return

计时(在我的机器上——编译器和硬件当然可能会影响计时(：

• 纯numpy：1个循环，3的最佳值：每个循环419毫秒
• 输入i：1循环的原始cython函数，每个循环3：428毫秒的最佳值
• func2:1个循环，每个循环的最佳时间为:336毫秒
• func3:1个循环，每个循环3:206毫秒的最佳值

Answer 2

很大程度上，你错了。Python在这些类型的循环中非常出色。

这是因为他们是重量级。正如Josh Adel所展示的，一个有效的纯C变体只能实现2到4倍的加速。

最好的办法是寻找几种类型的低效率：

• 许多切片(例如切片许多小行(

  这最好通过手动迭代并移动到Cython或Numba(或另一个JITing编译器(来解决

• 完整的Numpy数组计算，只是不够快

  Numba和Parakeet可以提供帮助，但您确实想要Theano或numexpr

• 更复杂，非JITable，不可表达的慢速

  Cython或C是您唯一的选择，因此请转到它们

另见答案的第一部分。