demo_posBPDN

 demo_posBPDN

 Solves the following basis pursuit denoising (BPDN) problem
 min_x  \tau \|x\|_1 + 1/2*||y-Ax||_2^2
 with a positivity constraint... i.e., x >= 0

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

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