Description of demo

0001 % demo_dynamicSeq
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 as new measurements w = Bx are received.
0007 %
0008 % Written by: Salman Asif, Georgia Tech
0009 % Email: sasif@gatech.edu
0010 % Created: August 2012
0011 %
0012 % Modified: May 2013 (using l1homotopy for update)
0013 %
0014 clear
0015 close all force
0016   
0017 %% Setup path
0018 mname = mfilename;
0019 mpath = mfilename('fullpath');
0020 mdir = mpath(1:end-length(mname));
0021 cd(mdir);
0022 
0023 addpath utils/
0024 addpath utils/qr/
0025 
0026 disp(['--------------------',datestr(now),'-------------------------'])
0027 
0028 % load RandomStates
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 = inf;       % additive Gaussian noise
0038   
0039 N = 512;   % signal length
0040 M = round(N/2);    % no. of measurements
0041 T = round(M/4);    % sparsity level
0042 Mn = round(M/10);  % new measurements
0043 
0044 % rank-1 update mode
0045 delx_mode = 'mil'; % mil or qr
0046 fprintf('Add sequential measurements..\n');
0047 str0 = sprintf('mType-%s, sType-%s, SNR = %d, (N,M,T) = %d, %d, %d; Mn = %d.', mType, sType, SNR, N, M, T,Mn);
0048 disp(str0);
0049 
0050  
0051 maxsim = 10;
0052 SIM_stack = cell(maxsim,1);
0053 SIM_memory = cell(maxsim,1);
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     % measurement matrix
0065     in = []; in.type = mType;
0066     A = genAmat(M,N,in);
0067     A = A*sqrt(M); % Because genAmat generates matrices with unit (expected) norm columns
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 after adding measurements (DynamicSeq)
0103     xh_old = xh; tau_old = tau;
0104     x_old = x; y_old = y; A_old = A;
0105     
0106     % Additional measurements
0107     in = []; in.type = mType;
0108     B = genAmat(Mn,N,in);
0109     B = B*sqrt(Mn); % Because genAmat generates matrices with unit (expected) norm columns
0110     
0111     % noise in the new measurements
0112     sigma = sqrt(norm(B*x)^2/10^(SNR/10)/Mn);    
0113     e = randn(Mn,1)*sigma;
0114     w = B*x+e;
0115     
0116     A = [A_old; B];
0117     y = [y_old; w];
0118     
0119     % parameter selection
0120     % tau = sigma*sqrt(log(N));
0121     tau = max(1e-4*max(abs(A'*y)),sigma*sqrt(log(N)));    
0122     
0123     homotopy_mode = 'dummy';
0124     switch homotopy_mode
0125         case 'dummy'
0126             % create a dummy variable...
0127             % use homotopy on the measurements...
0128             % in principle, we can start with any xh_old with this formulation
0129             % and any starting value of tau or W...
0130             W = tau; 
0131             W_old = tau_old;
0132             Atr = A'*(A*xh_old-y);
0133             u =  -W.*sign(xh_old)-Atr;
0134             pk_old = Atr+u;
0135         otherwise
0136             % Use direct information from the old solution
0137             W = tau; 
0138             W_old = tau_old;
0139             %             u = out.pk;
0140             %             yh = [y_old; B*xh_old];
0141             %             Atdyu = A'*(yh-y);
0142             Atr = A'*(A*xh_old-y);
0143             u =  -W.*sign(xh_old)-Atr;
0144             pk_old = Atr+u;
0145     end 
0146     
0147     in = out;
0148     gamma_old = in.gamma;
0149     switch delx_mode
0150         case 'mil';
0151             % in.delx_mode = 'mil';
0152             % The following gram matrix and its inverse can be used from the
0153             % previous homotopy. Too lazy to include that right now...
0154             % wt BPDN homotopy update
0155             AtAgx = A(:,gamma_old)'*A(:,gamma_old);
0156             iAtAgx = inv(AtAgx);
0157             in.iAtA = iAtAgx;
0158         case {'qr','chol'};
0159             % in.delx_mode = 'qr';
0160             [Q R] = qr(A(:,gamma_old),0);
0161             in.Q = Q; in.R = R;
0162         case 'qrM'
0163             % in.delx_mode = 'qrM';
0164             [Q0 R0] = qr(A(:,gamma_old));
0165             in.Q0 = Q0; in.R0 = R0;
0166     end
0167 
0168     in.xh_old = xh_old;
0169     in.pk_old = pk_old;
0170     in.u = u;
0171     in.W = W;
0172     in.W_old = W_old;
0173     in.delx_mode = delx_mode;
0174     in.debias = 0;
0175     in.verbose = 0;
0176     in.plots = 0;
0177     in.record = 1;
0178     in.err_fun = @(z) (norm(x-z)/norm(x))^2;
0179     tic
0180     out = l1homotopy(A,y,in);
0181     xh_dynSeq = out.x_out;
0182     gamma_dynSeq = out.gamma;
0183     iter_dynSeq = out.iter;
0184     time_dynSeq = toc;
0185          
0186     SIM_stack{sim} = [sim, tau, ...
0187         norm(x_old-xh)^2/norm(x)^2, sum(iter_bpdn,2), sum(time_bpdn,2), ...
0188         norm(x-xh_dynSeq)^2/norm(x)^2, sum(iter_dynSeq,2), sum(time_dynSeq,2)];  
0189     
0190     fprintf('sim %d. tau = %3.4g, (err,iter,time): l1homotopy from scratch-%3.4g,%3.4g,%3.4g; dynSeq l1homotopy-%3.4g,%3.4g,%3.4g. \n', ...
0191         SIM_stack{sim});
0192    
0193     %% plot recovered signals
0194     figure(1); plot([x xh xh_dynSeq]);
0195     title('Comparison betweeen the new and old homotopy code')
0196 end
0197 mS =  mean(cell2mat(SIM_stack),1);
0198 fprintf('Average results: maxsim %d. tau = %3.4g, (err,iter,time): l1homotopy from scratch-%3.4g,%3.4g,%3.4g; dynSeq l1homotopy-%3.4g,%3.4g,%3.4g. \n', maxsim, mS(2:end));
demo_dynamicSeq

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
