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