2
$\begingroup$

A basic deconvolution problem is defined by the cost function $C(f) = ||f \star h - g||^2_2$, where $f$ and $g$ are unknown latent and measured images, $h$ is known kernel. I want to derive the gradient $\frac{dC(f)}{df}$ to use later for the optimization purposes to obtain $f$. The approach is to express the convolution through a circulant operator $A$: $C(f) = ||A(h) \cdot f - g||^2_2$. It is very well known that diagonalization of $A$ gives:

$A = \frac{W^{\dagger}}{\sqrt{N}} \cdot D \cdot \frac{W}{\sqrt{N}}$, where $W$ and $W^{\dagger}$ are DFT and IDFT matrices, $D$ is diagonal containing Fourier transform of $h$.

Now $\frac{dC(f)}{df}=2 \cdot A^T \cdot (A \cdot f - g)=2 \cdot \frac{W}{\sqrt{N}} \cdot D \cdot \frac{W^{\dagger}}{\sqrt{N}} \cdot (A \cdot f - g)$. Numerically this can be efficiently computed as

$2 \cdot \texttt{FFT(FFT}(h) \cdot \texttt{IFFT}(A \cdot f - g))$.

However, looking at other references I found that one has to compute

$2 \cdot \texttt{IFFT(FFT}^*(h) \cdot \texttt{FFT}(A \cdot f - g))$.

One may actually notice that this is a complex conjugate of my result. However, I don't understand where my derivation went wrong.

$\endgroup$

1 Answer 1

1
$\begingroup$

A simple method to derive the above is using the matrix form of the convolution.
Since a convolution is a linear operator then it can be expressed using a matrix product:

$$ \boldsymbol{h} \ast \boldsymbol{f} = \boldsymbol{H} \boldsymbol{f} $$

Now, the objective function can be rewritten:

$$ C \left( \boldsymbol{f} \right) = \frac{1}{2} {\left\| \boldsymbol{H} \boldsymbol{f} - \boldsymbol{g} \right\|}_{2}^{2} $$

Then the gradient is given by:

$$ \nabla C \left( \boldsymbol{f} \right) = \boldsymbol{H}^{T} \left( \boldsymbol{H} \boldsymbol{f} - \boldsymbol{g} \right) $$

The question is, what's the operator behind the $\boldsymbol{H}^{T}$ matrix?
It turns out to be the correlation operator. Namely the convolution with the reflected kernel.

So it can be written:

$$ \nabla C \left( \boldsymbol{f} \right) = \boldsymbol{h} \star \left( \boldsymbol{h} \ast \boldsymbol{f} - \boldsymbol{g} \right) $$

You may have a look at The Gradient / Derivative of Least Squares of 2D Image Convolution.

Specifically, in your derivation, you did the whole analysis on the real numbers. Yet the diagonal of $\boldsymbol{D}$ are complex numbers (The DFT of the kernel). Hence the transpose Becomes a conjugation.

So, if $\boldsymbol{H} = \boldsymbol{W}^{H} \boldsymbol{D} \boldsymbol{W}$ then $\boldsymbol{H}^{T} = \boldsymbol{H}^{H} = \boldsymbol{W}^{H} \boldsymbol{D}^{H} \boldsymbol{W}$ which means it is conjugate in the complex domain which means reflection in time / spatial domain (Which is equivalent to correlation).

Remark: The analysis is done over discrete signals.
Remark: Diagonalization of a convolution matrix is done specifically for the case of periodic boundary conditions. Then the matrix becomes a circulant matrix which is diagonalizable by the DFT matrix.

$\endgroup$
9
  • $\begingroup$ Thanks! I think this is precisely the last formula in my question $2 \cdot IFFT(FFT^*(h) \cdot FFT(A \cdot f - g))$. But I also wanted to do detailed derivation and as you can see got $2 \cdot FFT(FFT(h) \cdot IFFT(A \cdot f - g))$, which is just complex conjugate of the first result. Since the derivative should be real then these two results should be equivalent. I checked it in Matlab. The first result is real but the latter is complex. However, their real parts are identical up to close to eps numbers . And the complex part of the latter result is not much larger then eps either. $\endgroup$
    – baronett
    Commented May 10 at 16:46
  • $\begingroup$ @baronett, I wrote why you missed the conjugate. See the passage below the link. $\endgroup$
    – Royi
    Commented May 10 at 16:48
  • $\begingroup$ @baronett, I'd also mention that many papers are sloppy about this and use this trick and neglect the fact it assumes circular convolution and not linear. In practice, to make it equivalent one must pad the input appropriately. $\endgroup$
    – Royi
    Commented May 10 at 16:52
  • $\begingroup$ alright, I see your point. But I don't quite understand why for complex matrices $(AB)^T$ is not equal to $B^T*A^T$ but to $B^H*A^H$? Just checked it Matlab, the transpose equality still holds. Also $H$ is real, shouldn't $H^T$ be equivalent to $H^H$? $\endgroup$
    – baronett
    Commented May 10 at 17:15
  • $\begingroup$ @baronett, Indeed for real matrix $\boldsymbol{H}^{T} = \boldsymbol{H}^{H}$ as I wrote. Yet when you use diagonalization of a real matrix, unless it is symmetric, you're moving to the complex field of matrices which then you must use the complex conjugate operator as I did above. $\endgroup$
    – Royi
    Commented May 10 at 17:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.