Numba代码比Cython替代方案快得多。

Question

我必须使用cython编写一个小模拟程序，通常可以使用numba进行加速。 但由于numba不支持我想要使用的scipy函数来修改函数，所以我不得不做出这个转变。

因此，我已经将我的模拟程序翻译成了cython，与numba相比，一切都变得非常缓慢。也许我的cython代码存在瓶颈，我没有看到。

我的cython代码：

import numpy as np
cimport numpy as cnp
cimport cython

cnp.import_array()

@cython.boundscheck(False)
@cython.wraparound(False)
def simulation(int N=10000, int N_save=10000):

    cdef cnp.ndarray[double] x = np.empty(N_save, dtype=np.double)
    cdef cnp.ndarray[double] r = np.random.standard_normal(size=N_save)
    
    cdef fs = int(N / N_save)
    cdef xold = 0
    x[0] = 0

    for i in range(1, N):

        if i%N_save == 0:
            r = np.random.standard_normal(size=N_save)

        xnew = xold + r[i%N_save]

        if (i % fs) == 0:
            x[int(i / fs)] = xnew

        xold = xnew

    return x

我用于编译的代码：

from setuptools import setup
from Cython.Build import cythonize
import numpy as np

setup(
    ext_modules=cythonize(
        "test.pyx",
        compiler_directives={'language_level' : "3"},
    ), 
    include_dirs=[np.get_include()]
)

用Numba加速的相同代码：

@nb.jit(nopython=True, fastmath=True)
def simulation_nb(N=10000, N_save=10000):

    x = np.zeros(N_save)
    r = np.random.standard_normal(size=N_save)
    
    fs = int(N / N_save)
    xold = 0
    x[0] = 0


    for i in range(1, N):

        if i%N_save == 0:
            r = np.random.standard_normal(size=N_save)

        xnew = xold + r[i%N_save]

        if (i % fs) == 0:
            x[int(i / fs)] = xnew

        xold = xnew

    return x

< p >使用模拟数据（N=1000000，N_save=10000）的基准测试结果如下：< /p> < ol >

Cython：183ms

Numba：31.4ms

原生Python：217ms

为什么我的Cython代码几乎和原生Python一样慢？

- Lars Dingeldein

---

当Cython运行较慢时，这可能是由于类型转换的原因，而且可能加剧了缺少类型注释的问题。另外，如果你在Cython中使用C数据结构，那么比起使用Python数据结构，速度会更快。 - Ameya

cython -a 会生成一个带注释的 HTML 文件，它将为您提供有关可能效率低下的提示。 - DavidW

Answer 1

你的主要问题在于没有为所有变量提供Cython编译器类型。这会导致你的Cython代码效率低下，因为它会与Python解释器进行许多交互。例如，变量i、fs、xnew和xold的类型对Cython编译器来说是未知的。正如评论中已经提到的那样，可以通过使用--annotated选项进行编译来轻松观察到这一点。

为了公平比较，让我们在我的机器上测试一下你的numba代码的时间：

In [7]: %timeit simulation_nb(N=1000000, N_save=10000)
24.7 ms ± 177 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)

通过使用类型内存视图（typed memoryviews）而不是 cnp.ndarrays（它们更快且语法更清晰），使用 C 除法代替调用 Python 的 int() 函数，以及为所有变量提供类型，您的 Cython 代码可以看起来像这样:

%%cython -a -c=-O3 -c=-march=native

cimport numpy as np
import numpy as np
cimport cython

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
def simulation_cy1(int N=10000, int N_save=10000):
    cdef double[::1] x = np.empty(N_save, dtype=np.float64)
    cdef double[::1] r = np.random.standard_normal(size=N_save)
    cdef int fs = int(N / N_save)
    cdef int i
    cdef double xnew, xold = 0
    x[0] = 0

    for i in range(1, N):
        if i%N_save == 0:
            r = np.random.standard_normal(size=N_save)
        xnew = xold + r[i%N_save]
        if (i % fs) == 0:
            x[i / fs] = xnew
        xold = xnew
    return np.asarray(x)

在我的计算机上，这个程序的表现与Numba类似。

In [10]: %timeit simulation_cy1(N=1000000, N_save=10000)
24.5 ms ± 171 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

当然，这可以进一步加速。Jérôme给出了一些提示，例如您可以替换取模运算符（如果您很幸运，C编译器已经优化了这些操作），并尝试通过使用C/C++函数来替换对np.random.standard_normal的调用。 C ++ STL提供了std::normal_distribution, 可以轻松地在Cython中使用。但是，老实说，我认为这不值得付出努力，因为在循环内部并不经常调用np.random.standard_normal，不会对性能产生显着影响。

- joni

---

有趣！由于结果与最优计算时间相矛盾，我认为Cython能够内联函数并进行与Numba相同的优化，但实际上不是这样的：在生成的C代码中，Cython仍然使用__Pyx_PyObject_Call调用standard_normal函数，因此它无法被内联。我发现表达式185e-6*1000000/10000=0.185存在错误（当我用我的大脑计算时，我错过了一个0 ^^''）：它是0.0185！因此，你的结果是有意义的。Cython生成了一个很好的代码，挑战了Numba，但经过优化的Numba代码更快 ;)。 - Jérôme Richard

standard_normal调用可以使用numpy随机C API，从而消除Python开销。请参见https://github.com/ev-br/mc_lib/blob/master/mc_lib/rndm.pyx - ev-br

Answer 2

Cython 无法像 Numba 那样内联 Numpy 函数。此外，Cython 代码未经过类型注释，导致其操作速度变慢，因为它在 CPython 对象上运行（请参阅 @joni 回答中的具体细节）。

对于小型数组，Numpy 函数通常具有显着的开销，但看起来这里的主要问题不是开销。该函数在我的机器上为生成 10000 个项需要 185 微秒 （非常缓慢）。这意味着代码必须在 Python 和 Cython 中至少需要 185e-6*1000000/10000=0.0185 秒（实际上在 Python 中需要 0.265 秒）。由于 Cython 的时间显著更长，这意味着大部分开销出现在其他地方。@joni 的回答表明是由于缺少类型信息而导致的。优化后，大部分时间都花费在 np.random.standard_normal 函数上，正如预期的那样。

Numba 可以内联函数，因为它有自己的 Numpy 实现，可以进行 JIT 编译。缺点是复杂函数编译可能相当耗时，并且许多函数尚未得到支持，因为它们都需要重新实现，这是一个庞大的工作量。在 np.random.standard_normal 的情况下，Numba 实现的速度明显比仅使用 Numpy 快得多，但只有在循环中才是如此。这肯定是因为计算缓慢的表达式可以在循环中预先计算一次。对汇编代码的分析表明，这个预计算似乎与内部函数 numba::cpython::randomimpl::_gauss_pair_impl::_3clocals_3e::compute_gauss_pair_248 有关，如果我理解正确，这是一个“循环拆分”优化的结果（所以编译器能够将计算循环拆分为两个部分，并且只需做一次第一部分，因为结果总是相同）。很难进一步描述哪一部分被优化了，因为生成的汇编代码相当庞大且非常难以理解。要记住的主要点是，所有这些都是由于 Numpy 计算函数的内联和积极优化的结合而实现的。

顺便说一下，模运算很慢。您应该像避开瘟疫一样避免它们。使用带有条件的计数器肯定会更快（特别是因为处理器可以预测条件，因此可以非常快速地丢弃它）。

以下是经过优化的代码：

# Use parallel=True for a parallel implementation (see later)
@nb.jit(nopython=True, fastmath=True)
def simulation_nb(N=10000, N_save=10000):
    np.random.seed(0)
    x = np.zeros(N_save)
    r = np.empty(N_save)
    for j in range(N_save):
        r[j] = np.random.standard_normal()
    
    fs = int(N / N_save)
    xold = 0
    x[0] = 0

    counter1 = 1
    counter2 = 1
    counter3 = 0

    for i in range(1, N):

        if counter1 == N_save:
            # Use nb.prange here for a parallel implementation (see later)
            for j in range(N_save):
                r[j] = np.random.standard_normal()
            counter1 = 0

        xnew = xold + r[counter1]

        if counter2 == fs:
            x[counter3] = xnew
            counter3 += 1
            counter2 = 0

        xold = xnew

        counter1 += 1
        counter2 += 1

    return x

这段代码比之前的代码快得多（17毫秒而不是28毫秒）。np.random.standard_normal函数大约占据了此函数时间的95%。在循环中逐个调用np.random.standard_normal可以避免昂贵的数组分配（这只有在Numba中才能通过内联实现）。

与@joni的Cython版本相比，此版本更快（17毫秒VS 21毫秒）。

提高这个Numba代码的速度的一种方法是并行生成随机数，但这并不简单，因为主循环是固有的顺序循环，并且N_save的默认值对于多个线程帮助不大。将N_save增加一个数量级应该足以满足大多数系统需求。请注意，Numba按线程设置随机种子，因此结果可能会略有不同。还要注意，据我所知，Cython不可能做到这一点，因为Numpy的np.random.standard_normal实现不是线程安全的（与Numba相反）。经过优化的并行实现在我的6核机器上使用相同参数需要4.5毫秒，在10倍较大的N_save下仅需3.8毫秒（由于并行开销）。

- Jérôme Richard