Description of Simulation_DynamicX

0001 % Simulation for BPDN homotopy update of DynamicX.
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 clear classes; close all; 
0007 
0008 % % load fixed random states
0009 % load RandomStates
0010 % rand('state', s_rand);
0011 % randn('state', s_randn);
0012 
0013 sim_runs = 10; % Run for large number of simulations to compare the performance.
0014 tau_table = [.5, .1, .05, .01];
0015 outer_iterations = length(tau_table);
0016 
0017 stack_time_solve = zeros(sim_runs*outer_iterations,4);
0018 stack_time_update = zeros(sim_runs*outer_iterations,3);
0019 stack_iter = zeros(sim_runs*outer_iterations,9);
0020 
0021 for out_iter = 1:outer_iterations
0022     for sim_iter = 1:sim_runs
0023         % n is the original signal length
0024         n = 1024;
0025 
0026         % k is number of observations to make
0027         k = 512;
0028 
0029         % number of spikes to put down
0030         % n_spikes = floor(.01*n);
0031         ratio_sp = 5;
0032         n_spikes = round(k/ratio_sp);
0033         T = n_spikes;
0034 
0035         % random +/- 1 signal
0036         x = zeros(n,1);
0037         q = randperm(n);
0038         x(q(1:n_spikes)) = sign(randn(n_spikes,1));
0039         % x(q(1:n_spikes)) = randn(n_spikes,1);
0040 
0041         % measurement matrix
0042         A = randn(k,n)/sqrt(k);
0043         % orthonormalize rows
0044         % A = orth(A')';
0045 
0046         % noisy observations
0047         sigma = 0.01; % how much noise to add should depend on power of Ax
0048         e = randn(k,1)*sigma;
0049         y = A*x + e;
0050 
0051         % regularization parameter
0052         tau = tau_table(out_iter)*max(abs(A'*y)); % l1_ls
0053         % if sigma>0
0054         %     tau = sigma * sqrt(log(n)*2); % BPDN
0055         %     %tau = max(abs(A'*(A*x-y))); % ideal ???
0056         % end
0057 
0058         % Homotopy method
0059         tic
0060         [xp, gamma_x, iter_xp, t_homotopy1] = BPDN_homotopy_function(A, y, tau,4*n); %BPDN
0061         t_homotopy1 = toc;
0062         tolA_h = tau*sum(abs(xp))+1/2*(norm(A*xp-y))^2;
0063 
0064         % Adjust tolerance for the objective value, relative to the
0065         % homotopy solution (which we assume to be accurate!).
0066         pdg_scale = 1e-6;
0067 
0068         tolA = tolA_h+pdg_scale;
0069         stopCri=4;
0070         debias = 0;
0071         x_BB_mono = xp;
0072         t_BB_mono = 0;
0073         first_tau_factor = 0.8*(max(abs(A'*y))/tau);
0074         steps = 5;
0075         
0076         tic;
0077         [x_BB_mono,x_debias_BB_mono,obj_BB_mono,...
0078             times_BB_mono,debias_start_BB_mono,mse]= ...
0079             GPSR_BB(y,A,tau,...
0080             'Debias',debias,...
0081             'Monotone',1,...
0082             'Initialization',0,...
0083             'Continuation',1,...
0084             'ContinuationSteps',steps,...
0085             'FirstTauFactor',first_tau_factor,...
0086             'StopCriterion',stopCri,...
0087             'Verbose', 0,...
0088             'True_x',xp,...
0089             'MaxiterA', 1000,...
0090             'ToleranceA',tolA,...
0091             'ToleranceD',0.001);
0092         t_BB_mono = toc;
0093 
0094         % FPC
0095         opts.mxitr = 4*n;
0096         opts.f_value_tol = tolA;
0097         M = [];
0098         tic
0099         [x_FPC, Out_FPC] = FPC_AS(n,A,y,tau,M,opts); % call FPC_AS
0100         t_FPC = toc;
0101         
0102         % Model the change in signal
0103         dx_model = 1;
0104         switch dx_model
0105             case  1
0106                 % Add some elements (T_in) or remove some existing elements (T_out)
0107                 dx  = zeros(n,1);
0108                 T_in = randint(1,1,round(T/20));
0109                 T_out = randint(1,1,round(T/20));
0110                 qin = randperm(n);
0111                 dx(qin(1:T_in)) = (randn(T_in,1));
0112                 qout = randperm(T);
0113                 oldx_indices = q(1:T);
0114                 dx(oldx_indices(qout(1:T_out))) = -x(oldx_indices(qout(1:T_out))).*rand(T_out,1);
0115                 dx(oldx_indices(qout(T_out+1:end))) = randn(T-T_out,1)*.1;
0116             case  2;
0117                 % Perturb the current locations by Gaussian noise
0118                 oldx_indices = q(1:T);
0119                 dx = zeros(n,1);
0120                 dx(oldx_indices) = .1*randn(T,1);
0121                 T_in = 0;
0122                 T_out = 0;
0123             otherwise
0124                 disp('Nooooo');
0125         end
0126         xt = x+dx;
0127         yt = y+A*dx;
0128 
0129         pk_old = A'*(A*xp-y);
0130         AtAgx = A(:,gamma_x)'*A(:,gamma_x);
0131         iAtAgx = inv(AtAgx);
0132 
0133         % Update using homotopy
0134         tic
0135         [xp_h, gamma_xh, iter_xp_update, t_homotopy_update2] = DynamicX_BPDN_function(A, AtAgx, iAtAgx, y, yt, xp, gamma_x, pk_old, tau,4*n);
0136         t_homotopy_update = toc;
0137         tolA_h = tau*sum(abs(xp_h))+1/2*(norm(A*xp_h-yt))^2;
0138 
0139         tic
0140         [xp2, gamma_x2, iter_xp2, t_homotopy2] = BPDN_homotopy_function(A, y, tau,4*n); %BPDN
0141         t_homotopy2 = toc;
0142 
0143         % % l1_ls
0144         % [x_l1_ls_update,status,history] = l1_ls(A,yt,2*tau,1e-4);
0145         % t_l1_ls_update = history(7,end);
0146         % tolA_ls = history(2,end)/2;
0147 
0148         % GPSR
0149         tolA = tolA_h+pdg_scale;
0150         stopCri=4;
0151         debias = 0;
0152         % x_BB_mono = xp;
0153         % t_BB_mono = 0;
0154         % PUTTIN CONTINUATION IN HERE MAKES THINGS WORSE
0155         tic;
0156         [x_BB_mono_update,x_debias_BB_mono,obj_BB_mono,...
0157             times_BB_mono,debias_start_BB_mono,mse]= ...
0158             GPSR_BB(yt,A,tau,...
0159             'Debias',debias,...
0160             'Monotone',1,...
0161             'Initialization',x_BB_mono,...
0162             'StopCriterion',stopCri,...
0163             'ToleranceA',tolA,...
0164             'True_x',xp_h,...
0165             'MaxiterA', 1000,...
0166             'Verbose', 0,...
0167             'ToleranceD',0.001);
0168         t_BB_mono_update = toc;
0169         
0170         % FPC
0171         opts.x0 = x_FPC;
0172         opts.f_value_tol = tolA;
0173         M = [];
0174         tic
0175         [x_FPC_update, Out_FPC_update] = FPC_AS(n,A,yt,tau,M,opts); % call FPC_AS
0176         t_FPC_update = toc;
0177         
0178         stack_time_solve((out_iter-1)*sim_runs+sim_iter,:) = [t_homotopy2 t_homotopy1 t_BB_mono t_FPC];
0179         stack_time_update((out_iter-1)*sim_runs+sim_iter,:) = [t_homotopy_update t_BB_mono_update t_FPC_update];
0180         % New mat files rangTau2 last entry is iteration count for GPSR
0181         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(xt-xp_h)/norm(xt) norm(xp_h-x_BB_mono_update)/norm(xp_h) norm(xp_h-x_FPC_update)/norm(xp_h) norm(xt-x_FPC_update)/norm(xt) Out_FPC_update.itr];
0182         [out_iter sim_iter]        
0183     end
0184 end
0185 
0186 % Generate table for performance
0187 disp(['       DynamicX          BPDN homotopy           GPSR            FPC_AS']);
0188 table_stats = zeros(outer_iterations,8);
0189 table_stats(:,3) = mean(reshape(stack_iter(:,1),sim_runs,outer_iterations));
0190 table_stats(:,1) = mean(reshape(stack_iter(:,2),sim_runs,outer_iterations));
0191 table_stats(:,5) = mean(reshape(stack_iter(:,3),sim_runs,outer_iterations));
0192 table_stats(:,7) = mean(reshape(stack_iter(:,4),sim_runs,outer_iterations));
0193 table_stats(:,4) = mean(reshape(stack_time_solve(:,1),sim_runs,outer_iterations));
0194 table_stats(:,2) = mean(reshape(stack_time_update(:,1),sim_runs,outer_iterations));
0195 table_stats(:,6) = mean(reshape(stack_time_update(:,2),sim_runs,outer_iterations));
0196 table_stats(:,8) = mean(reshape(stack_time_update(:,3),sim_runs,outer_iterations))
Simulation_DynamicX_BPDN

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
