How to use DFT to solve this minimization problem? 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.
 A: 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. 
A: 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.
