With symmetric padding, the inverse wavelet transform is not the adjoint of the wavelet transform

[zaccharieramzi]

I recently stumbled upon a bothering fact when using the pywt. When we use the default "symmetric" padding, the inverse wavelet transform is not the adjoint of the wavelet transform (and it's not the mathematical inverse of the wavelet transform, "only in one direction").

However, if you use the "zero" padding the adjoint property is verified (not the inverse property). And if you use the "periodization" padding then all the properties are verified.

I am not sure whether this is a theoretical issue or if it's down to implementation.

You can find a code below highlighting the said issue:

import numpy as np
import pywt

wavelet_type = "db4"
level = 4
mode = "symmetric"

coeffs_tpl = pywt.wavedecn(data=np.zeros((512, 512)), wavelet=wavelet_type, mode=mode, level=level)
coeffs_1d, coeff_slices, coeff_shapes = pywt.ravel_coeffs(coeffs_tpl)
coeffs_tpl_rec = pywt.unravel_coeffs(coeffs_1d, coeff_slices, coeff_shapes)
_ = pywt.waverecn(coeffs_tpl_rec, wavelet=wavelet_type, mode=mode)

def py_W(im):
    alpha = pywt.wavedecn(data=im, wavelet=wavelet_type, mode=mode, level=level)
    alpha, _, _ = pywt.ravel_coeffs(alpha)
    return alpha

def py_Ws(alpha):
    coeffs = pywt.unravel_coeffs(alpha, coeff_slices, coeff_shapes)
    im = pywt.waverecn(coeffs, wavelet=wavelet_type, mode=mode)
    return im

x_example = np.random.rand(*coeffs_1d.shape)
y_example = np.random.rand(512, 512)
print("Adjoint:")
x_Tadj_y = np.dot(x_example, np.conjugate(py_W(y_example)))
T_x_y = np.dot(py_Ws(x_example).flatten(), np.conjugate(y_example.flatten()))
print(np.allclose(x_Tadj_y, T_x_y))
print("\n Inverse from image to image:")
print(np.allclose(py_Ws(py_W(y_example)), y_example))
print("\n Inverse from coefficients to coefficients:")
print(np.allclose(py_W(py_Ws(x_example)), x_example))

---

EDIT: When reading the documentation about the different paddings, I could see that the problematic paddings compute a redundant representation of the original image (for example you have more coefficients than pixels in your original image). Therefore, it's normal that the inverse property is not verified from coefficients to coefficients (if taken at random without structure). However, this doesn't explain why the adjoint property is not verified.

Also I don't know whether this redundance is theoretical to ensure perfect reconstruction or is an implementation choice.

[grlee77]

Yes, you are correct that the matrix representation for these other boundary modes is not square (so not orthogonal). I do not think it is possible to perform perfect reconstruction for the general boundary modes without retaining the extra coefficients.

The following publication on arXiv discusses how to determine an adjoint in these cases (as well as for the case of biorthogonal wavelets): https://arxiv.org/pdf/1707.02018.pdf

I did not get a chance yet to read it in detail, but they provide an implementation for Matlab. https://github.com/jamesfolberth/ieee_adjoints

The wavelet adjoint portion of the code in that repository is MIT licensed, so it would a valid candidate to adapt for use in PyWavelets. Our 1D and 2D DWT give the same coefficients as Matlab's so it is probably reasonably straightforward to adapt to convert their approach to Python. I have not yet looked how straightforward it would be to extend the approach to the n-dimensional case.

[zaccharieramzi]

Ok thanks, makes sense. So, if I sum it up, this redundancy is theoretically needed. Do you have by any chance a link to a paper explaining why?

Also, since we are on the topic, from what I understand to compute the adjoint operator, I would need to build its matrix representation. I see you haven't read the article in details (neither have I), but does it seem reasonable ? For large images it doesn't sound so.

[grlee77]

Do you have by any chance a link to a paper explaining why?

Not offhand. I think this is covered in some textbooks, possibly in the one by Strang and Nguyen, but I don't recall for sure: http://www-math.mit.edu/~gs/books/wfb.html

from what I understand to compute the adjoint operator, I would need to build its matrix representation

I do not think you need to explicitly form the matrix representation in practice. See for example this section of code in their demo: https://github.com/jamesfolberth/ieee_adjoints_mit/blob/master/wavelet/fval_driver.m#L119-L146

From looking at that code, It looks like they use in practice is something like: 1.) extend the signal using the desired boundary mode by an amount equal to wavelet.dec_len - 1 (see issue #473 about adding wextend to the PyWavelets public API). 2.) apply the wavelet transform to the extended signal using mode='zero'.

That should be reasonably efficient. You have to make a copy of the data with the boundaries extended, but do not have to explicitly form a full DWT matrix.

[zaccharieramzi]

Ok makes sense, I was only looking at this: https://github.com/jamesfolberth/ieee_adjoints_mit/blob/master/wavelet/adjoints/analysis_mat_adjoint_2d.m

Regarding #473 , do you think it would make sense to also directly provide in the API the adjoint operator, rather than letting the user code it? I don't have time right now to implement wextend in pywt straight away, but if you think that it is not urgent I would be happy to do it when I find the time.

Other than that I think we can consider this issue closed (I guess the conversation will now move to #473 ). Thanks for taking the time!

[grlee77]

Regarding #473 , do you think it would make sense to also directly provide in the API the adjoint operator, rather than letting the user code it?

I would certainly be open to having an example demonstrating adjoint vs. inverse and how to do this in the demos folder of the repository. I am not sure yet whether it makes sense to include it directly in the public API.

It seems like it would be appropriate to include in the WaveletTransformBase operator class as implemented in ODL. The adjoint method as currently implemented there doesn't seem to be complete.

[zaccharieramzi]

I would certainly be open to having an example demonstrating adjoint vs. inverse and how to do this in the demos folder of the repository.

Ok, I will try working on that.

Also, I actually have one last question on this topic: why can only the periodic padding be used for an orthogonal transform ("periodization" mode)? Is there a theoretical reason behind it?

[grlee77]

Also, I actually have one last question on this topic: why can only the periodic padding be used for an orthogonal transform ("periodization" mode)? Is there a theoretical reason behind it?

From a linear algebra perspective, if you were to write out the matrix equivalent of the the DWT for these other modes, it is non-square and so cannot be orthogonal. The rectangular shape is due to the need to retain some extra boundary coefficients.

There are some circumtances where symmetric and antisymmetric boundary modes can be performed without retaining any extra coefficients, but I think they require that the filter itself is also either symmetric or antisymmetric which is not the case for the filter banks in PyWavelets (aside from Haar/db1). There is a paper by C. Brislawn on the various possible "nonexpansive symmetric extensions": https://doi.org/10.1006/acha.1996.0026

[robtovey]

As a follow-up to this, can I request/suggest that the documentation is updated to make it clear that some boundary modes do not behave like 'normal wavelets'? Maybe I'm an idealist, but my default expectation of wavelet transforms is that they are orthogonal. The documentation just refers to the transforms as '(Inverse) Discrete Wavelet' transforms. My understanding from here is that, depending on the boundary mode, sometimes the 'inverse' is not the adjoint and maybe not even the exact inverse? These properties depend almost entirely on the mode and not the wavelet itself? As the OP says, the documentation for boundary modes only discusses how DWT is effected. It would be nice if the guarantees of the DWT/IDWT pairs could be confirmed somewhere, at least to draw attention to the fact that the modes change more properties than described at the moment.