demo_BPDN

 demo_BPDN

 Solves the following basis pursuit denoising (BPDN) problem: 
 \min_x  \tau \|x\|_1 + \frac{1}{2}\|y-Ax\|_2^2 

 Written by: Salman Asif, Georgia Tech
 Email: sasif@gatech.edu
 Created: February 2013

0001 % demo_BPDN
0002 %
0003 % Solves the following basis pursuit denoising (BPDN) problem:
0004 % \min_x  \tau \|x\|_1 + \frac{1}{2}\|y-Ax\|_2^2
0005 %
0006 % Written by: Salman Asif, Georgia Tech
0007 % Email: sasif@gatech.edu
0008 % Created: February 2013
0009 
0010 clear
0011 close all force
0012 
0013 %% Setup path
0014 mname = mfilename;
0015 mpath = mfilename('fullpath');
0016 mdir = mpath(1:end-length(mname));
0017 cd(mdir);
0018 
0019 addpath utils/  
0020 
0021 disp(['--------------------',datestr(now),'-------------------------'])
0022 
0023 % load RandomStates
0024 %
0025 rseed = 2012;
0026 rand('state',rseed);
0027 randn('state',rseed);
0028 
0029 % simulation parameters
0030 mType = 'randn'; % {'randn','orth','rdct'};
0031 sType = 'randn'; % {'randn','sign','highD', 'blocks','pcwPoly'}
0032 SNR = 40;       % additive Gaussian noise
0033 
0034 N = 512;   % signal length
0035 M = round(N/2);    % no. of measurements
0036 T = round(M/3);    % sparsity level
0037 
0038 % rank-1 update mode
0039 delx_mode = 'mil'; % mil or qr
0040 lambda = 0;
0041 fprintf('Standard LASSO/BPDN homotopy..\n');
0042 str0 = sprintf('mType-%s, sType-%s, SNR = %d, (N,M,T) = %d, %d, %d.', mType, sType, SNR, N, M, T);
0043 disp(str0);
0044  
0045 maxsim = 10;
0046 SIM_stack = cell(maxsim,1);
0047 
0048 for sim = 1:maxsim
0049  
0050     % Generate a random signal
0051     in = []; in.type = sType; in.T = T; in.randgen = 1;
0052     x = genSignal(N,in); 
0053     [val ind] = sort(abs(x),'descend');
0054     ind_pos = ind(val>1e-1);
0055     gamma_orig = ind_pos(1:min(length(ind_pos),M-1));
0056      
0057     % measurement matrix
0058     in = []; in.type = mType;
0059     A = genAmat(M,N,in);
0060     
0061     % measurements
0062     sigma = sqrt(norm(A*x)^2/10^(SNR/10)/M);    
0063     e = randn(M,1)*sigma;
0064     y = A*x+e;
0065     
0066     %     % orthogonalize rows of A
0067     %     [Q, R] = qr(A',0);
0068     %     A = Q'; y = R' \ y;
0069     
0070     % parameter selection
0071     % tau = sigma*sqrt(log(N));
0072     tau = max(1e-4*max(abs(A'*y)),sigma*sqrt(log(N)));
0073     
0074     err_fun = @(z) (norm(x-z)/norm(x))^2;
0075     maxiter = 4*N;
0076     
0077     %% BPDN using l1homotopy
0078     in = [];   
0079     in.tau = tau;
0080     in.delx_mode = delx_mode;
0081     in.debias = 0;
0082     in.verbose = 0;
0083     in.plots = 0;
0084     in.record = 1;
0085     in.err_fun = err_fun;
0086     out = l1homotopy(A,y,in);
0087     xh = out.x_out;
0088     iter_bpdn = out.iter;
0089     time_bpdn = out.time;
0090     gamma_bpdn = out.gamma;
0091     err_bpdn = out.error_table;
0092         
0093     %% Old code
0094     in = [];
0095     in.tau = tau;
0096     in.delx_mode = delx_mode;
0097     in.debias = 0;
0098     in.verbose = 0;
0099     in.plots = 0;
0100     in.record = 1;
0101     in.err_fun = err_fun;
0102     out = l1homotopy_v1(A,y,in);
0103     xh_old = out.x_out;
0104     iter_old = out.iter;
0105     time_old = out.time;
0106     gamma_old = out.gamma;
0107     err_old = out.error_table;
0108     
0109     SIM_stack{sim} = [sim, tau, ...
0110         norm(x-xh)^2/norm(x)^2, sum(iter_bpdn,2), sum(time_bpdn,2), ...
0111         norm(x-xh_old)^2/norm(x)^2, sum(iter_old,2), sum(time_old,2)];  
0112     
0113     fprintf('sim %d. tau = %3.4g, (err,iter,time): l1homotopy-%3.4g,%3.4g,%3.4g; old homotopy-%3.4g,%3.4g,%3.4g. \n', ...
0114         SIM_stack{sim});
0115 
0116     %% plot recovered signals
0117     fig(1); plot([x xh xh_old]);    
0118     title('Comparison betweeen the new and old homotopy code')
0119 end
0120 mS =  mean(cell2mat(SIM_stack),1);
0121 fprintf('Average results: maxsim %d. tau = %3.4g, (err,iter,time): l1homotopy-%3.4g,%3.4g,%3.4g; old homotopy-%3.4g,%3.4g,%3.4g. \n', maxsim, mS(2:end));