Wrong wavelet adjoint operator

[zaccharieramzi]

The adjoint operator of the wavelet transform operator is not correct in all cases. Namely, when using pywt wavelets, we have redundant wavelet transforms (currently hardcoded to be so because of #79 ). In this case, wavelets are not orthogonal and therefore the adjoint is not the inverse anymore.

For the case of isap it is also wrong but it might be because of accuracy issues cf #53 .

This paper by Folberth and Becker shows how to compute the wavelet tranform adjoint. The code is available in Matlab. The operator wextend might become available in pywt under the name pad after this issue is taken care of. It is however more or less already available.

Also remember that we will have to compute 2 adjoints: one for the decomposition and one for the reconstruction (since the inverse is not the adjoint of the decomposition, there is no reason that the decomposition is going to be the adjoint of the inverse). The article only studies the synthesis case, that is the adjoint of the reconstruction operator.

Finally, we will need some unit tests to test that the adjoint property is indeed verified.

/!\ Important: this probably affects badly the code base. Indeed, people probably used the adj_op function in lieu of the reconstruction. This would need to be reviewed.

[LElgueddari]

Could you propose a PR for that at least to have orthogonnal wavelets with peridization or using zero-padding to have a well defined adjoint?

[zaccharieramzi]

Orthogonal wavelets already have a well-defined adjoint (or at least a very good approximation) in the inverse.

The problem with zero-padding is that in the case of pywavelet it is non-orthogonal. Luckily, it turns out that for zero-padding the inverse is also a rather good approximation of the adjoint (also detailed in the article).

[philouc]

Did you mean: When writing
"since the inverse is adjoint of the decomposition, there is no reason that the decomposition is going to be the adjoint of the inverse " -> since the inverse is NOT adjoint of the decomposition, there is no reason that the decomposition is going to be the adjoint of the inverse.

You're right, we should get unittests soon to check this property and Loubna's right too. Please send a PR.

[zaccharieramzi]

Did you mean: When writing ...

Yes edited.

RE the PR, I will do one soon on the padding for the pywt wavelets allowing us to switch easily to zero-padding (or periodization). However, the real adjoint PR is going to be tricky and I cannot guarantee to be able to do it right away.

[philouc]

For sure. The real adjoint PR is more demanding and is not the highest priority right now.

[zaccharieramzi]

RE this issue, and the fact that maybe isap wavelet adjoints were wrong because of accuracy issues:
it's more severe than this. The adjoint is currently computed using the synthesis operator, which in the undecimated case is not a good approximation.

To compute the adjoint you need to take the analysis filters and rotate them (using the filter_rot arg in https://github.com/CEA-COSMIC/ModOpt/blob/master/modopt/signal/wavelet.py#L172).

[chaithyagr]

Update here:
I have written a class WaveletUD2 based on convolution with impulse response functions and filter_rot=True for adj_op . (Using helper functions from Modopt) Thanks @zaccharieramzi for help with initial understanding.
I have got some runs with this and things look fine, I will write a test for this and get a PR for this.

[chaithyagr]

This issue needs a fix first in pysap where we expose a function that does the adjoint operator, which later can lead to possibly some changes in code in pysap-mri to interface right.
@LElgueddari , as I don't have the right permissions, can you please instance one part of this issue in pysap?
We can hold this issue till that if fixed and till we have updated adjoint operator exposed.

[chaithyagr]

@AGrigis , can you please instance a version to get the above comment done in pysap?
I see that this issue was transferred from pysap

[chaithyagr]

Going forward in future, it would also be helpful if PySAP exposed the binaries of Sparse2D at say .local/pysap/bin
Currently, we have ongoing PR in pysap-mri where we use these binaries with the help of Modopt to use the Undecimated wavelets. A dirty fix for this is to clone the repo and build locally for travis tests.

However, we would need this updated for CEA-COSMIC/pysap-mri#40