Example 1: Polynomial approximation + total variation filtering
Example: Simultaneous least-square polynomial approximation and total variation filtering.
Ivan Selesnick, Polytechnic Institute of NYU December 2011 Email: selesi@poly.edu Reference: Polynomial Smoothing of Time Series with Additive Step Discontinuities I. W. Selesnick, S. Arnold, and V. R. Dantham
Contents
Start
clear close all printme = @(filename) print('-dpdf', sprintf('figures/Example1_%s', filename));
Create test signals
N = 100; % N : length of data n = (1:N)'; s1 = n < 0.3*N; % s1 : step function s2 = 2-2*((n-N/2)/(N/2)).^2; s2 = s2 + n/N; % s2 : polynomial function randn('state',0); % Initialize randn so that example can be exactly reproduced noise = randn(N,1); sigma = 0.26; noise = noise / sqrt((1/N)*sum(noise.^2)) * sigma; % Gaussian (normal) noise
Polynomial approximation
Show polynomial approximation of noisy polynomial data
a = polyfit(n,s2+noise,2); s2pf = polyval(a, n); figure(1) clf plot(n, s2 + noise, '.k', n, s2pf); title('Polynomial approximation to noisy polynomial data') printme('polyfit')
TV filtering of noisy step data
lambda = 1.6; % lambda : regularization parameter for TV filter Nit = 100; % Nit : number of iterations for TV filter algorithm [x_tv, cost] = TVfilt(s1+noise, lambda, Nit, 1.5, 10); figure(1) plot(cost) title('TV filter - Cost function history') xlabel('Iteration')
Output of TV filter
figure(1) plot(n, s1 + noise, '.k', n, x_tv); title('(a) Total variation filtering of noisy step data'); printme('TV_filter')
Least-square filtering of noisy step data
lambda = 10; % lambda : regularzation parameter for least squares filter e = ones(N-1, 1); D = spdiags([e -e], [0 1], N-1, N); % D : first-order different matrix (sparse) H = speye(N) + lambda*(D'*D); % H : I + lambda D'*D (sparse) x_l2 = H \ (s1+noise); figure(1) plot(n, s1 + noise, '.k', n, x_l2); title('(b) Least-squares filtering of noisy step data'); printme('L2_filtered')
Create noisy signal
Create polynomial signal with additive step discontinuity
s = s1 + s2; % s : polynomial signal with additive step discontinuity y = s + noise; % y : noisy signal figure(1) plot(n, y, '.k') title('(a) Noisy data'); printme('data')
PATV: Polynomial approximation + total variation filter
Example
% parameters d = 2; % d : degree of approximation polynomial lambda = 3; % lambda : regularization parameter Nit = 100; % Nit : number of iterations mu0 = 20; % mu0 : ADMM parameter mu1 = 0.2; % mu1 : ADMM parameter [x, p, cost] = patv(y, d, lambda, Nit, mu0, mu1); % display cost function history figure(1) plot(cost) title('PATV algorithm - Cost function history'); xlabel('Iteration') printme('PATV_convergence')
Display calculated TV component
figure(1) plot(n, x, 'k') title('(b) Calculated TV component (PATV)'); printme('PATV_TV_component')
Display TV-compensated data
figure(2) clf plot(n, y - x, '.k') title('(c) TV-compensated data (PATV)'); printme('PATV_TV_compensated_data')
Display estimated signal
figure(1) plot(n, y,'.k', n, x+p, 'black') title('(d) Estimated signal (PATV)'); legend('Data','Estimated signal', 'Location','south') printme('PATV')
Display first-order difference of TV component. Note: there are 3 non-zero values
figure(2) stem(abs(diff(x)),'marker','none') title('|diff(x)|');
Lambda too small
When lambda is too small, the noise is not fully reduced
lambda = 0.2; [x, p, cost] = patv(y, d, lambda, Nit, mu0, mu1); figure(1) plot(n, y,'.k', n, x+p, 'black') title('(a) Estimated signal (PATV with \lambda too small)'); legend('Data','Estimated signal', 'Location','south') printme('PATV_lambda_small')
Lambda too large
When lambda is too large, the step discontinuity is under-estimated
lambda = 7; [x, p, cost] = patv(y, d, lambda, Nit, mu0, mu1); figure(2) plot(n, y,'.k', n, x+p, 'black') legend('Data','Estimated signal', 'Location','south') title('(b) Estimated signal (PATV with \lambda too large)'); printme('PATV_lambda_large')
C-PATV: Constrained formulation
r = sqrt(N) * sqrt(sum((1/N)*(noise.^2))); % r : constraint parameter Nit = 50; % Nit : number of iterations mu0 = 3.5; % mu0 : ADMM parameter mu1 = 0.5; % mu1 : ADMM parameter [x_constr, p_constr, cost, constr] = cpatv(y, d, r, Nit, mu0, mu1); % check constraint: e = y - x_constr - p_constr; sqrt((1/N)*sum ( e.^2)) % This value should be r/sqrt(N). % Cost function versus constraint function histories figure(1) plot(constr, cost,'.-'); xlabel('Constraint') ylabel('Cost') title('Cost vs. Constraint (C-PATV)') ax3 = axis; line([1 1]*r, ax3(3:4), 'linestyle', '--') axis(ax3) printme('CPATV_convergence')
ans = 0.2600
Display estimated signal
figure(1) plot(n, y,'.k', n, x_constr+p_constr, 'black') % , 'MarkerSize',MS) title('Estimated signal (C-PATV)'); legend('Data','Estimated signal', 'Location','south') printme('CPATV')
Enhanced PATV
Lp quasi-norm minimization with p < 1 leads to fewer extraneous steps in the estimated signal
lambda = 3; % lambda : regularization parameter Nit = 100; % Nit : number of iterations mu0 = 20; % mu0 : ADMM parameter mu1 = 0.2; % mu1 : ADMM parameter p = 0.7; % p : power (Lp quasi-norm) E = 0.02; % E : small number [x, p, cost] = patv_Lp(y, d, lambda, p, E, Nit, mu0, mu1); figure(1), clf plot(n, y,'.k', n, x+p, 'black') title('Estimated signal (Enhanced PATV)'); legend('Data','Estimated signal', 'Location','south') printme('enhanced_PATV')
It can be seen that enhanced PATV leads to fewer extraneous steps in the calculated TV component.
figure(1) plot(x) hold off title('Calculated TV component (Enhanced PATV)')
figure(2) stem(abs(diff(x)), 'marker','none') hold off title('|diff(x1)|')