demo_replaceDict

 demo_replaceDict

 Solves the following BPDN problem
 min_x  tau\|x\|_1 + 1/2*||y-Ax||_2^2
 
 and updates the solution as the A matrix is updated... 
 Infact, we can modify all the terms: y, A, and x simultaneously

 Applications:
       in dictionary learning where A is iteratively updated
       in tracking where y, A, and x, are updated..
 
 Written by: Salman Asif, Georgia Tech
 Email: sasif@gatech.edu
 Created: August 2012

0001 % demo_replaceDict
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 the A matrix is updated...
0007 % Infact, we can modify all the terms: y, A, and x simultaneously
0008 %
0009 % Applications:
0010 %       in dictionary learning where A is iteratively updated
0011 %       in tracking where y, A, and x, are updated..
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 
0032 disp(['--------------------',datestr(now),'-------------------------'])
0033 
0034 % load RandomStates
0035 %
0036 rseed = 2012;
0037 rand('state',rseed);
0038 randn('state',rseed);
0039 
0040 % simulation parameters
0041 mType = 'randn'; % {'randn','orth','rdct'};
0042 sType = 'randn'; % {'randn','sign','highD', 'blocks','pcwPoly'}
0043 SNR = inf;       % additive Gaussian noise
0044 
0045 % reweighted setup
0046 rwt = 5;        % number of reweighting iterations
0047 
0048 N = 512;   % signal length
0049 M = round(N/2);    % no. of measurements
0050 T = round(M/4);    % sparsity level
0051 Tn = 0; % round(T/10);  % update in the signal...
0052 
0053 % rank-1 update mode
0054 delx_mode = 'mil'; % mil or qr
0055 fprintf('Replace matrix A and/or signal x..\n');
0056 str0 = sprintf('mType-%s, sType-%s, SNR = %d, (N,M,T) = %d, %d, %d.', mType, sType, SNR, N, M, T);
0057 disp(str0);
0058  
0059 maxsim = 10;
0060 SIM_stack = cell(maxsim,1);
0061 SIM_memory = cell(maxsim,1);
0062 
0063  
0064 for sim = 1:maxsim
0065     
0066     % Generate a random signal
0067     in = []; in.type = sType; in.T = T; in.randgen = 1;
0068     x = genSignal(N,in);
0069     [val ind] = sort(abs(x),'descend');
0070     ind_pos = ind(val>1e-1);
0071     gamma_orig = ind_pos(1:min(length(ind_pos),M-1));
0072         
0073     % measurement matrix
0074     in = []; in.type = mType;
0075     A = genAmat(M,N,in);
0076     
0077     % measurements
0078     sigma = sqrt(norm(A*x)^2/10^(SNR/10)/M);    
0079     e = randn(M,1)*sigma;
0080     y = A*x+e;
0081     
0082     %     % orthogonalize rows of A
0083     %     [Q, R] = qr(A',0);
0084     %     A = Q'; y = R' \ y;
0085     
0086     % parameter selection
0087     % tau = sigma*sqrt(log(N));
0088     tau = max(1e-4*max(abs(A'*y)),sigma*sqrt(log(N)));
0089     
0090     maxiter = 4*N;
0091     err_fun = @(z) (norm(x-z)/norm(x))^2;
0092     
0093     %% BPDN using l1homotopy
0094     in = [];
0095     in.tau = tau;
0096     in.delx_mode = delx_mode;
0097     in.debias = 0;
0098     in.verbose = 0;
0099     in.plots = 0;
0100     in.record = 1;
0101     in.err_fun = err_fun;
0102     tic
0103     out = l1homotopy(A,y,in);
0104     xh = out.x_out;
0105     iter_bpdn = out.iter;
0106     time_bpdn = toc;
0107     gamma_bpdn = out.gamma;
0108     err_bpdn = out.error_table;
0109     
0110     %% Update the solution after updating the measurement matrix and/or the
0111     % sparse signal
0112     
0113     xh_old = xh; tau_old = tau;
0114     x_old = x; y_old = y; A_old = A;
0115     
0116     % Update the measurement matrix
0117     in = []; in.type = mType;
0118     B = genAmat(M,N,in);
0119     
0120     A = A_old+1/5*B;
0121     
0122     % Update the sparse signal.
0123     in = []; in.type = sType; in.T = Tn; in.randgen = 1;
0124     dx = genSignal(N,in);
0125     x = x_old + dx;
0126     % e = randn(M,1)*sigma;
0127     y = A*x+e;
0128            
0129     % parameter selection
0130     % tau = sigma*sqrt(log(N));
0131     tau = max(1e-4*max(abs(A'*y)),sigma*sqrt(log(N)));
0132     
0133     homotopy_mode = 'dummy';
0134     switch homotopy_mode
0135         case 'dummy'
0136             % create a dummy variable...
0137             % use homotopy on the measurements...
0138             % in principle, we can start with any xh_old with this formulation
0139             % and any starting value of tau or W...
0140             W = tau; 
0141             W_old = tau_old;
0142             % yh = A*xh_old;
0143             % u =  -sign(xh_old).*W_old;
0144             % Atdyu = A'*(yh-y)-u;
0145             yh = A*xh_old;
0146             Atr = A'*(A*xh_old-y);
0147             u =  -W.*sign(xh_old)-Atr;
0148             pk_old = Atr+u;
0149         otherwise            
0150             didp('Go back ... no escape');
0151     end 
0152     
0153     in = out;
0154     gamma_old = in.gamma;
0155     switch delx_mode
0156         case 'mil';
0157             % in.delx_mode = 'mil';
0158             % The following gram matrix and its inverse can be used from the
0159             % previous homotopy. Too lazy to include that right now...
0160             % wt BPDN homotopy update
0161             AtAgx = A(:,gamma_old)'*A(:,gamma_old);
0162             iAtAgx = inv(AtAgx);
0163             in.iAtA = iAtAgx;
0164         case {'qr','chol'};
0165             % in.delx_mode = 'qr';
0166             [Q R] = qr(A(:,gamma_old),0);
0167             in.Q = Q; in.R = R;
0168         case 'qrM'
0169             % in.delx_mode = 'qrM';
0170             [Q0 R0] = qr(A(:,gamma_old));
0171             in.Q0 = Q0; in.R0 = R0;
0172     end
0173 
0174     in.xh_old = xh_old;
0175     in.pk_old = pk_old;
0176     % in.Atdyu = Atdyu;
0177     in.u = u;
0178     in.W = W;
0179     in.W_old = W_old;
0180     in.delx_mode = delx_mode;
0181     in.debias = 0;
0182     in.verbose = 0;
0183     in.plots = 0;
0184     in.record = 1;
0185     in.err_fun = @(z) (norm(x-z)/norm(x))^2;
0186     tic
0187     out = l1homotopy(A,y,in);
0188     xh_dynDict = out.x_out;
0189     gamma_dynDict = out.gamma;
0190     iter_dynDict = out.iter;
0191     time_dynDict = toc; 
0192      
0193     SIM_stack{sim} = [sim, tau, ...
0194         norm(x_old-xh)^2/norm(x)^2, sum(iter_bpdn,2), sum(time_bpdn,2), ...
0195         norm(x-xh_dynDict)^2/norm(x)^2, sum(iter_dynDict,2), sum(time_dynDict,2)];  
0196     
0197     fprintf('sim %d. tau = %3.4g, (err,iter,time): l1homotopy from scratch-%3.4g,%3.4g,%3.4g; dynDict l1homotopy-%3.4g,%3.4g,%3.4g. \n', ...
0198         SIM_stack{sim});
0199    
0200     %% plot recovered signals
0201     figure(1); plot([x xh xh_dynDict]);
0202     title('Comparison betweeen the new and old homotopy code')
0203 end
0204 mS =  mean(cell2mat(SIM_stack),1);
0205 fprintf('Average results: maxsim %d. tau = %3.4g, (err,iter,time): l1homotopy from scratch-%3.4g,%3.4g,%3.4g; dynDict l1homotopy-%3.4g,%3.4g,%3.4g. \n', maxsim, mS(2:end));