Example1: Sparse Deconvolution

Sparse deconvolution by L1 norm minimization using majorization-minimization and fast solvers for banded systems.

Ivan Selesnick, Polytechnic Institute of New York University, October 2012

Contents

Start

clear
close all
printme = @(filename) print('-dpdf', sprintf('figures/Example1_%s', filename));

Create data

N = 300;                                % N : signal length
x = zeros(N,1);
x(50) = 1;
x(80) = 0.5;
x(100) = -0.3;

b = [1 0.8];                            % Define filter
r = 0.9; om = 0.95;
a = [1 -2*r*cos(om) r^2];
hx = filter(b,a,x);

h = filter(b, a, [1 zeros(1,N)]);       % Calculate impulse response.

randn('state',0);                       % Set seed to example reproducible
sigma = 0.2;
noise = sigma*randn(N,1);
y = hx + noise;

figure(1)
clf
subplot(2,1,1)
plot(x)
box off
title('Sparse Signal');

subplot(2,1,2)
plot(hx)
box off
title('Convolution');

printme('original')
3*sigma*sqrt(sum(abs(h).^2))
ans =

    2.0091

figure(2)
clf
subplot(2,1,1)
plot(y);
box off
title('Noisy data');

printme('data')

Verify banded matrices

Verify that the banded matrics (A,B) implement the difference equation. It should be consistent with the 'filter' commaned.

Nb = length(b);
Na = length(a);
e = ones(N, 1);
B = spdiags(b(ones(N,1), :), 0:-1:1-Nb, N, N);
A = spdiags(a(ones(N,1), :), 0:-1:1-Na, N, N);

H = @(x) A\(B*x);                                       % filter

err = H(noise) - filter(b,a,noise);
max(abs(err))

% This verifies that inv(A)*B*x is the same as filter(b,a,x)
ans =

   4.4409e-16

Sparse deconvolution

lam = 2.0;                                      % lam : regularization parameter
Nit = 100;                                      % Nit : number of iterations
[x_, cost] = deconvL1(y, lam, b, a, Nit);       % Run algorithm

figure(2)
clf
subplot(2,1,1)
plot(1:Nit, cost, '.-')
box off
title('Cost function history');
xlabel('Iteration')
xlim([0 20])

printme('cost')

figure(1)
clf
subplot(2,1,1)
plot(x_)
ylim([-0.5 1])
box off
title('Sparse deconvolution');
axnote(sprintf('\\lambda = %.2f', lam))

printme('deconv_L1')

Verify optimality conditions

Verify that the solution produced by 'deconvL1' satifies the optimality conditions. The points in the scatter plot should lie on the dashed lines.

HT = @(x) B'*(A'\x);            % filter transpose H'
z = HT(y - H(x_));

figure(1)
clf

m = max(abs(x_));
line([-m 0 0 m]*1.2, [-1 -1 1 1]*lam, 'color', [1 1 1]*0.5)

line(x_, z, 'marker', '.', 'linestyle', 'none')
title('Optimality conditions')

ylim([-1 1]*1.5*lam)
xlim([-1 1]*1.2*m)

ylabel('H^T(y-Hx)')
xlabel('x')
box off

set(gca,'fontname','symbol')
set(gca, 'ytick', [-1 0 1]*lam, 'yticklabel', {'-l', '0', 'l'});

printme('optim')

Optimality condition - bound

The signal H^T(y - Hx) should lie within the bounds [-lambda, lambda]

figure(1)
clf
subplot(2,1,1)
plot(z)
xlim([0 N])
box off
line([0 N],[1 1]*lam, 'linestyle', '--')
line([0 N],-[1 1]*lam, 'linestyle', '--')
title('H^T(y - Hx)');
ylim([-1.5 2]*lam)

set(gca,'fontname','symbol')
set(gca, 'ytick', [-1 0 1]*lam, 'yticklabel', {'-l', '0', 'l'});

printme('bound')

Selection of lambda

lambda may be chosen as the minimum of H^T(noise)

figure(2)
clf
subplot(2,1,1)
plot(HT(noise))
xlim([0 N])
box off
title('H^{T}w');
ylim([-1.5 2]*lam)

printme('HTnoise')
max(abs(HT(noise)))
ans =

    1.8953

Least squares solution

lam = 2;                                    % lam : regularization parameter
AAT = A*A';                                 % A*A' : banded matrix
g = (1/lam) * (A'\(B'*y));                  % (1/lam) H'*y
F = lam*AAT + B*B';                         % F : banded matrix
s = g - (B'*(F\(B*g)));

figure(2)
clf
subplot(2,1,1)
plot(s)
box off
ylim([-0.5 1])
title('Least square deconvolution');

printme('deconv_L2')