Description of demo

0001 % demo_streaming
0002 %
0003 % Solves the following streaming problem
0004 % min_x  \sum_i \|x_i-1\|_1 + 1/2*||B_i-1 x_i-1 + A_i x_i - y_i||_2^2
0005 %
0006 % by adding and removing one set of measurements at every time instant
0007 %
0008 % Applications:
0009 %       streaming signal recovery from overlapping measurements
0010 %
0011 % We can also add any other regularization operator in the reconstruction
0012 %
0013 % Written by: Salman Asif, Georgia Tech
0014 % Email: sasif@gatech.edu
0015 % Created: August 2012
0016 
0017 clear
0018 close all % force
0019 
0020 % Limit the number of computational threads (for profiling)
0021 % maxNumCompThreads(1);
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/qr/
0031 addpath solvers
0032 
0033 disp(['--------------------',datestr(now),'-------------------------'])
0034 
0035 % load RandomStates
0036 %
0037 rseed = 2012;
0038 rand('state',rseed);
0039 randn('state',rseed);
0040 
0041 % simulation parameters
0042 mType = 'sign'; % {'randn','orth','rdct'};
0043 sType = 'pcwreg'; % {'randn','sign','highD', 'blocks','pcwPoly'}
0044 SNR = inf;       % additive Gaussian noise
0045 
0046 N = 128;   % signal length
0047 M = round(N/4);    % no. of measurements
0048 T = round(M/4);    % sparsity level
0049 
0050 % Overlapping measurements
0051 LM = 2*N; 
0052 
0053 % streaming window
0054 P = 4; % size of the working window is P*N
0055 
0056 % remove the top-left edge of the system matrix
0057 LeftEdge_trunc = 1;
0058 % fprintf('Consider oldest set of coefficients to be fully known, and remove their contribution from the measurements... \n');
0059     
0060 % rank-1 update mode
0061 delx_mode = 'mil'; % mil or qr
0062 fprintf('Streaming measurements with a time-varying signal ... \n');
0063 str0 = sprintf('mType-%s, sType-%s, SNR = %d, (N,M,T) = %d, %d, %d, P = %d.', mType, sType, SNR, N, M, T, P);
0064 disp(str0);
0065 
0066 % Simulation parameters
0067 maxsim = 128;
0068 SIM_stack = cell(maxsim,1);
0069 SIM_memory = cell(maxsim,1);
0070  
0071 x_vec = zeros(N*maxsim,1);
0072 xh_vec = zeros(N*maxsim,1);
0073 sig_vec = zeros(N*maxsim,1);
0074 sigh_vec = zeros(N*maxsim,1);
0075 
0076 %% Signal and system Setup
0077 % Generate long sequence of signal and sensing matrices
0078 
0079 % Generate a random signal
0080 in = []; in.type = sType; in.T = maxsim*T; in.randgen = 0; 
0081 in.wType = 'daub79'; in.take_fwt = 1; in.J = 5; %log2(N)-3;
0082 genSignal_h = @(N) genSignal(N,in);
0083 [x_long sig wave_struct] = genSignal(maxsim*N,in);
0084 
0085 if isempty(wave_struct)
0086     W_h = @(z) z; iW_h = @(z) z; 
0087 else 
0088     W_h = wave_struct.W_h;
0089     iW_h = wave_struct.iW_h;    
0090 end
0091 
0092 in = []; in.type = mType;
0093 genAmat_h = @(M,N) genAmat(M,N,in);
0094 
0095 x = zeros(N*P+LM-N,1);
0096 A = zeros(M*P,N*P+LM-N);
0097 % A = sparse(M*P,N*(P+1));
0098 y = zeros(M*P,1);
0099  
0100 % Overlapping measurements (type imagesc(At) to view the structure)
0101 % LM = 2*N...
0102 s_ind = 0;
0103 % xt = [genSignal_h(N); genSignal_h(N)];
0104 xt = [W_h(sig(1:N)); W_h(sig(N+1:2*N))];
0105 s_ind = s_ind+2*N;
0106 
0107 At = genAmat_h(M,LM);
0108 sigma = sqrt(norm(At*xt)^2/10^(SNR/10)/M);
0109 et = randn(M,1)*sigma;
0110 yt = At*xt+et;
0111 
0112 x(1:2*N) = xt;
0113 A(1:M,1:2*N) = At;
0114 y(1:M) = yt;
0115 
0116 for ii = 2:P
0117     % xt = genSignal_h(N);
0118     xt = W_h(sig(s_ind+1:s_ind+N));
0119     s_ind = s_ind+N;
0120     
0121     x(ii*N+1:(ii+1)*N) = xt;
0122     
0123     At = genAmat_h(M,2*N);
0124     A((ii-1)*M+1:ii*M,(ii-1)*N+1:(ii+1)*N) = At;
0125     
0126     et = randn(M,1)*sigma;
0127     y((ii-1)*M+1:ii*M) = At*x((ii-1)*N+1:(ii+1)*N)+et;
0128 end
0129 
0130 % parameter selection
0131 % tau = sigma*sqrt(log(N));
0132 tau = max(1e-2*max(abs(A'*y)),sigma*sqrt(log(2*N)));
0133 
0134 
0135 maxiter = 2*P*N;
0136 err_fun = @(z) (norm(x-z)/norm(x))^2;
0137 
0138 %% Initialize by solving a rwt L1 problem
0139 in = [];
0140 in.tau = tau; W = tau;
0141 in.delx_mode = delx_mode;
0142 in.debias = 0;
0143 in.verbose = 0;
0144 in.plots = 0;
0145 in.record = 1;
0146 in.err_fun = err_fun;
0147 tic 
0148 for wt_itr = 1:5
0149     
0150     out = l1homotopy(A,y,in);
0151     xh = out.x_out;
0152     iter_bpdn = out.iter;
0153     time_bpdn = toc;
0154     gamma_bpdn = out.gamma;
0155     err_bpdn = out.error_table;
0156     
0157     % Update weights
0158     xh_old = xh;
0159     
0160     alpha = 1; epsilon = 1;
0161     beta = M*(norm(xh_old,2)/norm(xh_old,1))^2;
0162     W = tau/alpha./(beta*abs(xh_old)+epsilon);
0163     
0164     W_old = W;
0165     % yh = A*xh_old;
0166     % u = -sign(xh_old).*W_old; % out.pk;
0167     % Atdyu = A'*(yh-y)-u;
0168     yh = A*xh_old;
0169     Atr = A'*(A*xh_old-y);
0170     u =  -W.*sign(xh_old)-Atr;
0171     pk_old = Atr+u;
0172     
0173     in = out;
0174     in.xh_old = xh_old;
0175     in.pk_old = pk_old;
0176     in.u = u;
0177     in.W_old = W_old;
0178     in.W = W;
0179     % in.Atdyu = Atdyu;
0180 end 
0181 W = W_old;
0182 
0183 for sim = 1:maxsim
0184     
0185     %% Update the solution after updating the measurement matrix and/or the
0186     % sparse signal
0187     x_old = x;  
0188     y_old = y; A_old = A;
0189     
0190     % Time-varying signal setup
0191     % Shift the window (apply any prediction?)
0192     x(1:end-N) = x(N+1:end); 
0193     xt_old = x(end-N+1:end);
0194 
0195     % xt = genSignal_h(N);
0196     if s_ind+N > length(sig)
0197         break;
0198     end
0199     xt = W_h(sig(s_ind+1:s_ind+N));
0200     s_ind = s_ind+N;
0201     x(end-N+1:end) = xt; 
0202     
0203     % System matrix setup...
0204     % Shift up and left
0205     A(1:(P-1)*M,1:P*N) = A(M+1:end,N+1:end);     
0206     % new measurement matrix
0207     At = genAmat_h(M,2*N); 
0208     A(end-M+1:end,end-2*N+1:end) = At;
0209     
0210     % shift old measurements and add one new set of measurementts
0211     y(1:end-M) = y(M+1:end);
0212     e = randn(M,1)*sigma;
0213     yt = At*x(end-2*N+1:end)+e;
0214     y(end-M+1:end) = yt;
0215     
0216     % Update the signal estimate (for warm start)
0217     xh_old = xh;
0218     xh_old(1:end-N) = xh(N+1:end);
0219     xh_t = xh(end-N+1:end);
0220     xh_old(end-N+1:end) = 0*xh_t;   
0221     xh_old(abs(xh_old)<tau/sqrt(log(P*N))) = 0;
0222     
0223     % Remove the top-left edge of the system matrix
0224     if LeftEdge_trunc
0225         A_old = A; y_old = y;
0226         
0227         alpha0h = xh_old(1:N);
0228         xh_old = xh_old(N+1:end);
0229         y = y-A(:,1:N)*alpha0h;
0230         
0231         A = A(:,N+1:end);
0232         alpha0 = x(1:N);
0233         x = x(N+1:end);
0234     end
0235     
0236     % Update weights
0237     alpha = 1; epsilon = 1;
0238     beta = M*(norm(xh_old,2)/norm(xh_old,1))^1;
0239     W = tau/alpha./(beta*abs(xh_old)+epsilon);
0240     % W(end-N+1:1) = tau;
0241     W_old = W;
0242     
0243     homotopy_mode = 'dummy';
0244     switch homotopy_mode
0245         case 'dummy'
0246             % create a dummy variable...
0247             % use homotopy on the measurements...
0248             % in principle, we can start with any xh_old with this formulation
0249             % and any starting value of tau or W...
0250             gamma = find(xh_old);
0251             M_trunc = size(A,1); % P*(M-1);
0252             if length(gamma) > M_trunc
0253                 [xh_sort ind_sort] = sort(abs(xh_old),'descend');
0254                 xh_old(ind_sort(M_trunc+1:end)) = 0;
0255                 gamma = ind_sort(1:M_trunc);            
0256             end            
0257             % yh = A*xh_old;
0258             % u =  -sign(xh_old).*W_old;
0259             % Atdyu = A'*(yh-y)-u;
0260             
0261             yh = A*xh_old;
0262             Atr = A'*(A*xh_old-y);
0263             u =  -W.*sign(xh_old)-Atr;
0264             pk_old = Atr+u;
0265         otherwise
0266             didp('Go back ... no escape');
0267     end
0268     
0269     in = out;
0270     gamma_old = gamma;
0271     in.gamma = gamma_old;
0272     switch delx_mode
0273         case 'mil';
0274             % in.delx_mode = 'mil';
0275             % The following gram matrix and its inverse can be used from the
0276             % previous homotopy. Too lazy to include that right now...
0277             % wt BPDN homotopy update
0278             AtAgx = A(:,gamma_old)'*A(:,gamma_old);
0279             iAtAgx = inv(AtAgx);
0280             in.iAtA = iAtAgx;
0281         case {'qr','chol'};
0282             % in.delx_mode = 'qr';
0283             [Q R] = qr(A(:,gamma_old),0);
0284             in.Q = Q; in.R = R;
0285         case 'qrM'
0286             % in.delx_mode = 'qrM';
0287             [Q0 R0] = qr(A(:,gamma_old));
0288             in.Q0 = Q0; in.R0 = R0;
0289     end
0290     
0291     in.xh_old = xh_old;
0292     in.pk_old = pk_old;
0293     in.u = u;
0294     % in.Atdyu = Atdyu;
0295     in.W = W;
0296     in.W_old = W_old;
0297     in.delx_mode = delx_mode;
0298     in.debias = 0;
0299     in.verbose = 0;
0300     in.plots = 0;
0301     in.record = 1;
0302     in.err_fun = @(z) (norm(x-z)/norm(x))^2;
0303     tic
0304     out = l1homotopy(A,y,in);
0305     time_streamingRWT = toc;
0306     xh_streamingRWT = out.x_out;
0307     gamma_streamingRWT = out.gamma;
0308     iter_streamingRWT = out.iter;
0309     
0310     xh = xh_streamingRWT;
0311     
0312     % Check the solution using BPDN directly with l1homotopy
0313     %     in = []; x_old = x;
0314     %     in.W = W;
0315     %     in.delx_mode = delx_mode;
0316     %     in.debias = 0;
0317     %     in.verbose = 0;
0318     %     in.maxiter = maxiter;
0319     %     in.plots = 0;
0320     %     in.record = 1;
0321     %     in.err_fun = err_fun;
0322     %     tic
0323     %     out = l1homotopy(A,y,in);
0324     %     xh = out.x_out;
0325     %     iter_bpdn = out.iter;
0326     %     time_bpdn = toc;
0327     %     gamma_bpdn = out.gamma;
0328     %     err_bpdn = out.error_table;
0329     
0330     %% SpaRSA
0331     x_sparsa = xh_old; W_sparsa = W/tau;
0332     psi_function = @(x,tau) soft(x,tau*W_sparsa);
0333     phi_function = @(x) sum(abs(W_sparsa.*x));
0334     tic;
0335     [x_sparsa,x_debias_SpaRSA_m,obj_SpaRSA_m_cont,...
0336         times_SpaRSA_m_cont,debias_start_SpaRSA_m,mse_SpaRSA_m,taus_SpaRSA_m, numA, numAt]= ...
0337         SpaRSA_adpW(y,A,tau,...
0338         'Monotone',0,...
0339         'adp_wt',0,...
0340         'W_new',W_sparsa,...
0341         'Debias',0,...
0342         'Initialization',x_sparsa,...
0343         'StopCriterion',2,...
0344         'ToleranceA',1e-4,...
0345         'psi',psi_function,...
0346         'phi',phi_function,...
0347         'Safeguard',1,...
0348         'MaxiterA',maxiter,...
0349         'Verbose',0,...
0350         'True_x',x,...
0351         'Continuation',1,...
0352         'Continuationsteps',-1);
0353     
0354     time_sparsa = toc;
0355     iter_sparsa = (numA+numAt)/2;    
0356     error_sparsa = norm(x-x_sparsa)/norm(x);        
0357     
0358     %%
0359     SIM_stack{sim} = [sim, tau, ...
0360         norm(x-xh_streamingRWT)^2/norm(x)^2, sum(iter_streamingRWT,2), sum(time_streamingRWT,2), ... 
0361         norm(x-x_sparsa)^2/norm(x)^2, sum(iter_sparsa,2), sum(time_sparsa,2)];
0362     
0363     % l1homotopy-%3.4g,%3.4g,%3.4g;
0364     % norm(x-xh)^2/norm(x)^2, sum(iter_bpdn,2), sum(time_bpdn,2), ...
0365     
0366 %     fprintf('sim %d. tau = %3.4g, (err,iter,time): streamingRWT homotopy-%3.4g,%3.4g,%3.4g, SpaRSA-%3.4g,%3.4g,%3.4g. \n', ...
0367 %         SIM_stack{sim});
0368 %
0369     %% plot sparse coeffs.
0370     fig(1);     
0371     subplot(211);
0372     plot([x xh_streamingRWT x_sparsa]); 
0373     title('Comparison betweeen the original and reconstructed signal')
0374     drawnow;
0375     
0376     %% Reconstructed signal
0377     if LeftEdge_trunc
0378         xh = [alpha0h; xh_streamingRWT];
0379         x = [alpha0; x];
0380         A = A_old; y = y_old;
0381     else
0382         xh = xh_streamingRWT;
0383     end
0384     x_vec((sim-1)*N+1:sim*N,1) = x(1:N);
0385     xh_vec((sim-1)*N+1:sim*N,1) = xh(1:N); 
0386     sig_vec((sim-1)*N+1:sim*N,1) = iW_h(x(1:N));
0387     sigh_vec((sim-1)*N+1:sim*N,1) = iW_h(xh(1:N)); 
0388     
0389     %% plot recovered signals
0390     subplot(212);
0391     plot([sig_vec(1:s_ind(end)) sigh_vec(1:s_ind(end))]); 
0392         
0393 end
0394 mS =  mean(cell2mat(SIM_stack),1);
0395 fprintf('Average results: maxsim %d. tau = %3.4g, (err,iter,time): streamingRWT l1homotopy-%3.4g,%3.4g,%3.4g; SpaRSA-%3.4g,%3.4g,%3.4g. \n', maxsim, mS(2:end));
0396 % l1homotopy-%3.4g,%3.4g,%3.4g;
0397 
0398 % l1homotopy-%3.4g,%3.4g,%3.4g;
0399 L = N;
0400 err_l1homotopy = norm(sig_vec(L-N+1:N*sim)-sigh_vec(L-N+1:N*sim,1))^2/norm(sig_vec(L-N+1:N*sim))^2; 
0401 fprintf('Signal error: l1homotopy-%3.4g. \n',err_l1homotopy);
0402
demo_streaming

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
