update_dual

 update_dual.m

 This function computes the minimum step size in the dual 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_lambda_p - dual update direction
 ak 
 bk
 new_lambda - element entered in the support of lambda during primal update
 out_x - element of x shrunk to zero during primal update phase.

 Outputs: 
 i_theta - index corresponding to newly active dual constraint (new_x)
 out_lambda - element in lambda shrunk to zero
 theta - dual step size 
 chk_lambda - 1  an element is removed from support of lambda
              0  a new element enters the support of x
 
 Written by: Salman Asif, Georgia Tech
 Email: sasif@ece.gatech.edu

0001 % update_dual.m
0002 %
0003 % This function computes the minimum step size in the dual 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_lambda_p - dual update direction
0012 % ak
0013 % bk
0014 % new_lambda - element entered in the support of lambda during primal update
0015 % out_x - element of x shrunk to zero during primal update phase.
0016 %
0017 % Outputs:
0018 % i_theta - index corresponding to newly active dual constraint (new_x)
0019 % out_lambda - element in lambda shrunk to zero
0020 % theta - dual step size
0021 % chk_lambda - 1  an element is removed from support of lambda
0022 %              0  a new element enters the support of x
0023 %
0024 % Written by: Salman Asif, Georgia Tech
0025 % Email: sasif@ece.gatech.edu
0026 
0027 function [i_theta, out_lambda, theta, chk_lambda] = update_dual(gamma_x, gamma_lambda, z_lambda, lambda_k, del_lambda_p, ak, bk, new_lambda, out_x);
0028 
0029 N = length(lambda_k);
0030 
0031 % gamma_xc = setdiff([1:N]', [gamma_x; out_x]); WRONG
0032 % check out_x as well, that is if outgoing x switches sign in just one step
0033 temp_gamma = zeros(N,1);
0034 temp_gamma(gamma_x) = gamma_x;
0035 gamma_xc = find([1:N]' ~= temp_gamma);
0036 
0037 theta1_constr = (1-ak(gamma_xc))./bk(gamma_xc);
0038 theta1_pos_index = find(theta1_constr>0);
0039 theta1_pos = theta1_constr(theta1_pos_index);
0040 [theta1 i_theta1] = min(theta1_pos);
0041 if isempty(theta1)
0042     theta1 = inf;
0043 end
0044 theta2_constr = -(1+ak(gamma_xc))./bk(gamma_xc);
0045 theta2_pos_index = find(theta2_constr>0);
0046 theta2_pos = theta2_constr(theta2_pos_index);
0047 [theta2 i_theta2] = min(theta2_pos);
0048 if isempty(theta2)
0049     theta2 = inf;
0050 end
0051 
0052 if theta1 > theta2
0053     theta = theta2;
0054     i_theta = gamma_xc(theta2_pos_index(i_theta2));
0055 else
0056     theta = theta1;
0057     i_theta = gamma_xc(theta1_pos_index(i_theta1));
0058 end
0059 
0060 gamma_lambda_app = [gamma_lambda; new_lambda];
0061 theta3_constr = (-lambda_k(gamma_lambda_app)./del_lambda_p(gamma_lambda_app));
0062 theta3_pos_index = find(theta3_constr>0);
0063 [theta3 i_theta3] = min(theta3_constr(theta3_pos_index));
0064 out_lambda_index = gamma_lambda_app(theta3_pos_index(i_theta3));
0065 
0066 chk_lambda = 0;
0067 out_lambda = [];
0068 if theta3 > 0 & theta3<theta
0069     chk_lambda = 1;
0070     theta = theta3;
0071     out_lambda = out_lambda_index;
0072 end
0073 
0074 
0075 %%% THESE ARE PROBABLY UNNECESSARY
0076 %%% NEED TO REMOVE THEM.
0077 
0078 % This one is ONLY for those indices which are zero. And we don't know where
0079 % will its dlambda point in next steps, so after we calculate dlambda and its in opposite
0080 % direction to z_lambda, we will have to remove that index from the support.
0081 
0082 lambdak_1 = lambda_k+theta*del_lambda_p;
0083 lambdak_1(out_lambda) = 0;
0084 wrong_sign = find(sign(lambdak_1(gamma_lambda)).*z_lambda(gamma_lambda)==-1);
0085 if ~isempty(gamma_lambda(wrong_sign))
0086     chk_lambda = 1;
0087     theta = 0;
0088     out_lambda = gamma_lambda(wrong_sign(1));
0089 end
0090 
0091 % The following checks are just to deal with degenerate cases when more
0092 % than one elements want to enter or leave the support at any step
0093 % (e.g., Bernoulli matrix with small number of measurements)
0094 
0095 % This happens if more than one dual constraints become active in one step
0096 % so some of the new elements in support of x got missed, here we check if
0097 % they are still active.
0098 i_theta_temp = gamma_xc(find(abs(ak(gamma_xc)+theta*bk(gamma_xc))-1 >= 10*eps));
0099 if ~isempty(out_x)
0100     i_theta_more = i_theta_temp(find(i_theta_temp~=out_x));
0101 else
0102     i_theta_more = i_theta_temp;
0103 end
0104 if ~isempty(i_theta_more)
0105     theta = 0;
0106     i_theta = i_theta_more(1);
0107     out_lambda=[];
0108     chk_lambda=0;
0109 end