DS_homotopy_demo

 DS_homotopy_demo_test.m

 Solves the following Dantzig selector problem
 min_x  ||x||_1  subject to  ||A'(Ax-y)||_\infty <= epsilon

 using Primal Dual Pursuit Homotopy method.

 Written by: Salman Asif, Georgia Tech
 Email: sasif@gatech.edu

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

0001 % DS_homotopy_demo_test.m
0002 %
0003 % Solves the following Dantzig selector problem
0004 % min_x  ||x||_1  subject to  ||A'(Ax-y)||_\infty <= epsilon
0005 %
0006 % using Primal Dual Pursuit Homotopy method.
0007 %
0008 % Written by: Salman Asif, Georgia Tech
0009 % Email: sasif@gatech.edu
0010 %
0011 %-------------------------------------------+
0012 % Copyright (c) 2008.  Muhammad Salman Asif
0013 %-------------------------------------------+
0014 
0015 clear; close all; 
0016 
0017 rseed = 0;
0018 rand('state',rseed);
0019 randn('state',rseed);
0020 
0021 N = 512;   % signal length
0022 T = 20;    % sparsity level
0023 M = 250;    % no. of measurements
0024 
0025 % Generate a random signal
0026 x = zeros(N,1);
0027 q = randperm(N);
0028 x(q(1:T)) = sign(randn(T,1));
0029 
0030 % measurement matrix
0031 A = randn(M,N)/sqrt(M); % Gaussian
0032 
0033 % Hadamard
0034 % H = hadamard(N);
0035 % A = H(q(1:M),:)/sqrt(M);
0036 
0037 % Bernoulli
0038 % A = randsrc(M,N)/sqrt(M);
0039 
0040 % Random Projection
0041 % A = orth(A')';
0042 
0043 % AtA = A'*A;
0044 
0045 % measurements
0046 SNR = inf;
0047 s = A*x;
0048 sigma = sqrt((norm(s)^2/(10^(SNR/10)))/M);
0049 % sigma = 0.01;
0050 e = randn(M,1)*sigma;
0051 y = A*x+e;
0052 thresh = 1e-3; % sqrt(2*log(N))*sigma;
0053 
0054 % Initialization of primal and dual sign and support
0055 z_x = zeros(N,1);
0056 z_lambda = zeros(N,1);
0057 gamma_lambda = [];  % Dual support
0058 gamma_x = [];       % Primal support
0059 gamma_lambda_app = [];
0060 
0061 % Initial step
0062 Primal_constrk = -A'*y;
0063 [c i] = max(abs(Primal_constrk));
0064 
0065 gamma_lambdak = i;
0066 gamma_xk = gamma_lambdak;
0067 
0068 z_lambda(gamma_lambdak) = sign(Primal_constrk(gamma_lambdak));
0069 epsilon = c;
0070 Primal_constrk(gamma_lambdak) = sign(Primal_constrk(gamma_lambdak))*epsilon;
0071 xk_1 = zeros(N,1);
0072 
0073 lambdak_1 = zeros(N,1);
0074 lambdak_1(gamma_lambdak) = inv(A(:,gamma_lambdak)'*A(:,gamma_lambdak))*z_lambda(gamma_lambdak);
0075 
0076 Dual_constrk = A'*(A*lambdak_1);
0077 
0078 z_x(gamma_xk) = -sign(Dual_constrk(gamma_xk));
0079 Dual_constrk(gamma_xk) = sign(Dual_constrk(gamma_xk));
0080 
0081 z_xk = z_x;
0082 z_lambdak = z_lambda;
0083 
0084 % loop parameters
0085 done = 0;
0086 iter = 0;
0087 data_precision = eps;   % floating point precision
0088 
0089 old_delta = 0;
0090 out_lambda = [];
0091 out_x = [];
0092 count_delta_stop = 0;
0093 
0094 constraint_plots = 0;
0095 
0096 AtglAgx = A(:,gamma_lambdak)'*A(:,gamma_xk);
0097 iAtglAgx = inv(A(:,gamma_lambdak)'*A(:,gamma_xk));
0098 AtgxAgl = AtglAgx';
0099 iAtgxAgl = iAtglAgx';
0100 time_s = clock;
0101 
0102 itr_history = [iter epsilon 0 0 length(gamma_xk)]
0103 
0104 while ~done
0105     iter = iter+1;
0106     warning('off','MATLAB:divideByZero')
0107 
0108     %%%% IDEA implementation %%%%
0109     gamma_x = gamma_xk;
0110     gamma_lambda = gamma_lambdak;
0111     z_lambda = z_lambdak;
0112     z_x = z_xk;
0113     x_k = xk_1;
0114     lambda_k = lambdak_1;
0115 
0116     %%%%%%%%%%%%%%%%%%%%%
0117     %%%% update on x %%%%
0118     %%%%%%%%%%%%%%%%%%%%%
0119 
0120     % UPDATE INVERSE MATRIX
0121 
0122     % Update direction
0123     %del_x = -inv(A(:,gamma_lambda)'*A(:,gamma_x))*z_lambda(gamma_lambda);
0124     del_x = -iAtglAgx*z_lambda(gamma_lambda);
0125     del_x_vec = zeros(N,1);
0126     del_x_vec(gamma_x) = del_x;
0127 
0128     pk = Primal_constrk;
0129     % dk = AtA*del_x_vec;
0130     dk = A'*(A(:,gamma_x)*del_x);
0131     
0132     %%% CONTROL THE MACHINE PRECISION ERROR AT EVERY OPERATION: LIKE BELOW.
0133     pk_temp = Primal_constrk;
0134     % gammaL_temp = find(abs(abs(Primal_constrk)-epsilon)<min(epsilon,1e-12));
0135     % pk_temp(gammaL_temp) = sign(Primal_constrk(gammaL_temp))*epsilon;
0136     
0137     xk_temp = x_k;
0138     % gammaX_temp = find(abs(x_k)<2*eps);
0139     % xk_temp(gammaX_temp) = 0;
0140     %%%---
0141     
0142     % Step size computation
0143     [i_delta, out_x, delta, chk_x] = update_primal(gamma_x, gamma_lambda, z_x,  xk_temp, del_x_vec, pk_temp, dk, epsilon, out_lambda);
0144 
0145     if old_delta < 4*eps & delta < 4*eps
0146         count_delta_stop = count_delta_stop + 1;
0147     else
0148         count_delta_stop = 0;
0149     end
0150     if count_delta_stop >= 50
0151         disp('stuck in some corner');
0152         break;
0153     end
0154     old_delta = delta;
0155 
0156     xk_1 = x_k+delta*del_x_vec;
0157     Primal_constrk = pk+delta*dk;
0158     epsilon_old = epsilon;
0159     epsilon = epsilon-delta;
0160 
0161     if chk_x == 1
0162         % If an element is removed from gamma_x
0163         out_x = out_x(1);
0164         gl_old = gamma_lambda;
0165         gx_old = gamma_x;
0166 
0167         outx_index = find(gamma_x==out_x);
0168         gamma_x(outx_index) = gamma_x(end);
0169         gamma_x(end) = out_x;
0170         gamma_x = gamma_x(1:end-1);
0171 
0172         %CHECK THIS SINGULARITY CONDITION USING SCHUR COMPLEMENT IDEA !!!
0173         % X = [A B; C D];
0174         % detX = detA detS
0175         % S = D-C A^{-1} B
0176 
0177         % So the idea here is slightly different, we don't want the last
0178         % entry (Q22) in the inverse of modified Gram matrix to be zero, because
0179         % that will make inversion impossible.
0180 
0181         % This is the best thing, pick the index with maximum absolute
0182         % value in the column corresponding to out_x in iAtgxgl. This will
0183         % ensure that new matrix will be invertible.
0184         [max_val outl_index] = max(abs(iAtgxAgl(:,outx_index)));
0185         new_lambda = gamma_lambda(outl_index);
0186         some_lambda = new_lambda;
0187 
0188         outl_index = find(gamma_lambda==some_lambda);
0189         gamma_lambda(outl_index) = gamma_lambda(end);
0190         gamma_lambda(end) = some_lambda;
0191         gamma_lambdak = gamma_lambda;
0192         gamma_lambda = gamma_lambda(1:end-1);
0193 
0194         rowi = outx_index; % ith row of A is swapped with last row (out_x)
0195         colj = outl_index; % jth column of A is swapped with last column (out_lambda)
0196         AtgxAgl_ij = AtgxAgl;
0197         temp_row = AtgxAgl_ij(rowi,:);
0198         AtgxAgl_ij(rowi,:) = AtgxAgl_ij(end,:);
0199         AtgxAgl_ij(end,:) = temp_row;
0200         temp_col = AtgxAgl_ij(:,colj);
0201         AtgxAgl_ij(:,colj) = AtgxAgl_ij(:,end);
0202         AtgxAgl_ij(:,end) = temp_col;
0203         iAtgxAgl_ij = iAtgxAgl;
0204         temp_row = iAtgxAgl_ij(colj,:);
0205         iAtgxAgl_ij(colj,:) = iAtgxAgl_ij(end,:);
0206         iAtgxAgl_ij(end,:) = temp_row;
0207         temp_col = iAtgxAgl_ij(:,rowi);
0208         iAtgxAgl_ij(:,rowi) = iAtgxAgl_ij(:,end);
0209         iAtgxAgl_ij(:,end) = temp_col;
0210 
0211         AtgxAgl = AtgxAgl_ij(1:end-1,1:end-1);
0212         AtglAgx = AtgxAgl';
0213         iAtgxAgl = update_inverse(AtgxAgl_ij, iAtgxAgl_ij,2);
0214         iAtglAgx = iAtgxAgl';
0215         xk_1(out_x) = 0;
0216     else
0217         % If a new element is added to gamma_lambda
0218         gamma_lambdak = [gamma_lambda; i_delta];
0219         new_lambda = i_delta;
0220         lambda_k(new_lambda) = 0;
0221     end
0222 
0223     out_lambda = [];
0224     z_lambdak = zeros(N,1);
0225     z_lambdak(gamma_lambdak) = sign(Primal_constrk(gamma_lambdak));
0226     sgn_new_lambda = sign(Primal_constrk(new_lambda));
0227     Primal_constrk(gamma_lambdak) = sign(Primal_constrk(gamma_lambdak))*epsilon;
0228 
0229     if epsilon < thresh; %sqrt(2*log(N))*sigma; %1e-7 %|| iter > 5*T || (length(gamma_lambda) == K)
0230         delta_end = epsilon_old-thresh;
0231         Primal_constrk = pk+delta_end*dk;
0232         xk_1 = x_k + delta_end*del_x_vec;
0233         disp('done!');
0234         break;
0235     end
0236 
0237     if iter > 500*T
0238         disp('too many iterations ooooooooops');
0239     end
0240 
0241     %%%%%%%%%%%%%%%%%%%%%%%%
0242     %%% update on lambda %%%
0243     %%%%%%%%%%%%%%%%%%%%%%%%
0244 
0245     sgn_new_lambda = sign(Primal_constrk(new_lambda));
0246 
0247     AtgxAnl = A(:,gamma_x)'*A(:,new_lambda);
0248     %AtgxAnl = AtA(gamma_x,new_lambda);
0249 
0250     % Update direction
0251     %del_lambda = -inv(A(:,gamma_x)'*A(:,gamma_lambda))*(A(:,gamma_x)'*A(:,new_lambda));
0252     del_lambda = -iAtgxAgl*AtgxAnl;
0253     del_lambda_p = zeros(N,1);
0254     del_lambda_p(gamma_lambda) = del_lambda*sgn_new_lambda;
0255     del_lambda_p(new_lambda) = 1*sgn_new_lambda;
0256 
0257     ak = Dual_constrk;
0258     % bk = AtA*del_lambda_p;
0259     bk = A'*(A(:,[gamma_lambda; new_lambda])*del_lambda_p([gamma_lambda; new_lambda]));
0260 
0261     % check if the sign of update direction is correct
0262     if sign(bk(out_x)) == sign(ak(out_x)) & abs(bk(out_x)) >=eps
0263         bk = -bk;
0264         del_lambda_p = -del_lambda_p;
0265     end
0266 
0267     %%% CONTROL THE MACHINE PRECISION ERROR AT EVERY OPERATION: LIKE BELOW.
0268     ak_temp = Dual_constrk;
0269     % gammaX_temp = find(abs(abs(Dual_constrk)-1)<2*eps);
0270     % ak_temp(gammaX_temp) = sign(Dual_constrk(gammaX_temp));
0271     
0272     lambdak_temp = lambda_k;
0273     % gammaL_temp = find(abs(lambda_k)<1*eps);
0274     % lambdak_temp(gammaL_temp) = 0;
0275     %%%---
0276     
0277     % Step size computation
0278     [i_theta, out_lambda, theta, chk_lambda] = update_dual(gamma_x, gamma_lambda, z_lambda, lambdak_temp, del_lambda_p, ak_temp, bk, new_lambda, out_x);
0279 
0280     out_x = [];
0281 
0282     lambdak_1 = lambda_k + theta*del_lambda_p;
0283     Dual_constrk = ak+theta*bk;
0284 
0285     if chk_lambda == 1
0286         % If an element is removed from gamma_lambda
0287         outl_index = find(gamma_lambda==out_lambda);
0288         i_theta = gamma_x(end);
0289 
0290         if ~isempty(outl_index) % because otherwise new_lambda is removed.
0291             % Use Sherman-Woodburry-Morrison formula
0292             % (A+BC)^{-1} = A^{-1} - A^{-1}B(I+CA^{-1}B)^{-1}CA^{-1};
0293             iA = iAtgxAgl;
0294             B = AtgxAnl-A(:,gamma_x)'*A(:,out_lambda);
0295             C = zeros(1,length(gamma_lambda));
0296             C(outl_index) = 1;
0297             AtgxAgl(:,outl_index) = AtgxAnl;
0298             AtglAgx = AtgxAgl';
0299 
0300             iAB = iA*B;
0301             CiA = iA(outl_index,:);
0302             iAtgxAgl = iA-iAB*((CiA)/(1+iAB(outl_index)));
0303             iAtglAgx = iAtgxAgl';
0304 
0305             gamma_lambda(outl_index) = new_lambda;
0306         end
0307 
0308         z_xk = zeros(N,1);
0309         z_xk(gamma_x) = -sign(Dual_constrk(gamma_x));
0310         lambdak_1(out_lambda) = 0;
0311     else
0312         % If an element is added to gamma_x
0313 
0314         %%% adding columns
0315         % AtglAgx_mod = zeros(length(gamma_x));
0316         AtglAgx_mod = [AtglAgx A(:,gamma_lambda)'*A(:,i_theta); AtgxAnl' A(:,new_lambda)'*A(:,i_theta)];
0317         %AtglAgx_mod = [AtglAgx AtA(gamma_lambda,i_theta); AtgxAnl' AtA(new_lambda,i_theta)];
0318 
0319         %CHECK THIS SINGULARITY CONDITION USING SCHUR COMPLEMENT IDEA !!!
0320         % X = [A B; C D];
0321         % detX = detA detS
0322         % S = D-C A^{-1} B
0323         A11 = AtglAgx;
0324         A12 = AtglAgx_mod(1:end-1,end);
0325         A21 = AtglAgx_mod(end,1:end-1);
0326         A22 = AtglAgx_mod(end,end);
0327         S = A22 - A21*(iAtglAgx*A12);
0328         if abs(S) < 2*eps
0329             disp('matrix has become singular');
0330             % done = 1;
0331             % break;
0332         end
0333         AtglAgx = AtglAgx_mod;
0334         AtgxAgl = AtglAgx';
0335         iAtglAgx = update_inverse(AtglAgx, iAtglAgx,1);
0336         iAtgxAgl = iAtglAgx';
0337 
0338         out_lambda = [];
0339         gamma_lambda = [gamma_lambda; new_lambda];
0340 
0341         gamma_x = [gamma_x; i_theta];
0342         new_x = i_theta;
0343         z_xk = zeros(N,1);
0344         z_xk(gamma_x) = -sign(Dual_constrk(gamma_x));
0345         xk_1(new_x) = 0;
0346     end
0347     Dual_constrk(gamma_x) = sign(Dual_constrk(gamma_x));
0348     gamma_lambdak = gamma_lambda;
0349     gamma_xk = gamma_x;
0350    
0351     if constraint_plots
0352         fig = figure(1);
0353         subplot(211)
0354         hold off;
0355         plot(pk,'.r', 'MarkerSize',14);        
0356         hold on;
0357         plot(Primal_constrk, 'LineWidth',1);
0358         set(gca,'FontSize',16, 'XLim',[1 N] );
0359         title({'BPDN shrinkage constraints,'; ['n = ',num2str(N), ', m = ', num2str(M), ', T = ',num2str(T)]});
0360         plot(1:N, epsilon*ones(1,N),'--k','MarkerSize',12);
0361         plot(1:N, -epsilon*ones(1,N), '--k','MarkerSize',12);
0362         plot(1:N, epsilon_old*ones(1,N),'--m','MarkerSize',12);
0363         plot(1:N, -epsilon_old*ones(1,N), '--m','MarkerSize',12);
0364         
0365         
0366         subplot(212);
0367         hold off;
0368         plot(ak,'.r');
0369         hold on;
0370         plot(Dual_constrk);
0371         %legend('Old constraint', 'New constraint', 'Element change','Location','NorthEast');
0372         title('Dual constraints');
0373         plot(1:N, 1, '--k');
0374         plot(1:N, -1, '--k');
0375         axis([1 N  -1.5 1.5])
0376         
0377         if chk_x == 0
0378             subplot(211);
0379             plot(new_lambda, Primal_constrk(new_lambda),'or','MarkerSize',18,'LineWidth',2);
0380             text(new_lambda, Primal_constrk(new_lambda)*1.1, ['Incoming \gamma = ',num2str(new_lambda)],'FontSize',14);
0381         else
0382             subplot(212)
0383             plot(out_x, (ak(out_x)),'*k','MarkerSize',18,'LineWidth',2);
0384             text(out_x,ak(out_x)*1.1, ['Outgoing \gamma = ',num2str(out_x)],'FontSize',14);
0385         end
0386         
0387         if chk_lambda == 0
0388             subplot(212)
0389             plot(new_x, Dual_constrk(new_x),'or','MarkerSize',18,'LineWidth',2);
0390             text(new_x, Dual_constrk(new_x)*1.1, ['Incoming \gamma = ',num2str(new_x)],'FontSize',14);
0391         else
0392             subplot(211);
0393             plot(out_lambda, (Primal_constrk(out_lambda)),'*k','MarkerSize',18,'LineWidth',2);
0394             text(out_lambda,Primal_constrk(out_lambda)*1.1, ['Outgoing \gamma = ',num2str(out_x)],'FontSize',14);
0395         end
0396     end
0397     itr_history_temp = [iter epsilon delta theta length(gamma_xk)]
0398     itr_history = [itr_history; itr_history_temp];
0399 end
0400 time_e = clock;
0401 etime(time_e,time_s);
0402 
0403 figure(1); clf; 
0404 subplot(211); plot(pk,'.r'); hold on;
0405 plot(Primal_constrk); 
0406 plot(1:N, epsilon,'k'); plot(1:N, -epsilon, 'k'); plot(1:N, epsilon_old,'g'); plot(1:N, -epsilon_old, 'g');
0407 title('Primal constraints');
0408 
0409 subplot(212); plot(ak,'.r'); hold on; 
0410 plot(Dual_constrk); 
0411 plot(1:N, 1, '-k'); plot(1:N, -1, '-k'); axis([1 N  -1.5 1.5])
0412 title('Dual constraints');
0413 
0414 figure(2); clf; 
0415 plot(xk_1); hold on; plot(x,'.r');
0416 legend('estimation', 'original');
0417 title('Estimated signal')
0418 
0419 % %%
0420 % x0 = A'*inv(A*A')*y;
0421 % % large scale
0422 % Afun = @(z) A*z;
0423 % Atfun = @(z) A'*z;
0424 % tic
0425 % xp = l1eq_pd(x0, Afun, Atfun, y, 1e-8, 1000, 1e-8, 1000);
0426 % toc