Wrong fourier adjoint operator

[zaccharieramzi]

The ajdoint operator of the fourier operator is not correct. We need to normalize by the square root of the dimension of the data (or if we choose to use the numpy package -more complete- rather than the scipy one, we can just use the "norm" kwarg).

To see this, run the following experiment:

import numpy as np

from pysap.numerics.fourier import FFT2
from pysap.numerics.utils import convert_mask_to_locations

dim_1d = 512
dim = (dim_1d, dim_1d)

mask = np.ones(dim)
samples = convert_mask_to_locations(mask)

fourier_op = FFT2(samples=samples, shape=dim)

def fft(x):
    return fourier_op.op(x)

def ifft(x):
    return fourier_op.adj_op(x)

fft_example = np.random.rand(*dim) + 1j * np.random.rand(*dim)  # x
im_example = np.random.rand(*dim)  # y
print("Adjoint:")
x_Tadj_y = np.dot(fft_example.flatten(), np.conjugate(fft(im_example).flatten()))
T_x_y = np.dot(ifft(fft_example).flatten(), np.conjugate(im_example.flatten()))
print(np.allclose(x_Tadj_y, T_x_y))

I guess the easiest solution is to switch to numpy and use the "norm" kwarg.

We also need one (or more) unit test(s) to make sure it doesn't happen. (Not sure whether the one above can be used because of the randomness but sthg similar should do).

I haven't yet checked for NFFT2 but it's probably the same situation. If it is, we need to correct it as well (can probably be done in one PR).

/!\ Important: I am not sure how this affects all the existing code base (namely whether all the examples would still work, whether in some places the correction has already been made manually).

[LElgueddari]

Have you checked the unittest of this Fourier operator, it is just a matter of normalisation.. and this subject have already been discussed in the PR #56, and they asked to keep it like that so basically the normalization will all be in the IFFT operator. Moreover the NFFT and the FFT have been checked so as they give exactly the same result.

[zaccharieramzi]

Well but still it holds that the adjoint is wrong because it doesn't verify the definition.

If we want to keep the fourier operator non-unitary, then the adjoint needs to be computed differently than just the inverse fourier transform.

[LElgueddari]

No but the inverse Fourier transform is really the adjoint Fourier transform. We just need to pay attention when you create the fft_exemple.
And actually a unittest that tested this is already existing in pysap and running.
Moreover, the topic of symmetric Fourier operator has already been discussed in the PR #56 and it came out that we don't need to change this.
Don't forget that Fourier is linear which means that dividing the result of fft2 or the input of fft2 by the same number will give the extact same solution

[zaccharieramzi]

No it's not when the normalization is not done correctly, the code example I wrote proves that.

In the unit test you can by the way see that there is a renormalization happening, hiding the fact that the inverse is not the adjoint as coded in pysap (also btw the unit test is not complete since it only focuses on real images and real fourier coefficients).

[zaccharieramzi]

Re #56 , I don't see where the discussion occured on the matter of normalization.

I agree that fourier is linear, however multiplying the result of the fft2 and/or the ifft2 by a scalar changes whether one is the adjoint of the other (see for example this). Indeed you need the transformation to be unitary for that to be true.

[LElgueddari]

Here is the conversation I was talking about.

[zaccharieramzi]

Ok so I think the documentation is also wrong. Because we don't normalize in a symmetric way in our case.

Using the notations of the example I gave you can just do:

print(np.linalg.norm(fft(im_example)) / np.linalg.norm(im_example))
print(np.linalg.norm(ifft(fft_example)) / np.linalg.norm(fft_example))

The reason I would be against keeping the original fft convention is because then we would need 4 operators in the FFT2 class: the op, adj_op, inv_op, adj_inv_op. This is, I think, overkill given that we can keep only 2 operators if we stick to a symmetric normalization (which is not currently the case in pysap).

---

EDIT:
I had read the documentation from an earlier commit, the comment about the documentation is therefore wrong. The rest still stands.

[LElgueddari]

The PR is still open we can still change this one if @sfarrens et @AGrigis agreed on this?

[sfarrens]

No problem for me. Just let me know if you plan to commit any further changes. If not, I can probably merge this today.

[LElgueddari]

Yes I'm planing on addind some commmit related to this issue

[zaccharieramzi]

I am closing this following Loubna's PR merge #56 .