numpy.random.multinomial在版本1.16.6比后续版本快10倍。

Question

这里是代码和结果：

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

请问为什么 numpy 版本在 1.16.6 之后的版本会慢十倍？

我们升级了 pandas 到最新版本 1.3.4，因此需要使用 1.16.6 之后的 numpy 版本。

- Winston Guo

---

变更日志有什么新内容吗？ - Mad Physicist

我提交了这个问题，以获取Numpy团队的反馈，并询问他们是否可以改善这种情况。 - Jérôme Richard

Answer 1

TL;DR: 这是由于 numpy.random.multinomial 函数中的 额外检查 开销导致的 本地性能回归。由于所需检查的相对执行时间，非常小的数组 受到强烈影响。

---

底层实现

在Numpy代码的Git提交记录中进行二分查找，发现性能回归最初出现在2019年4月中旬。可以在提交dd77ce3cb中重现，但不能在7e8e19f9a中。其中有一些构建问题，但是通过一些快速修复，我们可以证明提交0f3dd0650是第一个导致问题的提交。该提交说明了它:

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

深入分析此提交发现它修改了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函数，因为您的输入数组非常小。在我的机器上，在一个简单的Numba函数中使用np.random.multinomial可以获得良好的性能。

- Jérôme Richard

Answer 2

我检查了一些生成器，发现时间上没有太大的变化。

我猜想这种差异可能是由于一些开销造成的，因为您只对单个值进行采样。这似乎是一个好的假设。当我增加生成的随机样本的大小到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个值比采样1000次1个值要快得多。

他们在1.16.6和1.17.0之间改变了很多代码，例如此提交，很难分析。抱歉帮不上忙 - 我建议在Numpy的github上提出问题。

- dankal444

---

1.17.0 版本是引入新 API 的版本，这可能与此有关。您知道 np.random.default_rng().multinomial 是否具有相同的开销吗？ - Mad Physicist

@MadPhysicist 是的，这个速度甚至更慢。将 np.random.default_rng() 移动到初始化可以提高性能，但仍然具有相同/更差的性能。 - dankal444