wtBPDN_homotopy_test

 BPDN_homotopy_Weighted

 Solves the following basis pursuit denoising (BPDN) problem
 min_x  \Sum \e_i |x_i| + 1/2*||y-Ax||_2^2

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

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

0001 % BPDN_homotopy_Weighted
0002 %
0003 % Solves the following basis pursuit denoising (BPDN) problem
0004 % min_x  \Sum \e_i |x_i| + 1/2*||y-Ax||_2^2
0005 %
0006 % Written by: Salman Asif, Georgia Tech
0007 % Email: sasif@ece.gatech.edu
0008 %
0009 %-------------------------------------------+
0010 % Copyright (c) 2010.  Muhammad Salman Asif
0011 %-------------------------------------------+
0012 
0013 clear
0014 clc
0015 
0016 % load RandomStates
0017 rseed = sum(100*clock);
0018 rseed = 0;
0019 % rand('state',rseed);
0020 % randn('state',rseed);
0021 RandStream.setDefaultStream(RandStream('mt19937ar','seed',rseed));
0022 
0023 N = 128;   % signal length
0024 T = 20;    % sparsity level
0025 M = 100;    % no. of measurements
0026 
0027 % Generate a random signal
0028 x = zeros(N,1);
0029 q = randperm(N);
0030 x(q(1:T)) = sign(randn(T,1));
0031 
0032 % measurement matrix
0033 A = randn(M,N)/sqrt(M); % Gaussian
0034 
0035 % Hadamard
0036 % H = hadamard(N);
0037 % A = H(q(1:M),:)/sqrt(M);
0038 
0039 % Bernoulli
0040 % A = randsrc(M,N)/sqrt(M);
0041 
0042 % Random Projection
0043 % A = orth(A')';
0044 
0045 % AtA = A'*A;
0046 
0047 % measurements
0048 sigma = 0.01;
0049 e = randn(M,1)*sigma;
0050 y = A*x+e;
0051 
0052 
0053 % Regularization parameters
0054 epsilon_vec = [ones(N/2,1); ones(N/2,1)*1e4]*1e-2;
0055 
0056 tau = 1e-2;
0057 W_new = ones(N,1);
0058 W_new(abs(x)>0) = 1e-5./abs(x(abs(x)>0));
0059 
0060 epsilon_vec = tau.*W_new;
0061 unique_eps = sort(unique(epsilon_vec),'descend');
0062 thresh = min(epsilon_vec);
0063 
0064 maxiter = 4*N;
0065 
0066 in = [];
0067 in.x_orig = x;
0068 in.epsilon_vec = epsilon_vec;
0069 in.maxiter = maxiter;
0070 in.delx_mode = 'mil';
0071 out = wtBPDN_function(A,y,in);
0072 xh_wt = out.x_out;
0073 gamma_wt = out.gamma;
0074 iter_wt = out.iter;
0075 
0076 %%
0077 
0078 % Initialization of primal sign and support
0079 z_x = zeros(N,1);
0080 gamma_x = [];       % Primal support
0081 
0082 % Initial step
0083 Primal_constrk = -A'*y; 
0084 constr_mask = abs(Primal_constrk)>epsilon_vec;
0085 [c i] = max(abs(Primal_constrk.*constr_mask));
0086 eps_iter = sum(unique_eps>c)+1;
0087 
0088 gamma_xk = i;
0089 
0090 z_x(gamma_xk) = -sign(Primal_constrk(gamma_xk));
0091 epsilon = c;
0092 Primal_constrk(gamma_xk) = sign(Primal_constrk(gamma_xk))*epsilon;
0093 xk_1 = zeros(N,1);
0094 
0095 z_xk = z_x;
0096 
0097 % loop parameters
0098 done = 0;
0099 iter = 0;
0100 data_precision = eps;   % floating point precision
0101 
0102 old_delta = 0;
0103 out_x = [];
0104 count_delta_stop = 0;
0105 
0106 % Turn the plots on and off
0107 constraint_plots = 1; % 1 - draw the plots after every homotopy iteration.
0108 
0109 AtAgx = A(:,gamma_xk)'*A(:,gamma_xk);
0110 iAtAgx = inv(AtAgx);
0111 
0112 while ~done
0113     iter = iter+1;
0114     % warning('off','MATLAB:divideByZero')
0115 
0116     gamma_x = gamma_xk;
0117     z_x = z_xk;
0118     x_k = xk_1;
0119 
0120     %%%%%%%%%%%%%%%%%%%%%
0121     %%%% update on x %%%%
0122     %%%%%%%%%%%%%%%%%%%%%
0123 
0124     % Update direction
0125 %     del_x = -inv(A(:,gamma_lambda)'*A(:,gamma_x))*z_lambda(gamma_lambda);
0126 %     diff_inv = [ max(max(abs(inv(A(:,gamma_lambda)'*A(:,gamma_x))-iAtglAgx)))]
0127 %     if diff_inv > 1e-8
0128 %         stp =1;
0129 %     end
0130     indicator_temp = epsilon>epsilon_vec(gamma_x);
0131     del_x = iAtAgx*(indicator_temp.*z_x(gamma_x));
0132     del_x_vec = zeros(N,1);
0133     del_x_vec(gamma_x) = del_x;
0134 
0135     pk = Primal_constrk;
0136     % dk = AtA*del_x_vec;
0137     dk = A'*(A(:,gamma_x)*del_x);
0138     
0139     %%% CONTROL THE MACHINE PRECISION ERROR AT EVERY OPERATION: LIKE BELOW.
0140     pk_temp = Primal_constrk;
0141     gammaL_temp = find(abs(abs(Primal_constrk)-epsilon)<min(epsilon,1e-12));
0142 %     pk_temp(gammaL_temp) = sign(Primal_constrk(gammaL_temp))*epsilon;
0143     
0144     xk_temp = x_k;
0145     gammaX_temp = find(abs(x_k)<1*eps);
0146     %%%---
0147     
0148     temp_gamma = zeros(N,1);
0149     temp_gamma(gamma_x) = gamma_x;
0150     gamma_xc = find([1:N]' ~= temp_gamma);
0151     
0152     % For different values of regularization parameters
0153     eps_temp = epsilon.*(epsilon>epsilon_vec(gamma_xc))+epsilon_vec(gamma_xc).*(epsilon<=epsilon_vec(gamma_xc));
0154     one_temp = epsilon>epsilon_vec(gamma_xc);
0155     if find(one_temp-1)
0156         stp = 1;
0157     end
0158     if find(eps_temp-epsilon)
0159         stp = 1;
0160     end
0161     delta1_constr = (eps_temp-pk(gamma_xc))./(one_temp+dk(gamma_xc));
0162     delta1_pos_ind = find(delta1_constr>0);
0163     delta1_pos = delta1_constr(delta1_pos_ind);
0164     [delta1 i_delta1] = min(delta1_pos);
0165     if isempty(delta1)
0166         delta1 = inf;
0167     end
0168     delta2_constr = (eps_temp+pk(gamma_xc))./(one_temp-dk(gamma_xc));
0169     delta2_pos_ind = find(delta2_constr>0);
0170     delta2_pos = delta2_constr(delta2_pos_ind);
0171     [delta2 i_delta2] = min(delta2_pos);
0172     if isempty(delta2)
0173         delta2 = inf;
0174     end
0175     
0176     if delta1>delta2
0177         delta = delta2;
0178         i_delta = gamma_xc(delta2_pos_ind(i_delta2));
0179     else
0180         delta = delta1;
0181         i_delta = gamma_xc(delta1_pos_ind(i_delta1));
0182     end
0183     
0184     delta3_constr = (-x_k(gamma_x)./del_x_vec(gamma_x));
0185     delta3_pos_index = find(delta3_constr>0);
0186     [delta3 i_delta3] = min(delta3_constr(delta3_pos_index));
0187     out_x_index = gamma_x(delta3_pos_index(i_delta3));
0188     
0189     chk_x = 0;
0190     out_x = [];
0191     if delta3 > 0 & delta3 <= delta
0192         chk_x = 1;
0193         delta = delta3;
0194         out_x = out_x_index;
0195     end
0196 
0197     [i_delta out_x delta -chk_x]
0198     if old_delta < 4*eps & delta < 4*eps
0199         count_delta_stop = count_delta_stop + 1;
0200     else
0201         count_delta_stop = 0;
0202     end
0203     if count_delta_stop >= 50
0204         disp('stuck in some corner');
0205         break;
0206     end
0207     old_delta = delta;
0208 
0209     xk_1 = x_k+delta*del_x_vec;
0210     Primal_constrk = pk+delta*dk;
0211     epsilon_old = epsilon;
0212     epsilon = epsilon-delta;
0213     
0214     if epsilon <= unique_eps(eps_iter)
0215         epsilon = unique_eps(eps_iter);
0216         delta_end = epsilon_old-epsilon;
0217         Primal_constrk = pk+delta_end*dk;
0218         epsilon_temp = epsilon.*(epsilon>epsilon_vec(gamma_xk))+epsilon_vec(gamma_xk).*(epsilon<=epsilon_vec(gamma_xk));
0219         Primal_constrk(gamma_x) = sign(Primal_constrk(gamma_x)).*epsilon_temp;
0220  
0221         xk_1 = x_k + delta_end*del_x_vec;
0222         eps_iter = eps_iter+1;
0223         if eps_iter > length(unique_eps)
0224             disp('done!');
0225             break;
0226         else
0227             disp('switch epsilon!');
0228             continue; 
0229         end
0230     end
0231     if epsilon <= min(epsilon_vec); %sqrt(2*log(N))*sigma; %1e-7 %|| iter > 5*T || (length(gamma_lambda) == K)
0232         delta_end = epsilon_old-thresh;
0233         Primal_constrk = pk+delta_end*dk;
0234         xk_1 = x_k + delta_end*del_x_vec;
0235         disp('done!');
0236         break;
0237     end
0238     
0239     if length(gamma_x)==M & chk_x ==0
0240         stp =1;
0241         disp('Cannot do it Sire'); % Commondos: are you out of your mind Sire!
0242         break;
0243     end
0244 
0245     if chk_x == 1
0246         % If an element is removed from gamma_x
0247         gx_old = gamma_x;
0248 
0249         outx_index = find(gamma_x==out_x);
0250         gamma_x(outx_index) = gamma_x(end);
0251         gamma_x(end) = out_x;
0252         gamma_xk = gamma_x(1:end-1);
0253 
0254         rowi = outx_index; % ith row of A is swapped with last row (out_x)
0255         colj = outx_index; % jth column of A is swapped with last column (out_lambda)
0256         AtAgx_ij = AtAgx;
0257         temp_row = AtAgx_ij(rowi,:);
0258         AtAgx_ij(rowi,:) = AtAgx_ij(end,:);
0259         AtAgx_ij(end,:) = temp_row;
0260         temp_col = AtAgx_ij(:,colj);
0261         AtAgx_ij(:,colj) = AtAgx_ij(:,end);
0262         AtAgx_ij(:,end) = temp_col;
0263         
0264         iAtAgx_ij = iAtAgx;
0265         temp_row = iAtAgx_ij(colj,:);
0266         iAtAgx_ij(colj,:) = iAtAgx_ij(end,:);
0267         iAtAgx_ij(end,:) = temp_row;
0268         temp_col = iAtAgx_ij(:,rowi);
0269         iAtAgx_ij(:,rowi) = iAtAgx_ij(:,end);
0270         iAtAgx_ij(:,end) = temp_col;
0271 
0272         AtAgx = AtAgx_ij(1:end-1,1:end-1);
0273         iAtAgx = update_inverse(AtAgx_ij, iAtAgx_ij,2);
0274         xk_1(out_x) = 0;
0275     else
0276         % If an element is added to gamma_x
0277         gamma_xk = [gamma_x; i_delta];
0278         i_theta = i_delta;
0279         new_x = i_delta;
0280         AtAnl = A(:,gamma_x)'*A(:,new_x);
0281         AtAgx_mod = [AtAgx A(:,gamma_x)'*A(:,i_theta); AtAnl' A(:,new_x)'*A(:,i_theta)];
0282 
0283         AtAgx = AtAgx_mod;
0284         iAtAgx = update_inverse(AtAgx, iAtAgx,1);
0285         xk_1(i_theta) = 0;
0286         
0287         gamma_x  = gamma_xk;
0288     end
0289 
0290     z_xk = zeros(N,1);
0291     z_xk(gamma_xk) = -sign(Primal_constrk(gamma_xk));
0292     epsilon_temp = epsilon.*(epsilon>epsilon_vec(gamma_x))+epsilon_vec(gamma_x).*(epsilon<=epsilon_vec(gamma_x));
0293     Primal_constrk(gamma_x) = sign(Primal_constrk(gamma_x)).*epsilon_temp;
0294     
0295     if iter > 500*T
0296         disp('too many iterations ooooooooops');
0297         break;
0298     end
0299 
0300     if constraint_plots
0301         fig = figure(1);
0302         subplot(2,1,1)
0303         hold off
0304         plot(pk,'.r', 'MarkerSize',14);
0305         hold on;
0306         plot(Primal_constrk, 'LineWidth',1);
0307 
0308         if chk_x == 0
0309             plot(new_x, Primal_constrk(new_x),'or','MarkerSize',18,'LineWidth',2);
0310             text(new_x, Primal_constrk(new_x)*1.1, ['Incoming \gamma = ',num2str(new_x)],'FontSize',14);
0311         else
0312             plot(out_x, Primal_constrk(out_x),'*k','MarkerSize',18,'LineWidth',2);
0313             text(out_x,Primal_constrk(out_x)*1.1, ['Outgoing \gamma = ',num2str(out_x)],'FontSize',14);
0314         end
0315         set(gca,'FontSize',16, 'XLim',[1 N] );
0316         title({'BPDN shrinkage constraints,'; ['n = ',num2str(N), ', m = ', num2str(M), ', T = ',num2str(T)]});
0317         plot(1:N, epsilon*ones(1,N),'--k','MarkerSize',12);
0318         plot(1:N, -epsilon*ones(1,N), '--k','MarkerSize',12);
0319         plot(1:N, epsilon_old*ones(1,N),'--m','MarkerSize',12);
0320         plot(1:N, -epsilon_old*ones(1,N), '--m','MarkerSize',12);
0321 
0322         figure(1);
0323         subplot(2,1,2)
0324         hold off
0325         plot(x_k,'.r','MarkerSize',14); hold on;
0326         plot(xk_1,'LineWidth',1);
0327         if chk_x == 1
0328             plot(out_x, 0,'ok', 'MarkerSize',18,'LineWidth',2);
0329         else
0330             plot(new_x, 0,'or', 'MarkerSize',18,'LineWidth',2);
0331         end
0332         set(gca,'FontSize',16,'XLim',[1 N]);
0333         title(['BPDN estimate at \tau = ',num2str(epsilon), ', iter = ', num2str(iter)]);
0334         
0335         if iter == 1
0336             disp('  ');
0337             disp('Every frame in the figure corresponds to a critical point on the homotopy path.')
0338             disp('Circle represents an incoming element, star represents an outgoing element.');
0339             disp(' ');
0340             disp('Put pause somewhere in the code to see this. ');
0341             disp('For now press some key to continue...');
0342             %pause
0343         end
0344 
0345         %drawnow;
0346         %print(gcf,'-dbmp','-r200',['BPDNPath_',num2str(iter)])
0347     end
0348 %     pause
0349 %     [iter epsilon delta]
0350 end
0351 
0352 fig = figure(1);
0353 subplot(2,1,1)
0354 hold off
0355 plot(pk,'.r', 'MarkerSize',14);
0356 hold on;
0357 plot(Primal_constrk, 'LineWidth',1);
0358 
0359 
0360 set(gca,'FontSize',16, 'XLim',[1 N] );
0361 title({'BPDN shrinkage constraints,'; ['n = ',num2str(N), ', m = ', num2str(M), ', T = ',num2str(T)]});
0362 plot(1:N, epsilon*ones(1,N),'--k','MarkerSize',12);
0363 plot(1:N, -epsilon*ones(1,N), '--k','MarkerSize',12);
0364 plot(1:N, epsilon_old*ones(1,N),'--m','MarkerSize',12);
0365 plot(1:N, -epsilon_old*ones(1,N), '--m','MarkerSize',12);
0366 
0367 figure(1);
0368 subplot(2,1,2)
0369 hold off
0370 plot(x_k,'.r','MarkerSize',14); hold on;
0371 plot(xk_1,'LineWidth',1);
0372 set(gca,'FontSize',16,'XLim',[1 N]);
0373 title(['BPDN estimate at \tau = ',num2str(thresh), ', iter = ', num2str(iter)]);
0374 
0375 %%
0376 digit = 6; if sigma > 0; digit = 7; end
0377 opts = [];
0378 opts.tol = 10^(-digit);
0379 W_new = ones(N,1);
0380 opts.weights = W_new;
0381 opts.print = 0;
0382 opts.maxit = maxiter;
0383 
0384 opts.nu = 0; opts.rho = tau;
0385 tic;
0386 [x_yall1,Out] = yall1(A,y,opts);
0387 
0388 %% CVX check
0389 cvx_begin
0390     cvx_precision best
0391     variable xc(N);
0392     minimize( sum(epsilon_vec.*abs(xc))+ 1/2*sum(square(A*xc-y)))
0393 cvx_end
0394 
0395 %%
0396 figure(101); plot([xc xk_1])
0397 figure(102); plot([xc-xk_1])