0
$\begingroup$

(Note: Corrected code is posted below the original code.)

For an exercise in optimization, I am interested in coding a simple example from scratch: a Tikhonov denoising routine using gradient descent, that matches exactly in output the LSQ-based Tikhonov denoising result.

For regularisation of a noised signal $b$ to recover a denoised version $x$, I apply the typical formulation of Tikhonov

$$ \|A x - b\|^2 + \|\lambda x\|^2 $$

With $A$ as the identity matrix and $\lambda x$ the weighted gradient matrix $|\lambda \nabla x|$.

In the LSQ case this can be solved using $\hat{A}$ and $\hat{b}$ where $ \hat{A} = \left[ I ; \sqrt{\lambda} \nabla x \right] $ , $b = \left[b ; 0 \right] $ and $x = A \backslash b$. The code below shows this result which smooths the signal.

With the gradient descent version, I now try to minimize $f$ where

$$f(x) = 0.5 \|x-b\|^2 + \|\lambda \nabla x\|^2$$

I do this by stepping in the direction that is the negative of the function's gradient, which I believe to be

$$f'(x) = x - b + 2\lambda \nabla ^2 x $$

with $\nabla ^2$ the Laplacian operator. But this code produces rubbish. Can anyone spot my mistake?

Code:

function tik_1d

signal = sin((1:512)/128*2*pi);
signal = signal + randn(size(signal))*0.05;
lambda = 15;
niter = 50;
signal_lsq = tik_lsq(signal, lambda);
signal_opt = tik_opt(signal, lambda, niter);
figure(1); set(gcf, 'color', 'w');
subplot(3, 1, 1); plot(signal); title('b');
subplot(3, 1, 2); plot(signal_lsq); title('LSQ x');
subplot(3, 1, 3); plot(signal_opt); title('GD x');

end

function signal_lsq = tik_lsq(signal, lambda)
    grad = gradient_matrix(signal)*lambda; %subfunction see below
    A = [eye(numel(signal)); grad];
    s = signal';
    b = [s; zeros(size(s))];
    signal_lsq = A\b;
end

function m = gradient_matrix(signal)
    n = numel(signal);
    diag1 = ones(n,1);
    diag2 = -diag1;
    m = spdiags([diag1 diag2], [0 1], n, n);
end

function signal_opt = tik_opt(signal, lambda, niter)
    e = zeros(niter,1);
    signal_opt = signal;
    tau = 0.1;
    for n = 1:niter
        e(n) = eval_f(signal_opt, signal, lambda); %subfunction see below
        diff = tau*grad_f(signal_opt, signal, lambda); %subfunction see below
        signal_opt = signal_opt - diff;
    end
    figure(2); plot(e);
end

function e = eval_f(signal_rec, signal, lambda)
    signal_grad = conv(signal_rec, [1 -1], 'same');
    e = 0.5*norm(signal_rec-signal) + norm(lambda*signal_grad);
end

function g = grad_f(signal_rec, signal, lambda)
    signal_lap = conv(signal_rec, [1 -2 1], 'same');
    g = (signal - signal_rec) + lambda*(signal_lap);
end

Output:

Tikhonov denoising results

UPDATE: Here is the code with Christian Clason's suggestions incorporated:

function tik_1d

signal = sin((1:512)/128*2*pi);
signal = signal + randn(size(signal))*0.05;
lambda = 15;
niter = 500;
signal_lsq = tik_lsq(signal, lambda);
signal_opt = tik_opt(signal, lambda, niter);
figure(1); set(gcf, 'color', 'w');
subplot(4, 1, 1); plot(signal); title('b'); ylim([-1 1]);
subplot(4, 1, 2); plot(signal_lsq); title('LSQ x'); ylim([-1 1]);
subplot(4, 1, 3); plot(signal_opt); title('GD x'); ylim([-1 1]);
subplot(4, 1, 4); plot(signal_opt - signal_lsq'); title('LSQ x - GD x'); ylim([-1 1]);

end

function signal_lsq = tik_lsq(signal, lambda)
    grad = gradient_matrix(signal)*sqrt(lambda); %subfunction see below
    A = [eye(numel(signal)); grad];
    s = signal';
    b = [s; zeros(size(s))];
    signal_lsq = A\b;
end

function m = gradient_matrix(signal)
    n = numel(signal);
    diag1 = ones(n,1);
    diag2 = -diag1;
    m = spdiags([diag1 diag2], [0 1], n, n);
end

function signal_opt = tik_opt(signal, lambda, niter)
    e = zeros(niter,1);
    signal_opt = signal;
    tau = 0.001;
    for n = 1:niter
        e(n) = eval_f(signal_opt, signal, lambda); %subfunction see below
        diff = tau*grad_f(signal_opt, signal, lambda); %subfunction see below
        signal_opt = signal_opt - diff;
    end
    figure(2); plot(e);
end

function e = eval_f(signal_rec, signal, lambda)
    signal_grad = conv(signal_rec, [1 -1], 'same');
    e = 0.5*norm(signal_rec-signal) + norm(lambda*signal_grad);
end

function g = grad_f(signal_rec, signal, lambda)
    signal_lap = conv(signal_rec, [1 -2 1], 'same');
    g = (signal - signal_rec) - lambda*(signal_lap);
end

Output:

enter image description here

$\endgroup$
  • $\begingroup$ Careful with the sign for the Laplacian: $\nabla^T \nabla = -\mathrm{div} \nabla = -\Delta$. Your step size might also be too large; either try smaller and smaller ones until it works, or (better) use an Armijo line search. Your effective $\lambda$ for LSQ and gradient descent are also off by a factor of two (there shouldn't be $1/2$ in $f(x)$ if there isn't in the LSQ). $\endgroup$ – Christian Clason Sep 14 '16 at 16:29
  • $\begingroup$ Have you checked each of your functions separately to make sure that they compute the respective terms correctly? (For example, the gradient of a linear function should be constant and the Laplacian zero; the gradient of a sine should be a cosine etc.) $\endgroup$ – Christian Clason Sep 14 '16 at 16:44
  • $\begingroup$ @ChristianClason Surgically dissected thank you! It was all of these together -- a wrong sign in the Laplace operator, a step size that was much too large, and a spurious factor of 2. I have attached code incorporating these changes and the new output, which looks great. $\endgroup$ – barnhillec Sep 15 '16 at 7:51
  • $\begingroup$ Glad to be of help! I'll put these points in an answer so it can be marked as "answered" and won't keep popping up on the front page. $\endgroup$ – Christian Clason Sep 15 '16 at 9:54
2
$\begingroup$

There are a few problems with your gradient descent:

  1. The gradient of $\|\nabla u\|^2$ is $$\nabla^* \nabla u = -\mathrm{div} \nabla u = -\Delta u,$$ i.e., the negative Laplacian.

  2. Gradient descent does not converge unconditionally, but only if the step sizes are chosen appropriately (in particular, not too large). Either you reduce the step size manually until you have convergence, or implement a proper line search, e.g., based on Armijo or Wolfe conditions.

  3. You're not minimizing the same functional as in the LSQR approach, since there's a factor of $1/2$ in front of the first term which isn't there in the LSQR functional (which effectively makes your $\lambda$ twice as large).

Finally, a nitpick: although this is consistent in both approaches, your gradient and matrix $\hat A$ correspond to $\lambda \|\nabla u\|^2$, not $\|\lambda \nabla u\|^2$. (The former is the way it's usually written.)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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