Adjoint or transpose for CDF 9/7 basis

[sparrowhawker]

Since the CDF 9/7 basis is not orthogonal, the inverse transform matrix is not the adjoint (or transpose) of the forward transform matrix. Would it be possible to modify the idwt() function, such that the use of idwt() with an adjoint=true flag simply returns the action of the transpose of idwt with the cdf 9/7 basis on a given vector? This would be very useful for change of basis problems where the adjoint of the inverse transform from cdf 9/7 to cartesian is required in applying the chain rule.

[gummif]

The lifting transforms are implemented as a sequence of update and predict steps. I'm not sure how you would implement the transpose transform. Do you have any sources for an algorithm for this?

[sparrowhawker]

I think it should be simple enough (in principle!) to implement if I can explain it right through the so called dot product test. I provide a reference at the end, but I can't vouch for its correctness as I don't quite understand the wavelet transform beyond looking at it from a linear algebra perspective. But here goes, I think I learnt a lot just trying to explain this!

In the case of DWT and IDWT with an orthogonal basis, we can think of the forward and inverse operations as a matrix A and A' acting on any random vectors x and y respectively.

Ax and A'y are implemented by dwt(x, wavelet(WT.haar)) and idwt(y, wavelet(WT.haar)). For two vectors x and y and a matrix A the following code illustrates the dot product test, to make sure that the correct adjoint (fancy word for a matrix transpose) is implemented for the forward transform.

using Wavelets
wt_haar = wavelet(WT.haar)
wt_db9 = wavelet(WT.db9)
wt_CDF = wavelet(WT.cdf97, WT.Lifting)

srand(121)
y = rand(128, 128)
x = rand(128, 128)

# define dot product test:
# essentially to check that the scalar y'Ax == (A'y)'x
# i.e., a test to check that <y, Ax> == <A'y, x>
# this checks that the linear operations of Ax and A'y
# are correctly implemented when A and A' are
# too large to explicitly form
dot_product_test(x,y,wt) = vecdot(y, dwt(x, wt)) - vecdot(idwt(y, wt), x)

# check with haar
dh = dot_product_test(x, y, wt_haar)
# check with db9
ddb = dot_product_test(x, y, wt_db9)
# check with cdf
dcdf = dot_product_test(x, y, wt_CDF)

For me, dh, ddb and dcdf are 4.973799150320701e-14, 4.085620730620576e-13, and -3.066038748151467,

Thus indicating that the Haar and Daubechies transforms satisfy the dot product test, but CDF9/7 does not. This is because the CDF transforms are not orthogonal, and the inverse transform is not simply the adjoint (transpose) of the forward transform.

This dot product test is simply a consequence of matrix algebra, and should hold for any linear operator A. For example, we could do a discrete Fourier transform using a matrix A with as many rows as frequencies (say m) to analyze and as many columns as time samples (say n). the first and last rows of A would be

exp(2pi*f1*1im*t1) ... exp(2pi*f1*1im*tn)
and
exp(2pi*fm*1im*t1) ... exp(2pi*fm*1im*tn)

the adjoint matrix, would simply be A' . However, we never actually do perform fft using this method with an explicit matrix multiply. In fact, if you run the dot product test using fft() and ifft() in Julia, you'll get a similar result as the dot product test above, with a normalization due to the scaling the forward transform with the number of time points. To illustrate this, using the same x and y as defined in the above code for the dot product test,

vecdot(y,fft(x)) - vecdot(ifft(y)*length(x[:]), x)
4.501998773775995e-11 - 2.9558577807620168e-12im

returns close to machine precision zero. So I would think of the DWT and IDWT operation in a similar fashion. The adjoint operation for IDWT with the biorthogonal basis can be given by the same code which produces the IDWT, with the rows and columns switched in the algorithm for IDWT.

Here is an example I found on github written in C but I can't vouch for its correctness:
https://github.com/ahay/src/blob/master/api/c/wavelet.c

Many thanks!

[gummif]

The lifting transform is implemented as a sequence of simple transforms, something like A = TRSU. So finding the adjoint simply means finding the adjoint for each step, A' = U'S'R'T'. I'm not sure when I will have time to work this out but you are welcome to take a look at it.