Overflow and admissible range for float 32 input parameters

[remy-abergel]

Hi,

Processing some unusual raw data related to Electron Paramagnetic Resonance imaging have led me to normalization issues that can be roughly summarized by the following code.

import numpy as np
import finufft

###########################
# generate dataset (fp32) #
###########################
M = 10000
N1, N2, N3 = 128, 128, 128
x = np.linspace(-np.pi, -np.pi, M, dtype=np.float32)
y = np.linspace(-np.pi, -np.pi, M, dtype=np.float32)
z = np.linspace(-np.pi, -np.pi, M, dtype=np.float32)

cof = float(1e13)
c_re = cof * np.linspace(0, 1, M, dtype=np.float32)
c_im = cof * np.linspace(0, 1, M, dtype=np.float32)
c = c_re + 1J * c_im

######################################
# method n. 1 : apply nufft3d1 as is #
######################################
f = finufft.nufft3d1(x, y, z, c, n_modes=(N1, N2, N3))
print("f.real: min = %.3e, max = %.3e\nf.imag: min = %.3e, max = %.3e" % (f.real.min(), f.real.max(), f.imag.min(), f.imag.max()))
# f.real: min = -inf, max = inf
# f.imag: min = -inf, max = inf

################################################################
# method n. 2 : normalize c input and rescale output afterward #
################################################################
c_nrm = c/cof
f = cof * finufft.nufft3d1(x, y, z, c_nrm, n_modes=(N1, N2, N3))
print("f.real: min = %.3e, max = %.3e\nf.imag: min = %.3e, max = %.3e" % (f.real.min(), f.real.max(), f.imag.min(), f.imag.max()))
# f.real: min = -5.000e+16, max = 5.000e+16
# f.imag: min = -5.000e+16, max = 5.000e+16

For sure, method n. 2 is easy to implement and may probably be preferred to method n. 1. However, provided that all the inputs x, y, z, c and the theoretical value of the real and imaginary parts of the output f fall into the admissible range for np.float32 arrays ([1e-38, 1e-38]) with reasonable security margins, should we expect method n. 1 to run without overflow issues? Otherwise, is there some guidelines about the admissible range for nufft parameters?

Many thanks

[lu1and10]

I guess it depends on your point distribution, M, N1, N2, N3, etc.
In your 3D case the dominate part seems to be $e^{3\beta} * cof$ ($\beta$ is the ES kernel parameter in the paper), I guess it's around 16 in your setup, $e^{48} * 1e13$ is already a quite large number, after the output from fft probably is inf.

[remy-abergel]

Thank you again for you reply. I am afraid I do not understand your point here.

My point is that the theoretical value of f = nufft3d1(x, y, z, c, n_modes=(N1, N2, N3)) does not overflow in float32 precision. (EDIT: Besides the user does not provide here inputs with 1e38 amplitude (but only 1e13)).

If nufft internally relies on fftw (or whatever) evaluation over arrays with maximal amplitude of 5e37, overflow is likely to happen. This kind of issue can be easily avoided by applying an appropriate normalization prior to nufft, I simply wonder if this normalization operation should be anticipated and handled by the user or handled internally during the nufft call.

My feeling is that the user cannot be aware that nufft calls over arrays with 1e13 maximal amplitude will rely on intermediate arrays with 5e37 maximal amplitude (and cause overflow).

[lu1and10]

Thank you again for you reply. I am afraid I do not understand your point here.

My point is that the theoretical value of f = nufft3d1(x, y, z, c, n_modes=(N1, N2, N3)) does not overflow in float32 precision. (EDIT: Besides the user does not provide here inputs with 1e38 amplitude (but only 1e13)).

If nufft internally relies on fftw (or whatever) evaluation over arrays with maximal amplitude of 5e37, overflow is likely to happen. This kind of issue can be easily avoided by applying an appropriate normalization prior to nufft, I simply wonder if this normalization operation should be anticipated and handled by the user or handled internally during the nufft call.

My feeling is that the user cannot be aware that nufft calls over arrays with 1e13 maximal amplitude will rely on intermediate arrays with 5e37 maximal amplitude (and cause overflow).

I'm not sure what is the easily appropriate normalization could be used handle to both overflow and underflow here, maybe you could propose a normalization scheme here and we could incorporate to the code. I guess there is a limitation on the range of an input, one numeric program can handle.
The reason you expect the input of range 1e13 should return output of range of similar order is because your input non-uniform points is actually uniform, if a user supply a different non-uniform distribution of points, the intermediate values could be much larger. I would prefer a good normalization scheme can handle underflow, overflow and different distributions of non-uniform points. Do you think normalize by the mean input value should be enough?

[remy-abergel]

I'm not sure what is the easily appropriate normalization could be used handle to both overflow and underflow here

You are right. Finding a good normalization seems easy in my example but things are probably more tricky in the general case. I will take a look at what can be an appropriate normalization, at least for my EPR datasets. I will let you know if I find something interesting.

The reason you expect the input of range 1e13 should return output of range of similar order is because your input non-uniform points is actually uniform, if a user supply a different non-uniform distribution of points, the intermediate values could be much larger.

In situations where the output is out of the range [1e-38, 1e38], I think that their is no reason to argue. The user must consider extended floating point precision (e.g. float64) or compute something else (e.g., in a totally different context, one may be satisfied with computing gammaln instead of gamma in order to avoid underflow/overflow issues in further computations).

The actual examples I am considering with EPR do not fall in this category. These datasets involve non-uniformly distributed frequency nodes (since EPR usually relies on regularly sampled 3D Radon transforms) and lead to outputs with significantly smaller amplitude than 1e38 (typically 1e10 if no normalization of the raw data is performed). I have the feeling that computation in float32 precision remains possible but I agree that finding the appropriate normalization is not easy.

If finding an appropriate normalization in the general case reveals to be too complicated, providing to the user some guidelines about the admissible range of parameters (depending on the available precision) would be nice anyway.

Anyway, the more I use finufft, the more I appreciate its speed and reduced memory consumption. Thank you.

[lu1and10]

If finding an appropriate normalization in the general case reveals to be too complicated, providing to the user some guidelines about the admissible range of parameters (depending on the available precision) would be nice anyway.

I'm not sure we should normalize respect to 1. user data or 2. just change the ES kernel from $e^{\beta\sqrt{x^2-1}}$ instead of $e^{\beta(\sqrt{x^2-1}-1)}$, the problem you have for overflow is coming from the $e^{\beta}$ term. I think we will discuss in our group meeting to see what we can do for this.

In summary, ideally you expect if input is in [1e-38, 1e+38] and output is also in [1e-38, 1e+38], nufft should not return inf/nan as a result, right?

[remy-abergel]

In summary, ideally you expect if input is in [1e-38, 1e+38] and output is also in [1e-38, 1e+38], nufft should not return inf/nan as a result, right?

Ideally yes, but honestly, I would feel scary to try to compute nufft (or even fft) of an input array with 1e38 maximal amplitude in float32 precision. However, for input and output falling in this range with a more reasonable security margin (to be defined), if something can be done, I think it would be nice to do it.

[ahbarnett]

Thanks for finding what is essentially a bug, and for the kind words. This is a failure mode I had never anticipated: because you're in 3D, the kernel max size gets cubed. for eps=1e-6 (about the min useful in single-prec), w=7, so kernel size is about 1e7, so its cube is 1e21, and mult by your strength 1e13 and M=1e4 points all falling on top of each other (in your data), and you get 1e38. Bingo! We will rescale the kernel internally, and this will fix it. Certainly this never bothered anyone in double-prec, but single I see is a problem. Thanks, Alex

…

On Mon, Jun 3, 2024 at 3:03 PM Rémy Abergel ***@***.***> wrote: In summary, ideally you expect if input is in [1e-38, 1e+38] and output is also in [1e-38, 1e+38], nufft should not return inf/nan as a result, right? Ideally yes, but honestly, I would feel scary to try to compute nufft (or even fft) of an input array with 1e38 maximal amplitude in float32 precision. However, for input *and* output falling in this range with a more reasonable security margin (to be defined), if something can be done, I think it would be nice to do it. — Reply to this email directly, view it on GitHub <#448 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACNZRSTDS7VFHI2BIHN5573ZFS4ZRAVCNFSM6AAAAABISVIWRSVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4TMNJTGA4DA> . You are receiving this because you were mentioned.Message ID: ***@***.*** com>

-- *-------------------------------------------------------------------~^`^~._.~' |\ Alex Barnett Center for Computational Mathematics, Flatiron Institute | \ http://users.flatironinstitute.org/~ahb 646-876-5942

[ahbarnett]

MCE in Matlab:

>> finufft3d1(single(1.0),single(1.0),single(1.0),single(1e17),+1,1e-6,1,1,1)
ans =
  single
           Inf +            0i
>> finufft3d1(single(1.0),single(1.0),single(1.0),single(1e16),+1,1e-6,1,1,1)
ans =
  single
  9.999985e+15 +            0i

[remy-abergel]

Great, thank you very much for your kind message and for the upcoming fix!
Rémy

[DiamonDinoia]

@remy-abergel this should be fixed in master now. Could you let us know if you still experience any issue?

you could:

git clone finufft
cd finufft
pip install ./python/finufft

this assumes you have a c/c++ compiler

[remy-abergel]

I tried on multiple real-life examples (that were previously failing without normalization) and it worked like a charm. Many thanks to all of you for this fix 🙏