demo_dynamicX

 demo_dynamicX

 Solves the following BPDN problem
 min_x  tau\|x\|_1 + 1/2*||y-Ax||_2^2
 
 and updates the solution using updated measurements 
 for a slighly modified x

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

0001 % demo_dynamicX
0002 %
0003 % Solves the following BPDN problem
0004 % min_x  tau\|x\|_1 + 1/2*||y-Ax||_2^2
0005 %
0006 % and updates the solution using updated measurements
0007 % for a slighly modified x
0008 %
0009 % Written by: Salman Asif, Georgia Tech
0010 % Email: sasif@gatech.edu
0011 % Created: February 2013
0012 
0013 clear
0014 close all force
0015   
0016 %% Setup path
0017 mname = mfilename;
0018 mpath = mfilename('fullpath');
0019 mdir = mpath(1:end-length(mname));
0020 cd(mdir);
0021 
0022 addpath utils/
0023 addpath utils/qr/
0024 
0025 disp(['--------------------',datestr(now),'-------------------------'])
0026 
0027 % load RandomStates
0028 %
0029 rseed = 2012;
0030 rand('state',rseed);
0031 randn('state',rseed);
0032 
0033 % simulation parameters
0034 mType = 'randn'; % {'randn','orth','rdct'};
0035 sType = 'randn'; % {'randn','sign','highD', 'blocks','pcwPoly'}
0036 SNR = inf;       % additive Gaussian noise
0037 
0038 % reweighted setup
0039 rwt = 5;        % number of reweighting iterations
0040 
0041 N = 512;   % signal length
0042 M = round(N/2);    % no. of measurements
0043 T = round(M/5);    % sparsity level
0044 Tn = round(T/20);  % innovations in the signal
0045 
0046 % rank-1 update mode
0047 delx_mode = 'mil'; % mil or qr
0048 fprintf('Update time-varying signal..\n');
0049 str0 = sprintf('mType-%s, sType-%s, SNR = %d, (N,M,T) = %d, %d, %d; Tn = %d.', mType, sType, SNR, N, M, T,Tn);
0050 disp(str0);
0051  
0052 maxsim = 10;
0053 SIM_stack = cell(maxsim,1);
0054 SIM_memory = cell(maxsim,1);
0055 
0056 for sim = 1:maxsim
0057     
0058     % Generate a random signal
0059     in = []; in.type = sType; in.T = T; in.randgen = 1;
0060     x = genSignal(N,in);
0061     [val ind] = sort(abs(x),'descend');
0062     ind_pos = ind(val>1e-1);
0063     gamma_orig = ind_pos(1:min(length(ind_pos),M-1));
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     %% Update the solution for a time-varying signal (DynamicX)
0103     xh_old = xh; tau_old = tau;
0104     x_old = x; y_old = y; 
0105     
0106     % Update the sparse signal
0107     in = []; in.type = sType; in.T = Tn; in.randgen = 1;
0108     dx = genSignal(N,in);
0109     x = x_old + dx;
0110     e = randn(M,1)*sigma;
0111     y = A*x+e;
0112     
0113     % parameter selection
0114     % tau = sigma*sqrt(log(N));
0115     tau = max(1e-4*max(abs(A'*y)),sigma*sqrt(log(N)));
0116 
0117     
0118     % create a dummy variable...
0119     % use homotopy on the measurements...
0120     % in principle, we can start with any xh_old with this formulation
0121     % and any starting value of tau or W...
0122     W = tau;
0123     Atr = A'*(A*xh_old-y);
0124     u =  -W.*sign(xh_old)-Atr;
0125     pk_old = Atr+u;
0126     
0127     in = out;
0128     in.xh_old = xh_old;
0129     in.pk_old = pk_old;
0130     in.u = u;
0131     in.W = W;
0132     in.delx_mode = delx_mode;
0133     in.debias = 0;
0134     in.verbose = 0;
0135     in.plots = 0;
0136     in.record = 1;
0137     in.err_fun = @(z) (norm(x-z)/norm(x))^2;
0138     tic
0139     out = l1homotopy(A,y,in);
0140     xh_dynX = out.x_out;
0141     gamma_dynX = out.gamma;
0142     iter_dynX = out.iter;
0143     time_dynX = toc;
0144          
0145     
0146     SIM_stack{sim} = [sim, tau, ...
0147         norm(x_old-xh)^2/norm(x)^2, sum(iter_bpdn,2), sum(time_bpdn,2), ...
0148         norm(x-xh_dynX)^2/norm(x)^2, sum(iter_dynX,2), sum(time_dynX,2)];  
0149     
0150     fprintf('sim %d. tau = %3.4g, (err,iter,time): l1homotopy from scratch-%3.4g,%3.4g,%3.4g; dynX l1homotopy-%3.4g,%3.4g,%3.4g. \n', ...
0151         SIM_stack{sim});
0152    
0153     %% plot recovered signals
0154     figure(1); plot([x xh xh_dynX]);
0155     title('Comparison betweeen the new and old homotopy code')
0156 end
0157 mS =  mean(cell2mat(SIM_stack),1);
0158 fprintf('Average results: maxsim %d. tau = %3.4g, (err,iter,time): l1homotopy from scratch-%3.4g,%3.4g,%3.4g; dynX l1homotopy-%3.4g,%3.4g,%3.4g. \n', maxsim, mS(2:end));