`poisson()` defaults to the approximate POISNORM method; patched `numpy.random.poisson` does not sample Poisson(lambda)

[vchamarthi]

Affects: mkl_random.MKLRandomState.poisson default, mkl_random.interfaces.numpy_random.poisson, and numpy.random.poisson after mkl_random.patch_numpy_random()

Issue

MKLRandomState.poisson defaults to method="POISNORM" (mklrand.pyx:5789), and the NumPy interface hardcodes it (interfaces/_numpy_random.py:552). oneMKL documents POISNORM for lambda >= 1 as an approximation: "method based on Poisson inverse CDF approximation by Gaussian inverse CDF".

The docstring gives the exact Poisson pmf and doesn't say the method is approximate. The NumPy interface refers users to numpy.random.poisson, which is exact.

The error is too small to detect at 10^6-10^7 draws but clear at 10^8:

import numpy as np
import mkl_random
from scipy import stats

lam, n = 10.0, 10**8
x = mkl_random.interfaces.numpy_random.RandomState(1).poisson(lam, size=n)
counts = np.bincount(x, minlength=4)[:4]
print(counts / (stats.poisson(lam).pmf(np.arange(4)) * n) - 1)
# k = 1: -0.031, k = 2: -0.016, k = 3: -0.0048   (1 SE: 0.0047, 0.0021, 0.0011)

How found

Chi-square goodness-of-fit against the exact Poisson pmf, computed in 60-digit arithmetic, with seeds 1, 2 and 3. The control is an exact inverse-CDF sampler under the same test. Under the null the test is calibrated at every N used (z ~ N(0, 1), uniform p-values). An injected error of -3% at P(X=1) and -1.6% at P(X=2) gives z = 0.3 at 10^7 and z = 11.2 at 10^8, as predicted by the test's noncentrality.

z for POISNORM (the default), per seed:

lambda	N = 10^6	N = 10^7	N = 10^8	control, N = 10^8
1.5	-0.3, -1.2, -0.5	1.0, -0.3, -0.7	3.5, 2.5, -0.1	-0.8, 0.2, -1.2
3	-0.4, -1.4, 0.8	-0.5, 1.1, 1.7	22.7, 24.8, 23.6	-1.1, -0.4, -1.0
10	1.8, -1.6, -0.9	4.3, 3.7, 0.7	26.1, 27.1, 24.6	-1.0, -0.8, -0.5
30	-0.3, 0.3, -0.8	-0.3, 0.6, -1.9	-0.2, 0.2, 0.8	-0.3, 0.4, 0.1
100	1.6, -0.4, -1.9	0.9, 1.0, -0.3	6.7, 5.2, 6.3	-0.9, -0.7, -0.8
1000	0.4, -0.9, -1.3	0.2, 0.6, -0.9	2.8, 2.0, 3.9	-0.4, -2.0, 0.3

Relative error of P(X = k) at lambda = 10, pooled over 3 * 10^8 draws:

k	0	1	2	3	4	5	6	7
rel. error	+0.004	-0.035	-0.015	-0.005	+0.000	+0.002	+0.001	+0.001
SE	+0.5	-12.8	-12.3	-7.6	+0.2	+5.1	+6.1	+3.4

The NumPy interface's output is bit-identical to method="POISNORM", and POISNORM's output is identical for oneMKL 2023.2 through 2026.1.

Environment: mkl_random 1.5.0 (py313h581662c_0, Intel channel), oneMKL 2026.1.0, NumPy 2.4.3, Linux x86-64, Xeon 6980P.

Performance

1 core, 10^7 draws, ns per draw:

lambda	POISNORM (current default)	PTPE	NumPy Generator	NumPy RandomState
10	0.88	13.24	48.63	54.18
100	0.88	8.02	27.36	33.79

The 31-62x speedup of the patched numpy.random.poisson over NumPy comes from the approximation. PTPE, which should be exact, is 3.4-3.7x faster than NumPy's Generator.

Possible fix

• Default the NumPy interface to an exact method, and keep POISNORM available through method=.
• PTPE cannot be that default yet: it is wrong for 27 <= lambda <= 100 in all oneMKL versions tested (reported separately against oneMKL). Until that's fixed, the NumPy interface needs another exact sampler for that range.
• State in the poisson docstring that POISNORM is approximate for lambda >= 1.
• Add a goodness-of-fit test for the default method at lambda = 3 and 10 with at least 10^8 draws.