DynamicX_DS_demo

 Test file for Dantzig selector (DS) update with dynamic change in x

0001 % Test file for Dantzig selector (DS) update with dynamic change in x
0002 
0003 % Author: Muhammad Salman Asif @ Georgia Tech
0004 % E-mail: sasif@ece.gatech.edu
0005 % Created: February 2009
0006 %
0007 %-------------------------------------------+
0008 % Copyright (c) 2009.  Muhammad Salman Asif
0009 %-------------------------------------------+
0010 
0011 clear; clc;
0012 
0013 addpath('../Pursuits_Homotopy')
0014 
0015 % % load fixed random states
0016 % load RandomStates
0017 % rand('state', s_rand);
0018 % randn('state', s_randn);
0019 
0020 % signal length
0021 N = 512;
0022 % number of spikes to put down
0023 T = 20;
0024 % number of observations to make
0025 M = 256;
0026 
0027 % random signal supported on q(1:T)
0028 x = zeros(N,1);
0029 q = randperm(N);
0030 x(q(1:T)) = sign(randn(T,1));
0031 
0032 
0033 % measurement matrix
0034 A = randn(M,N)/sqrt(M);
0035 
0036 % observations
0037 sigma = .01;
0038 e = randn(M,1)*sigma;
0039 y = A*x + e;
0040 
0041 % regularization parameter
0042 tau = 0.05*max(abs(A'*y)); %l1_ls
0043 if sigma>0
0044     tau = sigma * sqrt(log(N)*2); % BPDN
0045     %tau = max(abs(A'*(A*x-y))); % ideal ???
0046 end
0047 
0048 % Solve for measurements y
0049 [xp, lame, gamma_x, gamma_lambda, xp_iter, t1] = DS_homotopy_function(A, y, tau,4*N); %Lasso
0050 
0051 % Model the change in signal
0052 dx_model = 1; % 0 - Perturb only non-zero locations
0053               % 1 - Add some new elements (T_in) and remove some existing
0054               % elements (T_out)
0055 switch dx_model
0056     case  1
0057         % Add some elements (T_in) or remove some existing elements (T_out)
0058         dx  = zeros(N,1);
0059         dx(q(1:T)) = randn*.1;
0060         T_in = round(T/20);
0061         T_out = round(T/20);
0062         qin = randperm(N);
0063 
0064         dx(qin(1:T_in)) = sign(randn(T_in,1));
0065 
0066         qout = randperm(T);
0067         oldx_indices = q(1:T);
0068         dx(oldx_indices(qout(1:T_out))) = -x(oldx_indices(qout(1:T_out)));
0069 
0070     case  2;
0071         % Perturb the current locations by Gaussian noise
0072         oldx_indices = q(1:T);
0073         dx = zeros(N,1);
0074         dx(oldx_indices) = randn(T,1)*.2;
0075     otherwise
0076         disp('Nooooo');
0077 end
0078 
0079 xt = x+dx;
0080 yt = y+A*dx;
0081 
0082 pk_old = A'*(A*xp-y);
0083 ak_old = A'*(A*lame);
0084 
0085 Aglgx = A(:,gamma_lambda)'*A(:,gamma_x); 
0086 [Q_glgx R_glgx] = qr(Aglgx);
0087 
0088 [xp_h, lambda_h, gamma_xh, gamma_lh, xp_h_iter, th] = DynamicX_DS_function(A, Q_glgx, R_glgx, y, yt, xp, lame, gamma_x, gamma_lambda, pk_old, ak_old, tau, 4*N);
0089 
0090 [xp2, lame2, gamma_x2, gamma_lambda2, xp2_iter, t2] = DS_homotopy_function(A, yt, tau, 4*N);
0091 
0092 disp(' ');
0093 disp('old_homotopy , new_homotopy , update_homotopy')
0094 cputime_comparison = [t1 t2 th]
0095 iter_comparison = [xp_iter xp2_iter xp_h_iter]
0096 
0097 % figure(2); plot(xp_h-xp2); shg