Simulation_DynamicSeq_REC

 Simulation to test the performance of REC

0001 % Simulation to test the performance of REC
0002 
0003 clear
0004 close all
0005 
0006 % % load fixed random states
0007 % load RandomStates
0008 % rand('state', s_rand);
0009 % randn('state', s_randn);
0010 
0011 sim_runs = 10;
0012 additional_elements = 100;
0013 
0014 % No. of new measurements added at a time
0015 m_u = 10;
0016 inner_iterations = additional_elements/m_u;
0017 
0018 stack_time = zeros(sim_runs*inner_iterations,2);
0019 stack_iter = zeros(sim_runs*inner_iterations,5);
0020 
0021 
0022 for outer_iter = 1:sim_runs
0023     % % source length
0024     N = 150;
0025     % codeword length
0026     M = 2*N;
0027     M_st = M;
0028     
0029     % number of perturbations
0030     max_sparse_errors = .2; % e.g., 20% of M
0031     T = round(max_sparse_errors*M_st);
0032     
0033     T = round(N/2.5);
0034     
0035     % coding matrix
0036     Orth_mat = randn(M,M);
0037     A = randn(M,N)/sqrt(N);
0038     A = orth(A);
0039 
0040     % Annihilating projection matrix
0041     AtA = A'*A;
0042     iAtA = inv(AtA);
0043     AiAtA = A*iAtA;
0044     AiAtAAt = AiAtA*A';
0045     Q = eye(M)-AiAtAAt;
0046 
0047     % source word
0048     x = (randn(N,1));
0049 
0050     % channel: perturb T randomly chosen entries
0051     q = randperm(M);
0052     
0053     % Introduce sparse errors
0054     e = zeros(M,1);
0055     % e(q(1:T)) = randsrc(T,1); % Arbitrary sparse errors
0056     e(q(1:T)) = -A(q(1:T),:)*x; % Erasures!
0057 
0058     % Small noise
0059     x0 = randn(N,1);
0060     Ax0 = A*x0;
0061     sigma = median(abs(Ax0))/20; % control the power in small noise
0062     q_y = randn(M,1)*sigma;
0063 
0064     % Received codeword
0065     y = A*x+e+q_y;
0066 
0067     % Regularization parameter
0068     tau = 0.01;%*max(abs(Q*y)); % l1_ls
0069 %     if sigma>0
0070 %         tau = sigma * sqrt(log(N)*2); % BPDN
0071 %         % tau = max(abs(Q'*q_y)); % ideal ???
0072 %     end
0073 
0074     % Data recovery
0075     Qy = Q*y;
0076     [ep, gamma_e, ep_iter, t1] =  BPDN_homotopy_function(Q, Qy, tau, 4*M);
0077     % xp = inv(A'*A)*A'*(y-ep);
0078     xp = AiAtA'*(y-ep);
0079 
0080     tolA_h = tau*sum(abs(ep))+1/2*(norm(Q*(ep-y)))^2;
0081 
0082     for inn_iter = 1:inner_iterations
0083         [outer_iter inn_iter]
0084         if inn_iter > 1
0085 %             % Update the parameters from the previous run of inner iterations
0086 %             AtA = AtA_n;
0087 %             iAtA = iAtA_n;
0088 %             AiAtA = AiAtA_n;
0089 %             AiAtAAt = AiAtAAt_n;
0090 %             Q = Q_n;
0091 %             y = y_n;
0092 %             A = A_n;
0093 %             q_y = q_y_n;
0094 %             e = e_n;
0095 %             xp = xp_n;
0096 %             ep = ep_n;
0097 %             M = M_n;
0098 %             gamma_e = gamma_e_n;
0099 %             e_BB_mono = e_BB_mono_n;
0100             % Update the parameters for next run of inner iterations
0101             AtA = FtF;
0102             iAtA = iFtF;
0103             AiAtA = FiFtF;
0104             AiAtAAt = FiFtFFt;
0105             Q = P;
0106             y = s;
0107             A = F;
0108             q_y = q_yw;
0109             e = c;
0110             xp = xp_h;
0111             ep = cp_h;
0112             M = M+m_u;
0113             gamma_e = gamma_h;
0114         end
0115         
0116         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0117         % Setup for adding m_u new observations %
0118         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0119 
0120         % No. of new observations
0121         % m_u =10;
0122 
0123         b = randn(m_u,N)/sqrt(N); % new rows in coding matrix
0124         max_sparse_errors_new = max_sparse_errors;
0125         max_sparse_errors_new = .2;
0126         T_new = round((m_u*max_sparse_errors_new)*(rand>=.5)+rand*.5); % new gross/sparse errors
0127         d = zeros(m_u,1);
0128         q_new = randperm(m_u);
0129         
0130         % Introduce gross errors in the new elements of the codeowrd
0131         % d(q_new(1:T_new)) = randsrc(T_new,1); % Arbitrary sparse errors
0132         d(q_new(1:T_new)) = -b(q_new(1:T_new),:)*x; % Erasures!
0133 
0134         % small noise in the new measurements
0135         q_w = randn(m_u,1)*sigma;
0136 
0137         w = b*x+d+q_w;
0138         F = [A; b];
0139         s = [y;w];
0140         c = [e; d];
0141         q_yw = [q_y; q_w];
0142 
0143 
0144         iAtAbt = iAtA*b';
0145         biAtAbt = b*iAtAbt;
0146         S_biAtAbt = inv(eye(m_u)+biAtAbt);
0147         AiAtAbt = AiAtA*b';
0148         iAtAbt_S = (iAtAbt*S_biAtAbt);
0149         AiAtAbt_S = A*iAtAbt_S;
0150 
0151         FtF = AtA+b'*b; %F'*F;
0152         iFtF = iAtA - iAtAbt_S*iAtAbt'; % inv(FtF);
0153         FiFtF = [AiAtA - AiAtAbt_S*iAtAbt'; iAtAbt'-biAtAbt*S_biAtAbt*iAtAbt']; % F*iFtF;
0154         FiFtFFt = [AiAtAAt - AiAtAbt_S*AiAtAbt' AiAtAbt-AiAtAbt_S*biAtAbt; AiAtAbt'-biAtAbt*AiAtAbt_S' biAtAbt-biAtAbt*S_biAtAbt*biAtAbt]; %FiFtF*F';
0155         P = eye(M+m_u)-FiFtFFt;
0156 
0157         dp = w-b*xp;
0158         z_d = sign(dp);
0159         cp_h = [ep; dp];
0160         gamma_n = M+find(abs(dp)>2*eps);
0161         gamma_n_old = gamma_n;
0162         gamma_h = [gamma_e; gamma_n];
0163         epsilon = 0;
0164         e0 = 0;
0165 
0166         QgQ = Q(gamma_e,gamma_e);
0167         PgP = P(gamma_h,gamma_h);
0168         uQ = AiAtAbt_S(gamma_e,:);
0169         vQ = AiAtAbt(gamma_e,:);
0170         QgQ_update = QgQ + uQ*vQ';
0171         PgP_update = [[QgQ_update; P(gamma_n,gamma_e)] P(gamma_h,gamma_n)];
0172         PgP = PgP_update;
0173 
0174         iQgQ = inv(QgQ);
0175         iQgQ_update = iQgQ - (iQgQ*uQ)*(inv(eye(m_u)+vQ'*iQgQ*uQ))*(vQ'*iQgQ);
0176 
0177         P11 = QgQ_update;
0178         P12 = P(gamma_e,gamma_n);
0179         P21 = P12';
0180         P22 = P(gamma_n,gamma_n);
0181         S_P = inv(P22-P21*iQgQ_update*P12);
0182         iP11_P12 = iQgQ_update*P12;
0183         iPgP_update = [iQgQ_update+(iP11_P12*S_P)*iP11_P12' -iP11_P12*S_P; -S_P*iP11_P12' S_P];
0184         iPgP = iPgP_update;
0185 
0186         pk_old = P*(cp_h-s);
0187         [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);
0188         xp_h = FiFtF'*(s-cp_h);
0189 
0190         % % Check the solution using homotopy from scratch.
0191         [cp2, gamma_c2, cp2_iter, t2] = BPDN_homotopy_function(P, P*s, tau, 4*M);
0192         
0193         stack_time((outer_iter-1)*inner_iterations+inn_iter,:) = [t2 th];
0194         stack_iter((outer_iter-1)*inner_iterations+inn_iter,:) = [cp2_iter cp_h_iter T_new norm(x-xp_h)/norm(x) norm(cp_h-cp2)];
0195 %         iter_count = (outer_iter-1)*inner_iterations+inn_iter
0196     end
0197 end
0198 
0199 ['BPDN      REC']
0200 Average_time = mean(stack_time,1)
0201 Average_iterations = mean(stack_iter(:,1:2))