astro-informatics / s2fft

Differentiable and accelerated spherical transforms with JAX
https://astro-informatics.github.io/s2fft
MIT License
136 stars 9 forks source link

Wigner FFT + IFFT is not the identity for 'DH' sampling and low bandwidths #242

Open PhilippMisofCH opened 1 hour ago

PhilippMisofCH commented 1 hour ago

First of all, thanks again for this nice library. While working with it, I've noticed some unexpected results in my project, which eventually boil down to the fact that applying the inverse Wigner FFT on the Wigner FFT does not yield the original input, not even at low bandwidth (L=8) and dh sampling (This is actually one of the only cases that I've tested thoroughly). In contrast, an equivalent setup for the $S^2$ transforms seems to work as expected.

So I wondered if this is indeed a problem, or it is to be expected.

I've also attached a small example where I test this out:

L = 8
sampling = 'dh'
a = s2fft.sampling.s2_samples.phis_equiang(L, sampling=sampling)
b = s2fft.sampling.s2_samples.thetas(L, sampling=sampling)
g = s2fft.sampling.s2_samples.phis_equiang(L, sampling=sampling)
f_grid = jnp.meshgrid(g, b, a, sparse=True, indexing='ij')
f_s2 = jnp.squeeze(jnp.sin(f_grid[1]) * jnp.cos(f_grid[2]))
flm = s2fft.transforms.spherical.forward_jax(f_s2, L, sampling=sampling, reality=True)
f_s2_back = s2fft.transforms.spherical.inverse_jax(flm, L, sampling=sampling)
assert jnp.allclose(f_s2, f_s2_back, atol=1e-4), ("Spherical forward + inverse transform is"
                                                  " not the identity")

f_so3 = jnp.sin(f_grid[0]) * jnp.cos(f_grid[1]) * jnp.sin(f_grid[2])
flmn = s2fft.transforms.wigner.forward_jax(f_so3, L, L, sampling=sampling, reality=True)
f_back = s2fft.transforms.wigner.inverse_jax(flmn, L, L, sampling=sampling)
assert jnp.allclose(f_so3, f_back, atol=1e-2), ("Wigner forward + inverse transform is not "
                                                "the identity")

The first assertion passes, but the second one fails.

I've also looked around at other reported issues and found e.g.

It would be great if you could tell me if I misunderstood something or if this is indeed a bug / numerical issue. Thanks in advance!

PhilippMisofCH commented 1 hour ago

Adding on this: At first sight it looks like one of the notebooks is demonstrating explicitly that this works, but here, one starts in the Fourier space, applies the inverse transformation and then goes back to the Fourier space, where all values agree.

Similarly, doing the following also passes the assertion:

rng = np.random.default_rng(83459)
flmn = s2fft.utils.signal_generator.generate_flmn(rng, L, L, reality=True)
f_so3 = s2fft.transforms.wigner.inverse_jax(flmn, L, L, sampling=sampling)
flmn = s2fft.transforms.wigner.forward_jax(f_so3, L, L, sampling=sampling, reality=True)
f_so3_back = s2fft.transforms.wigner.inverse_jax(flmn, L, L, sampling=sampling)
assert jnp.allclose(f_so3, f_so3_back, atol=1e-8), ("Wigner forward + inverse transform is not "
                                                "the identity")

But for an arbitrary (smooth) signal generated on $SO(3)$ this does not work as shown before.