ENH: reduce the overhead of the checks in the multinomial function

[zephyr111]

Proposed new feature or change:

Hello,

We have discovered a performance regression of the numpy.random.multinomial function for small arrays apparently caused by the commit 0f3dd06 which add new checks. The checks are expensive regarding their purpose (especially for small arrays). This can be reproduced with Python 3.9.9 on a x86-64 machine using both Windows and Linux. Is there a way to reduce the overhead of such checks?

For more information, please read this Stack-Overflow post.

Thank you.

[mattip]

Saving a click: this microbenchmark

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

gives these results:

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

The commit was part of the enormous PR #13163 to integrate randomgen into numpy.

[mattip]

A very good and extensive analysis of the overhead is given by in the linked StackOverflow post linked to above, the analysis points to this line in mtrand.pyx causing 75% of the overhead

numpy/numpy/random/mtrand.pyx

Line 4249 in 7acb0fd

check_array_constraint(parr, 'pvals', CONS_BOUNDED_0_1)

and this line as responsible for 20% more

numpy/numpy/random/mtrand.pyx

Line 4278 in 7acb0fd

check_constraint(ni, 'n', CONS_NON_NEGATIVE)

The analysis also points out that calling multinomial 1,000,000 times for a single value in a tight loop is not a good design pattern.

[mattip]

Changing the code to use the new Generator

rng = np.random.default_rng()
res = timeit('rng.multinomial(1, [0.1, 0.2, 0.3, 0.4])', 
             number=1000000, globals=globals())

produces the same result since the similar code exists in _generator.pyx. The function check_array_constraint contains calls to np.any, np.all, and more which are not nicely optimized by Cython. They call back to the Python C-API to look up the np.any and np.all functions in numpy and call them with PyOjbect*.

We could refactor check_array_constraint faster by avoiding the C-API. I think using a 1d memory view and explicitly looping might work, but would make the code much more unwieldy.

[bashtage]

The check_constraint can probably be improved since this is for scalars.

numpy/numpy/random/_common.pyx

Lines 400 to 402 in 01079e2

	if cons == CONS_NON_NEGATIVE:
	if not np.isnan(val) and np.signbit(val):
	raise ValueError(name + " < 0")

I suspect the np. calls are not optimized by Cython.

[bashtage]

Annotated Cython shows this is exactly the case:

There are a couple of places where Python API are used (the if and the first elif) that could be pretty easily avoided.

[bashtage]

Fixing the performance of the array checker in a tidy manner is a lot harder since the type of the array being checked is not always the same, and so would need paths for different possible types.

[wguocbp]

Just installed the latest numpy and it works like a charm, thank you all for the fix and the unexpected fast turnaround, really impressed!

[mattip]

Closing since we improved what we can. @zephyr111 please open a new issue or reopen this one if you have an idea for more preformance improvements.