demo_dynamicRWT

 demo_dynamicRWT

 Solves the following dynamic BPDN problem
 min_x  \|Wx\|_1 + 1/2*||y-Ax||_2^2  

 which updates W and the solution as the signal changes 

   
 Applications:
       tracking where y, A, and x, are updated..
       predict an estimate of the solution and
       update weights according to the predicted solution

 We can also add a difference operator ||x-xh||_2^2 in the reconstruction

 Written by: Salman Asif, Georgia Tech
 Email: sasif@gatech.edu
 Created: August 2012

0001 % demo_dynamicRWT
0002 %
0003 % Solves the following dynamic BPDN problem
0004 % min_x  \|Wx\|_1 + 1/2*||y-Ax||_2^2
0005 %
0006 % which updates W and the solution as the signal changes
0007 %
0008 %
0009 % Applications:
0010 %       tracking where y, A, and x, are updated..
0011 %       predict an estimate of the solution and
0012 %       update weights according to the predicted solution
0013 %
0014 % We can also add a difference operator ||x-xh||_2^2 in the reconstruction
0015 %
0016 % Written by: Salman Asif, Georgia Tech
0017 % Email: sasif@gatech.edu
0018 % Created: August 2012
0019 
0020 clear
0021 close all force
0022  
0023 %% Setup path
0024 mname = mfilename;
0025 mpath = mfilename('fullpath');
0026 mdir = mpath(1:end-length(mname));
0027 cd(mdir);
0028 
0029 addpath utils/ 
0030 addpath utils/utils_wavelet/
0031 addpath solvers/
0032 % addpath WeightedBPDN/src/
0033 
0034 disp(['--------------------',datestr(now),'-------------------------'])
0035 
0036 % load RandomStates
0037 %
0038 rseed = 2012;
0039 rand('state',rseed);
0040 randn('state',rseed);
0041 
0042 % simulation parameters
0043 mType = 'randn'; % {'randn','orth','rdct'};
0044 sType = 'pcwpoly'; % {'randn','sign','highD', 'blocks','pcwPoly'}
0045 SNR = 40;       % additive Gaussian noise
0046 
0047 N = 256;   % signal length
0048 M = round(N/4);    % no. of measurements
0049 T = round(M/3);    % sparsity level
0050 
0051 % signal dynamics
0052 cshift = -2;
0053 rshift_h = @(z) (rand-0.5)/20;
0054 
0055 % reg. parameter for temporal difference..
0056 lambda_diff = sqrt(1);
0057 
0058 % rank-1 update mode
0059 delx_mode = 'mil'; % mil or qr
0060 lambda = 0;
0061 fprintf('Tracking a dynamical signal using reweighting..\n');
0062 str0 = sprintf('mType-%s, sType-%s, SNR = %d, (N,M,T) = %d, %d, %d.', mType, sType, SNR, N, M, T);
0063 disp(str0);
0064  
0065 maxsim = 200;
0066 SIM_stack = cell(maxsim,1);
0067 SIM_memory = cell(maxsim,1);
0068 
0069 x_vec = zeros(N*maxsim,1);
0070 xh_vec = zeros(N*maxsim,1);
0071 
0072 %% Setup signal and system
0073 
0074 % Generate a random signal
0075 in = []; in.type = sType; in.T = T; in.randgen = 0; in.wType = 'daub79'; in.take_fwt = 1;
0076 [x sig wave_struct] = genSignal(N,in);
0077 [val ind] = sort(abs(x),'descend');
0078 ind_pos = ind(val>1e-1);
0079 gamma_orig = ind_pos(1:min(length(ind_pos),M-1));
0080 
0081 % Setup dynamical model
0082 x_orig = x; % Initial signal
0083 
0084 % At every time instance, add to the original/previous signal
0085 % an integer circular shift that is known
0086 % a random drift that is unknown
0087 if isempty(wave_struct)
0088     % cshfits
0089     F_h = @(x,cshift,rshift) interp1(1:N,circshift(x,cshift),[1:N]+rshift,'linear','extrap')';
0090     W_h = @(z) z; iW_h = @(z) z; 
0091 else
0092       
0093     W_h = wave_struct.W_h;
0094     iW_h = wave_struct.iW_h;    
0095     
0096     F_h = @(x,cshift,rshift) W_h(interp1(1:N,circshift(iW_h(x),cshift),[1:N]+rshift,'linear','extrap')');
0097 end
0098 
0099 % measurement matrix
0100 in = []; in.type = mType;
0101 A = genAmat(M,N,in);
0102 
0103 % measurements
0104 sigma = sqrt(norm(A*x)^2/10^(SNR/10)/M);
0105 e = randn(M,1)*sigma;
0106 y = A*x+e;
0107 
0108 %     % orthogonalize rows of A
0109 %     [Q, R] = qr(A',0);
0110 %     A = Q'; y = R' \ y;
0111 
0112 % parameter selection
0113 % tau = sigma*sqrt(log(N));
0114 tau = max(1e-4*max(abs(A'*y)),sigma*sqrt(log(N)));
0115 
0116 maxiter = 4*N;
0117 err_fun = @(z) (norm(x-z)/norm(x))^2;
0118 
0119 %% initialize...!
0120 %
0121 % Assume a previous signal is known..
0122 % [val_sort ind_sort] = sort(abs(x),'descend');
0123 % xh_old = x;
0124 % xh_old(ind_sort(M+1:end)) = 0;
0125 %
0126 % lambda_diff = sqrt(1e-0);
0127 %
0128 % A = [A; lambda_diff*eye(N)];
0129 % y = [A(1:M,:)*x+e; lambda_diff*xh_old];
0130 
0131 
0132 % rwt L1 using l1homotopy
0133 in = [];
0134 in.tau = tau; W = tau;
0135 in.delx_mode = delx_mode;
0136 in.debias = 0;
0137 in.verbose = 0;
0138 in.plots = 0;
0139 in.record = 1;
0140 in.err_fun = err_fun;
0141 tic 
0142 for wt_itr = 1:5 
0143     
0144     W_old = W;
0145     
0146     out = l1homotopy(A,y,in);
0147     xh = out.x_out;
0148     iter_bpdn = out.iter;
0149     time_bpdn = toc;
0150     gamma_bpdn = out.gamma;
0151     err_bpdn = out.error_table;
0152     
0153     % Update weights
0154     xh_old = xh;
0155     
0156     alpha = 1; epsilon = 1;
0157     beta = M*(norm(xh_old,2)/norm(xh_old,1))^2;
0158     W = tau/alpha./(beta*abs(xh_old)+epsilon);
0159     
0160     % yh = A*xh_old;
0161     % u = -sign(xh_old).*W_old; % out.pk;
0162     % Atdyu = A'*(yh-y)-u;
0163     
0164     yh = A*xh_old;
0165     Atr = A'*(A*xh_old-y);
0166     u =  -W.*sign(xh_old)-Atr;
0167     pk_old = Atr+u;
0168     
0169     in = out;
0170     in.xh_old = xh_old;
0171     in.pk_old = pk_old;
0172     in.u = u;
0173     in.W_old = W_old;
0174     in.W = W;
0175     % in.Atdyu = Atdyu;
0176 end
0177 W = W_old;
0178  
0179 
0180 % in = [];
0181 % in.tau = tau;
0182 % in.delx_mode = delx_mode;
0183 % in.debias = 0;
0184 % in.verbose = 0;
0185 % in.plots = 0;
0186 % in.record = 1;
0187 % in.err_fun = err_fun;
0188 % tic
0189 % out = l1homotopy(A,y,in);
0190 % xh = out.x_out;
0191 % iter_bpdn = out.iter;
0192 % time_bpdn = toc;
0193 % gamma_bpdn = out.gamma;
0194 % err_bpdn = out.error_table;
0195 %
0196 % W = tau;
0197 %
0198 % x_sparsa = 0;
0199 % tic;
0200 % [x_sparsa,x_debias_SpaRSA_m,obj_SpaRSA_m_cont,...
0201 %     times_SpaRSA_m_cont,debias_start_SpaRSA_m,mse_SpaRSA_m,taus_SpaRSA_m, numA, numAt]= ...
0202 %     SpaRSA_adpW(y,A,tau,...
0203 %     'Monotone',0,...
0204 %     'adp_wt',0,...
0205 %     'Debias',0,...
0206 %     'StopCriterion',2,...
0207 %     'ToleranceA',1e-4,...
0208 %     'Safeguard',1,...
0209 %     'MaxiterA',maxiter,...
0210 %     'Verbose',0,...
0211 %     'True_x',x,...
0212 %     'Continuation',1,...
0213 %     'Continuationsteps',-1);
0214 %
0215 % time_sparsa = toc;
0216 % iter_sparsa = (numA+numAt)/2;
0217 % error_sparsa = norm(x-x_sparsa)/norm(x);
0218 
0219 
0220 for sim = 1:maxsim
0221     
0222     % Update the solution after updating the measurement matrix and/or the
0223     % sparse signal
0224     x_old = x; sig_old = sig;
0225     y_old = y; A_old = A(1:M,:);
0226     
0227     % Time-varying signal setup
0228     dyn_ref = 'previous';
0229     switch dyn_ref
0230         case 'previous'
0231             % In this model, x may become dense after some time
0232             
0233             % Update the sparse signal.
0234             rshift = rshift_h(1);
0235             x = F_h(x_old,cshift, rshift);
0236             % Threshold the original signal...
0237             % (to stop it from becoming dense) ???
0238             % [x_sort ind_sort] = sort(abs(x),'descend');
0239             % x(ind_sort(M+1:end)) = 0;
0240             
0241             % Update the signal estimate
0242             xh_old = F_h(xh,cshift,0);
0243             % xh_old(abs(xh_old)<max(xh_old)/1e3) = 0;
0244             xh_old(abs(xh_old)<tau/sqrt(log(N))) = 0;
0245             
0246         case 'initial'
0247             % In this model, x will probably remain sparse for all times
0248             
0249             % Update the sparse signal.
0250             cshift = -sim; rshift = rshift_h(1);
0251             x = F_h(x_orig,cshift, rshift);            
0252             
0253             % Update the signal estimate
0254             xh_old = F_h(xh,cshift-sim+1,0);
0255             % xh_old(abs(xh_old)<max(xh_old)/1e3) = 0;
0256             xh_old(abs(xh_old)<tau/sqrt(log(N))) = 0;
0257     end
0258     
0259     % Update weights
0260     W_old = W;
0261     
0262     alpha = 1; epsilon = 1;
0263     beta = M*(norm(xh_old,2)/norm(xh_old,1))^1;
0264     W = tau/alpha./(beta*abs(xh_old)+epsilon);
0265     
0266     % Update the measurements
0267     in = []; in.type = mType;
0268     B = genAmat(M,N,in);
0269     e = randn(M,1)*sigma;
0270 
0271     % add difference regularization
0272     A = [0*A_old+1*B; lambda_diff*eye(N)];
0273     y = [A(1:M,:)*x+e; lambda_diff*xh_old];
0274     
0275     % To remove diff. regularization, set lambda_diff = 0;
0276     if lambda_diff == 0
0277         A = A(1:M,:);
0278         y = y(1:M,:);
0279     end
0280     
0281     % parameter selection (you may change it at every iteration...)
0282     %         % tau = sigma*sqrt(log(N));
0283     %         tau = max(1e-3*max(abs(A'*y)),sigma*sqrt(log(N)));
0284     
0285     homotopy_mode = 'dummy';
0286     switch homotopy_mode
0287         case 'dummy'
0288             % create a dummy variable...
0289             % use homotopy on the measurements...
0290             % in principle, we can start with any xh_old with this formulation
0291             % and any starting value of tau or W...
0292             W_old = W;
0293             gamma = find(xh_old);
0294             M_trunc = size(A,1); % M-1;
0295             if length(gamma) > M_trunc
0296                 [xh_sort ind_sort] = sort(abs(xh_old),'descend');
0297                 xh_old(ind_sort(M_trunc+1:end)) = 0;
0298                 gamma = ind_sort(1:M_trunc);            
0299             end            
0300             % yh = A*xh_old;
0301             % u =  -sign(xh_old).*W_old;
0302             % Atdyu = A'*(yh-y)-u;
0303             yh = A*xh_old;
0304             Atr = A'*(A*xh_old-y);
0305             u =  -W.*sign(xh_old)-Atr;
0306             pk_old = Atr+u;
0307         otherwise
0308             didp('Go back ... no escape');
0309     end
0310     
0311     tic
0312     in = out;
0313     gamma_old = gamma;
0314     in.gamma = gamma_old;
0315     switch delx_mode
0316         case 'mil';
0317             % in.delx_mode = 'mil';
0318             % The following gram matrix and its inverse can be used from the
0319             % previous homotopy. Too lazy to include that right now...
0320             % wt BPDN homotopy update
0321             AtAgx = A(:,gamma_old)'*A(:,gamma_old);
0322             iAtAgx = inv(AtAgx);
0323             in.iAtA = iAtAgx;
0324         case {'qr','chol'};
0325             % in.delx_mode = 'qr';
0326             [Q R] = qr(A(:,gamma_old),0);
0327             in.Q = Q; in.R = R;
0328         case 'qrM'
0329             % in.delx_mode = 'qrM';
0330             [Q0 R0] = qr(A(:,gamma_old));
0331             in.Q0 = Q0; in.R0 = R0;
0332     end
0333     
0334     in.xh_old = xh_old;
0335     in.pk_old = pk_old;
0336     in.u = u;
0337     % in.Atdyu = Atdyu;
0338     in.W = W;
0339     in.W_old = W_old;
0340     in.delx_mode = delx_mode;
0341     in.debias = 0;
0342     in.verbose = 0;
0343     in.plots = 0;
0344     in.record = 0;
0345     in.err_fun = @(z) (norm(x-z)/norm(x))^2;
0346     out = l1homotopy(A,y,in);
0347     xh_dynRWT = out.x_out;
0348     gamma_dynRWT = out.gamma;
0349     iter_dynRWT = out.iter;
0350     time_dynRWT = toc;
0351     
0352     xh = xh_dynRWT;
0353     
0354     %% Check the solution using BPDN directly with l1homotopy
0355     %     in = []; x_old = x;
0356     %     in.W = W;
0357     %     in.delx_mode = delx_mode;
0358     %     in.debias = 0;
0359     %     in.verbose = 0;
0360     %     in.maxiter = maxiter;
0361     %     in.plots = 0;
0362     %     in.record = 1;
0363     %     in.err_fun = err_fun;
0364     %     tic
0365     %     out = l1homotopy(A,y,in);
0366     %     xh = out.x_out;
0367     %     iter_bpdn = out.iter;
0368     %     time_bpdn = toc;
0369     %     gamma_bpdn = out.gamma;
0370     %     err_bpdn = out.error_table;
0371     
0372     %% SpaRSA
0373     x_sparsa = xh_old; W_sparsa = W/tau;
0374     psi_function = @(x,tau) soft(x,tau*W_sparsa);
0375     phi_function = @(x) sum(abs(W_sparsa.*x));
0376     tic;
0377     [x_sparsa,x_debias_SpaRSA_m,obj_SpaRSA_m_cont,...
0378         times_SpaRSA_m_cont,debias_start_SpaRSA_m,mse_SpaRSA_m,taus_SpaRSA_m, numA, numAt]= ...
0379         SpaRSA_adpW(y,A,tau,...
0380         'Monotone',0,...
0381         'adp_wt',0,...
0382         'W_new',W_sparsa,...
0383         'Debias',0,...
0384         'Initialization',x_sparsa,...
0385         'StopCriterion',2,...
0386         'ToleranceA',1e-4,...
0387         'psi',psi_function,...
0388         'phi',phi_function,...
0389         'Safeguard',1,...
0390         'MaxiterA',maxiter,...
0391         'Verbose',0,...
0392         'True_x',x,...
0393         'Continuation',1,...
0394         'Continuationsteps',-1);
0395     
0396     time_sparsa = toc;
0397     iter_sparsa = (numA+numAt)/2;    
0398     error_sparsa = norm(x-x_sparsa)/norm(x);     
0399     
0400     %%
0401     SIM_stack{sim} = [sim, tau, ...
0402         norm(x-xh_dynRWT)^2/norm(x)^2, sum(iter_dynRWT,2), sum(time_dynRWT,2),...
0403         norm(x-x_sparsa)^2/norm(x)^2, sum(iter_sparsa,2), sum(time_sparsa,2)];
0404     
0405     fprintf('sim %d. tau = %3.4g, (err,iter,time): dynRWT l1homotopy-%3.4g,%3.4g,%3.4g; SpaRSA-%3.4g,%3.4g,%3.4g. \n', ...
0406         SIM_stack{sim});
0407     
0408     %% plot recovered signals
0409     figure(1); 
0410     subplot(211); 
0411     plot([x xh xh_dynRWT]);
0412     subplot(212); 
0413     plot([iW_h(x) iW_h(xh)]);
0414     title('Comparison betweeen the original and reconstructed signal')
0415     
0416     %% Reconstructed signal
0417     x_vec((sim-1)*N+1:sim*N,1) = x;
0418     xh_vec((sim-1)*N+1:sim*N,1) = xh_dynRWT; 
0419 end
0420 % mS =  mean(cell2mat(SIM_stack),1);
0421 mS =  mean(cell2mat(SIM_stack(ceil(maxsim*0.2):maxsim,:)),1);
0422 fprintf('Average results for last 80 percent iterations: maxsim %d. tau = %3.4g, (err,iter,time): dynRWT l1homotopy-%3.4g,%3.4g,%3.4g; SpaRSA-%3.4g,%3.4g,%3.4g. \n', maxsim, mS(2:end));