Simulation_DynamicSeq_BPDN

PURPOSE ^

Simulation for dynamic update with sequential measurements.

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

 Simulation for dynamic update with sequential measurements.
 Compare warm start of GPSR and FPC with the homotopy update of BPDN.
 Need GPSR and FPC packages in your matlab path to run this simulation.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 % Simulation for dynamic update with sequential measurements.
0002 % Compare warm start of GPSR and FPC with the homotopy update of BPDN.
0003 % Need GPSR and FPC packages in your matlab path to run this simulation.
0004 
0005 clear
0006 close all
0007 clear classes;
0008 
0009 sim_runs = 20; % Use a large number of simulation run to compare the performance
0010 tau_table = [.5 .1 .05 .01];
0011 outer_iterations = length(tau_table);
0012 stack_time_solve = zeros(sim_runs*outer_iterations,4);
0013 stack_time_update = zeros(sim_runs*outer_iterations,3);
0014 stack_iter = zeros(sim_runs*outer_iterations,7);
0015 
0016 for out_iter = 1:outer_iterations
0017     for sim_iter = 1:sim_runs
0018         % n is the original signal length
0019         n = 1024;
0020         % m is number of observations to make
0021         m = 512;
0022 
0023         % number of spikes to put down
0024         % n_spikes = floor(.01*n);
0025         ratio_sp = 5;
0026         n_spikes = round(m/ratio_sp);
0027         T = n_spikes;
0028 
0029         % random +/- 1 signal
0030         x = zeros(n,1);
0031         q = randperm(n);
0032         x(q(1:n_spikes)) = sign(randn(n_spikes,1));
0033         % x(q(1:n_spikes)) = randn(n_spikes,1);
0034 
0035         % measurement matrix
0036         G = randn(m,n)/sqrt(m);
0037         % orthonormalize rows
0038         % G = orth(G')';
0039         m_u =1;
0040         A = G(1:end-m_u,:);
0041         b = G(end-m_u+1:end,:);
0042 
0043         hR = @(x) A*x;
0044         hRt = @(x) A'*x;
0045 
0046         % noisy observations
0047         sigma = 0.01; % Add noise depending on the amplitude of y
0048         e = sigma*randn(m-m_u,1);
0049         y = A*x + e;
0050 
0051         % regularization parameter
0052         tau = tau_table(out_iter)*max(abs(A'*y));
0053 
0054         % Homotopy method
0055         tic
0056         [xp, gamma_x, iter_xp, t_homotopy1] = BPDN_homotopy_function(A, y, tau,4*n); %Lasso
0057         t_homotopy1 = toc;
0058         tolA_h = tau*sum(abs(xp))+1/2*(norm(A*xp-y))^2;
0059 
0060         % GPSR
0061         pdg_scale = 1e-12;
0062         tolA = tolA_h*(1+pdg_scale);
0063         stopCri=4;
0064         debias = 0;
0065 
0066         first_tau_factor = 0.8*(max(abs(A'*y))/tau);
0067         steps = 5;
0068         tic;
0069         x_BB_mono = xp;
0070         [x_BB_mono,x_debias_BB_mono,obj_BB_mono,...
0071             times_BB_mono,debias_start_BB_mono,mse]= ...
0072             GPSR_BB(y,hR,tau,...
0073             'Debias',debias,...
0074             'Monotone',1,...
0075             'AT',hRt,...
0076             'Initialization',0,...
0077             'Continuation',1,...
0078             'ContinuationSteps',steps,...
0079             'FirstTauFactor',first_tau_factor,...
0080             'StopCriterion',stopCri,...
0081             'ToleranceA',tolA,...
0082             'Verbose', 0, ...
0083             'ToleranceD',0.001);
0084         t_BB_mono = toc;
0085         
0086         % FPC
0087         opts.gtol = 1e-6;
0088         %opts.f_value_tol = tolA;
0089         opts.PrintOptions = 0;
0090         M = [];
0091         tic;
0092         [x_FPC, Out_FPC] = FPC_AS(n,A,y,tau,M,opts); % call FPC_AS
0093         t_FPC = toc;
0094         
0095         % New measurement
0096         chk_e0 = 1;  % 0 or 1: This selects if we want to take U = (A'A+b'b) (1) or (A'A+e0*b'b) (0)...!
0097         % In simulation U = (A'A+b'b) works better
0098 
0099         % Gram matrix update
0100         AtAgx = A(:,gamma_x)'*A(:,gamma_x)+chk_e0*b(:,gamma_x)'*b(:,gamma_x);
0101         iAtAgx = inv(AtAgx);
0102         pk_old = A'*(A*xp-y);
0103         w = b*x+randn(m_u,1)*sigma;
0104 
0105         % Update using homotopy
0106         tic
0107         [xp_h, gamma_xh, iter_xp_update, t_homotopy_update] =  DynamicSeq_BPDN_function(A, b, AtAgx, iAtAgx, y, w, xp, gamma_x, pk_old, tau, chk_e0, 4*n);
0108         t_homotopy_update = toc;
0109         
0110         
0111         
0112         % Compute using homotopy from scratch
0113         y2 = [y; w];
0114         G = [A; b];
0115         tic
0116         [xp2, gamma_x2, iter_xp2, t_homotopy2] = BPDN_homotopy_function(G, y2, tau,4*n); %Lasso
0117         t_homotopy2 = toc;
0118         
0119         tolA_h = tau*sum(abs(xp2))+1/2*(norm(G*xp2-y2))^2;
0120 
0121         % GPSR
0122         tolA = tolA_h*(1+pdg_scale);
0123         stopCri=4;
0124         debias = 0;
0125 
0126         % PUTTIN CONTINUATION IN HERE MAKES THINGS WORSE
0127         tic;
0128         [x_BB_mono_update,x_debias_BB_mono,obj_BB_mono,...
0129             times_BB_mono,debias_start_BB_mono,mse]= ...
0130             GPSR_BB(y2,G,tau,...
0131             'Debias',debias,...
0132             'Monotone',1,...
0133             'Initialization',x_BB_mono,...
0134             'StopCriterion',stopCri,...
0135             'Verbose', 0, ...
0136             'ToleranceA',tolA,...
0137             'ToleranceD',0.001);
0138         t_BB_mono_update = toc;
0139 
0140         % FPC
0141         opts.x0 = x_FPC;
0142 %       opts.f_value_tol = tolA;
0143         M = [];
0144         tic
0145         [x_FPC_update, Out_FPC_update] = FPC_AS(n,G,y2,tau,M,opts); % call FPC_AS
0146         t_FPC_update = toc;
0147 %         stack_time_solve((out_iter-1)*sim_runs+sim_iter,:) = [t_homotopy1 t_BB_mono];
0148 %         stack_time_update((out_iter-1)*sim_runs+sim_iter,:) = [t_homotopy_update t_BB_mono_update];
0149 %         stack_iter((out_iter-1)*sim_runs+sim_iter,:) = [iter_xp iter_xp_update norm(x-xp) norm(x-xp_h) norm(x-x_BB_mono_update) length(times_BB_mono)];
0150 
0151         stack_time_solve((out_iter-1)*sim_runs+sim_iter,:) = [t_homotopy2 t_homotopy1 t_BB_mono t_FPC];
0152         stack_time_update((out_iter-1)*sim_runs+sim_iter,:) = [t_homotopy_update t_BB_mono_update t_FPC_update];
0153         stack_iter((out_iter-1)*sim_runs+sim_iter,:) = [iter_xp2 iter_xp_update , length(times_BB_mono), (Out_FPC_update.nProdA+Out_FPC_update.nProdAt)/2 norm(x-xp_h)/norm(x) norm(xp_h-x_BB_mono_update)/norm(xp_h) norm(xp_h-x_FPC_update)/norm(xp_h)];
0154         [out_iter sim_iter]
0155     end
0156 end
0157 
0158 % Generate table for performance
0159 disp(['      DynamicSeq         BPDN homotopy           GPSR            FPC_AS']);
0160 
0161 table_stats = zeros(4,8);
0162 table_stats(:,3) = mean(reshape(stack_iter(:,1),sim_runs,outer_iterations));
0163 table_stats(:,1) = mean(reshape(stack_iter(:,2),sim_runs,outer_iterations));
0164 table_stats(:,5) = mean(reshape(stack_iter(:,3),sim_runs,outer_iterations));
0165 table_stats(:,7) = mean(reshape(stack_iter(:,4),sim_runs,outer_iterations));
0166 table_stats(:,4) = mean(reshape(stack_time_solve(:,1),sim_runs,outer_iterations));
0167 table_stats(:,2) = mean(reshape(stack_time_update(:,1),sim_runs,outer_iterations));
0168 table_stats(:,6) = mean(reshape(stack_time_update(:,2),sim_runs,outer_iterations));
0169 table_stats(:,8) = mean(reshape(stack_time_update(:,3),sim_runs,outer_iterations));
0170 table_stats
Generated on Mon 10-Jun-2013 23:03:23 by m2html © 2005