Solved – Variational method — implementation of function gradient for image denoise

gradient descentimage processingMATLABoptimization

I am studying Variational Methods and decided to code a simple matlab implementation of the famous image denoise example:

$E(u) = \displaystyle\frac{1}{2}\sum_{i=1}^N (u_i – f_i)^2 + \displaystyle \frac{\lambda}{2}\sum_{i=1}^N | \nabla u_i|^2$

where the first term (data) measures how close the denoised image $u$ is to my original image $f$ and the second term is the regularization which preserves spatial smoothness.

Nevertheless, I believe there is something wrong with the way I calculate the gradient. Basically, by taking the $\frac{dE}{du}$, it should be just the sum of data and regularization terms, isn't it right?

options = optimset('GradObj','on','MaxIter',100000);
% image
image = randi(100,10,10);

% initialGuess
u = ones(size(image));

% weighting on the regularization
lambda = 0.125;

% optimization
[optTheta] = fminunc(@cost_function, u, options, lambda, image);

where my cost_function.m is implemented as follows:

function [val,grad] = cost_function(u, lambda, f)

    data_term = u - f;
    smooth_term = gradient(u);

    val = 0.5*norm(data_term) + 0.5*lambda*norm(smooth_term);
    grad = data_term + lambda*smooth_term;

end

The problem is in the gradient computation for the optimization algorithm. If I add the parameter 'DerivativeCheck','on' in the options, then matlab complains that the gradient is not correct:

____________________________________________________________
DerivativeCheck Information

Objective function derivatives:
Maximum relative difference between user-supplied 
and finite-difference derivatives = 0.999536.
  User-supplied derivative element (58,1):     -80.0002
  Finite-difference derivative element (58,1): -0.0371037
____________________________________________________________
Error using validateFirstDerivatives (line 96)
DerivativeCheck failed:
User-supplied and forward finite-difference derivatives do not match within 1e-06 relative tolerance.

Best Answer

As noted in the other answer, your problem is not using both outputs from gradient. However, here I want to note important issues with the formulation.

First, your discrete problem can be written as a standard linear least squares problem. That is, your objective function is (half of) $$E[u] = \|u-f\|^2 + \lambda(\|\Delta_xu\|^2+\|\Delta_yu\|^2)$$ where now $u$ and $f$ are vectors created by flattening the corresponding image arrays (i.e. u=U(:) and f=F(:)). This is a sum of squares, which is equivalent to the least-squares residual $$ \begin{bmatrix} I \\ \sqrt{\lambda}D \end{bmatrix} u \approx \begin{bmatrix} f \\ 0 \end{bmatrix} \,,\, D = \begin{bmatrix} D_x \\ D_y \end{bmatrix} $$ where $I$ is the identity matrix, $D$ is a discrete gradient operator. The least squares solution of the system $Au=b$ can then be computed as u=A\b, using Matlab's backslash operator.

Second, for this problem the discrete gradient would typically be computed using forward differences (rather than the centered differences used by gradient). This is because the least squares solution of the above problem is $$(I+\lambda{D^T}D)u=f$$ With 1st order differences, the regularization matrix $D^TD$ is then a 2nd order discretization of the Laplacian $\nabla^2u$. The equation above is then a discrete Poisson equation, and moreover is consistent with the continuous (variational) form of the optimization problem.

If all of the above is done with sparse matrices, then the solution is very efficient, and with no iteration required. As a side note, this can all be done in Matlab using gridfit (with the "springs" option).

For completeness, below is an example Matlab implementation.

[ny,nx]=size(F);
Dfnc=@(sz)spdiags(ones(sz(1),1)*[-1,1],0:1,sz(1)-1,sz(1));
Dy=kron(eye(nx),Dfnc(ny,nx)); Dx=kron(Dfnc(nx,ny),eye(ny)); D=[Dx;Dy];
% Fx=diff(F,[],2); Fy=diff(F,[],1); Derr=norm(D*F(:)-[Fx(:);Fy(:)]) % Derr=0
I=speye(ny*nx); A=[I;lam*D]; b=[F(:);zeros(size(D,1),1)];
U=zeros(ny,nx); U(:)=A\b;
Related Question