Simulation_DynamicX_BPDN_Pathological

PURPOSE ^

Simulation of homotopy update for signals with large dynamic range.

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

 Simulation of homotopy update for signals with large dynamic range.
 Comparison with FPC_AS and simple BPDN homotopy
 Need FPC package and accompanying problems to test parts of this file.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 % Simulation of homotopy update for signals with large dynamic range.
0002 % Comparison with FPC_AS and simple BPDN homotopy
0003 % Need FPC package and accompanying problems to test parts of this file.
0004 
0005 clear
0006 % close all
0007 % clear classes;
0008 
0009 % % load fixed random states
0010 % load RandomStates
0011 % rand('state', s_rand);
0012 % randn('state', s_randn);
0013 
0014 sim_runs = 100;
0015 % tau_table = [1e-1 1e-3 1e-5];
0016 outer_iterations = 1;
0017 
0018 stack_time_solve = zeros(sim_runs*outer_iterations,2);
0019 stack_time_update = zeros(sim_runs*outer_iterations,2);
0020 stack_iter = zeros(sim_runs*outer_iterations,5);
0021 
0022 % Test difficult problems, adopted from FPC_AS package.
0023 FPC_problems = 2;
0024 
0025 if FPC_problems == 1
0026     % Following list of problems can be used from FPC_AS package
0027     Problist = {'Ameth6Xmeth2n1024m512k154seed200', 'Ameth6Xmeth6n1024m512k154seed200', ...
0028         'CaltechTest1', 'CaltechTest2', 'CaltechTest3', 'CaltechTest4'};
0029     di = 6; % choose the problem type
0030     load( Problist{di}, 'b','Omega','n','xs');
0031 
0032     m = length(b);
0033     x = xs;
0034     q = (find(abs(xs)>0));
0035     T = length(q);
0036     large_mag = max(abs(x(q)));
0037     small_mag = min(abs(x(q)));
0038     DCT_mat = dct(eye(n));
0039     A = DCT_mat(Omega,:);
0040 else
0041     n = 512;
0042     m = 228;
0043     T = 30;
0044     q = randperm(n);
0045     T1 = 15;
0046     T2 = T-T1;
0047     x = zeros(n,1);
0048     large_mag = 1e6;
0049     small_mag = 1e2;
0050     x(q(1:T1)) = randsrc(T1,1)*large_mag;
0051     x(q(T1+1:T1+T2)) = randsrc(T2,1)*small_mag;
0052 
0053     % measurement matrix
0054     A = randn(m,n)/sqrt(m);
0055     % orthonormalize rows
0056     % A = orth(A')';
0057 
0058 end
0059 
0060 % noisy observations
0061 sigma = 0.0; % how much noise to add should depend on power of Ax
0062 e = randn(m,1)*sigma;
0063 y = A*x + e;
0064 
0065 % regularization parameter
0066 tau = .01*small_mag/large_mag*max(abs(A'*y)); 
0067 %tau = tau_table;
0068 % if sigma>0
0069 %     tau = sigma * sqrt(log(n)*2); % BPDN
0070 %     %tau = max(abs(A'*(A*x-y))); % ideal ???
0071 % end
0072 
0073 % Homotopy method
0074 tic
0075 [xp, gamma_x, iter_xp, t_homotopy1] = BPDN_homotopy_function(A, y, tau,4*n); %BPDN
0076 t_homotopy1 = toc;
0077 tolA_h = tau*sum(abs(xp))+1/2*(norm(A*xp-y))^2;
0078 
0079 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0080 % The relaxation in the tolerance of objective value for FPC
0081 % The performance of FPC is best for the small values of these parameters
0082 % and that increases the iteration count.
0083 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0084 pdg_scale = 1e-3; 
0085 opts.gtol = 1e-8;
0086 
0087 tolA = tolA_h*(1+pdg_scale);
0088 
0089 M = [];
0090 opts.f_value_tol = tolA;
0091 
0092 tic
0093 [x_FPC, Out_FPC] = FPC_AS(n,A,y,tau,M,opts); % call FPC_AS
0094 t_FPC = toc
0095 [t_homotopy1 t_FPC]
0096 
0097 for out_iter = 1:outer_iterations
0098     for sim_iter = 1:sim_runs
0099 
0100         % Model the change in signal
0101         dx_model = 1;
0102         switch dx_model
0103             case  1
0104                 % Add some elements (T_in) or remove some existing elements (T_out)
0105                 dx  = zeros(n,1);
0106                 T_in = randint(1,1,round(T/10))+1;
0107                 T_out = randint(1,1,round(T/10));
0108                 new_large_small = round(rand(T_in,1));
0109                 qin = randperm(n);
0110                 dx(qin(1:T_in)) = randn(T_in,1).*(new_large_small.*large_mag + (1-new_large_small).*small_mag);
0111                 
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)*.01;
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('DIY');
0125         end
0126         xt = x+dx;
0127         yt = y+A*dx;
0128 
0129         pk_old = A'*(A*xp-y);
0130         pk_old(gamma_x) = sign(pk_old(gamma_x))*tau;
0131         AtAgx = A(:,gamma_x)'*A(:,gamma_x);
0132         iAtAgx = inv(AtAgx);
0133 
0134         % Update using homotopy
0135         tic;
0136         [xp_h, gamma_xh, iter_xp_update, t_homotopy_update] = DynamicX_BPDN_function(A, AtAgx, iAtAgx, y, yt, xp, gamma_x, pk_old, tau,2*n);
0137         t_homotopy_update = toc;
0138         
0139         tolA_h = tau*sum(abs(xp_h))+1/2*(norm(A*xp_h-yt))^2;
0140 
0141         % Homotopy method from scratch
0142         tic
0143         [xp2, gamma_x2, iter_xp2, t_homotopy2] = BPDN_homotopy_function(A, yt, tau,2*n); %BPDN
0144         t_homotopy2 = toc;
0145         tolA_h = tau*sum(abs(xp))+1/2*(norm(A*xp-y))^2;
0146         
0147         % % l1_ls
0148         % [x_l1_ls_update,status,history] = l1_ls(A,yt,2*tau,1e-4);
0149         % t_l1_ls_update = history(7,end);
0150         % tolA_ls = history(2,end)/2;
0151 
0152         tolA = tolA_h*(1+pdg_scale);
0153         
0154         opts.x0 = x_FPC;
0155         opts.f_value_tol = tolA;
0156         M = [];
0157 
0158         tic;
0159         [x_FPC_update, Out_FPC_update] = FPC_AS(n,A,yt,tau,M,opts); % call FPC_AS
0160         t_FPC_update = toc;
0161         
0162         stack_time_solve((out_iter-1)*sim_runs+sim_iter,:) = [t_homotopy2 t_FPC];
0163         stack_time_update((out_iter-1)*sim_runs+sim_iter,:) = [t_homotopy_update t_FPC_update];
0164         stack_iter((out_iter-1)*sim_runs+sim_iter,:) = [iter_xp2 iter_xp_update (Out_FPC_update.nProdA+Out_FPC_update.nProdAt)/2 norm(xt-xp_h) norm(xt-x_FPC_update)];
0165         [out_iter sim_iter]
0166     end
0167 end
0168 
0169 % Generate table for performance
0170 disp(['       DynamicX          BPDN homotopy          FPC_AS']);
0171 disp(['   nProdAtA    CPU     nProdAtA     CPU    nProdAtA     CPU'])
0172 table_stats = zeros(length(outer_iterations),6);
0173 table_stats(:,3) = mean(reshape(stack_iter(:,1),sim_runs,outer_iterations));
0174 table_stats(:,1) = mean(reshape(stack_iter(:,2),sim_runs,outer_iterations));
0175 table_stats(:,5) = mean(reshape(stack_iter(:,3),sim_runs,outer_iterations));
0176 table_stats(:,4) = mean(reshape(stack_time_solve(:,1),sim_runs,outer_iterations));
0177 table_stats(:,2) = mean(reshape(stack_time_update(:,1),sim_runs,outer_iterations));
0178 table_stats(:,6) = mean(reshape(stack_time_update(:,2),sim_runs,outer_iterations));
0179 table_stats
Generated on Mon 10-Jun-2013 23:03:23 by m2html © 2005