wtBPDN_Update_function_v2

 BPDN_homotopy_ReWeighted.m

 Solves the following basis pursuit denoising (BPDN) problem
 min_x  \|x\|_1 + 1/2*||y-AWx||_2^2

 and dynamically update the weights w_i

 Inputs:
 A - m x n measurement matrix
 y - measurement vector
 in - input structure
   tau - value of regularization parameter
   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
   yhat - dummy vector for old measuremnents
   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@ece.gatech.edu
 Created: March 16, 2011

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