BUG: fix negative samples generated by Wald distribution #29609

mmaaz-git · 2025-08-21T06:59:25Z

In numpy.random, when generating Wald samples, if the mean and scale have a large discrepancy, the current implementation sometimes generates negative samples. An explicit example of this has been added to test_generator_mt19937.py, mean=1e8, scale=2.25 (which is specific to the seed used in that test file).

The key line implicated in distributions.c:

X = mean + mu_2l * (Y - sqrt(4 * scale * Y + Y * Y));

which has numerical issues when Y >> scale.

I have replaced this with the equivalent rationalized form:

d = Y + sqrt(Y) * sqrt(Y + 4 * scale);
X = mean - (2 * mean * Y) / d;

This can be seen my noting the following algebraic identity:

Y - sqrt(4 * scale * Y + Y * Y))
= - 4 * scale * Y / (Y + sqrt(Y^2 + 4 * scale * Y))

As mu_2l = mean / (2 * scale), we then have:

X = mean - 2 * mean * Y / (Y + sqrt(Y^2 + 4 * scale * Y))

Introduce a new variable for the denominator and factoring gives the final expressions:

d = Y + sqrt(Y) * sqrt(Y + 4 * scale)
X = mean - (2 * mean * Y) / d

And now the test passes, i.e., all samples are non-negative.

I have not made this change to the legacy implementation (nor the legacy test) as I figure they are deprecated, although it also suffers from the same issue.

Disclaimer: I found this failing property during a research project I am working on using LLM agents to write property-based tests (i.e., in Hypothesis). This property, that all Wald samples should be non-negative, was proposed and written "autonomously" by the agent. I spent several hours analyzing the failure and becoming acquainted with the NumPy codebase myself in order to verify it was a real issue, running it with the seed in the test suite, etc. I also pinpointed the cause of the problem and wrote this fix and the PR myself. I take full responsibility for it.

The Hypothesis test that caught the issue was simply:

import numpy.random
from hypothesis import given, strategies as st, settings
from numpy.random import MT19937, Generator

@given(
    st.floats(min_value=1e8, max_value=1e15),
    st.floats(min_value=0.1, max_value=10.0)
)
@settings(max_examples=50)
def test_wald_negative_values(mean, scale):
    """Wald distribution should never produce negative values"""
    random = Generator(MT19937(1234567890))
    samples = random.wald(mean, scale, size=1000)
    assert all(s >= 0 for s in samples), f"Found negative values with mean={mean}, scale={scale}"

Note that even after this PR, this property fails, albeit the minimal example is now mean=280432326968266.0, scale=0.25 (i.e., even more extreme values). But, at such values, it seems like we hit the limits of floating point computations.

WarrenWeckesser · 2025-08-21T14:59:24Z

Thanks @mmaaz-git, I confirmed that the problem exists on the main development branch.

At first glance, your reformulation looks mathematically correct and avoids the problem. I probably won't get to a closer look until this weekend (anyone else can jump in, of course).

One change I suggest before I click the "Approve workflows to run": you say the new test fails for you without the modified formulas, but in my build of numpy, your new test passes using the main development branch. I was able to get the test to fail by increasing mean from 1e8 to 1e9. Then I get 3 negative values. Could you update the test with a change like that, to increase the chances that it fails on more platforms?

mmaaz-git · 2025-08-21T16:56:19Z

I've changed the test case to mean=1e9 and can confirm I also get 3 negative values on my machine.

WarrenWeckesser

For anyone interested, the method that is used to generate variates is described on the wikipedia page
https://en.wikipedia.org/wiki/Inverse_Gaussian_distribution#Sampling_from_an_inverse-Gaussian_distribution

Our current implementation uses the following formulation to compute x. I changed the notation slightly to be closer to the wikipedia article and avoid reusing variable names, so mean is mu and scale is lam (i.e. lambda).

n = Normal(0, 1)  # Standard normal variate
mu_2l = mu / (2 * lam)
y = mu * n**2
x = mu + mu_2l * (y - sqrt(4 * lam * y + y * y))

That is just the first step, where x is computed.

The problem is the loss of precision in y - np.sqrt(4 * lam * y + y * y). When y is sufficiently large, the two terms in that subtraction are relatively close, which results in a catastrophic loss of precision.

In this pull request, the formulation is changed to

n = Normal(0, 1)  # Standard normal variate
y = mu * n**2
d = y + sqrt(y) * sqrt(y + 4 * lam)
x = mu - (2 * mu * y) / d

In effect, @mmaaz-git has "rationalized the numerator" and replaced the problematic subtraction with an addition in the denominator.

The following plot of the relative error in x as a function of the normal variate n, with mean = mu = 1e8 and scale = lam = 2.25, shows the effectiveness of this reformulation. The relative error was determined by computing a high precision result using mpmath.

@mmaaz-git, one problem with your version is that the normal variate (and therefore y) could be 0, and your formula generates nan in that case. You'll need to check for the normal variate being 0; just return mu in that case.

We can also make a few small changes to the expression to cut down on the number of floating point operations, and end up with something like the following (where I also show the check for the normal variate being 0):

n = Normal(0, 1)  # Standard normal variate
if n == 0:
    return mu
y = mu * n**2
d = 1 + sqrt(1 + 4 * lam / y)
x = mu * (1 - 2 / d)

(No need to change the actual variable names; keep mean, scale etc.)

mmaaz-git · 2025-08-24T05:15:38Z

Thanks for the detailed review! Really satisfying to see the relative error plot, and great points about the normal rv being 0, and how to simplify some of the expressions. I have made the requested changes.

djeada · 2025-08-24T21:00:47Z

numpy/random/src/distributions/distributions.c


-  mu_2l = mean / (2 * scale);
  Y = random_standard_normal(bitgen_state);
+  if (Y == 0) {


does the Y need to be exactly 0 or close to zero within some epsilon like fabs(Y) < 1e-12 ?

That's a good question. I tried it with forcing Y=1e-10 or 1e-15, etc, and it doesn't seem to cause any issues.

For what it's worth, I realized that even without the check for if (Y==0), it's still fine when Y=0, because you get d=inf (per the IEEE standard), and then X = mean * (1 - d/inf) = mean.

I'm not sure what is the usual numpy practice here -- can we rely on inf as an intermediate step?

Maybe mostly to wake up this PR again. I think relying on the inf gracefully overflowing seems OK if @WarrenWeckesser agrees (I don't think the floating point error that gets set here matters).

But maybe it doesn't matter and this is just ready, @WarrenWeckesser are you happy with this?

I'll get back to this tomorrow, if not sooner.

For what it's worth, I realized that even without the check for if (Y==0), it's still fine when Y=0, because you get d=inf (per the IEEE standard), and then X = mean * (1 - d/inf) = mean.

Yes, you're right--the inf propagates correctly. (We know scale is not 0, so we won't ever compute 0/0 resulting in nan.) And this is C, so we don't have to worry about exceptions, and we are free to ignore the floating point status flags.

So my initial suggestion to check for Y == 0 was unnecessary--sorry for the extra work! I'm OK now with undoing that change, and letting Y = 0 propagate unchecked.

If we did want to avoid dividing by 0, then instead of using a small threshold as @djeada suggested, we could check for mean * Y * Y being 0, since that is what is in the denominator of the expression for d. If mean and the initial Y are sufficiently small, mean * Y * Y could underflow to 0. So instead of a small threshold on Y, the product could be computed and checked for 0. But, unless one of the resident C gurus says we really should avoid dividing by 0, I don't think it is necessary.

Thanks for looking into it! I've removed the check now.

The previously suggested change is not necessary.

In numpy.random, when generating Wald samples, if the mean and scale have a large discrepancy, the current implementation suffers from catastrophic cancellation. An explicit example of this has been added to test_generator_mt19937.py. The key line implicated in distributions.c: `X = mean + mu_2l * (Y - sqrt(4 * scale * Y + Y * Y));` which has numerical issues when Y >> scale. I have replaced this with the equivalent rationalized form: ``` d = Y + sqrt(Y) * sqrt(Y + 4 * scale); X = mean - (2 * mean * Y) / d; ``` And now the test passes.

mmaaz-git · 2025-10-10T19:01:54Z

Had to rebase and re-sync submodules

mmaaz-git · 2025-10-10T20:24:12Z

Looks like everything is passing! What does it mean to be a backport candidate? Do I need to make any additional changes?

WarrenWeckesser · 2025-10-10T21:57:45Z

Thanks @mmaaz-git, merged.

This change is now in the main development branch, so it will part of NumPy 2.4.0.

"Backport candidate" means the change might also be included in the next point release of the 2.3 series, which would be 2.3.4 (if there is one).

In numpy.random, when generating Wald samples, if the mean and scale have a large discrepancy, the current implementation suffers from catastrophic cancellation. An explicit example of this has been added to test_generator_mt19937.py. The key line implicated in distributions.c: `X = mean + mu_2l * (Y - sqrt(4 * scale * Y + Y * Y));` which has numerical issues when Y >> scale. I have replaced this with the equivalent rationalized form: ``` d = Y + sqrt(Y) * sqrt(Y + 4 * scale); X = mean - (2 * mean * Y) / d; ``` And now the test passes.

BUG: fix negative samples generated by Wald distribution (#29609)

seberg · 2025-10-11T05:41:13Z

Should be good, I wonder a bit if we shoukd be very conservative about possibly changing random streams in bug fix releases?
This isn't the old generator though of course.

In numpy.random, when generating Wald samples, if the mean and scale have a large discrepancy, the current implementation suffers from catastrophic cancellation. An explicit example of this has been added to test_generator_mt19937.py. The key line implicated in distributions.c: `X = mean + mu_2l * (Y - sqrt(4 * scale * Y + Y * Y));` which has numerical issues when Y >> scale. I have replaced this with the equivalent rationalized form: ``` d = Y + sqrt(Y) * sqrt(Y + 4 * scale); X = mean - (2 * mean * Y) / d; ``` And now the test passes.

github-actions bot added the 00 - Bug label Aug 21, 2025

WarrenWeckesser added the component: numpy.random label Aug 21, 2025

melissawm added this to NumPy first-time contributor PRs Aug 21, 2025

melissawm moved this to Awaiting a code review in NumPy first-time contributor PRs Aug 22, 2025

WarrenWeckesser previously requested changes Aug 24, 2025

View reviewed changes

djeada reviewed Aug 24, 2025

View reviewed changes

mmaaz-git requested a review from WarrenWeckesser August 28, 2025 03:47

charris added the 09 - Backport-Candidate PRs tagged should be backported label Sep 1, 2025

mmaaz-git added 5 commits October 10, 2025 14:55

Removed mu_2l variable

139fa1f

changed test case

70ac43a

check for 0 normal rv and simplify expressions

ad4344c

remove check for Y==0

73545ea

mmaaz-git force-pushed the stable-wald branch from 28c4589 to 73545ea Compare October 10, 2025 19:00

WarrenWeckesser merged commit 0ce0ef0 into numpy:main Oct 10, 2025
77 checks passed

github-project-automation bot moved this from Awaiting a code review to Completed in NumPy first-time contributor PRs Oct 10, 2025

charris removed the 09 - Backport-Candidate PRs tagged should be backported label Oct 10, 2025

charris added a commit that referenced this pull request Oct 10, 2025

Merge pull request #29926 from charris/backport-29609

e2f0383

BUG: fix negative samples generated by Wald distribution (#29609)

Uh oh!

BUG: fix negative samples generated by Wald distribution #29609

BUG: fix negative samples generated by Wald distribution #29609

Uh oh!

Conversation

mmaaz-git commented Aug 21, 2025

Uh oh!

WarrenWeckesser commented Aug 21, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

mmaaz-git commented Aug 21, 2025

Uh oh!

WarrenWeckesser left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mmaaz-git commented Aug 24, 2025

Uh oh!

djeada Aug 24, 2025

Choose a reason for hiding this comment

Uh oh!

mmaaz-git Aug 25, 2025

Choose a reason for hiding this comment

Uh oh!

seberg Oct 9, 2025

Choose a reason for hiding this comment

Uh oh!

WarrenWeckesser Oct 9, 2025

Choose a reason for hiding this comment

Uh oh!

WarrenWeckesser Oct 10, 2025

Choose a reason for hiding this comment

Uh oh!

mmaaz-git Oct 10, 2025

Choose a reason for hiding this comment

Uh oh!

mmaaz-git commented Oct 10, 2025

Uh oh!

mmaaz-git commented Oct 10, 2025

Uh oh!

Uh oh!

WarrenWeckesser commented Oct 10, 2025

Uh oh!

seberg commented Oct 11, 2025

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

5 participants

WarrenWeckesser commented Aug 21, 2025 •

edited

Loading

WarrenWeckesser left a comment •

edited

Loading