l1Decode_homotopy_qr

 l1Decode_homotopy_fast.m

 Test out l1 decoding with homotopy update scheme
 Keeps adding new measurements until the original signal is recovered
 exactly.

 In this script, we start the decoding with $m=n$ measurements and
 sequentially add m_u new measurements and solve using homotopy.

 Created by Salman Asif @ Georgia Tech

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

0001 % l1Decode_homotopy_fast.m
0002 %
0003 % Test out l1 decoding with homotopy update scheme
0004 % Keeps adding new measurements until the original signal is recovered
0005 % exactly.
0006 %
0007 % In this script, we start the decoding with $m=n$ measurements and
0008 % sequentially add m_u new measurements and solve using homotopy.
0009 %
0010 % Created by Salman Asif @ Georgia Tech
0011 %
0012 %-------------------------------------------+
0013 % Copyright (c) 2009.  Muhammad Salman Asif
0014 %-------------------------------------------+
0015 
0016 
0017 % Modification history :
0018 % 09/29/08
0019 % 10/02/08 : multiple measurements
0020 % 02/15/09 : fast update of inverse matrix
0021 % 06/21/09 : Solve using qr factorization
0022 
0023 close all
0024 clear
0025 
0026 % % load fixed random states
0027 % load RandomStates
0028 % rand('state', s_rand);
0029 % randn('state', s_randn);
0030 
0031 % source length
0032 N = 256;
0033 
0034 % codeword length
0035 M = 4*N;
0036 
0037 % number of perturbations
0038 T = round(.3*N);
0039 
0040 % coding matrix
0041 Orth_mat = randn(M,M);
0042 Orth_mat = randsrc(M,M);
0043 % Orth_mat = orth(Orth_mat);
0044 G_mat = Orth_mat(:,1:N);
0045 
0046 % Homotopy setup
0047 m_st = N;   % number of measurements to start with
0048 A = G_mat(1:m_st,:);
0049 B_stack = G_mat(m_st+1:end,:);
0050 
0051 % source word
0052 x = randn(N,1);
0053 
0054 % channel: perturb T randomly chosen entries
0055 q = randperm(m_st-1);
0056 c = zeros(m_st,1);
0057 c(q(1:T)) = randn(T,1);
0058 y = A*x-c;
0059 err_loc = (rand(M-N,1)>=.9);
0060 err_new = err_loc.*randn(M-N,1);
0061 w_vec = B_stack*x-err_new;
0062 figure(1); clf; plot([y c]);
0063 title('Codeword and sparse errors');
0064 legend('Codeword', 'Sparse'); shg
0065 
0066 
0067 % initial estimate of x
0068 % l1_min or direct solution (LS) depending on n_st
0069 % iA = inv(A);
0070 % x0 = iA*y;
0071 [Qg Rg] = qr(A');
0072 x0 = Qg*(Rg'\y);
0073 e0 = A*x0-y;
0074 lambda_0 = zeros(m_st,1);
0075 
0076 A_old = A;
0077 y_old = y;
0078 x_old = x0;
0079 c_old = e0;
0080 lambda_old = lambda_0;
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 m_u = 4; % number of new measurements added at a time
0090 
0091 total_iter = 0;
0092 
0093 for iter_mat = 1:m_u:M-m_st
0094     m_new = m_old + m_u;
0095     x_k = x_old;
0096     lambda_k = lambda_old;
0097     gamma_k = gamma_old;
0098     %     gamma_kc = gamma_kc_old;
0099 
0100     B = B_stack(iter_mat:iter_mat+m_u-1,:);
0101     G = [A_old' B'];
0102     w = w_vec(iter_mat:iter_mat+m_u-1);
0103     q = [y_old; w];
0104     c_original = [c; err_new(1:iter_mat+m_u-1)];
0105     c_old = [c_old; B*x_k-w];
0106     
0107     % Error estimate for new measurements.
0108     if abs(abs(B*x_k-w))>=10*eps
0109         nu = sign(B*x_k-w);
0110         i_nu = [m_new-m_u+1:m_new]';
0111         gammak_new = [gamma_k; i_nu];
0112     else
0113         stp = 1; % Something is wrong! (degenracy)
0114         break;
0115     end
0116 
0117     lambda_old = [lambda_k; nu]; % xi_old
0118 
0119     done = 0;
0120     epsilon_old = 0;
0121     iter_in = 0;
0122     while ~done
0123         iter_in = iter_in+1;
0124         % [iter_mat iter_in]
0125     
0126         gamma_k = gammak_new;
0127         x_k = x_old;
0128         c_k = c_old;
0129         lambda_k = lambda_old; % xi_k
0130 
0131         gamma_kc = gamma_kc_new;
0132 
0133         % INVERSE UPDATE
0134         %iGg = inv(G(:,gamma_kc));
0135         % iGg = iGg_new;
0136         % del_lambda = -iGg*(G(:,i_nu)*sign(lambda_k(i_nu))); % del_xi
0137         
0138         del_lambda = -Rg\(Qg'*(G(:,i_nu)*sign(lambda_k(i_nu))));
0139         
0140         del_lambda_vec = zeros(m_new,1); % del_xi_vec
0141         del_lambda_vec(gamma_kc) = del_lambda;
0142 
0143         % find incoming elements in support of error vector
0144         constr1 = (1-lambda_k(gamma_kc))./del_lambda_vec(gamma_kc);
0145         constr1_pos = constr1(constr1>2*eps);
0146         delta1 = min(constr1_pos);
0147         idelta1 = gamma_kc(find(constr1==delta1));
0148 
0149         constr2 = -(1+lambda_k(gamma_kc))./del_lambda_vec(gamma_kc);
0150         constr2_pos = constr2(constr2>2*eps);
0151         delta2 = min(constr2_pos);
0152         idelta2 = gamma_kc(find(constr2==delta2));
0153 
0154         if delta1>delta2
0155             delta = delta2;
0156             idelta = idelta2;
0157         else
0158             delta = delta1;
0159             idelta = idelta1;
0160         end
0161 
0162         if epsilon_old+delta <=1
0163             lambdak_1 = lambda_k + delta*del_lambda_vec;
0164             lambdak_1(idelta) = sign(lambdak_1(idelta));
0165             gammak_temp = [gamma_k; idelta];
0166 
0167             % outgoing element from the support of
0168             uz_vec = zeros(N,1);
0169             idelta_gamma_kc = find(idelta==gamma_kc);
0170             uz_vec(idelta_gamma_kc) = sign(lambdak_1(idelta));
0171 
0172             %delx = iGg'*uz_vec;
0173             
0174             delx = Qg*(Rg'\uz_vec);
0175             
0176             delc= G'*delx;
0177             constr3 = -(c_k(gamma_k)./delc(gamma_k));
0178             constr3_pos = constr3(constr3>2*eps);
0179             delta3 = min(constr3_pos);
0180             idelta3 = gamma_k(find(constr3==delta3));
0181             xk_1 = x_k + delta3*delx;
0182             ck_1 = c_k + delta3*delc;
0183             ck_1(idelta3) = 0;
0184 
0185             out_lambda = idelta3;
0186             
0187             % See if signal is decoded exactly...
0188             if [m_new-N~=length(find(abs(ck_1)>=1e4*eps))]
0189                 gamma_final = setdiff(gammak_temp,find(abs(G'*xk_1-q)<=100*eps));
0190                 x_old = xk_1;
0191                 c_old = ck_1;
0192                 lambda_old = lambdak_1;
0193                 stp = 1; %% WHAT TO DO NOW ???
0194                 all_done = 1; % Take a break!
0195                 break; break; break;
0196             end
0197             gammak_new = setdiff(gammak_temp,idelta3);
0198             gamma_kc_new = gamma_kc;
0199             gamma_kc_new(idelta_gamma_kc) = idelta3;
0200 
0201             % UPDATE INVERSE OF G(gamma_kc)
0202             % Gg = Gg + (g_new - g_old)*1_idelta_gamma_kc';
0203             % diff_G_in_out = G(:,idelta3)-G(:,idelta);
0204             % iGg_new = iGg - (iGg*diff_G_in_out/(1+iGg(idelta_gamma_kc,:)*diff_G_in_out))*(iGg(idelta_gamma_kc,:));
0205             
0206             % UPDATE QR FACTORS FOR G(gamma_kc)
0207             diff_G_in_out = G(:,idelta3)-G(:,idelta);
0208             loc_vector = zeros(length(diff_G_in_out),1);
0209             loc_vector(idelta_gamma_kc) = 1;
0210             Q0 = Qg; R0 = Rg;
0211             [Qg Rg] = qrupdate(Q0,R0,diff_G_in_out,loc_vector);
0212             
0213             x_old = xk_1;
0214             c_old = ck_1;
0215             lambda_old = lambdak_1;
0216             gamma_old = gammak_new;
0217 
0218             if ~isempty(find(idelta3 == i_nu, 1))
0219                 lambda_old(idelta3) = (epsilon_old+delta)*lambdak_1(idelta3);
0220                 i_nu = setdiff(i_nu,idelta3);
0221                 if isempty(i_nu)
0222                     epsilon = 1;
0223                 else
0224                     epsilon = delta + epsilon_old;
0225                 end
0226             else
0227                 epsilon = delta + epsilon_old;
0228             end
0229 
0230             if epsilon >=1
0231                 done = 1;
0232                 
0233                 A_old = G';
0234                 y_old = q;
0235                 m_old = m_new;
0236             end
0237             epsilon_old = epsilon;
0238 
0239         else
0240             new_delta = 1-epsilon_old;
0241             lambdak_1 = lambda_k + new_delta*del_lambda_vec;
0242             lambda_old = [lambdak_1];
0243             gamma_old = gammak_new;
0244 
0245             done = 1;
0246             
0247             A_old = G';
0248             y_old = q;
0249             m_old = m_new;
0250         end
0251     end
0252     figure(2);  clf; 
0253     subplot(211);plot(x_old); hold on; plot(x,'r.');
0254     title(['Estimated signal with ', num2str(iter_mat+m_u-1+m_st), ' measurements in the presence of ',num2str(T),  '  errors']);
0255     legend('Estimation', 'Original');
0256     subplot(212);plot(c_old); hold on; plot(c_original,'r.');
0257     title(['Estimated and original sparse error vectors']);
0258     legend('Estimation', 'Original'); drawnow
0259     
0260     total_iter = total_iter+iter_in;
0261     disp(['No. of errors = ',num2str(nnz(c_original)), ', Current no. of measurements = ',num2str(iter_mat+m_u-1+m_st), ', Homotopy iterations = ', num2str(iter_in), ',  Running sum of iterations = ', num2str(total_iter)]);
0262             
0263     if all_done == 1
0264         % Need L1-magic package for this verification.
0265         % x02 = G'\q;
0266         %  [xp2 lamf1 lamf2] = l1decode_pd2(x02, G', [], q, 1e-12, 30);
0267         break;
0268     end
0269 end
0270 x_out = x_old;
0271 figure(2);  clf;
0272 subplot(211);plot(x_old); hold on; plot(x,'r.');
0273 title(['Estimated signal with ', num2str(iter_mat+m_u-1+m_st), ' measurements in the presence of ',num2str(T),  '  errors']);
0274 legend('Estimation', 'Original');
0275 subplot(212);plot(c_old); hold on; plot(c_original,'r.');
0276 title(['Estimated and original sparse error vectors']);
0277 legend('Estimation', 'Original'); drawnow