How to use DFT to solve this minimization problem?

Question

This is a problem when I'm reading a paper.

Equation:

$min\{\sum_p(S_p-I_p)^2+\beta((\partial_xS_p-h_p)^2+(\partial_yS_p-v_p)^2) \} $

where $S,I,h,v$ are all $M*N$ matrices and p stands for every element in the Matrix. $I,h,v$ are known.

The paper just mentioned "we diagonalize derivative operators after Fast Fourier Transform for speedup" and get the solution

$S=\mathscr{F}^{-1}\left(\frac{\mathscr{F}(I)+\beta(\mathscr{F}(\partial_x)^*\mathscr{F}(h)+\mathscr{F}(\partial_y)^*\mathscr{F}(v))}{\mathscr{F}(1)+\beta(\mathscr{F}(\partial_x)^*\mathscr{F}(\partial_x)+\mathscr{F}(\partial_y)^*\mathscr{F}(\partial_y)}\right)$

where $\mathscr{F}(1)$ stands for FFT of delta function. Plus, multiplication and division are all component-wise. "*" means conjugation.

I've looked up some books but still can't get how this happen. I don't find any connection between DFT and minimization or quadratic function.

Here's the paper, equation is on page 4.

Answer 1

You need to minimize the objective function

$$f:S \mapsto \| S-I\|_2^2+\beta(\| \partial_xS-h\|^2_2+\|\partial_yS-v\|_2^2) $$

where $\| \cdot \|_2$ is the "entrywise" $\ell^2$ norm. This is done as usual: the global minimum $S$ must be a critical point and then the derivative must be zero at $S$. The differential of the objective function is given by

$$Df(S)(M) = 2\langle S-I,M\rangle + \beta\left( 2 \langle\partial_x S -h,\partial_x M\rangle + 2 \langle\partial_y S - v,\partial_y M \rangle \right)$$

where $\langle A,B \rangle = \sum_p A_p B_p$ (the "entrywise" inner product). Now, that for any $M$ the last expression vanishes implies that

$$ S-I + \beta( \partial_x^T \partial_x S-\partial_x^Th + \partial_y^T \partial_y S -\partial_y^Tv) =0.$$

($\cdot^T$ is the adjoint operator). Now you can take the DFT on both sides,

$$ \mathscr{F}(S) - \mathscr{F}(I) + \beta \left( \mathscr{F}(\partial_x^T\partial_x S) - \mathscr{F}(\partial_x^Th)+ \mathscr{F}(\partial_y^T\partial_y S) - \mathscr{F} (\partial_y^Tv)\right)=0,$$ and then use the fact that the DFT diagonalizes the gradient $\mathscr{F}(\partial_. S) = \mathscr{F}(\partial_\cdot ) \mathscr{F} (S)$ and some algebra to find that $$ (\mathscr{F}(1) + \beta (\mathscr{F}(\partial_x)^*\mathscr{F}(\partial_x) + \mathscr{F}(\partial_y)^*\mathscr{F}(\partial_y) ))\mathscr{F}(S)= \mathscr{F}(I) + \beta \left( \mathscr{F}(\partial_x)^*\mathscr{F}(h) + \mathscr{F}(\partial_x)^*\mathscr{F} (v)\right)$$ and solving for $S$ this gives the formula asked. (Remark: We could have written the objective function directly in terms of $\mathscr{F}(S)$ using Parseval and then optimize).

As for the speedup part, we would have to invert a very large ($NM \times NM$) matrix representing the operator $$ 1 + \beta( \partial_x^T \partial_x + \partial_y^T \partial_y )$$ if we do it directly in image space.

Answer 2

FFT is used only to "speed-up". If no need to speed-up this formula looks like standard minimal least square answer.

Which works like this - assume you need to minimize in x : |Mx-v|^2 for some rectangular matrix M. The answer is well-known x= inverse(M*M^t)*M V. Remark: The matrix M * inverse(M*M^t) is called Moore-Penrose pseaudoinverse.

Now what is M in yours example ?

M = 1 + d/dx + d/y

So M*M^t - will be what stands in denominator, while Mv - this what stands in numerator.

If you want to speed-up - put "F" everywhere.

I am not sure my answer is clear. Please comment if you need more details.