Description of DynamicSeq_REC

0001 % DynamicSeq_REC_function.m
0002 % Robust ell_1 decoding with homotopy update scheme for new measurements.
0003 %
0004 % [cp_h gamma_h iter th] = DynamicSeq_REC_function(P, s, iPgP, cp_old, gamma_h, gamma_n, pk, tau, M, m_u, maxiter);
0005 %
0006 % Outputs:
0007 % cp_h: sparse error estimate
0008 % gamma_h: support of sparse errors
0009 % iter: number of homotopy steps taken
0010 %
0011 % Inputs
0012 % P: I-F*inv(F'F)*F'
0013 % s: [y; w]
0014 % iPgP: inv(PgP) where PgP = P(gamma_xh,gamma_xh)
0015 % cp_old: Old estimate of error (includes the estimate for new measurements)
0016 % gamma_h: support of cp_old
0017 % gamma_n: support of estimated errors in new measurements
0018 % pk: P*(cp_old-s)
0019 % tau: regularization parameter
0020 % M: Number of old measurements
0021 % m_u: Number of new measurements
0022 % maxiter: Maximum allowed iterations
0023 %
0024 % Created by Salman Asif @ Georgia Tech
0025 % February 2009
0026 %
0027 %-------------------------------------------+
0028 % Copyright (c) 2009.  Muhammad Salman Asif
0029 %-------------------------------------------+
0030 
0031 
0032 function [cp_h gamma_h iter th] = REC_homotopy_function(P, s, iPgP, cp_old, gamma_h, gamma_n, pk, tau, M, m_u, maxiter);
0033 
0034 t0 = cputime;
0035 cp_h = cp_old;
0036 epsilon = 0;
0037 e0 = 0;
0038 pk_old = pk; % last indices will be zero at the start.
0039 
0040 pk_old([gamma_h]) = sign(pk_old([gamma_h]))*tau;
0041 pk_old(gamma_n) = pk_old(gamma_n)*e0*tau;
0042 
0043 done = 0;
0044 iter = 0;
0045 tu0 = cputime;
0046 outc = [];
0047 
0048 while iter < maxiter
0049     iter = iter +1 ;
0050 
0051     PgP = P(gamma_h,gamma_h);
0052 
0053     c_k = cp_h;
0054     len_e = length(gamma_h)-length(gamma_n);
0055     z_d = sign(c_k(gamma_n));
0056     igamma_n_gamma_h = zeros(length(gamma_n),1);
0057     for ii = 1:length(gamma_n)
0058         igamma_n_gamma_h(ii) = find(gamma_n(ii) == gamma_h);
0059     end
0060     temp_vec = zeros(length(gamma_h),1);
0061     temp_vec(igamma_n_gamma_h) = z_d;
0062     delc = -iPgP*temp_vec;
0063     % figure(102); imagesc(iPgP-inv(P(gamma_h,gamma_h))); colorbar; shg
0064     % delc = -inv(P(gamma_h,gamma_h))*temp_vec;
0065     
0066     delc_vec = zeros(M+m_u,1);
0067     delc_vec(gamma_h) = delc;
0068     pk = pk_old;
0069     dk = P(:,gamma_h)*delc;
0070     outc = [];
0071     temp_gamma = zeros(M+m_u,1);
0072     temp_gamma(gamma_h) = gamma_h;
0073     gamma_hc = find([1:M+m_u]' ~= temp_gamma);
0074 
0075     b_constr1 = (tau-pk(gamma_hc))./dk(gamma_hc);
0076     b_constr2 = (tau+pk(gamma_hc))./-dk(gamma_hc);
0077     b_constr3 = (-c_k(gamma_h)./delc_vec(gamma_h));
0078     itheta_1 = find(b_constr1>2*eps);
0079     itheta_2 = find(b_constr2>2*eps);
0080     itheta_3 = find(b_constr3>2*eps);
0081     theta1 = min(b_constr1(itheta_1));
0082     theta2 = min(b_constr2(itheta_2));
0083     theta3 = min(b_constr3(itheta_3));
0084     if theta1>theta2
0085         theta = theta2;
0086         itheta = gamma_hc(find(b_constr2==theta2));
0087         flag = 1;
0088     else
0089         theta = theta1;
0090         itheta = gamma_hc(find(b_constr1==theta1));
0091         flag = 1;
0092     end
0093 
0094     if theta3 < theta
0095         theta = theta3;
0096         itheta = gamma_h(find(b_constr3==theta3));
0097         outc = itheta;
0098         flag = 0;
0099     end
0100     
0101     if length(gamma_h)>=M+m_u-128 & flag ==1
0102         stp = 1;
0103     end
0104 
0105     epsilon = theta/tau + e0;
0106 
0107     if epsilon < 0
0108         epsilon= inf;
0109     end
0110     gamma_h_old = gamma_h;
0111     gamma_n_old = gamma_n;
0112     cp_old = c_k;
0113     cp_h = c_k+theta*delc_vec;
0114 
0115     if epsilon >=1
0116         theta_end = tau*(1-e0);
0117         cp_h = cp_old + theta_end*delc_vec;
0118         gamma_h = gamma_h_old;
0119         pk_old = pk+theta_end*dk;
0120         break;
0121     end
0122 
0123     if flag == 1
0124         gamma_h = [gamma_h; itheta];
0125         cp_h(itheta) = 0;
0126         PgP_mod = P(gamma_h,gamma_h);
0127 
0128         iPgP = update_inverse(PgP_mod, iPgP,1);
0129         PgP = PgP_mod;
0130 
0131     else
0132         outc_index = find(gamma_h==itheta);
0133         gamma_h(outc_index) = gamma_h(end);
0134         gamma_h(end) = itheta;
0135         gamma_h = gamma_h(1:end-1);
0136 
0137         rowi = outc_index; % ith row of P is swapped with last row (out_c)
0138         colj = outc_index; % jth column of P is swapped with last column (out_c)
0139         PgP_ij = PgP;
0140         temp_row = PgP_ij(rowi,:);
0141         PgP_ij(rowi,:) = PgP_ij(end,:);
0142         PgP_ij(end,:) = temp_row;
0143         temp_col = PgP_ij(:,colj);
0144         PgP_ij(:,colj) = PgP_ij(:,end);
0145         PgP_ij(:,end) = temp_col;
0146         iPgP_ij = iPgP;
0147         temp_row = iPgP_ij(colj,:);
0148         iPgP_ij(colj,:) = iPgP_ij(end,:);
0149         iPgP_ij(end,:) = temp_row;
0150         temp_col = iPgP_ij(:,rowi);
0151         iPgP_ij(:,rowi) = iPgP_ij(:,end);
0152         iPgP_ij(:,end) = temp_col;
0153 
0154         PgP = PgP_ij(1:end-1,1:end-1);
0155         iPgP = update_inverse(PgP_ij, iPgP_ij,2);
0156 
0157         cp_h(itheta) = 0;
0158         if length(find(itheta == gamma_n))
0159             outn_index = find(gamma_n==itheta);
0160             gamma_n(outn_index) = gamma_n(end);
0161             gamma_n(end) = itheta;
0162             gamma_n = gamma_n(1:end-1);
0163             % gamma_n = setdiff(gamma_n,itheta);
0164             if isempty(gamma_n)
0165                 e0_old = epsilon;
0166                 epsilon = 1;
0167                 pk_old = pk+theta*dk;
0168                 break;
0169             end
0170         end
0171     end
0172 
0173     e0 = (theta/tau+e0);
0174 
0175     pk_old = pk+theta*dk;
0176     pk_old([gamma_h; itheta]) = sign(pk_old([gamma_h; itheta]))*tau;
0177     pk_old(gamma_n_old) = pk_old(gamma_n_old)*e0;
0178 
0179 end
0180 th= cputime-t0;
DynamicSeq_REC_function

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
