Description of demo

0001 % demo_wtBPDN
0002 %
0003 % Solves the following basis pursuit denoising (BPDN) problem
0004 % min_x  \Sum \w_i |x_i| + 1/2*||y-Ax||_2^2
0005 %
0006 % and dynamically update the weights w_i
0007 %
0008 % Written by: Salman Asif, Georgia Tech
0009 % Email: sasif@gatech.edu
0010 % Created: March 16, 2011
0011 
0012 % Comparison between IRW-H and ARW-H,
0013 % INIT: solution of standard BPDN problem (un-weighted L1)
0014 % also included OMP, oracle-support LS, and oracle-weighted L1
0015 
0016 clear
0017 close all force
0018 
0019 %% parfor setup
0020 % numCores = str2double(getenv('NUMBER_OF_PROCESSORS'));
0021 % mpSize = numCores-1;
0022 % if matlabpool('size') ~= mpSize
0023 %     if matlabpool('size')~=0
0024 %         matlabpool close;
0025 %     else
0026 %         matlabpool('open', mpSize);
0027 %     end
0028 % end
0029 
0030 %% Setup path
0031 mname = mfilename;
0032 mpath = mfilename('fullpath');
0033 mdir = mpath(1:end-length(mname));
0034 cd(mdir);
0035 
0036 addpath ../Pursuits_Homotopy/
0037 addpath ../utils/ 
0038 addpath ../utils/utils_Wavelet/
0039 addpath src/
0040 
0041 disp(['--------------------',datestr(now),'-------------------------'])
0042 
0043 % load RandomStates
0044 %
0045 % rseed = 0;
0046 % rseed = sum(100*clock);
0047 % RandStream.setDefaultStream(RandStream('mt19937ar','seed',rseed));
0048 
0049 rseed = 2012;
0050 rand('state',rseed);
0051 randn('state',rseed);
0052 
0053 % simulation parameters
0054 mType = 'randn'; % {'randn','orth','rdct'};
0055 sType = 'blocks'; % {'randn','sign','highD', 'blocks','pcwPoly'}
0056 SNR = 40;
0057 rwt_mode = 5;   % mode for selecting parameters in iteratirve reweighting
0058 lambda = 0;     % mode for selecting regularization parameter for BPDN
0059 
0060 N = 512;   % signal length
0061 R = 4;
0062 M = round(N/R);    % no. of measurements
0063 T = round(M/5);    % sparsity level
0064 
0065 % reweighted setup
0066 rwt = 5;        % number of reweighting iterations
0067 rwt_adp = 0;    % number of reweighting iterations after adaptive reweighting
0068 
0069 % rank-1 update mode
0070 delx_mode = 'mil'; % mil or qr
0071 
0072 str0 = sprintf('mType-%s, sType-%s, SNR = %d, (N,M,T) = %d, %d, %d, rwt_mode-%d, lambda%3.4g.', mType, sType, SNR, N, M, T, rwt_mode, lambda);
0073 disp(str0);
0074 
0075 
0076 %% setup waitbar
0077 % h = waitbar(0,'1','Name','Adaptive wtBPDN...',...
0078 %     'CreateCancelBtn',...
0079 %     'setappdata(gcbf,''canceling'',1)');
0080 % setappdata(h,'canceling',0)
0081 % reverseStr = '';
0082 
0083 maxsim = 10;
0084 SIM_stack = cell(maxsim,1);
0085 SIM_memory = cell(maxsim,1);
0086                                 
0087 
0088 %% Setting up figures for proper display...
0089 fig(1); 
0090 Position = get(gcf,'Position');
0091 ScreenSize = get(0,'ScreenSize');
0092 Position(1) = 1;
0093 Position(2) = ScreenSize(4)/3; 
0094 Position(3:4) = ScreenSize(3:4)/2;
0095 set(gcf,'Position',Position);
0096 fig(2); 
0097 Position = get(gcf,'Position');
0098 ScreenSize = get(0,'ScreenSize');
0099 Position(1) = ScreenSize(3)/2;
0100 Position(2) = ScreenSize(4)/3; 
0101 Position(3:4) = ScreenSize(3:4)/2;
0102 set(gcf,'Position',Position);
0103 for sim = 1:maxsim
0104     %     % Check for Cancel button press
0105     %     if getappdata(h,'canceling')
0106     %         break
0107     %     end
0108     %     % Report current estimate in the waitbar's message field
0109     %     waitbar(sim/maxsim,h,sprintf('sim %d of %d',sim,maxsim))
0110     
0111     % Generate a random signal
0112     in = []; in.type = sType; in.T = T; in.randgen = 1;
0113     switch sType
0114         case 'blocks'
0115             in.take_fwt = 1;
0116             in.wType = 'haar';
0117         case 'heavisine'
0118             in.take_fwt = 1;
0119             in.wType = 'daub4';
0120         case 'pcwpoly';
0121             in.take_fwt = 1;
0122             in.wType = 'daub8';
0123     end
0124             
0125     x = genSignal(N,in); 
0126     [val ind] = sort(abs(x),'descend');
0127     ind_pos = ind(val>1e-1);
0128     gamma_orig = ind_pos(1:min(length(ind_pos),M-1));
0129     
0130 
0131     % measurement matrix
0132     in = []; in.type = mType;
0133     A = genAmat(M,N,in);
0134     
0135     % measurements
0136     sigma = sqrt(norm(A*x)^2/10^(SNR/10)/M);
0137     %sigma = .05;
0138     e = randn(M,1)*sigma;
0139     y = A*x+e;
0140     
0141     %     % orthogonalize rows of A
0142     %     [Q, R] = qr(A',0);
0143     %     A = Q'; y = R' \ y;
0144     
0145     % parameter selection
0146     if lambda > 0
0147         tau = lambda*max(abs(A'*y));
0148     else
0149         % tau = sigma*sqrt(log(N));
0150         tau = max(1e-4*max(abs(A'*y)),sigma*sqrt(log(N)));
0151     end
0152     
0153     maxiter = 4*N;
0154     
0155     %% Adaptive support & weight selection (ARW-H in the paper)
0156     in = [];
0157     in.A = A; in.y = y; in.x = x;
0158     in.x_init = zeros(N,1); in.max_rwt = rwt_adp;
0159     in.tau = tau;
0160     in.rwt_mode = rwt_mode;
0161     in.delx_mode = delx_mode;
0162     in.debias = 0;
0163     in.verbose = 0;
0164     if in.verbose
0165         fprintf('sType-%s, mType-%s, (N,M,T)-%d,%d,%d ',sType,mType,N,M,T);    
0166     end
0167     out = script_rwtBPDN_adaptive(in);
0168     xh_adp = out.x_out;
0169     iter_adp = out.iter;
0170     time_adp = out.time;
0171     gamma_adp = out.gamma;
0172     W_adp = out.W_new;
0173     err_adp = out.err;
0174     supp_diff_adp = out.supp_diff;
0175     
0176     
0177     %% Iterative reweighting (IRW-H in the paper)
0178     in = [];
0179     in.A = A; in.y = y; in.x = x;
0180     in.x_init = zeros(N,1); in.max_rwt = rwt;
0181     in.tau = tau;
0182     in.rwt_mode = rwt_mode;
0183     in.delx_mode = delx_mode;
0184     out = script_rwtBPDN_iterative(in);
0185     xh_rwt = out.x_out;
0186     xh_rwt_init = out.x_init;
0187     iter_rwt = out.iter;
0188     time_rwt = out.time;
0189     gamma_rwt = out.gamma;
0190     err_rwt = out.err;
0191     supp_diff_rwt = out.supp_diff;
0192     W_rwt = out.W_new;
0193     
0194     %% Oracle rwt BPDN
0195     alpha = 5; beta = 10;
0196     W_new = ones(N,1); W_new(gamma_orig) = 1/alpha./(beta*abs(x(gamma_orig)));    
0197     % epsilon = 0.1; W_new = 1/alpha./(beta*(abs(x))+epsilon);
0198     
0199     % To check the accuracy of rwtBPDN_adaptive set
0200     % W_new = W_adp/tau;
0201     
0202     in = [];
0203     in.tau = tau;
0204     in.maxiter = maxiter;
0205     in.x_orig = x;
0206     in.record = 1;
0207     in.delx_mode = delx_mode;
0208     %     in.Te = rwt_itr;
0209     AW = A*diag(1./W_new);
0210     out_new = BPDN_homotopy_function(AW, y, in); %BPDN
0211     xh_orac = out_new.x_out.*(1./W_new);
0212     gamma_orac = out_new.gamma;
0213     iter_orac = out_new.iter;
0214     time_orac = out_new.time;
0215     
0216     %% oracle LS
0217     x_LS = zeros(N,1);
0218     x_LS(gamma_orig) = A(:,gamma_orig)\y;
0219                                         
0220     %% OMP
0221     in = [];
0222     in.Te = round(1.2*T);
0223     out = OMP_function(y,A,in);
0224     x_omp = out.x_out;
0225     iter_omp = out.iter;
0226     gamma_omp = out.gamma;
0227     
0228     
0229     % fprintf('sim %d -- iter: %d; %d; %d; %d. Error: %3.4g; %3.4g; %3.4g; %3.4g; %3.4g. Diff-supp = %d; %d; %d; %d.\n', sim, iter_orac, iter_adp, iter_rwt, iter_omp, norm(x-x_LS)/norm(x), norm(x-xh_orac)/norm(x), norm(x-xh_adp)/norm(x), norm(x-xh_rwt)/norm(x), norm(x-x_omp)/norm(x), length(setdiff(gamma_orac,gamma_orig)), length(setdiff(gamma_adp,gamma_orig)),length(setdiff(gamma_rwt,gamma_orig)),length(setdiff(gamma_omp,gamma_orig)));
0230     % msg = sprintf('sim %d -- iter: %d; %d; %d; %d. Error: %3.4g; %3.4g; %3.4g; %3.4g; %3.4g. Diff-supp = %d; %d; %d; %d.\n', sim, iter_orac, iter_adp, iter_rwt, iter_omp, norm(x-x_LS)/norm(x), norm(x-xh_orac)/norm(x), norm(x-xh_adp)/norm(x), norm(x-xh_rwt)/norm(x), norm(x-x_omp)/norm(x), length(setdiff(gamma_orac,gamma_orig)), length(setdiff(gamma_adp,gamma_orig)),length(setdiff(gamma_rwt,gamma_orig)),length(setdiff(gamma_omp,gamma_orig)));
0231     % fprintf([reverseStr, msg]);
0232     % reverseStr = repmat(sprintf('\b'), 1, length(msg));
0233     
0234     % SIM_stack(sim,:) = [sim, iter_orac, iter_adp(1), sum(iter_adp,2), iter_rwt(1), sum(iter_rwt,2), iter_omp, norm(x-x_LS)/norm(x),  norm(x-xh_orac)/norm(x), err_adp(1), norm(x-xh_adp)/norm(x), err_rwt(1), norm(x-xh_rwt)/norm(x), norm(x-x_omp)/norm(x), length(setdiff(gamma_orac,gamma_orig)), length(setdiff(gamma_adp,gamma_orig)),length(setdiff(gamma_rwt,gamma_orig)),length(setdiff(gamma_omp,gamma_orig))];
0235     % fprintf('sim %d -- iter: %d; %d; %d; %d. Error: %3.4g; %3.4g; %3.4g; %3.4g; %3.4g. Diff-supp = %d; %d; %d; %d.\n', sim, iter_orac, iter_adp, iter_rwt, iter_omp, norm(x-x_LS)/norm(x), norm(x-xh_orac)/norm(x), norm(x-xh_adp)/norm(x), norm(x-xh_rwt)/norm(x), norm(x-x_omp)/norm(x), length(setdiff(gamma_orac,gamma_orig)), length(setdiff(gamma_adp,gamma_orig)),length(setdiff(gamma_rwt,gamma_orig)),length(setdiff(gamma_omp,gamma_orig)));
0236     
0237     SIM_stack{sim} = [sim, tau, norm(x-x_LS)^2/norm(x)^2, norm(x-xh_orac)^2/norm(x)^2, iter_orac, ...
0238         norm(x-xh_adp)^2/norm(x)^2, sum(iter_adp,2), sum(time_adp,2), ...
0239         norm(x-xh_rwt)^2/norm(x)^2, sum(iter_rwt,2), sum(time_rwt,2), ...
0240         norm(x-x_omp)^2/norm(x)^2, iter_omp];  
0241     
0242     % adaptive rwt
0243     exp = 1;
0244     SIM_memory{sim}{exp,1} = 'ARW-H';
0245     SIM_memory{sim}{exp,2} = iter_adp;
0246     SIM_memory{sim}{exp,3} = err_adp;
0247     SIM_memory{sim}{exp,4} = supp_diff_adp;
0248     SIM_memory{sim}{exp,5} = time_adp;
0249     % iterative rwt
0250     exp = exp+1;
0251     SIM_memory{sim}{exp,1} = 'IRW-H';
0252     SIM_memory{sim}{exp,2} = iter_rwt;
0253     SIM_memory{sim}{exp,3} = err_rwt;
0254     SIM_memory{sim}{exp,4} = supp_diff_rwt;
0255     SIM_memory{sim}{exp,5} = time_rwt;
0256 
0257     mS =  SIM_stack{sim};
0258     fprintf('sim %d. tau = %3.4g, oracLS=%3.4g. (err,iter,time): oracWT-%3.4g,%3.4g; ARW-H-%3.4g,%3.4g,%3.4g; IRW-H-%3.4g,%3.4g,%3.4g; OMP-%3.4g,%3.4g. \n', sim, mS(2:end));
0259 
0260     %% plot recovered signals
0261     fig(1); clf
0262     subplot(221); plot([x xh_adp xh_rwt xh_rwt_init]); title('ORIG .. ARW-H .. IRW-H .. INIT estimates');
0263     subplot(222); plot([W_adp W_rwt]); title('Weights: ARW-H .. IRW-H'); 
0264     subplot(223); plot(xh_adp,x,'bd'); title('Original vs ARW-H estimate'); axis square
0265     subplot(224); plot(xh_rwt_init, x,'k.'); 
0266     hold on; plot(xh_rwt, x,'rs');    title('Original vs INIT and IRW-H estimate'); axis square;
0267     fig(2); clf
0268     hold on; plot(xh_rwt_init, x,'k.');    
0269     hold on; plot(xh_rwt, x,'rs');        
0270     plot(xh_adp,x,'bd');
0271     XLim = get(gca,'XLim');
0272     XLim(1) = sign(XLim(1))*max(abs(XLim(:)));
0273     XLim(2) = sign(XLim(2))*max(abs(XLim(:)));
0274     set(gca,'XLim',XLim);
0275     set(gca,'YLim',XLim);    
0276     axis square; shg    
0277     title('Comparison betweeen INIT .. IRW-H .. ARW-H and the original signal')
0278     legend('INIT','IRW-H','ARW-H','Location','NorthWest');
0279     
0280     if sim == 1
0281         fprintf('\n Hit any key to continue... resize figures (if you want) \n')    
0282         pause;   
0283     end
0284 end
0285 
0286 mS =  mean(cell2mat(SIM_stack),1);
0287 fprintf('Average results: maxsim %d. tau = %3.4g, oracLS=%3.4g. (err,iter,time): oracWT-%3.4g,%3.4g; ARW-H-%3.4g,%3.4g,%3.4g; IRW-H-%3.4g,%3.4g,%3.4g; OMP-%3.4g,%3.4g. \n', maxsim, mS(2:end));
0288 
0289 %% Plot results...
0290 EXP = {'ARW-H','IRW-H'};
0291 exp_plot = [1 2];
0292 fig(101);
0293 for exp_no = 1:length(exp_plot)
0294     exp = exp_plot(exp_no);
0295     sim_err = [];
0296     sim_itr = [];
0297     sim_time = [];
0298     sim_supp_diff = [];
0299     for ii = 1:maxsim
0300         sim_err = [sim_err; [SIM_memory{ii}{exp,3}]];
0301         sim_itr = [sim_itr; [SIM_memory{ii}{exp,2}]];
0302         sim_time = [sim_time; [SIM_memory{ii}{exp,5}]];
0303         sim_supp_diff = [sim_supp_diff; [SIM_memory{ii}{exp,4}]];
0304     end
0305     subplot(2,length(exp_plot),exp_no); plot(sim_err); title(sprintf('rel. error for %s',char(EXP{exp})));
0306     subplot(2,length(exp_plot),length(exp_plot)+exp_no); plot(sim_itr); title(sprintf('iter. count for %s',char(EXP{exp})));
0307     % subplot(3,length(exp_plot),2*length(exp_plot)+exp_no); plot(sim_supp_diff); title(sprintf('supp diff for %s',char(EXP{exp})));
0308 end
0309 
0310 % delete(h)
0311 % DELETE the waitbar; don't try to CLOSE it.
0312 % close all force % After Ctrl+C
demo_wtBPDN

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
