%% 1D inverse problems regularization and the SVD
% In this tutorial we explore some one-dimensional inverse problems, and look at the singular values and regularization. We use Per Christian Hanson's Regularization Tools (reg tools). Follow this link for the code and manual <http://people.compute.dtu.dk/pcha/Regutools/index.html> The manual itself is a nice tutorial for numerical treatment of 1D inverse problems

addpath RegTools  % use the subdirectory where you downloaded reg tools

%
% The example deriv2 is esentially and inverse problem where the forward problem involves integrating twice, so the discrete inverse problem is differentiating twice numerically. Which we expect t be unstable. First set up the problem. In general the examples in reg tools generate a matrix A, a data vector b and a solution x. We can then play with solving Ax=b and see what happens when we add noise to b etc.

n=100 %Size of problem
example=3;
[A,b,x]=deriv2(n,example); % The second derivative inverse problem
size(A)
norm(A*x-b) % we see Ax is close to b
plot(x);hold on 
plot(b)
hold off
title(['x and b for example number ',num2str(example)])

%
% The matrix A is n by n, and the condition number measures who hard it is to invert

cond(A)

%
% In tis case about 10^4 so resonably hard.

plot(x)
title(['The original x for deriv2 with n=',num2str(n),' for example number ',num2str(example)])
x0=A\b;
plot(x0)
title('Soution x0 calculated with backslash')
norm(A*x0-b)% error in solution
plot(x-x0)
title('Difference between actual and calculated solution')

%
% Suggestion: try different values of n and example
% Clearly this is a fairly good solution but we did not deliberately add any noise to b (ie we did an 'inverse crime').
% Now we add some gaussian pseudo random noise with standard deviation std*norm(b)

std= 0.001;
bnoise = b + std*randn(size(b))*norm(b);
plot(b); hold on
plot(bnoise)
hold off


x1=A\bnoise;
plot(x1)
title(['Soution x1 with ',num2str(std*100),' percent noise in b'])
norm(A*x1-b)% error in solution
plot(x-x1)
title('Difference between actual and calculated solution')

%
% Suggestion: Try changing std to see what level of noise gives a recogisable solution.
% Now lets try one of the simplest tricks for solving an inverse problem: Tikonov regularization.

TikFac=5e-3
xtik=(A'*A+ TikFac^2 *eye(n))\(A'*bnoise);
plot(xtik)
title(['Tikhonov reg soln  with Tikhonov factor ',num2str(TikFac)])
norm(A*xtik-b)% error in solution
plot(x-xtik)
title('Difference between actual and Tikhonov solution')

%
% Suggestion: Try changing the Tikhonov factor TikFac. Too small and the solution is jagged, to big and it is over smooth.

[U,S,V]=svd(A); %Calculate singular values and vectors
semilogy(diag(S)) % plot 
title('Singular values of deriv2 on a log linear scale')
figure, loglog(diag(S)) %plot
title('Singular values of deriv2 on a log log scale')