DynamicSeq_BPDN_function

 DynamicLasso_update_function.m
 
 Solves  
 min_x  \epsilon ||x||_1 + 1/2*||Ax-y||_2^2 + 1/2* ||bx-w||_2^2

 using homotopy update method

 [xp_h, gamma_xh, iter] = DynamicSeq_BPDN_function(A, b, AtAgx, iAtAgx, y, w, xp, gamma_x, pk, tau, chk_e0, maxiter)

 A: mxn matrix
 b: new rows in measurement matrix
 AtAgx: AtA submatrix at row and column indices at gamma_x.
 iAtAgx: inv(AtAgx)
 y: old set of measurements
 w: new measurements
 xp: old solution
 gamma_x: old support
 pk: old constraints
 tau: threshold parameter. 
 chk_e0: 0 or 1: This selects if we want to take U = (A'A+b'b) (1) or (A'A+e0*b'b) (0)...!
 In simulation U = (A'A+b'b) works better

0001 % DynamicLasso_update_function.m
0002 %
0003 % Solves
0004 % min_x  \epsilon ||x||_1 + 1/2*||Ax-y||_2^2 + 1/2* ||bx-w||_2^2
0005 %
0006 % using homotopy update method
0007 %
0008 % [xp_h, gamma_xh, iter] = DynamicSeq_BPDN_function(A, b, AtAgx, iAtAgx, y, w, xp, gamma_x, pk, tau, chk_e0, maxiter)
0009 %
0010 % A: mxn matrix
0011 % b: new rows in measurement matrix
0012 % AtAgx: AtA submatrix at row and column indices at gamma_x.
0013 % iAtAgx: inv(AtAgx)
0014 % y: old set of measurements
0015 % w: new measurements
0016 % xp: old solution
0017 % gamma_x: old support
0018 % pk: old constraints
0019 % tau: threshold parameter.
0020 % chk_e0: 0 or 1: This selects if we want to take U = (A'A+b'b) (1) or (A'A+e0*b'b) (0)...!
0021 % In simulation U = (A'A+b'b) works better
0022 
0023 % BPDN solution update with dynamic measurements
0024 % Author: Salman Asif, Georgia Tech
0025 % Email: sasif@ece.gatech.edu
0026 % Created: February 2008
0027 %
0028 %-------------------------------------------+
0029 % Copyright (c) 2008.  Muhammad Salman Asif
0030 %-------------------------------------------+
0031 
0032 function [xp_h, gamma_xh, iter, th] = DynamicSeq_BPDN_function(A, b, AtAgx, iAtAgx, y, w, xp, gamma_x, pk, tau, chk_e0, maxiter)
0033 
0034 t0 = cputime;
0035 N = length(xp);
0036 xp_h = xp;
0037 gamma_xh = gamma_x;
0038 itheta = [];
0039 e0 = 0;
0040 
0041 % Gram matrix update
0042 % AtAgx = A(:,gamma_xh)'*A(:,gamma_xh)+chk_e0*b(:,gamma_xh)'*b(:,gamma_xh);
0043 pk_old = pk;
0044 
0045 notdone = 1;
0046 iter = 0;
0047 
0048 while iter < maxiter
0049     iter = iter +1 ;
0050 
0051     e0h = e0*(1-chk_e0)+1*chk_e0; 
0052     % this way we can set e0h to 1 or e0 which changes the definition of theta and U accordingly
0053 
0054     %%% UPDATE
0055     %     U = A(:,gamma_xh)'*A(:,gamma_xh)+e0h*b(:,gamma_xh)'*b(:,gamma_xh);
0056     %     iU = inv(U);
0057     U = AtAgx;
0058     iU = iAtAgx;
0059 
0060     d = b(:,gamma_xh)'*(b(:,gamma_xh)*xp_h(gamma_xh)-w);
0061     u = (b(:,gamma_xh)*iU*b(:,gamma_xh)');
0062     dx = -iU*d;
0063     pk = pk_old; %A'*(A*xp_h-y)+e0*b'*(b*xp_h-w);
0064     %dk = b'*(b*xp_h-w)+(Atf(A(:,gamma_xh)*dx)+b'*(e0h*(b(:,gamma_xh)*dx)));
0065     delx_vec = zeros(N,1);
0066     delx_vec(gamma_xh) = dx;
0067     Adelx = A(:,gamma_xh)*dx;
0068     AtAdelx = (A')*Adelx;
0069     dk = b'*(b*xp_h-w)+AtAdelx+b'*(e0h*(b*delx_vec));
0070 
0071     temp_gamma = zeros(N,1);
0072     temp_gamma(gamma_xh) = gamma_xh;
0073     gamma_xc = find([1:N]' ~= temp_gamma);
0074     % gamma_xc = setdiff([1:N],[gamma_xh]);
0075 
0076     b_constr1 = (tau-pk(gamma_xc))./dk(gamma_xc);
0077     b_constr2 = (tau+pk(gamma_xc))./-dk(gamma_xc);
0078     b_constr3 = (-xp_h(gamma_xh)./delx_vec(gamma_xh));
0079     itheta_1 = find(b_constr1>2*eps);
0080     itheta_2 = find(b_constr2>2*eps);
0081     itheta_3 = find(b_constr3>2*eps);
0082     theta1 = min(b_constr1(itheta_1));
0083     theta2 = min(b_constr2(itheta_2));
0084     theta3 = min(b_constr3(itheta_3));
0085     if isempty(theta1)
0086         theta1 = inf;
0087     end
0088     if isempty(theta2)
0089         theta2 = inf;
0090     end
0091     if isempty(theta3)
0092         theta3 = inf;
0093     end
0094 
0095     if theta1>theta2
0096         theta = theta2;
0097         itheta = gamma_xc(find(b_constr2==theta2));
0098         flag = 1;
0099     else
0100         theta = theta1;
0101         itheta = gamma_xc(find(b_constr1==theta1));
0102         flag = 1;
0103     end
0104     if theta3 < theta
0105         theta = theta3;
0106         itheta = gamma_xh(find(b_constr3==theta3));
0107         flag = 0;
0108     end
0109     epsilon = e0+(theta)/(1-theta*u);
0110     epsilon = (theta+e0-e0h*u*theta)/(1-theta*u);
0111     
0112     if epsilon < 0
0113         epsilon= inf;
0114     end
0115 
0116     if epsilon >=1
0117         theta_end = (1-e0)/(1+u-e0h*u);
0118         xp_h = xp_h + theta_end*delx_vec;
0119         th = cputime-t0;
0120         break;
0121     end
0122 
0123     xp_h_old = xp_h;
0124     gamma_xh_old = gamma_xh;
0125     xp_h = xp_h+theta*delx_vec;
0126     
0127     if flag == 1
0128         if e0h == e0
0129             AtgxAnx = A(:,gamma_xh)'*A(:,itheta)+b(:,gamma_xh)'*(epsilon*b(:,itheta));
0130             % btbgx = b(:,gamma_xh)'*b(:,gamma_xh);
0131             AtAgx_new = AtAgx + ((epsilon-e0)*b(:,gamma_xh)')*b(:,gamma_xh);
0132             iAtAgx_bt = iAtAgx*b(:,gamma_xh)';
0133             % iAtAgx_new = iAtAgx - (epsilon-e0)*(iAtAgx_bt)*(inv(1+(epsilon-e0)*b(:,gamma_xh)*iAtAgx_bt)*iAtAgx_bt');
0134             iAtAgx_new = iAtAgx - (((epsilon-e0)/(1+(epsilon-e0)*(b(:,gamma_xh)*iAtAgx_bt)))*iAtAgx_bt)*iAtAgx_bt';
0135             AtAgx_mod = [AtAgx_new AtgxAnx; AtgxAnx' A(:,itheta)'*A(:,itheta)+epsilon*b(:,itheta)'*b(:,itheta)];
0136 
0137             iAtAgx = update_inverse(AtAgx_mod, iAtAgx_new,1);
0138             AtAgx = AtAgx_mod;
0139             
0140             AtAgx_mod = [];
0141             AtAgx_new = [];
0142             iAtAgx_new = [];
0143             
0144             gamma_xh = [gamma_xh; itheta];
0145             xp_h(itheta) = 0;
0146         else
0147             AtgxAnx = A(:,gamma_xh)'*A(:,itheta)+b(:,gamma_xh)'*b(:,itheta);
0148             % btbgx = b(:,gamma_xh)'*b(:,gamma_xh);
0149             iAtAgx_bt = iAtAgx*b(:,gamma_xh)';
0150             AtAgx_mod = [AtAgx AtgxAnx; AtgxAnx' A(:,itheta)'*A(:,itheta)+b(:,itheta)'*b(:,itheta)];
0151 
0152             iAtAgx = update_inverse(AtAgx_mod, iAtAgx,1);
0153             AtAgx = AtAgx_mod;
0154             
0155             AtAgx_mod = [];
0156             gamma_xh = [gamma_xh; itheta];
0157             xp_h(itheta) = 0;
0158         end
0159     else
0160         outx_index = find(gamma_xh==itheta);
0161         gamma_xh(outx_index) = gamma_xh(end);
0162         gamma_xh(end) = itheta;
0163         gamma_xh = gamma_xh(1:end-1);
0164 
0165         if e0h == e0
0166             % btbgx = b(:,gamma_xh_old)'*b(:,gamma_xh_old);
0167             AtAgx_new = AtAgx + ((epsilon-e0)*b(:,gamma_xh_old)')*b(:,gamma_xh_old);
0168             iAtAgx_bt = iAtAgx*b(:,gamma_xh_old)';
0169             % iAtAgx_new = iAtAgx - (epsilon-e0)*(iAtAgx_bt)*(inv(1+(epsilon-e0)*b(:,gamma_xh_old)*iAtAgx_bt)*iAtAgx_bt');
0170             iAtAgx_new = iAtAgx - (((epsilon-e0)/(1+(epsilon-e0)*(b(:,gamma_xh_old)*iAtAgx_bt)))*iAtAgx_bt)*iAtAgx_bt';
0171 
0172             rowi = outx_index; % ith row of A is swapped with last row (out_x)
0173             colj = outx_index; % jth column of A is swapped with last column (out_x)
0174             AtAgx_ij = AtAgx_new;
0175             temp_row = AtAgx_ij(rowi,:);
0176             AtAgx_ij(rowi,:) = AtAgx_ij(end,:);
0177             AtAgx_ij(end,:) = temp_row;
0178             temp_col = AtAgx_ij(:,colj);
0179             AtAgx_ij(:,colj) = AtAgx_ij(:,end);
0180             AtAgx_ij(:,end) = temp_col;
0181             iAtAgx_ij = iAtAgx_new;
0182             temp_row = iAtAgx_ij(colj,:);
0183             iAtAgx_ij(colj,:) = iAtAgx_ij(end,:);
0184             iAtAgx_ij(end,:) = temp_row;
0185             temp_col = iAtAgx_ij(:,rowi);
0186             iAtAgx_ij(:,rowi) = iAtAgx_ij(:,end);
0187             iAtAgx_ij(:,end) = temp_col;
0188 
0189             AtAgx = AtAgx_ij(1:end-1,1:end-1);
0190             iAtAgx = update_inverse(AtAgx_ij, iAtAgx_ij,2);
0191             
0192             AtAgx_mod = [];
0193             AtAgx_new = [];
0194             iAtAgx_new = [];
0195             AtAgx_ij = [];
0196             iAtAgx_ij = [];
0197         else
0198 
0199             rowi = outx_index; % ith row of A is swapped with last row (out_x)
0200             colj = outx_index; % jth column of A is swapped with last column (out_x)
0201             AtAgx_ij = AtAgx;
0202             temp_row = AtAgx_ij(rowi,:);
0203             AtAgx_ij(rowi,:) = AtAgx_ij(end,:);
0204             AtAgx_ij(end,:) = temp_row;
0205             temp_col = AtAgx_ij(:,colj);
0206             AtAgx_ij(:,colj) = AtAgx_ij(:,end);
0207             AtAgx_ij(:,end) = temp_col;
0208             iAtAgx_ij = iAtAgx;
0209             temp_row = iAtAgx_ij(colj,:);
0210             iAtAgx_ij(colj,:) = iAtAgx_ij(end,:);
0211             iAtAgx_ij(end,:) = temp_row;
0212             temp_col = iAtAgx_ij(:,rowi);
0213             iAtAgx_ij(:,rowi) = iAtAgx_ij(:,end);
0214             iAtAgx_ij(:,end) = temp_col;
0215 
0216             AtAgx = AtAgx_ij(1:end-1,1:end-1);
0217             iAtAgx = update_inverse(AtAgx_ij, iAtAgx_ij,2);
0218             
0219             AtAgx_ij = [];
0220             iAtAgx_ij = [];
0221         end
0222 
0223         xp_h(itheta) = 0;
0224     end
0225 
0226     e0 = epsilon;
0227     pk_old = pk+theta*dk;
0228     pk_old([gamma_xh; itheta]) =  sign(pk_old([gamma_xh; itheta]))*tau;
0229 end
0230 
0231