wtBPDN_Update_function_v1

 Iterative reweighting updates for BPDN

 v1 (each weight treated as a separate regularization parameter...)
 This version is reported in the paper.
 
 Solves the following basis pursuit denoising (BPDN) problem
 min_x  \Sum \w_i |x_i| + 1/2*||y-Ax||_2^2

 and dynamically update the weights w_i

 Inputs:
 A - m x n measurement matrix
 y - measurement vector
 in - input structure
   x_old - old estimate
   gamma - old signal support
   AtA - A(:,gamma_xh)'*A(:,gamma_xh);
   iAtA - inv(AtA)
    or QR/Cholesky factors
   pk_old - old constrtaints set
   W_old - old values of regularization parameter
   W_new - new values of regularization parameter
   maxiter - maximum number of homotopy iterations

 Outputs:
 out - output structure
   x_out - output for BPDN
   gamma - support of the solution
   iter - number of homotopy iterations taken by the solver
   time - time taken by the solver

 Written by: Salman Asif, Georgia Tech
 Email: sasif@gatech.edu
 Created: March 16, 2011

0001 % Iterative reweighting updates for BPDN
0002 %
0003 % v1 (each weight treated as a separate regularization parameter...)
0004 % This version is reported in the paper.
0005 %
0006 % Solves the following basis pursuit denoising (BPDN) problem
0007 % min_x  \Sum \w_i |x_i| + 1/2*||y-Ax||_2^2
0008 %
0009 % and dynamically update the weights w_i
0010 %
0011 % Inputs:
0012 % A - m x n measurement matrix
0013 % y - measurement vector
0014 % in - input structure
0015 %   x_old - old estimate
0016 %   gamma - old signal support
0017 %   AtA - A(:,gamma_xh)'*A(:,gamma_xh);
0018 %   iAtA - inv(AtA)
0019 %    or QR/Cholesky factors
0020 %   pk_old - old constrtaints set
0021 %   W_old - old values of regularization parameter
0022 %   W_new - new values of regularization parameter
0023 %   maxiter - maximum number of homotopy iterations
0024 %
0025 % Outputs:
0026 % out - output structure
0027 %   x_out - output for BPDN
0028 %   gamma - support of the solution
0029 %   iter - number of homotopy iterations taken by the solver
0030 %   time - time taken by the solver
0031 %
0032 % Written by: Salman Asif, Georgia Tech
0033 % Email: sasif@gatech.edu
0034 % Created: March 16, 2011
0035 
0036 function out = wtBPDN_Update_function_v1(A, y, in);
0037 
0038 % x_old, gamma_xh, AtAgx, iAtAgx, pk_old, W_old, W_new, maxiter
0039 x_old = in.x_old;
0040 gamma_xh = in.gamma;
0041 pk_old = in.pk_old;
0042 W_old = in.W_old;
0043 W_new = in.W_new;
0044 maxiter = in.maxiter;
0045 
0046 t0 = cputime;
0047 N = size(A,2);
0048 M = size(A,1);
0049 
0050 % pk_old = A'*(A*x_k-y);
0051 % pk_old(gamma_xh) = sign(pk_old(gamma_xh)).*W_old(gamma_xh);
0052 xk_1 = x_old;
0053 
0054 epsilon = 0;
0055 Wk = (1-epsilon)*W_old + epsilon*W_new;
0056 dW = W_new-W_old;
0057 z_x = -sign(pk_old);
0058 
0059 done = 0;
0060 iter = 0;
0061 idelta = gamma_xh(1); flag = 1;
0062 itr_history = [];
0063 
0064 % initialize delx
0065 in_delx = []; 
0066 delx_mode = in.delx_mode;
0067 rhs = -dW.*z_x; 
0068 opts = in; % update_delx is now using opts as input parameter
0069 update_mode = 'init1';
0070 update_delx;
0071 
0072 
0073 while ~done
0074     iter = iter+1;
0075     x_k = xk_1;
0076     
0077     delx_vec = zeros(N,1);
0078     delx_vec(gamma_xh) = delx;
0079     
0080     if sign(delx_vec(idelta)) == sign(pk_old(idelta)) && iter > 1
0081         delta = 0; flag = 0;
0082     else
0083         pk = pk_old;
0084         dk = A'*(A(:,gamma_xh)*delx);
0085         
0086         if flag == 0 && dW(idelta) < sign(pk(idelta))*dk(idelta)
0087             delta = 0; flag = 1;
0088         else
0089             temp_gamma = zeros(N,1);
0090             temp_gamma(gamma_xh) = gamma_xh;
0091             gamma_xc = find([1:N]' ~= temp_gamma);
0092             
0093             %% Calculate homotopy step size
0094             in.delta_flag = 2;
0095             in.pk = pk; in.dk = dk;
0096             in.ak = Wk; in.bk = dW;
0097             in.gamma = gamma_xh; in.gamma_c = gamma_xc;
0098             in.delx_vec = delx_vec; in.x = xk_1;
0099             out = compute_delta(in);
0100             delta = out.delta; idelta = out.idelta;
0101             flag = out.flag;
0102         end
0103     end
0104     %% Update support
0105     e0 = epsilon;
0106     epsilon = e0 + delta;
0107     
0108     if epsilon > 1
0109         delta_end = 1-e0;
0110         xk_1 = x_k + delta_end*delx_vec;
0111         break;
0112     end
0113     
0114     xk_1 = x_k + delta*delx_vec;
0115     gamma_x_old = gamma_xh;
0116     
0117     % disp(sprintf(['iter = %d, delta = %3.4g, idelta = %d, flag = %d'], iter, delta, idelta, flag));
0118     itr_history = [itr_history; idelta delta flag];
0119     
0120     % update support
0121     update_supp; 
0122 
0123     pk_old = pk+delta*dk;
0124     Wk_old = Wk;
0125     Wk = (1-epsilon)*W_old + epsilon*W_new;
0126     pk_old([gamma_xh; idelta]) = sign(pk_old([gamma_xh; idelta])).*Wk([gamma_xh; idelta]);
0127     
0128     % update delx
0129     z_x = -sign(pk_old);
0130     rhs = -dW.*z_x;
0131     update_mode = 'update';
0132     update_delx;
0133     %     AtAgx = (A(:,gamma_xh)'*A(:,gamma_xh));
0134     %     delx2 = AtAgx\rhs(gamma_xh); % -AtAgx\(dW(gamma_xh).*z_x);
0135     %     figure(112); plot([delx delx2]);
0136     %     if norm(delx-delx2) > 1e-5
0137     %         stp = 1;
0138     %     end
0139     
0140     constr_violation = nnz((abs(pk_old(gamma_xc))-Wk(gamma_xc))>1e-10);
0141     sign_violation = nnz((sign(pk_old(gamma_xh))+sign(xk_1(gamma_xh)))>1);
0142     if constr_violation
0143         chk = gamma_xc((abs(pk_old(gamma_xc))-Wk(gamma_xc))>1e-10);
0144         stp = 1;
0145         fprintf('problem... with constraint violation -- %s \n', mfilename);
0146         fprintf('Refactorize the matrix... recompute delx \n');
0147         % some times it comes here due to bad conditioning of AtAgx.
0148         update_mode = 'init0';
0149         update_delx;
0150     end
0151     if sign_violation>=1
0152         chk = gamma_xh(sign(pk_old(gamma_xh))+sign(xk_1(gamma_xh))>1);
0153         stp = 1;
0154         fprintf('problem... sign mismatch -- %s \n',mfilename);
0155         fprintf('Refactorize the matrix... recompute delx \n');
0156         update_mode = 'init0';
0157         update_delx;
0158     end
0159 end
0160 
0161 total_iter = iter;
0162 total_time = cputime-t0;
0163 
0164 out = [];
0165 out.x_out = xk_1;
0166 out.gamma = gamma_xh;
0167 out.iter = total_iter;
0168 out.time = total_time;