BPDN_homotopy_demo

 BPDN_homotopy_demo.m
 
 Solves the following basis pursuit denoising problem
 min_x  \epsilon ||x||_1 + 1/2*||y-Ax||_2^2

 This script is a modified version of PD_pursuit method which
 uses only Primal update in PD_pursuit method.

 Here support of primal and dual vectors remain same at every step.
 The update direction here can be considered as equivalent to dual in Dantzig
 selector.

0001 % BPDN_homotopy_demo.m
0002 %
0003 % Solves the following basis pursuit denoising problem
0004 % min_x  \epsilon ||x||_1 + 1/2*||y-Ax||_2^2
0005 %
0006 % This script is a modified version of PD_pursuit method which
0007 % uses only Primal update in PD_pursuit method.
0008 %
0009 % Here support of primal and dual vectors remain same at every step.
0010 % The update direction here can be considered as equivalent to dual in Dantzig
0011 % selector.
0012 
0013 % In S-step solution property (for noiseless case) this becomes more clear when,
0014 % if primal update direction and dual vector match at every step
0015 % (maybe opposite in sign, depending on the formulation) then
0016 % homotopy path taken by BPDN and Dantzig selector is exactly same.
0017 %
0018 % Written by: Salman Asif, Georgia Tech
0019 % Email: sasif@gatech.edu
0020 %
0021 %-------------------------------------------+
0022 % Copyright (c) 2008.  Muhammad Salman Asif
0023 %-------------------------------------------+
0024 
0025 clear; close all
0026 
0027 rseed = 0;
0028 rand('state',rseed);
0029 randn('state',rseed);
0030 
0031 N = 512;   % signal length
0032 T = 50;    % sparsity level
0033 M = 250;    % no. of measurements
0034 
0035 % Generate a random signal
0036 x = zeros(N,1);
0037 q = randperm(N);
0038 x(q(1:T)) = sign(randn(T,1));
0039 
0040 % measurement matrix
0041 A = randn(M,N)/sqrt(M); % Gaussian
0042 
0043 % Hadamard
0044 % H = hadamard(N);
0045 % A = H(q(1:M),:)/sqrt(M);
0046 
0047 % Bernoulli
0048 % A = randsrc(M,N)/sqrt(M);
0049 
0050 % Random Projection
0051 % A = orth(A')';
0052 
0053 % AtA = A'*A;
0054 
0055 % measurements
0056 SNR = inf;
0057 s = A*x;
0058 sigma = sqrt((norm(s)^2/(10^(SNR/10)))/M);
0059 % sigma = 0.01;
0060 e = randn(M,1)*sigma;
0061 y = A*x+e;
0062 thresh = 1e-3; %sqrt(2*log(N))*sigma;
0063 
0064 % Initialization of primal and dual sign and support
0065 z_x = zeros(N,1);
0066 z_lambda = zeros(N,1);
0067 gamma_lambda = [];  % Dual support
0068 gamma_x = [];       % Primal support
0069 
0070 % Initial step
0071 Primal_constrk = -A'*y; 
0072 [c i] = max(abs(Primal_constrk));
0073 
0074 gamma_lambdak = i;
0075 gamma_xk = gamma_lambdak;
0076 
0077 z_lambda(gamma_lambdak) = sign(Primal_constrk(gamma_lambdak));
0078 epsilon = c;
0079 Primal_constrk(gamma_lambdak) = sign(Primal_constrk(gamma_lambdak))*epsilon;
0080 xk_1 = zeros(N,1);
0081 
0082 lambdak_1 = zeros(N,1);
0083 lambdak_1(gamma_lambdak) = inv(A(:,gamma_lambdak)'*A(:,gamma_lambdak))*z_lambda(gamma_lambdak);
0084 
0085 Dual_constrk = A'*(A*lambdak_1);
0086 
0087 z_x(gamma_xk) = -sign(Dual_constrk(gamma_xk));
0088 Dual_constrk(gamma_xk) = sign(Dual_constrk(gamma_xk));
0089 
0090 z_xk = z_x;
0091 z_lambdak = z_lambda;
0092 
0093 % loop parameters
0094 done = 0;
0095 iter = 0;
0096 data_precision = eps;   % floating point precision
0097 
0098 old_delta = 0;
0099 out_lambda = [];
0100 out_x = [];
0101 count_delta_stop = 0;
0102 
0103 % Turn the plots on and off
0104 constraint_plots = 1; % 1 - draw the plots after every homotopy iteration.
0105 
0106 AtglAgx = A(:,gamma_lambdak)'*A(:,gamma_xk);
0107 iAtglAgx = inv(A(:,gamma_lambdak)'*A(:,gamma_xk));
0108 AtgxAgl = AtglAgx';
0109 iAtgxAgl = iAtglAgx';
0110 
0111 while ~done
0112     iter = iter+1;
0113     % warning('off','MATLAB:divideByZero')
0114 
0115     gamma_x = gamma_xk;
0116     gamma_lambda = gamma_lambdak;
0117     z_lambda = z_lambdak;
0118     z_x = z_xk;
0119     x_k = xk_1;
0120     lambda_k = lambdak_1;
0121 
0122     %%%%%%%%%%%%%%%%%%%%%
0123     %%%% update on x %%%%
0124     %%%%%%%%%%%%%%%%%%%%%
0125 
0126     % Update direction
0127 %     del_x = -inv(A(:,gamma_lambda)'*A(:,gamma_x))*z_lambda(gamma_lambda);
0128 %     diff_inv = [ max(max(abs(inv(A(:,gamma_lambda)'*A(:,gamma_x))-iAtglAgx)))]
0129 %     if diff_inv > 1e-8
0130 %         stp =1;
0131 %     end
0132     del_x = -iAtglAgx*z_lambda(gamma_lambda);
0133     del_x_vec = zeros(N,1);
0134     del_x_vec(gamma_x) = del_x;
0135 
0136     pk = Primal_constrk;
0137     % dk = AtA*del_x_vec;
0138     dk = A'*(A(:,gamma_x)*del_x);
0139     
0140     %%% CONTROL THE MACHINE PRECISION ERROR AT EVERY OPERATION: LIKE BELOW.
0141     pk_temp = Primal_constrk;
0142     gammaL_temp = find(abs(abs(Primal_constrk)-epsilon)<min(epsilon,1e-12));
0143     pk_temp(gammaL_temp) = sign(Primal_constrk(gammaL_temp))*epsilon;
0144     
0145     xk_temp = x_k;
0146     gammaX_temp = find(abs(x_k)<1*eps);
0147     xk_temp(gammaX_temp) = 0;
0148     %%%---
0149     
0150     % Compute the step size
0151     [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);
0152     
0153     if old_delta < 4*eps & delta < 4*eps
0154         count_delta_stop = count_delta_stop + 1;
0155     else
0156         count_delta_stop = 0;
0157     end
0158     if count_delta_stop >= 50
0159         disp('stuck in some corner');
0160         break;
0161     end
0162     old_delta = delta;
0163 
0164     xk_1 = x_k+delta*del_x_vec;
0165     Primal_constrk = pk+delta*dk;
0166     epsilon_old = epsilon;
0167     epsilon = epsilon-delta;
0168     
0169     if epsilon < thresh; %sqrt(2*log(N))*sigma; %1e-7 %|| iter > 5*T || (length(gamma_lambda) == K)
0170         delta_end = epsilon_old-thresh;
0171         Primal_constrk = pk+delta_end*dk;
0172         xk_1 = x_k + delta_end*del_x_vec;
0173         disp('done!');
0174         break;
0175     end
0176     
0177     if length(gamma_x)==M & chk_x ==0
0178         stp =1;
0179         disp('Cannot do it Sire'); % Commondos: are you out of your mind Sire!
0180         break;
0181     end
0182 
0183     if chk_x == 1
0184         % If an element is removed from gamma_x
0185         gl_old = gamma_lambda;
0186         gx_old = gamma_x;
0187 
0188         outx_index = find(gamma_x==out_x);
0189         gamma_x(outx_index) = gamma_x(end);
0190         gamma_x(end) = out_x;
0191         gamma_x = gamma_x(1:end-1);
0192 
0193         outl_index = outx_index;
0194         gamma_lambda = gamma_x;
0195         gamma_lambdak = gamma_lambda;
0196         gamma_xk = gamma_x;
0197 
0198         rowi = outx_index; % ith row of A is swapped with last row (out_x)
0199         colj = outl_index; % jth column of A is swapped with last column (out_lambda)
0200         AtgxAgl_ij = AtgxAgl;
0201         temp_row = AtgxAgl_ij(rowi,:);
0202         AtgxAgl_ij(rowi,:) = AtgxAgl_ij(end,:);
0203         AtgxAgl_ij(end,:) = temp_row;
0204         temp_col = AtgxAgl_ij(:,colj);
0205         AtgxAgl_ij(:,colj) = AtgxAgl_ij(:,end);
0206         AtgxAgl_ij(:,end) = temp_col;
0207         iAtgxAgl_ij = iAtgxAgl;
0208         temp_row = iAtgxAgl_ij(colj,:);
0209         iAtgxAgl_ij(colj,:) = iAtgxAgl_ij(end,:);
0210         iAtgxAgl_ij(end,:) = temp_row;
0211         temp_col = iAtgxAgl_ij(:,rowi);
0212         iAtgxAgl_ij(:,rowi) = iAtgxAgl_ij(:,end);
0213         iAtgxAgl_ij(:,end) = temp_col;
0214 
0215         AtgxAgl = AtgxAgl_ij(1:end-1,1:end-1);
0216         AtglAgx = AtgxAgl;
0217         iAtgxAgl = update_inverse(AtgxAgl_ij, iAtgxAgl_ij,2);
0218         iAtglAgx = iAtgxAgl;
0219         xk_1(out_x) = 0;
0220     else
0221         % If an element is added to gamma_x
0222         gamma_lambdak = [gamma_lambda; i_delta];
0223         gamma_xk = gamma_lambdak;
0224         i_theta = i_delta;
0225         new_lambda = i_delta;
0226         AtgxAnl = A(:,gamma_x)'*A(:,new_lambda);
0227         AtglAgx_mod = [AtglAgx A(:,gamma_lambda)'*A(:,i_theta); AtgxAnl' A(:,new_lambda)'*A(:,i_theta)];
0228 
0229         AtglAgx = AtglAgx_mod;
0230         AtgxAgl = AtglAgx;
0231         iAtglAgx = update_inverse(AtglAgx, iAtglAgx,1);
0232         iAtgxAgl = iAtglAgx;
0233         xk_1(i_theta) = 0;
0234     end
0235 
0236     z_lambdak = zeros(N,1);
0237     z_lambdak(gamma_lambdak) = sign(Primal_constrk(gamma_lambdak));
0238     z_xk = -z_lambdak;
0239     Primal_constrk(gamma_lambdak) = sign(Primal_constrk(gamma_lambdak))*epsilon;
0240     
0241     if iter > 500*T
0242         disp('too many iterations ooooooooops');
0243         break;
0244     end
0245 
0246     if constraint_plots
0247         fig = figure(1);
0248         subplot(2,1,1)
0249         hold off
0250         plot(pk,'.r', 'MarkerSize',14);
0251         hold on;
0252         plot(Primal_constrk, 'LineWidth',1);
0253 
0254         if chk_x == 0
0255             plot(new_lambda, Primal_constrk(new_lambda),'or','MarkerSize',18,'LineWidth',2);
0256             text(new_lambda, Primal_constrk(new_lambda)*1.1, ['Incoming \gamma = ',num2str(new_lambda)],'FontSize',14);
0257         else
0258             plot(out_x, Primal_constrk(out_x),'*k','MarkerSize',18,'LineWidth',2);
0259             text(out_x,Primal_constrk(out_x)*1.1, ['Outgoing \gamma = ',num2str(out_x)],'FontSize',14);
0260         end
0261         set(gca,'FontSize',16, 'XLim',[1 N] );
0262         title({'BPDN shrinkage constraints,'; ['n = ',num2str(N), ', m = ', num2str(M), ', T = ',num2str(T)]});
0263         plot(1:N, epsilon*ones(1,N),'--k','MarkerSize',12);
0264         plot(1:N, -epsilon*ones(1,N), '--k','MarkerSize',12);
0265         plot(1:N, epsilon_old*ones(1,N),'--m','MarkerSize',12);
0266         plot(1:N, -epsilon_old*ones(1,N), '--m','MarkerSize',12);
0267 
0268         subplot(2,1,2)
0269         hold off
0270         plot(x_k,'.r','MarkerSize',14); hold on;
0271         plot(xk_1,'LineWidth',1);
0272         if chk_x == 1
0273             plot(out_x, 0,'ok', 'MarkerSize',18,'LineWidth',2);
0274         else
0275             plot(new_lambda, 0,'or', 'MarkerSize',18,'LineWidth',2);
0276         end
0277         set(gca,'FontSize',16,'XLim',[1 N]);
0278         title(['BPDN estimate at \tau = ',num2str(epsilon), ', iter = ', num2str(iter)]);
0279         
0280         if iter == 1
0281             disp('  ');
0282             disp('Every frame in the figure corresponds to a critical point on the homotopy path.')
0283             disp('Circle represents an incoming element, star represents an outgoing element.');
0284             disp(' ');
0285             disp('Put pause somewhere in the code to see this. ');
0286             disp('For now press some key to continue...');
0287             pause
0288         end
0289 
0290         %drawnow;
0291         %print(gcf,'-dbmp','-r200',['BPDNPath_',num2str(iter)])
0292     end
0293     % pause
0294     [iter epsilon delta]
0295 end
0296 
0297 fig = figure(1);
0298 subplot(2,1,1)
0299 hold off
0300 plot(pk,'.r', 'MarkerSize',14);
0301 hold on;
0302 plot(Primal_constrk, 'LineWidth',1);
0303 
0304 
0305 set(gca,'FontSize',16, 'XLim',[1 N] );
0306 title({'BPDN shrinkage constraints,'; ['n = ',num2str(N), ', m = ', num2str(M), ', T = ',num2str(T)]});
0307 plot(1:N, epsilon*ones(1,N),'--k','MarkerSize',12);
0308 plot(1:N, -epsilon*ones(1,N), '--k','MarkerSize',12);
0309 plot(1:N, epsilon_old*ones(1,N),'--m','MarkerSize',12);
0310 plot(1:N, -epsilon_old*ones(1,N), '--m','MarkerSize',12);
0311 
0312 figure(1);
0313 subplot(2,1,2)
0314 hold off
0315 plot(x,'.r','MarkerSize',14); hold on;
0316 plot(xk_1,'LineWidth',1);
0317 set(gca,'FontSize',16,'XLim',[1 N]);
0318 title(['BPDN estimate at \tau = ',num2str(thresh), ', iter = ', num2str(iter)]);