# How to code gradient descent-based Tikhonov denoising that exactly matches LSQ Tikhonov denoise?

(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)
s = signal';
b = [s; zeros(size(s))];
signal_lsq = A\b;
end

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');
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:

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)
s = signal';
b = [s; zeros(size(s))];
signal_lsq = A\b;
end

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');
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:

• 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). – Christian Clason Sep 14 '16 at 16:29
• 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.) – Christian Clason Sep 14 '16 at 16:44
• @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. – barnhillec Sep 15 '16 at 7:51
• 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. – Christian Clason Sep 15 '16 at 9:54

1. The gradient of $\|\nabla u\|^2$ is $$\nabla^* \nabla u = -\mathrm{div} \nabla u = -\Delta u,$$ i.e., the negative Laplacian.
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.)