1.16.6 版本的 numpy.random.multinomial 比更高版本快 10 倍答案

【问题标题】：numpy.random.multinomial at version 1.16.6 is 10x faster than later version1.16.6 版本的 numpy.random.multinomial 比更高版本快 10 倍
【发布时间】：2022-01-22 11:19:45
【问题描述】：

这里是代码和结果：

python -c "import numpy as np; from timeit import timeit; print('numpy version {}: {:.1f} seconds'.format(np.__version__, timeit('np.random.multinomial(1, [0.1, 0.2, 0.3, 0.4])', number=1000000, globals=globals())))"

numpy version 1.16.6:  1.5 seconds # 10x faster
numpy version 1.18.1: 15.5 seconds
numpy version 1.19.0: 17.4 seconds
numpy version 1.21.4: 15.1 seconds

注意，固定随机种子，不同numpy版本的输出是一样的

python -c "import numpy as np; np.random.seed(0); print(np.__version__); print(np.random.multinomial(1, [0.1, 0.2, 0.3, 0.4], size=10000))" /tmp/tt

关于为什么 1.16.6 之后的 numpy 版本慢 10 倍的任何建议？

我们已经将pandas升级到最新版本1.3.4，1.16.6之后需要numpy版本

【问题讨论】：

更新日志是否显示了什么？
我填写了this issue 以获得 Numpy 团队的一些反馈，并询问他们是否可以改善这种情况。

标签： numpy performance multinomial

【解决方案1】：

我检查了一些引擎盖下的发电机，发现时间没有太大变化。

我猜差异可能是由于一些开销，因为您只对单个值进行采样。这似乎是一个很好的假设。当我将生成的随机样本的大小增加到 1000 时，1.16.6 和 1.19.2（我当前的 Numpy 版本）之间的差异减少到 ~20%。

python -c "import numpy as np; from timeit import timeit; print('numpy version {}: {:.1f} seconds'.format(np.__version__, timeit('np.random.
multinomial(1, [0.1, 0.2, 0.3, 0.4], size=1000)', number=10000, globals=globals())))"

numpy version 1.16.6: 1.1 seconds
numpy version 1.19.2: 1.3 seconds

请注意，两个版本都有这个开销，只是新版本的开销要大得多。在这两个版本中，对 1000 个值进行一次采样比对 1 个值进行 1000 次采样要快得多。

他们在 1.16.6 和 1.17.0 之间更改了很多代码，例如 this commit，很难分析。抱歉不能更好地帮助您 - 我建议在 Numpy 的 github 上提出问题。

【讨论】：

版本 1.17.0 是引入新 API 的地方，这可能与它有关。你知道np.random.default_rng().multinomial是否有同样的开销？
@MadPhysicist 是的，这甚至要慢得多。将 np.random.default_rng() 移动到初始化有帮助，但这仍然具有相同/稍差的性能

【解决方案2】：

TL;DR：这是local performance regression 由numpy.random.multinomial 函数中的附加检查 开销引起的。 非常小的数组由于所需检查的相对执行时间而受到很大影响。

引擎盖下

对Numpy code 的 Git 提交进行二分搜索显示，性能回归首次出现在 2019 年 4 月中旬。它可以在提交 dd77ce3cb 中重现，但不能在 7e8e19f9a 中重现。中间的提交存在一些构建问题，但通过一些快速修复，我们可以证明提交 0f3dd0650 是第一个导致问题的。提交说它：

扩展多项式以允许广播
修复 NumPy 中遗漏的 zipf 更改
启用 0 作为超几何的有效输入

对 this commit 的深入分析表明，它修改了 Cython 文件 mtrand.pyx 中定义的 multinomial 函数，以执行以下两项额外检查：

def multinomial(self, np.npy_intp n, object pvals, size=None):
    cdef np.npy_intp d, i, sz, offset
    cdef np.ndarray parr, mnarr
    cdef double *pix
    cdef int64_t *mnix
    cdef int64_t ni

    d = len(pvals)
    parr = <np.ndarray>np.PyArray_FROM_OTF(pvals, np.NPY_DOUBLE, np.NPY_ALIGNED)
    pix = <double*>np.PyArray_DATA(parr)
    check_array_constraint(parr, 'pvals', CONS_BOUNDED_0_1)   # <==========[HERE]
    if kahan_sum(pix, d-1) > (1.0 + 1e-12):
        raise ValueError("sum(pvals[:-1]) > 1.0")

    if size is None:
        shape = (d,)
    else:
        try:
            shape = (operator.index(size), d)
        except:
            shape = tuple(size) + (d,)

    multin = np.zeros(shape, dtype=np.int64)
    mnarr = <np.ndarray>multin
    mnix = <int64_t*>np.PyArray_DATA(mnarr)
    sz = np.PyArray_SIZE(mnarr)
    ni = n
    check_constraint(ni, 'n', CONS_NON_NEGATIVE)              # <==========[HERE]
    offset = 0
    with self.lock, nogil:
        for i in range(sz // d):
            random_multinomial(self._brng, ni, &mnix[offset], pix, d, self._binomial)
            offset += d

    return multin

这两项检查是代码稳健所必需的。但是，考虑到它们的用途，它们目前相当昂贵。

确实，在我的机器上，第一次检查负责 ~75% 的开销，第二次检查负责 ~20%。检查需要几微秒，但由于您的输入非常小，与计算时间相比，开销很大。

解决此问题的一种解决方法是为此编写一个特定的 Numba 函数，因为您的输入数组非常小。在我的机器上，np.random.multinomial 在一个简单的 Numba 函数中会产生良好的性能。

【讨论】：