DynamicX_BPDN_Visual_demo

 Test file for LASSO (BPDN) update with dynamic change in x

0001 % Test file for LASSO (BPDN) update with dynamic change in x
0002 
0003 % Author: Muhammad Salman Asif @ Georgia Tech
0004 % E-mail: sasif@ece.gatech.edu
0005 % Created: November 2008
0006 % Modified: February 2009
0007 %
0008 %-------------------------------------------+
0009 % Copyright (c) 2009.  Muhammad Salman Asif
0010 %-------------------------------------------+
0011 
0012 
0013 clear
0014 
0015 % % load fixed random states
0016 % load RandomStates
0017 % rand('state', rand_state);
0018 % randn('state', randn_state);
0019 
0020 % signal length
0021 N = 128;
0022 % number of spikes to put down
0023 T = 20;
0024 % number of observations to make
0025 M = 64;
0026 
0027 % random signal supported on q(1:T)
0028 x = zeros(N,1);
0029 q = randperm(N);
0030 
0031 %N_S = 1: Gaussian signal, N_S = 0: Random signed sequence
0032 N_S = 0;
0033 if N_S == 1
0034     x(q(1:T)) = randn(T,1);
0035 else
0036     x(q(1:T)) = sign(randn(T,1));
0037 end
0038 
0039 % measurement matrix
0040 %disp('Creating measurment matrix...');
0041 G = randn(M,N)/sqrt(M);
0042 A = (orth(G'))';
0043 % A = G;
0044 
0045 % observations
0046 sigma = .01;
0047 e = randn(M,1)*sigma;
0048 y = A*x + e;
0049 
0050 % regularization parameter
0051 tau = 1e-2;
0052 
0053 % Solve for measurements y
0054 t0 = clock;
0055 [xp, gamma_x, xp_iter, xp_time] = BPDN_homotopy_function(A, y, tau,4*N); %BPDN
0056 t1 = clock;
0057 t0_1 = etime(t1,t0);
0058 
0059 % Model the change in signal
0060 dx_model = 1; % 0 - Perturb only non-zero locations
0061               % 1 - Add some new elements (T_in) and remove some existing
0062               % elements (T_out)
0063 switch dx_model
0064     case  1
0065         % Add some elements (T_in) or remove some existing elements (T_out)
0066         dx  = zeros(N,1);
0067         dx(q(1:T)) = randn*.1;
0068         T_in = round(T/20);
0069         T_out = round(T/20);
0070         qin = randperm(N);
0071 
0072         dx(qin(1:T_in)) = (randn(T_in,1));
0073 
0074         qout = randperm(T);
0075         oldx_indices = q(1:T);
0076         dx(oldx_indices(qout(1:T_out))) = -x(oldx_indices(qout(1:T_out)));
0077 
0078     case  2;
0079         % Perturb the current locations by Gaussian noise
0080         oldx_indices = q(1:T);
0081         dx = zeros(N,1);
0082         dx(oldx_indices) = randn(T,1)*.2;
0083     otherwise
0084         disp('Nooooo');
0085 end
0086 
0087 xt = x+dx;
0088 yt = y+A*dx;
0089 
0090 epsilon = 0;
0091 e0 = 0;
0092 xp2_h = xp;
0093 gamma_xh = gamma_x;
0094 itheta = [];
0095 
0096 done = 0;
0097 iter = 0;
0098 pk_old = A'*(A*xp-y);
0099 AtAgx = A(:,gamma_xh)'*A(:,gamma_xh);
0100 iAtAgx = inv(AtAgx);
0101 
0102 tu0 = clock;
0103 while ~done
0104     iter = iter +1 ;
0105 
0106     %%% NEED UPDATE
0107     % iAtAgx = inv(A(:,gamma_xh)'*A(:,gamma_xh));
0108 
0109     d = A(:,gamma_xh)'*(-y+yt);
0110     delx = iAtAgx*d;
0111     pk = pk_old;
0112     dk = A'*(A(:,gamma_xh)*delx+y-yt);
0113     delx_vec = zeros(N,1);
0114     delx_vec(gamma_xh) = delx;
0115 
0116     temp_gamma = zeros(N,1);
0117     temp_gamma(gamma_xh) = gamma_xh;
0118     gamma_xc = find([1:N]' ~= temp_gamma);
0119     % gamma_xc = setdiff([1:N],[gamma_xh]);
0120 
0121     b_constr1 = (tau-pk(gamma_xc))./dk(gamma_xc);
0122     b_constr2 = (tau+pk(gamma_xc))./-dk(gamma_xc);
0123     b_constr3 = (-xp2_h(gamma_xh)./delx_vec(gamma_xh));
0124     itheta_1 = find(b_constr1>2*eps);
0125     itheta_2 = find(b_constr2>2*eps);
0126     itheta_3 = find(b_constr3>2*eps);
0127     theta1 = min(b_constr1(itheta_1));
0128     theta2 = min(b_constr2(itheta_2));
0129     theta3 = min(b_constr3(itheta_3));
0130     if theta1>theta2
0131         theta = theta2;
0132         itheta = gamma_xc(find(b_constr2==theta2));
0133         flag = 1;
0134     else
0135         theta = theta1;
0136         itheta = gamma_xc(find(b_constr1==theta1));
0137         flag = 1;
0138     end
0139     if theta3 < theta
0140         theta = theta3;
0141         itheta = gamma_xh(find(b_constr3==theta3));
0142         flag = 0;
0143 
0144     end
0145     %     if flag == 0
0146     %         [epsilon+theta theta itheta flag xp2_h(itheta) delx_vec(itheta)]
0147     %     else
0148     %         [epsilon+theta theta itheta flag pk(itheta) dk(itheta)]
0149     %     end
0150     
0151     epsilon = (theta+e0);
0152     if epsilon < 0
0153         epsilon= inf;
0154     end
0155     gamma_xh_old = gamma_xh;
0156     xp2_old = xp2_h;
0157     xp2_h = xp2_h+theta*delx_vec;
0158     if flag == 1
0159         AtgxAnx = A(:,gamma_xh)'*A(:,itheta);
0160         AtAgx_mod = [AtAgx AtgxAnx; AtgxAnx' A(:,itheta)'*A(:,itheta)];
0161 
0162         iAtAgx = update_inverse(AtAgx_mod, iAtAgx,1);
0163         AtAgx = AtAgx_mod;
0164 
0165         gamma_xh = [gamma_xh; itheta];
0166         xp2_h(itheta) = 0;
0167     else
0168         outx_index = find(gamma_xh==itheta);
0169         gamma_xh(outx_index) = gamma_xh(end);
0170         gamma_xh(end) = itheta;
0171         gamma_xh = gamma_xh(1:end-1);
0172 
0173         rowi = outx_index; % ith row of A is swapped with last row (out_x)
0174         colj = outx_index; % jth column of A is swapped with last column (out_x)
0175         AtAgx_ij = AtAgx;
0176         temp_row = AtAgx_ij(rowi,:);
0177         AtAgx_ij(rowi,:) = AtAgx_ij(end,:);
0178         AtAgx_ij(end,:) = temp_row;
0179         temp_col = AtAgx_ij(:,colj);
0180         AtAgx_ij(:,colj) = AtAgx_ij(:,end);
0181         AtAgx_ij(:,end) = temp_col;
0182         iAtAgx_ij = iAtAgx;
0183         temp_row = iAtAgx_ij(colj,:);
0184         iAtAgx_ij(colj,:) = iAtAgx_ij(end,:);
0185         iAtAgx_ij(end,:) = temp_row;
0186         temp_col = iAtAgx_ij(:,rowi);
0187         iAtAgx_ij(:,rowi) = iAtAgx_ij(:,end);
0188         iAtAgx_ij(:,end) = temp_col;
0189 
0190         AtAgx = AtAgx_ij(1:end-1,1:end-1);
0191         iAtAgx = update_inverse(AtAgx_ij, iAtAgx_ij,2);
0192 
0193         xp2_h(itheta) = 0;
0194     end
0195 
0196     %[epsilon theta itheta flag]
0197     %     itheta_stack = [itheta_stack; itheta];
0198     %     flag_stack = [flag_stack; flag];
0199 
0200     if epsilon >=1
0201         theta_end = (1-e0);
0202         xp2_h = xp2_old + theta_end*delx_vec;
0203         gamma_xh = gamma_xh_old;
0204         epsilon = 1;
0205         disp('done!');
0206         break;
0207     end
0208 
0209     e0 = (theta+e0);
0210 
0211     pk_old = pk+theta*dk;
0212     pk_old([gamma_xh; itheta]) = sign(pk_old([gamma_xh; itheta]))*tau;
0213 
0214     figure(2);
0215     subplot(2,1,1);
0216     hold off;
0217     plot(pk,'.r','MarkerSize',14);
0218     hold on;
0219     plot(pk_old, 'LineWidth',1);
0220 
0221     if flag == 1
0222         plot(itheta, pk_old(itheta),'or','MarkerSize',18, 'LineWidth',2);
0223         text(itheta,pk_old(itheta)*1.2, ['Incoming \gamma = ',num2str(itheta)],'FontSize',14);
0224     else
0225         plot(itheta, pk_old(itheta),'*k','MarkerSize',18, 'LineWidth',2);
0226         text(itheta,pk_old(itheta)*1.2, ['Outgoing \gamma = ',num2str(itheta)],'FontSize',14);
0227     end
0228     set(gca,'FontSize',16, 'XLim',[1 N] );
0229     title('BPDN shrinkage constraints');
0230     %legend('Old constraint', 'New constraint', 'Element change','Location','NorthEast');
0231     title(['BPDN constraints, n = ',num2str(N), ', m = ', num2str(M), ', T = ', num2str(T), ', Tin = ', num2str(T_in)]);
0232     plot(1:N, tau*ones(1,N),'--k');
0233     plot(1:N, -tau*ones(1,N), '--k');
0234     axis([1 N -tau*1.5 tau*1.5])
0235 
0236     figure(2);
0237     subplot(2,1,2)
0238     hold off
0239     plot(xp2_old,'.r','MarkerSize',14); hold on;
0240     plot(xp2_h,'LineWidth',1);
0241     set(gca,'FontSize',16,'XLim',[1 N]);
0242     title({['BPDN estimate at \epsilon = ',num2str(epsilon),', Iteration count = ',num2str(iter)]});
0243     if flag == 0
0244         plot(itheta, 0,'ok', 'MarkerSize',18,'LineWidth',2);
0245     else
0246         plot(itheta, 0,'or', 'MarkerSize',18,'LineWidth',2);
0247     end
0248     drawnow;
0249     % pause;
0250 
0251     %print(gcf,'-dbmp','-r200',['BPDNPath_',num2str(iter)])
0252 end
0253 figure(2);
0254 subplot(2,1,1);
0255 hold off;
0256 plot(pk,'.r','MarkerSize',14);
0257 hold on;
0258 plot(pk_old, 'LineWidth',1);
0259 
0260 
0261 set(gca,'FontSize',16, 'XLim',[1 N] );
0262 title('BPDN shrinkage constraints');
0263 %legend('Old constraint', 'New constraint', 'Element change','Location','NorthEast');
0264 title(['BPDN constraints, n = ',num2str(N), ', m = ', num2str(M), ', T = ', num2str(T), ', Tin = ', num2str(T_in)]);
0265 plot(1:N, tau*ones(1,N),'--k');
0266 plot(1:N, -tau*ones(1,N), '--k');
0267 axis([1 N -tau*1.5 tau*1.5])
0268 
0269 figure(2);
0270 subplot(2,1,2)
0271 hold off
0272 plot(xp2_old,'.r','MarkerSize',14); hold on;
0273 plot(xp2_h,'LineWidth',1);
0274 set(gca,'FontSize',16,'XLim',[1 N]);
0275 title({['BPDN estimate at \epsilon = ',num2str(epsilon),', Iteration count = ',num2str(iter)]});
0276 drawnow;