Description of demo

0001 % demo_rwtL1
0002 %
0003 % Solves the following reweighted BPDN problem
0004 % min_x  \sum w_i |x_i| + 1/2*||y-Ax||_2^2
0005 %
0006 % Written by: Salman Asif, Georgia Tech
0007 % Email: sasif@gatech.edu
0008 % Created: August 2012
0009 
0010 clear
0011 close all force
0012  
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 addpath utils/qr/
0022 addpath utils/utils_Wavelet
0023 
0024 disp(['--------------------',datestr(now),'-------------------------'])
0025 
0026 % load RandomStates
0027 %
0028 rseed = 2012;
0029 rand('state',rseed);
0030 randn('state',rseed);
0031 
0032 % simulation parameters
0033 mType = 'randn'; % {'randn','orth','rdct'};
0034 sType = 'randn'; % {'randn','sign','highD', 'blocks','pcwPoly'}
0035 SNR = 40;       % additive Gaussian noise
0036 
0037 % reweighted setup
0038 rwt = 5;        % number of reweighting iterations
0039 
0040 N = 256;   % signal length
0041 M = round(N/2);    % no. of measurements
0042 T = round(M/3);    % sparsity level
0043 
0044 % rank-1 update mode
0045 delx_mode = 'mil'; % mil or qr
0046 fprintf('Iterative reweighting..\n');
0047 str0 = sprintf('mType-%s, sType-%s, SNR = %d, (N,M,T) = %d, %d, %d.', mType, sType, SNR, N, M, T);
0048 disp(str0);
0049  
0050 maxsim = 10;
0051 SIM_stack = cell(maxsim,1);
0052 SIM_memory = cell(maxsim,1);
0053                                 
0054  
0055 for sim = 1:maxsim
0056    
0057     % Generate a random signal
0058     in = []; in.type = sType; in.T = T; in.randgen = 1;
0059     x = genSignal(N,in); 
0060     [val ind] = sort(abs(x),'descend');
0061     ind_pos = ind(val>1e-1);
0062     gamma_orig = ind_pos(1:min(length(ind_pos),M-1));
0063     
0064 
0065     % measurement matrix
0066     in = []; in.type = mType;
0067     A = genAmat(M,N,in);
0068     
0069     % measurements
0070     sigma = sqrt(norm(A*x)^2/10^(SNR/10)/M);    
0071     e = randn(M,1)*sigma;
0072     y = A*x+e;
0073     
0074     %     % orthogonalize rows of A
0075     %     [Q, R] = qr(A',0);
0076     %     A = Q'; y = R' \ y;
0077     
0078     % parameter selection
0079     % tau = sigma*sqrt(log(N));
0080     tau = max(1e-4*max(abs(A'*y)),sigma*sqrt(log(N)));
0081     
0082     maxiter = 4*N;
0083     err_fun = @(z) (norm(x-z)/norm(x))^2;
0084     
0085     %% BPDN using l1homotopy
0086     in = [];   
0087     in.tau = tau;
0088     in.delx_mode = delx_mode;
0089     in.debias = 0;
0090     in.verbose = 0;
0091     in.plots = 0;
0092     in.record = 1;
0093     in.err_fun = err_fun;
0094     tic
0095     out = l1homotopy(A,y,in);
0096     xh = out.x_out;
0097     iter_bpdn = out.iter;
0098     time_bpdn = toc;
0099     gamma_bpdn = out.gamma;
0100     err_bpdn = out.error_table;
0101         
0102     %% Iterative reweighting...
0103     xh_old = xh; W_old = tau; 
0104     error_rwt = []; iter_rwt = []; time_rwt = []; 
0105     for rwt_itr = 1:rwt
0106         
0107         % Update weights
0108         alpha = 1; epsilon = 1;
0109         beta = M*(norm(xh_old,2)/norm(xh_old,1))^2;
0110         W = tau/alpha./(beta*abs(xh_old)+epsilon);
0111         
0112         homotopy_mode = 'dummy';
0113         switch homotopy_mode
0114             case 'dummy'                
0115                 % create a dummy variable...
0116                 % use homotopy on the measurements...
0117                 % in principle, we can start with any xh_old with this formulation
0118                 % and any starting value of W_old (for instance W)...
0119                 W_old = W;
0120                 % yh = A*xh_old;
0121                 % u = -sign(xh_old).*W_old; % out.pk;
0122                 % Atdyu = A'*(yh-y)-u;
0123                 yh = A*xh_old;
0124                 Atr = A'*(A*xh_old-y);
0125                 u =  -W.*sign(xh_old)-Atr;
0126                 pk_old = Atr+u;
0127             otherwise
0128                 % use homotopy on weights
0129                 u = out.pk;
0130                 Atdyu = 0;
0131         end
0132         
0133         in = out;
0134         in.xh_old = xh_old;
0135         in.pk_old = pk_old;
0136         in.u = u;
0137         in.W_old = W_old;
0138         in.W = W;
0139         % in.Atdyu = Atdyu;
0140         in.delx_mode = delx_mode;
0141         in.debias = 0;
0142         in.verbose = 0;
0143         in.plots = 0;
0144         in.record = 1;
0145         in.err_fun = @(z) (norm(x-z)/norm(x))^2;
0146         tic
0147         out = l1homotopy(A,y,in);
0148         xh_old = out.x_out;
0149         W_old = W;
0150         gamma_rwt = out.gamma;
0151         error_rwt = [error_rwt err_fun(xh_old)];
0152         iter_rwt = [iter_rwt out.iter];
0153         time_rwt = [time_rwt toc];
0154     end
0155         
0156 %     %% Check the solution using BPDN directly with l1homotopy
0157 %     in = []; x_old = x;
0158 %     in.W = W;
0159 %     in.delx_mode = delx_mode;
0160 %     in.debias = 0;
0161 %     in.verbose = 0;
0162 %     in.plots = 0;
0163 %     in.record = 1;
0164 %     in.err_fun = err_fun;
0165 %     tic
0166 %     out = l1homotopy(A,y,in);
0167 %     xh = out.x_out;
0168 %     iter_bpdn = out.iter;
0169 %     time_bpdn = toc;
0170 %     gamma_bpdn = out.gamma;
0171 %     err_bpdn = out.error_table;
0172     
0173     %%
0174     SIM_stack{sim} = [sim, tau, ...
0175         norm(x-xh)^2/norm(x)^2, sum(iter_bpdn,2), sum(time_bpdn,2), ...
0176         norm(x-xh_old)^2/norm(x)^2, sum(iter_rwt,2), sum(time_rwt,2)];  
0177     
0178     fprintf('sim %d. tau = %3.4g, (err,iter,time): l1homotopy from scratch-%3.4g,%3.4g,%3.4g; rwt l1homotopy-%3.4g,%3.4g,%3.4g. \n', ...
0179         SIM_stack{sim});
0180 
0181     %% plot recovered signals
0182     figure(1); plot([x xh xh_old]);    
0183     title('Comparison betweeen the new and old homotopy code')
0184 end
0185 mS =  mean(cell2mat(SIM_stack),1);
0186 fprintf('Average results: maxsim %d. tau = %3.4g, (err,iter,time): l1homotopy from scratch-%3.4g,%3.4g,%3.4g; rwt l1homotopy-%3.4g,%3.4g,%3.4g. \n', maxsim, mS(2:end));
demo_rwtL1

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
