Inverse Problems Lecture 3/2017: deconvolution with truncated SVD, part 2/2

and

This is screen capture of Matlab programming I did when teaching my course Inverse Problems at University of Helsinki.

The lecture was given at January 25, 2017.

Here is the final code:

% Simple illustration of deconvolution in 1D, the method we use is
% Truncated Singular Value Decomposition
%
% Samuli Siltanen January 2017

% Parameters for controlling the plot appearance
fsize = 16;
lwidth = 2;

%% Simulate the measurement
% Build a Point Spread Function (PSF)
M = 17;
psf = ones(1,2*M+1); % This makes sure psf has a unique centre
psf = psf/sum(psf); % Normalization of the psf

% Construct "unknown signal" f
N = 400;
x = linspace(0,1,N);
% f = zeros(N,1);
% f(1:(end/2)) = 1;
f = sin(2*pi*x);
f = f(:); % Force vector f to be vertical

% Construct the convolution matrix
A = convmtx(psf,N);
A = A(:,(M+1):(end-M));

% Simulate the "measurement"
m = A*f;

% Simulate "noisy measurement"
sigma = .1;
mn = A*f + sigma*randn(N,1);

%% Compute truncated SVD solution

% Determine the SVD of matrix A
[U,D,V] = svd(A);
svals = diag(D);

% Compute reconstruction
r_alpha = 15;
Dp_alpha = zeros(size(A));
for iii = 1:r_alpha
Dp_alpha(iii,iii) = 1/svals(iii);
end
f0 = V*Dp_alpha*(U.')*m;
fn = V*Dp_alpha*(U.')*mn;

%% plots
% Take a look at the matrix and its singular values
figure(2)
clf
subplot(1,2,1)
spy(A)
title('Nonzero elements in A','fontsize',fsize)
subplot(1,2,2)
semilogy(svals,'k.')
hold on
semilogy(svals(1:r_alpha),'r.')
title('Singular values of A (log plot)','fontsize',fsize)

% Take a look at few first singular vectors
figure(3)
clf
subplot(5,1,1)
plot(V(:,1))
title('Singular vector 1','fontsize',fsize)
subplot(5,1,2)
plot(V(:,2))
title('Singular vector 2','fontsize',fsize)
subplot(5,1,3)
plot(V(:,3))
title('Singular vector 3','fontsize',fsize)
subplot(5,1,4)
plot(V(:,4))
title('Singular vector 4','fontsize',fsize)
subplot(5,1,5)
plot(V(:,5))
title('Singular vector 5','fontsize',fsize)

% Take a look
figure(1)
clf
subplot(3,1,1)
plot(f,'k','linewidth',lwidth)
hold on
plot(mn,'r','linewidth',lwidth)
set(gca,'ytick',[0 max(f)],'fontsize',fsize)
subplot(3,1,2)
plot(f,'k','linewidth',lwidth)
hold on
plot(f0,'b','linewidth',lwidth)
set(gca,'ytick',[0 max(f)],'fontsize',fsize)
subplot(3,1,3)
plot(f,'k','linewidth',lwidth)
hold on
plot(fn,'b','linewidth',lwidth)
set(gca,'ytick',[0 max(f)],'fontsize',fsize)

[saygnileri]

Hi Prof. Siltanen, thank you very much for this amazing video. I plan to follow all the courses and now I am at the begining. I have a question regarding computing tsvd. Dp_alpha(iii, iii) = 1/svals(iii);
Why is it not like this Dp_alpha(iii, iii) = svals(iii); ?