DynamicSeq_REC_demo

 DynamicSeq_REC_demo.m
 Robust ell_1 decoding with homotopy update scheme for new measurements.

 Created by Salman Asif @ Georgia Tech
 February 2009

 Model changes in the error patterns with new measurements and observe the
 effect on the solution.

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

0001 % DynamicSeq_REC_demo.m
0002 % Robust ell_1 decoding with homotopy update scheme for new measurements.
0003 %
0004 % Created by Salman Asif @ Georgia Tech
0005 % February 2009
0006 %
0007 % Model changes in the error patterns with new measurements and observe the
0008 % effect on the solution.
0009 %
0010 %-------------------------------------------+
0011 % Copyright (c) 2009.  Muhammad Salman Asif
0012 %-------------------------------------------+
0013 
0014 clear
0015 close all
0016 
0017 % % load fixed random states
0018 % load RandomStates
0019 % rand('state', s_rand);
0020 % randn('state', s_randn);
0021   
0022 % % source length
0023 N = 256;
0024 % codeword length
0025 M = 3*N;
0026 
0027 % number of perturbations
0028 T = round(.2*M);
0029 
0030 % coding matrix
0031 Orth_mat = randn(M,M);
0032 G = Orth_mat(:,1:N);
0033 A = randn(M,N)/sqrt(M);
0034 
0035 % Annihilating projection matrix
0036 AtA = A'*A;
0037 iAtA = inv(AtA);
0038 AiAtA = A*iAtA;
0039 AiAtAAt = AiAtA*A';
0040 Q = eye(M)-AiAtAAt;
0041 
0042 % source word
0043 x = randn(N,1);
0044 
0045 % channel: perturb T randomly chosen entries
0046 q = randperm(M);
0047 e = zeros(M,1);
0048 e(q(1:T)) = randsrc(T,1);
0049 
0050 % Small noise
0051 x0 = randn(N,1);
0052 Ax0 = A*x0;
0053 sigma = median(abs(Ax0))/20; % control the power in small noise
0054 q_y = randn(M,1)*sigma;
0055 
0056 % Received codeword
0057 y = A*x+e+q_y;
0058 figure(1); clf; plot([y e q_y]); %hold on; plot(e,'g-.'); plot(q_y,'r--');
0059 title('Codeword, sparse errors and small noise');
0060 legend('Codeword', 'Sparse', 'Noise'); shg
0061 
0062 % Regularization parameter
0063 tau = 0.01*max(abs(Q*y)); % l1_ls
0064 % if sigma>0
0065 %     tau = sigma * sqrt(log(N)*2); % BPDN
0066 %     % tau = max(abs(Q'*q_y)); % ideal ???
0067 % end
0068 
0069 % Data recovery
0070 t1_s = cputime;
0071 Qy = Q*y;
0072 [ep, gamma_e, ep_iter] =  BPDN_homotopy_function(Q, Qy, tau, 4*M);
0073 t1_e = cputime;
0074 t1 = t1_e-t1_s;
0075 % xp = inv(A'*A)*A'*(y-ep);
0076 xp = AiAtA'*(y-ep);
0077 
0078 figure(2); clf; plot(ep); hold on; plot(e,'ro');
0079 title(['Estimated and original sparse errors with ', num2str(M),' measurements']);
0080 legend('Estimated', 'Original'); shg
0081 
0082 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0083 % Setup for adding m_u new observations %
0084 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0085 
0086 % No. of new observations
0087 m_u =10;
0088 
0089 b = randn(m_u,N)/sqrt(M); % new rows in coding matrix
0090 T_new = randint(1,1,m_u/5); % new gross/sparse errors
0091 d = zeros(m_u,1);
0092 q_new = randperm(m_u);
0093 d(q_new(1:T_new)) = randsrc(T_new,1);
0094 q_w = randn(m_u,1)*sigma;
0095 
0096 w = b*x+d+q_w;
0097 F = [A; b];
0098 s = [y;w];
0099 c = [e; d];
0100 q_yw = [q_y; q_w];
0101 
0102 iAtAbt = iAtA*b';
0103 biAtAbt = b*iAtAbt;
0104 S_biAtAbt = inv(eye(m_u)+biAtAbt);
0105 AiAtAbt = AiAtA*b';
0106 iAtAbt_S = (iAtAbt*S_biAtAbt);
0107 AiAtAbt_S = A*iAtAbt_S;
0108 
0109 FtF = AtA+b'*b; %F'*F;
0110 iFtF = iAtA - iAtAbt_S*iAtAbt'; % inv(FtF);
0111 FiFtF = [AiAtA - AiAtAbt_S*iAtAbt'; iAtAbt'-biAtAbt*S_biAtAbt*iAtAbt']; % F*iFtF;
0112 FiFtFFt = [AiAtAAt - AiAtAbt_S*AiAtAbt' AiAtAbt-AiAtAbt_S*biAtAbt; AiAtAbt'-biAtAbt*AiAtAbt_S' biAtAbt-biAtAbt*S_biAtAbt*biAtAbt]; %FiFtF*F';
0113 P = eye(M+m_u)-FiFtFFt;
0114 
0115 dp = w-b*xp;
0116 z_d = sign(dp);
0117 cp_h = [ep; dp];
0118 gamma_n = M+find(abs(dp)>2*eps);
0119 gamma_n_old = gamma_n;
0120 gamma_h = [gamma_e; gamma_n];
0121 epsilon = 0;
0122 e0 = 0;
0123 
0124 QgQ = Q(gamma_e,gamma_e);
0125 PgP = P(gamma_h,gamma_h); 
0126 uQ = AiAtAbt_S(gamma_e,:);
0127 vQ = AiAtAbt(gamma_e,:);
0128 QgQ_update = QgQ + uQ*vQ';
0129 PgP_update = [[QgQ_update; P(gamma_n,gamma_e)] P(gamma_h,gamma_n)];
0130 PgP = PgP_update;
0131 
0132 iQgQ = inv(QgQ);
0133 iQgQ_update = iQgQ - (iQgQ*uQ)*(inv(eye(m_u)+vQ'*iQgQ*uQ))*(vQ'*iQgQ);
0134 
0135 P11 = QgQ_update;
0136 P12 = P(gamma_e,gamma_n);
0137 P21 = P12';
0138 P22 = P(gamma_n,gamma_n);
0139 S_P = inv(P22-P21*iQgQ_update*P12);
0140 iP11_P12 = iQgQ_update*P12;
0141 iPgP_update = [iQgQ_update+(iP11_P12*S_P)*iP11_P12' -iP11_P12*S_P; -S_P*iP11_P12' S_P];
0142 iPgP = iPgP_update;
0143 
0144 pk_old = P*(cp_h-s); % last indices will be zero.
0145 
0146 % REC homotopy
0147 [cp_h gamma_h cp_h_iter th] = DynamicSeq_REC_function(P, s, iPgP, cp_h, gamma_h, gamma_n, pk_old, tau, M, m_u, 4*M);
0148 figure(3); clf; plot(cp_h); hold on; plot(c,'ro');
0149 title(['Estimated and original sparse errors after new ', num2str(m_u),' measurements']);
0150 legend('Estimated', 'Original'); shg
0151 
0152 % % Check the solution using cvx
0153 % cvx_begin
0154 %     cvx_precision high
0155 %     variables cp_cvx(M+m_u) zp_cvx(M+m_u) ;
0156 %     minimize( tau*(norm(cp_cvx,1))+ .5*(sum(square(zp_cvx))))
0157 %     subject to
0158 %     zp_cvx == P*(cp_cvx-s)
0159 % cvx_end
0160 %
0161 % % figure(103); clf; plot(cp_h); hold on; plot(cp2,'kx'); plot([cp_cvx; dp_cvx],'r.')
0162 % figure(10); hold on; plot(cp_cvx-cp_h); shg
0163 
0164 % Check the solution using homotopy from scratch.
0165 [cp2, gamma_c2, cp2_iter, t2] = BPDN_homotopy_function(P, P*s, tau, 4*M);
0166 
0167 
0168 pk2 = P*(cp2-s);
0169 pk_h = P*(cp_h-s);
0170 
0171 % [old_homotopy new_homotopy update_homotopy]
0172 disp(' ');
0173 disp('old_homotopy , new_homotopy , update_homotopy')
0174 time_table = [t1 t2 th]
0175 iter_table = [ep_iter cp2_iter cp_h_iter]
0176 
0177 figure(4); subplot(2,1,1); hold on; plot(cp_h-cp2);
0178 title('Difference between the solution of BPDN homotopy and REC update')
0179 subplot(2,1,2); hold on; plot(pk2-pk_h); drawnow
0180 title('Difference between the primal constraints of BPDN homotopy and REC update')