DynamicX_BPDN_function

 DynamicX_BPDN_function.m
 
 Solves  
 min_x  \tau ||x||_1 + 1/2*||Ax-(1-\epsilon)y-(\epsilon)yt||_2^2

 using homotopy update method while changing epsilon from 0 to 1.

 [xp_h, gamma_xh, iter, cpu] = DynamicX_BPDN_function(A, AtAgx, iAtAgx, y, yt, xp, gamma_x, pk, tau, maxiter)

 A: mxn matrix
 AtAgx: AtA submatrix at row and column indices at gamma_x.
 iAtAgx: inv(AtAgx)
 y: old set of measurements
 yt: new measurements
 xp: old solution
 gamma_x: old support
 pk: old constraints
 tau: threshold parameter. 
 cpu: CPU time

0001 % DynamicX_BPDN_function.m
0002 %
0003 % Solves
0004 % min_x  \tau ||x||_1 + 1/2*||Ax-(1-\epsilon)y-(\epsilon)yt||_2^2
0005 %
0006 % using homotopy update method while changing epsilon from 0 to 1.
0007 %
0008 % [xp_h, gamma_xh, iter, cpu] = DynamicX_BPDN_function(A, AtAgx, iAtAgx, y, yt, xp, gamma_x, pk, tau, maxiter)
0009 %
0010 % A: mxn matrix
0011 % AtAgx: AtA submatrix at row and column indices at gamma_x.
0012 % iAtAgx: inv(AtAgx)
0013 % y: old set of measurements
0014 % yt: new measurements
0015 % xp: old solution
0016 % gamma_x: old support
0017 % pk: old constraints
0018 % tau: threshold parameter.
0019 % cpu: CPU time
0020 
0021 % LASSO update with dynamic change in x
0022 % Author: Muhammad Salman Asif, Georgia Tech
0023 % Email: sasif@ece.gatech.edu
0024 % Created: July 2008
0025 %
0026 %-------------------------------------------+
0027 % Copyright (c) 2009.  Muhammad Salman Asif
0028 %-------------------------------------------+
0029 
0030 function [xp_h, gamma_xh, iter,cpu] = DynamicX_BPDN_function(A, AtAgx, iAtAgx, y, yt, xp, gamma_x, pk, tau, maxiter)
0031 t0 = cputime;
0032 N = length(xp);
0033 xp_h = xp;
0034 gamma_xh = gamma_x;
0035 pk_old = pk;
0036 
0037 iter = 0;
0038 done = 0;
0039 
0040 % if ~isa(A, 'function_handle')
0041 %   Atf = @(x) A'*x;
0042 %   Af = @(x) A*x;
0043 %   AtAf = @(x)Atf(Af(x));
0044 % end
0045 warning('off','MATLAB:divideByZero');
0046 
0047 e0 = 0;
0048 itheta = [];
0049 d = A'*(y-yt);
0050 
0051 while iter < maxiter
0052     iter = iter +1 ;
0053 
0054     delx = -iAtAgx*d(gamma_xh);
0055     delx_vec = zeros(N,1);
0056     delx_vec(gamma_xh) = delx;
0057 
0058     pk = pk_old;
0059     %dk = A'*(A*delx_vec+y-yt);
0060     % Agx = A(:,gamma_xh);
0061     Ag_delx = A(:,gamma_xh)*delx;
0062     AtAg_delx = A'*Ag_delx;
0063     dk = AtAg_delx+d;
0064     %dk = AtAf(delx_vec)+d;
0065     
0066     temp_gamma = zeros(N,1);
0067     temp_gamma(gamma_xh) = gamma_xh;
0068     gamma_xc = find([1:N]' ~= temp_gamma);
0069     % gamma_xc = setdiff([1:N],[gamma_xh]);
0070 
0071     b_constr1 = (tau-pk(gamma_xc))./dk(gamma_xc);
0072     b_constr2 = (tau+pk(gamma_xc))./-dk(gamma_xc);
0073     b_constr3 = (-xp_h(gamma_xh)./delx_vec(gamma_xh));
0074     itheta_1 = find(b_constr1>2*eps);
0075     itheta_2 = find(b_constr2>2*eps);
0076     itheta_3 = find(b_constr3>2*eps);
0077     theta1 = min(b_constr1(itheta_1));
0078     theta2 = min(b_constr2(itheta_2));
0079     theta3 = min(b_constr3(itheta_3));
0080     if theta1>theta2
0081         theta = theta2;
0082         itheta = gamma_xc(find(b_constr2==theta2));
0083         flag = 1;
0084     else
0085         theta = theta1;
0086         itheta = gamma_xc(find(b_constr1==theta1));
0087         flag = 1;
0088     end
0089     if theta3 < theta
0090         theta = theta3;
0091         itheta = gamma_xh(find(b_constr3==theta3));
0092         flag = 0;
0093     end
0094 
0095     epsilon = (theta+e0);
0096     if epsilon < 0
0097         epsilon= inf;
0098     end
0099     
0100     if epsilon >=1
0101         theta_end = (1-e0);
0102         xp_h = xp_h + theta_end*delx_vec;
0103         cpu = cputime-t0;
0104         break;
0105     end
0106     
0107     gamma_xh_old = gamma_xh;
0108     xp_old = xp_h;
0109     xp_h = xp_h+theta*delx_vec;
0110     
0111     if flag == 1
0112         AtgxAnx = A(:,gamma_xh)'*A(:,itheta);
0113         AtAgx_mod = [AtAgx AtgxAnx; AtgxAnx' A(:,itheta)'*A(:,itheta)];
0114 
0115         iAtAgx = update_inverse(AtAgx_mod, iAtAgx,1);
0116         AtAgx = AtAgx_mod;
0117         
0118         AtAgx_mod = [];
0119         gamma_xh = [gamma_xh; itheta];
0120         xp_h(itheta) = 0;
0121     else
0122         g_old = gamma_xh;
0123         outx_index = find(gamma_xh==itheta);
0124         len_gamma = length(gamma_xh);
0125         gamma_xh(outx_index) = gamma_xh(len_gamma);
0126         gamma_xh = gamma_xh(1:len_gamma-1);
0127         
0128         rowi = outx_index; % ith row of A is swapped with last row (out_x)
0129         colj = outx_index; % jth column of A is swapped with last column (out_lambda)
0130         
0131         AtAgx_ij = AtAgx;
0132         temp_row = AtAgx_ij(rowi,:);
0133         AtAgx_ij(rowi,:) = AtAgx_ij(len_gamma,:);
0134         AtAgx_ij(len_gamma,:) = temp_row;
0135         temp_col = AtAgx_ij(:,colj);
0136         AtAgx_ij(:,colj) = AtAgx_ij(:,len_gamma);
0137         AtAgx_ij(:,len_gamma) = temp_col;
0138         
0139         iAtAgx_ij = iAtAgx;
0140         temp_row = iAtAgx_ij(colj,:);
0141         iAtAgx_ij(colj,:) = iAtAgx_ij(len_gamma,:);
0142         iAtAgx_ij(len_gamma,:) = temp_row;
0143         temp_col = iAtAgx_ij(:,rowi);
0144         iAtAgx_ij(:,rowi) = iAtAgx_ij(:,len_gamma);
0145         iAtAgx_ij(:,len_gamma) = temp_col;
0146         
0147         AtAgx = AtAgx_ij(1:len_gamma-1,1:len_gamma-1);
0148         iAtAgx = update_inverse(AtAgx_ij, iAtAgx_ij,2);
0149         
0150         xp_h(itheta) = 0;
0151     end
0152 
0153     e0 = (theta+e0);
0154 
0155     pk_old = pk+theta*dk;
0156     pk_old([gamma_xh; itheta]) = sign(pk_old([gamma_xh; itheta]))*tau;
0157     %pk_old(gamma_xh) = sign(pk_old(gamma_xh))*tau;
0158 end
0159 cpu = cputime-t0;
0160 
0161 %
0162 % function iAtB = update_inverse(AtB, iAtB_old,flag);
0163 %
0164 % n = size(AtB,1);
0165 %
0166 % %A11 = AtB(1:n-1,1:n-1);
0167 % A12 = AtB(1:n-1,n);
0168 % A21 = AtB(n,1:n-1);
0169 % A22 = AtB(n,n);
0170 %
0171 % % add columns
0172 % if flag == 1
0173 %     iA11 = iAtB_old;
0174 %     iA11A12 = iA11*A12;
0175 %     A21iA11 = A21*iA11;
0176 %     S = A22-A21*iA11A12;
0177 %     Q11_right = iA11A12*(A21iA11/S);
0178 % %     Q11 = iA11+ Q11_right;
0179 % %     Q12 = -iA11A12/S;
0180 % %     Q21 = -A21iA11/S;
0181 % %     Q22 = 1/S;
0182 %
0183 %     iAtB = zeros(n);
0184 %     %iAtB = [Q11 Q12; Q21 Q22];
0185 %     iAtB(1:n-1,1:n-1) = iA11+ Q11_right;
0186 %     iAtB(1:n-1,n) = -iA11A12/S;
0187 %     iAtB(n,1:n-1) =  -A21iA11/S;
0188 %     iAtB(n,n) = 1/S;
0189 % %delete columns
0190 % else if flag == 2
0191 %         Q11 = iAtB_old(1:n-1,1:n-1);
0192 %         Q12 = iAtB_old(1:n-1,n);
0193 %         Q21 = iAtB_old(n,1:n-1);
0194 %         Q22 = iAtB_old(n,n);
0195 %         Q12Q21_Q22 = Q12*(Q21/Q22);
0196 %         iAtB = Q11 - Q12Q21_Q22;
0197 %         %iAtB = iA11;
0198 %     end
0199 % end