Simulation_l1Decode

 Test out l1 decoding with homotopy update scheme

 Created by Salman Asif @ Georgia Tech

0001 % Test out l1 decoding with homotopy update scheme
0002 %
0003 % Created by Salman Asif @ Georgia Tech
0004 
0005 close all
0006 clear
0007 
0008 randn('state',0);
0009 rand('state',0);
0010 
0011 % source length
0012 N = 128;
0013 
0014 % Number of errors introduced at random
0015 T_list = [10 20 30 40 50];
0016 P_list = [1 2 4 8];
0017 
0018 Stack_meas = [];
0019 Stack_iter = [];
0020 
0021 no_exp = 20; % number of simulations
0022 
0023 for P = P_list
0024     % Number of new elements added at a time
0025     No_new_element = P;
0026     
0027     table_iter = zeros(no_exp, length(T_list));
0028     table_error = zeros(no_exp, length(T_list));
0029     
0030     fprintf('N = %d, P = %d, sim_count = %d...\n',N,No_new_element,no_exp);
0031     
0032     for iT = 1:length(T_list)
0033         for inn_iter = 1:no_exp
0034             
0035             % Maximum allowed length of message
0036             M = 20*N;
0037             
0038             % number of perturbations
0039             T = T_list(iT); %round(.3*N);
0040             
0041             % coding matrix
0042             Orth_mat = randn(M,M);
0043             % Orth_mat = randsrc(M,M);
0044             % [Orth_mat R_mat] = qr(Orth_mat,0);
0045             G_mat = Orth_mat(:,1:N);
0046             
0047             % Homotopy setup
0048             m_st = N;   % number of measurements to start with
0049             A = G_mat(1:m_st,:);
0050             B_stack = G_mat(m_st+1:end,:);
0051             
0052             % source word
0053             x = randn(N,1);
0054             
0055             % channel: perturb T randomly chosen entries
0056             q = randperm(m_st-1);
0057             c = zeros(m_st,1);
0058             c(q(1:T)) = randn(T,1);
0059             
0060             y = A*x-c;
0061             
0062             % Introduce sparse errors
0063             err_loc = (rand(M-N,1)>1);
0064             err_new = err_loc.*randn(M-N,1);
0065             w_vec = B_stack*x+err_new;
0066             
0067             % initial estimate of x
0068             % l1_min or direct inversion depending on n_st
0069             iA = inv(A);
0070             x0 = iA*y;
0071             e0 = A*x0-y;
0072             lambda_0 = zeros(m_st,1);
0073             
0074             A_old = A;
0075             y_old = y;
0076             x_old = x0;
0077             c_old = e0;
0078             lambda_old = lambda_0;
0079             
0080             m_u = No_new_element;
0081             gamma_old = []; %find(A*x0-y) or find(abs(lambda_0)==1);
0082             m_old = m_st;
0083             all_done = 0;
0084             cond_list = [];
0085             cond_iter = 1;
0086             gamma_kc_new = [1:N]';
0087             iGg_new = iA';
0088             
0089             iter_total = 0;
0090             iter_in = 0;
0091             for iter_mat = 1:m_u:M-m_st
0092                 m_new = m_old + m_u;
0093                 x_k = x_old;
0094                 lambda_k = lambda_old;
0095                 gamma_k = gamma_old;
0096                 %     gamma_kc = gamma_kc_old;
0097                 
0098                 B = B_stack(iter_mat:iter_mat+m_u-1,:);
0099                 G = [A_old' B'];
0100                 w = w_vec(iter_mat:iter_mat+m_u-1);
0101                 q = [y_old; w];
0102                 
0103                 c_old = [c_old; B*x_k-w];
0104                 
0105                 if abs(abs(B*x_k-w))>=10*eps
0106                     nu = sign(B*x_k-w);
0107                     i_nu = [m_new-m_u+1:m_new]';
0108                     gammak_new = [gamma_k; i_nu];
0109                 else
0110                     stp = 1;
0111                 end
0112                 
0113                 lambda_old = [lambda_k; nu]; % xi_old
0114                 
0115                 done = 0;
0116                 epsilon_old = 0;
0117                 iter_in = 0;
0118                 while ~done
0119                     iter_in = iter_in+1;
0120                     % [iter_mat iter_in]
0121                     
0122                     gamma_k = gammak_new;
0123                     x_k = x_old;
0124                     c_k = c_old;
0125                     lambda_k = lambda_old; % xi_k
0126                     
0127                     gamma_kc = gamma_kc_new;
0128                     
0129                     % INVERSE UPDATE
0130                     %iGg = inv(G(:,gamma_kc));
0131                     iGg = iGg_new;
0132                     del_lambda = -iGg*(G(:,i_nu)*sign(lambda_k(i_nu))); % del_xi
0133                     
0134                     del_lambda_vec = zeros(m_new,1); % del_xi_vec
0135                     del_lambda_vec(gamma_kc) = del_lambda;
0136                     
0137                     % find incoming elements in support of error vector
0138                     constr1 = (1-lambda_k(gamma_kc))./del_lambda_vec(gamma_kc);
0139                     constr1_pos = constr1(constr1>2*eps);
0140                     delta1 = min(constr1_pos);
0141                     idelta1 = gamma_kc(find(constr1==delta1));
0142                     
0143                     constr2 = -(1+lambda_k(gamma_kc))./del_lambda_vec(gamma_kc);
0144                     constr2_pos = constr2(constr2>2*eps);
0145                     delta2 = min(constr2_pos);
0146                     idelta2 = gamma_kc(find(constr2==delta2));
0147                     
0148                     if delta1>delta2
0149                         delta = delta2;
0150                         idelta = idelta2;
0151                     else
0152                         delta = delta1;
0153                         idelta = idelta1;
0154                     end
0155                     
0156                     if epsilon_old+delta <=1
0157                         lambdak_1 = lambda_k + delta*del_lambda_vec;
0158                         lambdak_1(idelta) = sign(lambdak_1(idelta));
0159                         gammak_temp = [gamma_k; idelta];
0160                         
0161                         % outgoing element from the support of
0162                         uz_vec = zeros(N,1);
0163                         idelta_gamma_kc = find(idelta==gamma_kc);
0164                         uz_vec(idelta_gamma_kc) = sign(lambdak_1(idelta));
0165                         delx = iGg'*uz_vec;
0166                         
0167                         delc= G'*delx;
0168                         constr3 = -(c_k(gamma_k)./delc(gamma_k));
0169                         constr3_pos = constr3(constr3>2*eps);
0170                         delta3 = min(constr3_pos);
0171                         idelta3 = gamma_k(find(constr3==delta3));
0172                         xk_1 = x_k + delta3*delx;
0173                         ck_1 = c_k + delta3*delc;
0174                         ck_1(idelta3) = 0;
0175                         
0176                         out_lambda = idelta3;
0177                         
0178                         if [m_new-N~=length(find(abs(G'*xk_1-q)>=1e5*eps))]
0179                             gamma_final = setdiff(gammak_temp,find(abs(G'*xk_1-q)<=100*eps));
0180                             x_old = xk_1;
0181                             c_old = ck_1;
0182                             lambda_old = lambdak_1;
0183                             stp = 1; %% WHAT TO DO NOW ???
0184                             all_done = 1;
0185                             iter_total = iter_total+iter_in; 
0186                             break; break; break;
0187                         end 
0188                         gammak_new = setdiff(gammak_temp,idelta3);
0189                         gamma_kc_new = gamma_kc;
0190                         gamma_kc_new(idelta_gamma_kc) = idelta3;
0191                         
0192                         % UPDATE INVERSE OF G(gamma_kc)
0193                         % Gg = Gg + (g_new - g_old)'*1_idelta_gamma_kc;
0194                         diff_G_in_out = G(:,idelta3)-G(:,idelta);
0195                         iGg_new = iGg - (iGg*diff_G_in_out/(1+iGg(idelta_gamma_kc,:)*diff_G_in_out))*(iGg(idelta_gamma_kc,:));
0196                         
0197                         x_old = xk_1;
0198                         c_old = ck_1;
0199                         lambda_old = lambdak_1;
0200                         gamma_old = gammak_new;
0201                         
0202                         if length(find(idelta3 == i_nu))
0203                             lambda_old(idelta3) = (epsilon_old+delta)*lambdak_1(idelta3);
0204                             i_nu = setdiff(i_nu,idelta3);
0205                             if isempty(i_nu)
0206                                 epsilon = 1;
0207                             else
0208                                 epsilon = delta + epsilon_old;
0209                             end
0210                         else
0211                             epsilon = delta + epsilon_old;
0212                         end
0213                         
0214                         if epsilon >=1
0215                             done = 1;
0216                             
0217                             A_old = G';
0218                             y_old = q;
0219                             m_old = m_new;
0220                             iter_total = iter_total+iter_in;
0221                         end
0222                         epsilon_old = epsilon;
0223                         
0224                     else
0225                         new_delta = 1-epsilon_old;
0226                         lambdak_1 = lambda_k + new_delta*del_lambda_vec;
0227                         lambda_old = [lambdak_1];
0228                         gamma_old = gammak_new;
0229                         
0230                         iter_total = iter_total+iter_in;
0231                         done = 1;
0232                         
0233                         A_old = G';
0234                         y_old = q;
0235                         m_old = m_new;
0236                     end
0237                 end
0238                 % Current
0239                 % disp(['No. of errors = ',num2str(T_list(iT)), ', Simulation run = ', num2str(inn_iter),', Homotopy iterations = ', num2str(iter_total), ', No. of measurements = ',num2str(iter_mat+m_u-1+m_st),', Error in recovery = ',num2str(norm(x-xk_1))]);
0240                 
0241                 if all_done == 1
0242                     break;
0243                 end
0244             end
0245             % Total number of homotopy iterations.
0246             table_iterT(inn_iter,iT) = iter_total;
0247             % Total number of measurements required for perfect recovery!
0248             table_meas(inn_iter,iT) = iter_mat+m_u-1+m_st;
0249         end
0250     end
0251     % With few measurements smaller number of steps taken to update per
0252     % measuremnt and as we approach final solution number of steps per
0253     % measurement increase slightly.
0254     disp('Final results');
0255     Error_count = [T_list] % Total number of sparse errors
0256     Total_measurements = mean(table_meas,1) % Total number of measurements required for perfect recovery!
0257     Total_iterations = mean(table_iterT,1)  % Total number of homotopy iterations.
0258     Average_iterations = Total_iterations./(Total_measurements-N) % Avg homotopy iterations per new measurement
0259     
0260     Stack_meas = [Stack_meas; Total_measurements];
0261     Stack_iter = [Stack_iter; Total_iterations];    
0262 end
0263 
0264 %%
0265 plot_results = 1;
0266 if plot_results
0267 exp_name = 'l1decoding_N64_Exp20';
0268 exp_name = 'l1decoding_N128_Exp20';
0269 
0270 load(exp_name);
0271 
0272 addpath ../
0273 addpath ../utils/utils_fig 
0274 
0275 % figure setup
0276 marker = {'bx','r+','k*','gd'};
0277 linewidth = 1.5;
0278 
0279 
0280 axis_prop = {};
0281 axis_prop{1,4} = 'FontSize'; axis_prop{2,4} = 14;
0282 axis_prop{1,5} = 'FontWeight'; axis_prop{2,5} = 'normal';
0283 axis_prop{1,8} = 'XGrid'; axis_prop{2,8} = 'on';
0284 axis_prop{1,9} = 'GridLineStyle'; axis_prop{2,9} = '--';
0285 
0286 ydata = Stack_iter./(Stack_meas-N); 
0287 figure(1)
0288 clf; hold on;
0289 set(gca,'FontSize',14,'FontWeight','normal');
0290 plot(T_list, ydata(1,:),'-bx', 'LineWidth',1.5,'MarkerSize',10);
0291 plot(T_list, ydata(2,:),'-.r+', 'LineWidth',1.5,'MarkerSize',10);
0292 plot(T_list, ydata(3,:),':k*', 'LineWidth',1.5,'MarkerSize',10);
0293 plot(T_list, ydata(4,:),'--md', 'LineWidth',1.5,'MarkerSize',10);
0294 xlabel('Number of errors');
0295 ylabel('Number of homotopy iterations');
0296 % title('Average number of homotopy per measurement');
0297 legend('P=1','P=2','P=4','P=8','Location','NorthWest');
0298 axis tight;
0299 
0300 img_name = sprintf('%s-AvgIter',exp_name);
0301 set(gcf, 'Color', 'w');
0302 set(gcf,'Position',[300 0 500 600]);
0303 
0304 eval(sprintf('export_fig %s.png',img_name));
0305 eval(sprintf('export_fig %s.pdf',img_name));
0306 
0307 ydata = Stack_iter; 
0308 figure(2)
0309 clf; hold on;
0310 set(gca,'FontSize',14,'FontWeight','normal');
0311 plot(T_list, ydata(1,:),'-bx', 'LineWidth',1.5,'MarkerSize',10);
0312 plot(T_list, ydata(2,:),'-.r+', 'LineWidth',1.5,'MarkerSize',10);
0313 plot(T_list, ydata(3,:),':k*', 'LineWidth',1.5,'MarkerSize',10);
0314 plot(T_list, ydata(4,:),'--md', 'LineWidth',1.5,'MarkerSize',10);
0315 xlabel('Number of errors');
0316 ylabel('Number of homotopy iterations');
0317 % title('Total number of homotopy iteration');
0318 legend('P=1','P=2','P=4','P=8','Location','NorthWest');
0319 axis tight;
0320 
0321 img_name = sprintf('%s-TotalIter',exp_name);
0322 set(gcf, 'Color', 'w');
0323 set(gcf,'Position',[300 0 500 600]);
0324 
0325 eval(sprintf('export_fig %s.png',img_name));
0326 eval(sprintf('export_fig %s.pdf',img_name));
0327 end