Adding a numerically stable inverse CDF of standard Gaussian

1

Hi,

The current icdf in torch.distributions.Normal seems to be numerically unstable compared to scipy.stats.norm.

I implemented a torch version by translating the scipy source code to python and it works the same as the scipy version. I wonder if we could integrate the new code to the current torch.distributions.Normal. Thank you!

2

Could you explain the numerical instability a bit? I.e. where is it observed and how scipy deals with it?
I’m sure a PR would be welcome, if the instability is reproducible or breaks certain use cases.

3

Numerical instability issues occur when the input to the inverse CDF is very close to 0 or 1.
scipy handles these cases with approximations and here is the screenshot of the comments in scipy’s source code:

image
image954×524 34.8 KB

One example is when y=1-1e-8:

  • torch.distributions.Normal(0, 1).icdf(torch.tensor(y)) returns tensor(inf) (because torch.erfinv returns tensor(inf));
  • scipy.stats.norm.ppf(y) returns 5.612001243305505.
1 Like
4

Thanks for the update!
Would you like to create a feature request on GitHub and post the explanation there, so that we could discuss it there with the module owners and track it?
Btw. thanks for digging into this. :slight_smile:

1 Like
5

that’s kinda related to float32 precision, as 1e-8 < 2**-23 = 1.19e-7

and
torch.erfinv(torch.tensor([1-2**-23])) → 3.7439

issue is this:

np.float32(1)-np.float32(1e-8)
1.0
np.frexp(np.float32(1)-np.float32(1e-8))
(0.5, 1)

2 Likes
6

I see, this is great, thank you @googlebot !!
Hi @ptrblck, I think there’s no need for a feature request since using DoubleTensor fixes the issue.

Thank you both for your help and suggestions~!