update_primal

 update_primal.m

 This function computes the minimum step size in the primal update direction and
 finds change in the primal or dual support with that step.

 Inputs:
 gamma_x - current support of x
 gamma_lambda - current support of lambda
 z_x - sign sequence of x
 z_lambda - sign sequence of lambda
 del_x_vec - primal update direction
 pk
 dk
 epsilon - current value of epsilon
 out_lambda - element removed from support of lambda in previous step (if any)

 Outputs:
 i_delta - index corresponding to newly active primal constraint (new_lambda)
 out_x - element in x shrunk to zero
 delta - primal step size
 chk_x - 1  an element is removed from support of x
         0  a new element enters the support of lambda

 Written by: Salman Asif, Georgia Tech
 Email: sasif@ece.gatech.edu

0001 % update_primal.m
0002 %
0003 % This function computes the minimum step size in the primal update direction and
0004 % finds change in the primal or dual support with that step.
0005 %
0006 % Inputs:
0007 % gamma_x - current support of x
0008 % gamma_lambda - current support of lambda
0009 % z_x - sign sequence of x
0010 % z_lambda - sign sequence of lambda
0011 % del_x_vec - primal update direction
0012 % pk
0013 % dk
0014 % epsilon - current value of epsilon
0015 % out_lambda - element removed from support of lambda in previous step (if any)
0016 %
0017 % Outputs:
0018 % i_delta - index corresponding to newly active primal constraint (new_lambda)
0019 % out_x - element in x shrunk to zero
0020 % delta - primal step size
0021 % chk_x - 1  an element is removed from support of x
0022 %         0  a new element enters the support of lambda
0023 %
0024 % Written by: Salman Asif, Georgia Tech
0025 % Email: sasif@ece.gatech.edu
0026 
0027 function [i_delta, out_x, delta, chk_x] = update_primal(gamma_x, gamma_lambda, z_x, x_k, del_x_vec, pk, dk, epsilon, out_lambda);
0028 
0029 N = length(x_k);
0030 
0031 % gamma_lc = setdiff([1:N]', [gamma_lambda; out_lambda]); WRONG
0032 % check out_lambda as well, that is if outgoing lambda switches sign in just one step
0033 temp_gamma = zeros(N,1);
0034 temp_gamma(gamma_lambda) = gamma_lambda;
0035 gamma_lc = find([1:N]' ~= temp_gamma);
0036 
0037 delta1_constr = (epsilon-pk(gamma_lc))./(1+dk(gamma_lc));
0038 delta1_pos_ind = find(delta1_constr>0);
0039 delta1_pos = delta1_constr(delta1_pos_ind);
0040 [delta1 i_delta1] = min(delta1_pos);
0041 if isempty(delta1)
0042     delta1 = inf;
0043 end
0044 delta2_constr = (epsilon+pk(gamma_lc))./(1-dk(gamma_lc));
0045 delta2_pos_ind = find(delta2_constr>0);
0046 delta2_pos = delta2_constr(delta2_pos_ind);
0047 [delta2 i_delta2] = min(delta2_pos);
0048 if isempty(delta2)
0049     delta2 = inf;
0050 end
0051 
0052 if delta1>delta2
0053     delta = delta2;
0054     i_delta = gamma_lc(delta2_pos_ind(i_delta2));
0055 else
0056     delta = delta1;
0057     i_delta = gamma_lc(delta1_pos_ind(i_delta1));
0058 end
0059 
0060 delta3_constr = (-x_k(gamma_x)./del_x_vec(gamma_x));
0061 delta3_pos_index = find(delta3_constr>0);
0062 [delta3 i_delta3] = min(delta3_constr(delta3_pos_index));
0063 out_x_index = gamma_x(delta3_pos_index(i_delta3));
0064 
0065 chk_x = 0;
0066 out_x = [];
0067 if delta3 > 0 & delta3 <= delta
0068     chk_x = 1;
0069     delta = delta3;
0070     out_x = out_x_index;
0071 end
0072 
0073 
0074 % If program goes in any of the following loops, there is some problem!!!
0075 %
0076 %%% THESE ARE PROBABLY UNNECESSARY
0077 %%% NEED TO REMOVE THEM.
0078 %
0079 % The following checks are just to deal with degenerate cases when more
0080 % than one elements want to enter or leave the support at any step
0081 % (e.g., Bernoulli matrix with small number of measurements)
0082 
0083 % This one is ONLY for those indices which are zero. And we don't know where
0084 % will its dx point in next steps, so after we calculate dx and its in opposite
0085 % direction to z_x, we will have to remove that index from the support.
0086 xk_1 = x_k+delta*del_x_vec;
0087 xk_1(out_x) = 0;
0088 wrong_sign = find(sign(xk_1(gamma_x)).*z_x(gamma_x)==-1);
0089 if ~isempty(gamma_x(wrong_sign))
0090     disp('Sign mismatch!'); 
0091     % keyboard;
0092     
0093     chk_x = 1;
0094     delta = 0;
0095     % can also choose specific element which became non-zero first but all
0096     % that matters here is AtA(gx,gl) doesn't become singular.
0097     % [val_wrong_x ind_wrong_x] =  sort(abs(del_x_vec(gamma_x(wrong_sign))),'descend');
0098     out_x = gamma_x(wrong_sign(1));
0099 end
0100 
0101 % If more than one primal constraints became active in previous iteration i.e.,
0102 % more than one elements wanted to enter the support and we added only one.
0103 % So here we need to check if those remaining elements are still active.
0104 i_delta_temp = gamma_lc(abs(pk(gamma_lc)+delta*dk(gamma_lc))-(epsilon-delta) >= 10*eps);
0105 if ~isempty(i_delta_temp)
0106     if ~isempty(out_lambda)
0107         i_delta_more = i_delta_temp;%(find(i_delta_temp~=out_lambda));
0108     else
0109         i_delta_more = i_delta_temp;
0110     end
0111     if length(i_delta_more)>=1 & ~sum((i_delta_temp==i_delta))
0112         
0113         disp('Degenerate indices in previous iteration!'); 
0114         % keyboard;
0115         
0116         % ideal way would be to check that incoming element doesn't make AtA
0117         % singular!
0118         [v_temp i_temp] = max(-pk(i_delta_more)./dk(i_delta_more));
0119         i_delta = i_delta_more(i_temp);
0120         delta = 0;
0121         chk_x = 0;
0122         out_x = [];
0123     end
0124 end