BUG: random: dirichlet(alpha) can return nans in some cases.

Describe the issue:

When all the values in alpha are less than 0.1, and alpha ends in two or more zeros, the components of the variates returned by dirichlet(alpha) corresponding to those final zeros will be nan.

For example,

In [18]: rng.dirichlet([0.09, 0.085, 0, 0, 0])
Out[18]: array([0.91378938, 0.08621062,        nan,        nan,        nan])

When all the values in alpha are less than 0.1, dirichlet uses the algorithm that is based on the beta distribution. The problem occurs because dirichlet ends up calling the C function random_beta with both parameters a and b being 0, which results in random_beta returning nan. Currently, the public API for beta requires both a and b to be positive; this is checked before the beta method calls the C function random_beta. The dirichlet code calls random_beta directly, so that validation is bypassed.

It looks like random_beta handles one parameter being 0 in a manner consistent with the reasoning that allows dirichlet to have zeros in alpha. That's why nans are produced only when there are two or more zeros at the end of alpha, because that is the only case where dirichlet will call random_beta with both parameters being 0.

It shouldn't be too difficult to fix to dirichlet to handle the case where all the values in alpha are less than 0.1, and two or more values at the end of alpha are 0.

But there is a remaining question that was not brought up in #22547 or #23440: how should an input that is all zeros be handled? Some options (ordered by my preference):

Raise a ValueError (i.e. disallow alpha being all zeros).
Return a vector of zeros.
Return a random unit vector (i.e. a vector with len(alpha) - 1 zeros and single 1 at a random position in the vector). (On second thought, there probably isn't any reasonable justification for this.)

Runtime information:

In [4]: import sys, numpy; print(numpy.__version__); print(sys.version)
1.25.0rc1+530.g2e668061db
3.11.4 (main, Jul  3 2023, 14:49:40) [GCC 11.3.0]

In [5]: print(numpy.show_runtime())
[{'numpy_version': '1.25.0rc1+530.g2e668061db',
  'python': '3.11.4 (main, Jul  3 2023, 14:49:40) [GCC 11.3.0]',
  'uname': uname_result(system='Linux', node='pop-os', release='6.2.6-76060206-generic', version='#202303130630~1689015125~22.04~ab2190e SMP PREEMPT_DYNAMIC Mon J', machine='x86_64')},
 {'simd_extensions': {'baseline': [], 'found': [], 'not_found': []}},
 {'architecture': 'Zen',
  'filepath': '/usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so',
  'internal_api': 'openblas',
  'num_threads': 24,
  'prefix': 'libopenblas',
  'threading_layer': 'pthreads',
  'user_api': 'blas',
  'version': '0.3.20'}]
None

Context for the issue:

No response

BUG: random: Fix generation of nan by dirichlet.

I proposed a fix in #24220.

In that pull request, I disallowed an input that is all zeros. However, I'm now wondering if instead of an error, the return value in that case should just be a vector of zeros. Any opinions?

In that pull request, I disallowed an input that is all zeros. However, I'm now wondering if instead of an error, the return value in that case should just be a vector of zeros. Any opinions?

I think that returning a vector of zeros would be more consistent with #23440

The only thing I don't like about returning a vector of zeros is that it breaks the invariant that the sum of the vector is 1. But if the invariant is stated more carefully as "the sum of the variate components corresponding to nonzero alpha components will be 1", then it works.

Unless anyone brings up other points to consider, I'll update the PR later today or tomorrow to return all zeros instead of raising an exception.

BUG: random: Fix generation of nan by dirichlet.

I was looking at the PR going to suggest to simplify the branch a bit and to note that subnormals may need to be accounted for. Which led me to this:

rng.dirichlet([np.finfo(np.double).tiny/2, np.finfo(np.double).tiny/4], size=100)

which returns a mix of nan, nan, 0., 1. and 1., 0. (currently).

To me this looks like it isn't correct to return all zeros? Dirichlet doesn't have all dimensions result converge to return towards 0 in the limit of alphas -> 0. The ratio of alphas matters so the limit is very much undefined to me.

Even more obvious: A single dimensional dirichlet always returns 1 right now.

(The subnormal problem is that within random_beta you do 1/alpha, I think. For most subnormal numbers, 1./subnormal will overflow to inf giving you nan as well though (not sure about the exact parameters where this happens in practice).

One thing I could imagine being useful (but not sure!) would be to rescale alpha ahead of time: I.e. something like alphas /= alphas.sum(). NVM, sorry, I don't think that is how dirichlet works, unless things can be rephrased a lot more

... that subnormals may need to be accounted for.

Yeah, and to be clear, the handling of subnormals by the C function random_beta means the beta method can generate nan. This is an existing issue that is not the result of any recent work on dirichlet:

In [3]: np.__version__
Out[3]: '1.25.0rc1+551.gc41f9193c7'

In [4]: rng = np.random.default_rng()

In [5]: tiny = np.finfo(np.float64).tiny

In [6]: rng.beta(tiny/4, tiny/16, size=100)
Out[6]: 
array([ 1.,  0., nan,  0., nan,  1.,  1., nan,  1.,  1.,  1., nan,  0.,
        0., nan, nan,  1., nan,  1.,  1., nan, nan, nan,  1.,  1.,  0.,
        0.,  0., nan, nan,  1.,  1.,  1., nan, nan,  1., nan,  1., nan,
       nan,  1.,  1., nan,  0.,  1.,  1., nan,  0.,  1.,  1.,  0.,  0.,
        1., nan,  1.,  0.,  1.,  1.,  0.,  0.,  1., nan,  0.,  1.,  1.,
        0., nan,  1.,  1.,  1.,  0.,  1., nan,  1., nan, nan, nan,  1.,
        1., nan,  1., nan,  1.,  0., nan,  0., nan, nan,  1.,  0.,  1.,
       nan,  1.,  1.,  1.,  1.,  1.,  1.,  1., nan])

Right, I first thought you could just use if value < tiny rather thatn == 0. But if all values are very small that isn't right.

Another quick note: beta can generate nan even with normal values if they are very close to the smallest normal:

In [81]: np.__version__
Out[81]: '1.25.0rc1+563.g1c51f77fa5'

In [82]: tiny = np.finfo(np.float64).tiny

In [83]: rng = np.random.default_rng()

In [84]: x = rng.beta(2*tiny, 1.5*tiny, size=500000000)

In [85]: np.isnan(x).sum()
Out[85]: 454

In [97]: batch = 1

In [98]: while True:
    ...:     x = rng.beta(3*tiny, 3*tiny, size=100000000)
    ...:     n = np.isnan(x).sum()
    ...:     if n > 0:
    ...:         print(f'Found {n} nan in batch {batch}')
    ...:         break
    ...:     batch += 1
    ...: 
Found 1 nan in batch 35

I created a separate issue about beta: #24266

I'll submit a PR with a fix for that issue.

BUG: random: Fix generation of nan by dirichlet.

Navigation Menu

Repository navigation

BUG: random: `dirichlet(alpha)` can return nans in some cases. #24210

Describe the issue:

Runtime information:

Context for the issue:

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

BUG: random: dirichlet(alpha) can return nans in some cases. #24210

Description

Describe the issue:

Runtime information:

Context for the issue:

Activity

WarrenWeckesser commented on Jul 19, 2023

MatteoRaso commented on Jul 21, 2023

WarrenWeckesser commented on Jul 21, 2023

seberg commented on Jul 22, 2023

WarrenWeckesser commented on Jul 23, 2023

seberg commented on Jul 23, 2023

WarrenWeckesser commented on Jul 24, 2023

WarrenWeckesser commented on Jul 26, 2023

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions

BUG: random: `dirichlet(alpha)` can return nans in some cases. #24210