DynamicSeq_DS_demo

 Dantzig selector sequential measurement update homotopy comparison
 Author: Salman Asif
 Created: February 2009:
 Modified: June 2009

-------------------------------------------+
 Copyright (c) 2009.  Muhammad Salman Asif 
-------------------------------------------+

0001 % Dantzig selector sequential measurement update homotopy comparison
0002 % Author: Salman Asif
0003 % Created: February 2009:
0004 % Modified: June 2009
0005 %
0006 %-------------------------------------------+
0007 % Copyright (c) 2009.  Muhammad Salman Asif
0008 %-------------------------------------------+
0009 
0010 clear; clc
0011 
0012 % % load fixed random states
0013 % load RandomStates
0014 % randn('state',s_randn);
0015 % rand('state',s_rand);
0016 
0017 % signal length
0018 N = 256;
0019 % number of spikes to put down
0020 T = 20;
0021 % number of observations to make
0022 K = 100;
0023 
0024 % random +/- 1 signal
0025 x = zeros(N,1);
0026 q = randperm(N);
0027 x(q(1:T)) = randn(T,1);%sign(randn(T,1)).*randint(T,1,256);
0028 
0029 % measurement matrix
0030 G = randn(K,N)/sqrt(K);
0031 % G = orth(G')';
0032 A = G(1:end-1,:);
0033 b = G(end,:);
0034 
0035 x(q(1:T)) = randsrc(T,1);
0036 
0037 % observations
0038 sigma = .01;
0039 e = randn(K-1,1)*sigma;
0040 y = A*x + e;
0041 
0042 % regularization parameter
0043 tau = 0.005*max(abs(A'*y)); %l1_ls
0044 if sigma>0
0045     tau = sigma * sqrt(log(N)*2); % BPDN or DS
0046     %tau = max(abs(A'*(A*x-y))); % ideal ???
0047 end
0048 
0049 % Initial update for the DS
0050 [xp, lame, gamma_x, gamma_lambda, xp_iter, t1] = DS_homotopy_function(A, y, tau, 4*N);
0051 
0052 
0053 Aglgx = A(:,gamma_lambda)'*A(:,gamma_x); 
0054 [Q_glgx R_glgx] = qr(Aglgx);
0055 
0056 pk = A'*(A*xp-y);
0057 ak = A'*(A*lame);
0058 w = b*x+randn*sigma;
0059 
0060 [xp_h, lambda_h, gamma_xh, gamma_lh, xp_h_iter, th] = DynamicSeq_DS_function(A, b, Q_glgx, R_glgx, y, w, xp, lame, gamma_x, gamma_lambda, pk, ak, tau, 4*N);
0061 
0062 y2 = [y; w];
0063 G = [A; b];
0064 [xp2, lame2, gamma_x2, gamma_lambda2, xp2_iter, t2] = DS_homotopy_function(G, y2, tau, 4*N);
0065 
0066 disp(' ');
0067 disp('Results for ');
0068 disp('old_homotopy , new_homotopy , update_homotopy')
0069 cputime_comparison = [t1 t2 th]
0070 Iter_comparison = [xp_iter xp2_iter xp_h_iter]