l1Decode_homotopy_fast

 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 
0022 close all
0023 clear; clc
0024 
0025 % % load fixed random states
0026 % load RandomStates
0027 % rand('state', s_rand);
0028 % randn('state', s_randn);
0029 
0030 % source length
0031 N = 256;
0032 
0033 % codeword length
0034 M = 4*N;
0035 
0036 % number of sparse errors
0037 T = round(.2*N);
0038 
0039 % Coding matrix
0040 Orth_mat = randn(M,M);
0041 Orth_mat = randsrc(M,M);
0042 % Orth_mat = orth(Orth_mat);
0043 G_mat = Orth_mat(:,1:N);
0044 
0045 % Homotopy setup
0046 m_st = N;   % number of measurements to start with
0047 A = G_mat(1:m_st,:);
0048 B_stack = G_mat(m_st+1:end,:);
0049 
0050 % Source word
0051 x = randn(N,1);
0052 
0053 % Channel: perturb T randomly chosen entries
0054 q = randperm(m_st-1);
0055 c = zeros(m_st,1);
0056 c(q(1:T)) = randn(T,1);
0057 y = A*x-c;
0058 err_loc = (rand(M-N,1)>=.9);
0059 err_new = err_loc.*randn(M-N,1);
0060 w_vec = B_stack*x-err_new;
0061 
0062 
0063 % initial estimate of x
0064 % l1_min or direct solution (LS) depending on n_st
0065 iA = inv(A);
0066 x0 = iA*y;
0067 e0 = A*x0-y;
0068 lambda_0 = zeros(m_st,1);
0069 
0070 A_old = A;
0071 y_old = y;
0072 x_old = x0;
0073 c_old = e0;
0074 lambda_old = lambda_0;
0075 gamma_old = []; %find(A*x0-y) or find(abs(lambda_0)==1);
0076 m_old = m_st;
0077 all_done = 0;
0078 cond_list = [];
0079 cond_iter = 1;
0080 gamma_kc_new = [1:N]';
0081 iGg_new = iA';
0082 
0083 % number of new measurements added at a time.
0084 m_u = 8; 
0085 
0086 total_iter = 0;
0087 
0088 for iter_mat = 1:m_u:M-m_st
0089     m_new = m_old + m_u;
0090     x_k = x_old;
0091     lambda_k = lambda_old;
0092     gamma_k = gamma_old;
0093     %     gamma_kc = gamma_kc_old;
0094 
0095     B = B_stack(iter_mat:iter_mat+m_u-1,:);
0096     G = [A_old' B'];
0097     w = w_vec(iter_mat:iter_mat+m_u-1);
0098     q = [y_old; w];
0099     
0100     c_original = [c; err_new(1:iter_mat+m_u-1)];
0101     c_old = [c_old; B*x_k-w];
0102 
0103     if abs(abs(B*x_k-w))>=10*eps
0104         nu = sign(B*x_k-w);
0105         i_nu = [m_new-m_u+1:m_new]';
0106         gammak_new = [gamma_k; i_nu];
0107     else
0108         stp = 1; % Something is wrong! (degenracy)
0109         break;
0110     end
0111 
0112     lambda_old = [lambda_k; nu]; % xi_old
0113 
0114     done = 0;
0115     epsilon_old = 0;
0116     iter_in = 0;
0117     while ~done
0118         iter_in = iter_in+1;
0119         %[iter_mat iter_in]
0120     
0121         gamma_k = gammak_new;
0122         x_k = x_old;
0123         c_k = c_old;
0124         lambda_k = lambda_old; % xi_k
0125 
0126         gamma_kc = gamma_kc_new;
0127 
0128         % INVERSE UPDATE
0129         %iGg = inv(G(:,gamma_kc));
0130         iGg = iGg_new;
0131         del_lambda = -iGg*(G(:,i_nu)*sign(lambda_k(i_nu))); % del_xi
0132 
0133         del_lambda_vec = zeros(m_new,1); % del_xi_vec
0134         del_lambda_vec(gamma_kc) = del_lambda;
0135 
0136         % find incoming elements in support of error vector
0137         constr1 = (1-lambda_k(gamma_kc))./del_lambda_vec(gamma_kc);
0138         constr1_pos = constr1(constr1>2*eps);
0139         delta1 = min(constr1_pos);
0140         idelta1 = gamma_kc(find(constr1==delta1));
0141 
0142         constr2 = -(1+lambda_k(gamma_kc))./del_lambda_vec(gamma_kc);
0143         constr2_pos = constr2(constr2>2*eps);
0144         delta2 = min(constr2_pos);
0145         idelta2 = gamma_kc(find(constr2==delta2));
0146 
0147         if delta1>delta2
0148             delta = delta2;
0149             idelta = idelta2;
0150         else
0151             delta = delta1;
0152             idelta = idelta1;
0153         end
0154 
0155         if epsilon_old+delta <=1
0156             lambdak_1 = lambda_k + delta*del_lambda_vec;
0157             lambdak_1(idelta) = sign(lambdak_1(idelta));
0158             gammak_temp = [gamma_k; idelta];
0159 
0160             % outgoing element from the support of
0161             uz_vec = zeros(N,1);
0162             idelta_gamma_kc = find(idelta==gamma_kc);
0163             uz_vec(idelta_gamma_kc) = sign(lambdak_1(idelta));
0164             delx = iGg'*uz_vec;
0165 
0166             delc= G'*delx;
0167             constr3 = -(c_k(gamma_k)./delc(gamma_k));
0168             constr3_pos = constr3(constr3>2*eps);
0169             delta3 = min(constr3_pos);
0170             idelta3 = gamma_k(find(constr3==delta3));
0171             xk_1 = x_k + delta3*delx;
0172             ck_1 = c_k + delta3*delc;
0173             ck_1(idelta3) = 0;
0174 
0175             out_lambda = idelta3;
0176 
0177             % See if signal is decoded exactly...
0178             if [m_new-N~=length(find(abs(ck_1)>=1e-10))]
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                 break; break; break;
0186             end
0187             gammak_new = setdiff(gammak_temp,idelta3);
0188             gamma_kc_new = gamma_kc;
0189             gamma_kc_new(idelta_gamma_kc) = idelta3;
0190 
0191             % UPDATE INVERSE OF G(gamma_kc)
0192             % Gg = Gg + (g_new - g_old)*1_idelta_gamma_kc';
0193             diff_G_in_out = G(:,idelta3)-G(:,idelta);
0194             iGg_new = iGg - (iGg*diff_G_in_out/(1+iGg(idelta_gamma_kc,:)*diff_G_in_out))*(iGg(idelta_gamma_kc,:));
0195 
0196             x_old = xk_1;
0197             c_old = ck_1;
0198             lambda_old = lambdak_1;
0199             gamma_old = gammak_new;
0200 
0201             if ~isempty(find(idelta3 == i_nu, 1))
0202                 lambda_old(idelta3) = (epsilon_old+delta)*lambdak_1(idelta3);
0203                 i_nu = setdiff(i_nu,idelta3);
0204                 if isempty(i_nu)
0205                     epsilon = 1;
0206                 else
0207                     epsilon = delta + epsilon_old;
0208                 end
0209             else
0210                 epsilon = delta + epsilon_old;
0211             end
0212 
0213             if epsilon >=1
0214                 done = 1;
0215 
0216                 A_old = G';
0217                 y_old = q;
0218                 m_old = m_new;
0219             end
0220             epsilon_old = epsilon;
0221 
0222         else
0223             new_delta = 1-epsilon_old;
0224             lambdak_1 = lambda_k + new_delta*del_lambda_vec;
0225             lambda_old = [lambdak_1];
0226             gamma_old = gammak_new;
0227 
0228             done = 1;
0229 
0230             A_old = G';
0231             y_old = q;
0232             m_old = m_new;
0233         end
0234     end
0235     figure(2);  clf;
0236     subplot(211);plot(x_old); hold on; plot(x,'r.');
0237     title(['Estimated signal with ', num2str(iter_mat+m_u-1+m_st), ' measurements in the presence of ',num2str(T),  '  errors']);
0238     legend('Estimation', 'Original');
0239     subplot(212);plot(c_old); hold on; plot(c_original,'r.');
0240     title(['Estimated and original sparse error vectors']);
0241     legend('Estimation', 'Original'); drawnow
0242 
0243     total_iter = total_iter+iter_in;
0244     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)]);
0245             
0246     if all_done == 1
0247         % Need L1-magic package for this verification.
0248         % x02 = G'\q; %inv(G*G')*G*q;
0249         % [xp2 lamf1 lamf2] = l1decode_pd2(x02, G', [], q, 1e-12, 30);
0250         break;
0251     end
0252 end
0253 x_out = x_old;
0254 figure(2);  clf;
0255 subplot(211);plot(x_old); hold on; plot(x,'r.');
0256 title(['Estimated signal with ', num2str(iter_mat+m_u-1+m_st), ' measurements in the presence of ',num2str(T),  '  errors']);
0257 legend('Estimation', 'Original');
0258 subplot(212);plot(c_old); hold on; plot(c_original,'r.');
0259 title(['Estimated and original sparse error vectors']);
0260 legend('Estimation', 'Original'); drawnow