DynamicX_DS_function

 DynamicX_update_function.m

 Solves
 min_x  ||x||_1 s.t. \|A'*(Ax-(1-\epsilon)y-(\epsilon)yt)\|_\infty <= \tau

 using homotopy update method

 [xp_h, lambda_h, gamma_xh, gamma_lh, iter, th] = DynamicX_DS_function(A, Q_glgx, R_glgx, y, yt, xp, lame, gamma_x, gamma_lambda, pk, ak, tau, maxiter)

 A: mxn matrix
 [Q_glgx, R_glgx] = qr factors for A_glgx
 y: old set of measurements
 yt: new measurements
 xp: old primal solution
 lame: old dual solution
 gamma_x: old primal support
 gamma_lambda: old dual support
 pk: old primal constraints
 ak: old dual constraints
 tau: threshold parameter.

0001 % DynamicX_update_function.m
0002 %
0003 % Solves
0004 % min_x  ||x||_1 s.t. \|A'*(Ax-(1-\epsilon)y-(\epsilon)yt)\|_\infty <= \tau
0005 %
0006 % using homotopy update method
0007 %
0008 % [xp_h, lambda_h, gamma_xh, gamma_lh, iter, th] = DynamicX_DS_function(A, Q_glgx, R_glgx, y, yt, xp, lame, gamma_x, gamma_lambda, pk, ak, tau, maxiter)
0009 %
0010 % A: mxn matrix
0011 % [Q_glgx, R_glgx] = qr factors for A_glgx
0012 % y: old set of measurements
0013 % yt: new measurements
0014 % xp: old primal solution
0015 % lame: old dual solution
0016 % gamma_x: old primal support
0017 % gamma_lambda: old dual support
0018 % pk: old primal constraints
0019 % ak: old dual constraints
0020 % tau: threshold parameter.
0021 
0022 % DS update with dynamic change in x
0023 % Author: Muhammad Salman Asif, Georgia Tech
0024 % Email: sasif@ece.gatech.edu
0025 % Created: February 2009
0026 % Modified; June 2009
0027 %
0028 %-------------------------------------------+
0029 % Copyright (c) 2009.  Muhammad Salman Asif
0030 %-------------------------------------------+
0031 
0032 function [xp_h, lambda_h, gamma_xh, gamma_lh, iter, th] = DynamicX_DS_function(A, Q_glgx, R_glgx, y, yt, xp, lame, gamma_x, gamma_lambda, pk, ak, tau, maxiter)
0033 
0034 t0 = cputime;
0035 N = length(xp);
0036 e0 = 0;
0037 gamma_xh = gamma_x;
0038 gamma_lh = gamma_lambda;
0039 xp_h = xp;
0040 lambda_h = lame;
0041 idelta = [];
0042 itheta = [];
0043 epsilon = 0;
0044 pk_old = pk;
0045 pk_old(gamma_lh) = sign(pk_old(gamma_lh))*tau;
0046 ak_old = ak;
0047 ak_old(gamma_xh) = sign(ak_old(gamma_xh));
0048 out_lambda = [];
0049 out_x = [];
0050 in_x = [];
0051 in_lambda = [];
0052 
0053 iter = 0;
0054 
0055 d = A'*(y-yt);
0056 
0057 % figure(10); clf; plot(A'*(A*xp-y)); hold on
0058 % figure(11); clf; plot(A'*A*lame); hold on
0059 
0060 while iter < maxiter
0061     iter = iter +1 ;
0062     in_lambda = [];
0063     out_lambda = [];
0064     in_x = [];
0065     out_x = [];
0066     
0067     %%%%%%%%%%%%%%%%%%%%%
0068     %%% Primal update %%%
0069     %%%%%%%%%%%%%%%%%%%%%
0070 
0071     %delx = -inv(A(:,gamma_lh)'*A(:,gamma_xh))*d(gamma_lh);
0072     delx = -R_glgx\(Q_glgx'*d(gamma_lh));
0073     delx_vec = zeros(N,1);
0074     delx_vec(gamma_xh) = delx;
0075 
0076     pk = pk_old;
0077     %dk = A'*(A*delx_vec+y-yt);
0078     dk = A'*(A(:,gamma_xh)*delx)+d;
0079     %dk = AtAf(delx_vec)+d;
0080 
0081     %%% CONTROL THE MACHINE PRECISION ERROR AT EVERY OPERATION: LIKE BELOW.
0082     pk_temp = pk;
0083     gammaL_temp = find(abs(abs(pk)-tau)<min(epsilon,1e-12));
0084     pk_temp(gammaL_temp) = sign(pk(gammaL_temp))*tau;
0085     
0086     temp_gamma = zeros(N,1);
0087     temp_gamma(gamma_lh) = gamma_lh;
0088     gamma_lc = find([1:N]' ~= temp_gamma);
0089     % gamma_lc = setdiff([1:N],[gamma_lh]);
0090 
0091     b_constr1 = (tau-pk_temp(gamma_lc))./dk(gamma_lc);
0092     b_constr2 = (tau+pk_temp(gamma_lc))./-dk(gamma_lc);
0093     b_constr3 = (-xp_h(gamma_xh)./delx_vec(gamma_xh));
0094     itheta_1 = find(b_constr1>2*eps);
0095     itheta_2 = find(b_constr2>2*eps);
0096     itheta_3 = find(b_constr3>2*eps);
0097     theta1 = min(b_constr1(itheta_1));
0098     theta2 = min(b_constr2(itheta_2));
0099     theta3 = min(b_constr3(itheta_3));
0100     if theta1>theta2
0101         theta = theta2;
0102         itheta = gamma_lc(find(b_constr2==theta2));
0103         flag_theta = 1;
0104     else
0105         theta = theta1;
0106         itheta = gamma_lc(find(b_constr1==theta1));
0107         flag_theta = 1;
0108     end
0109     if theta3 < theta
0110         theta = theta3;
0111         itheta = gamma_xh(find(b_constr3==theta3));
0112         flag_theta = 0;
0113     end
0114     % [theta itheta flag_theta]
0115     epsilon = (theta+e0);
0116     if epsilon < 0
0117         epsilon= inf;
0118     end
0119     gamma_lh_old = gamma_lh;
0120     gamma_xh_old = gamma_xh;
0121     xp_old = xp_h;
0122     xp_h = xp_h+theta*delx_vec;
0123     if flag_theta == 1
0124         % a new element is added to the dual support
0125         new_lambda = itheta;
0126         gamma_lh_app = [gamma_lh; itheta];
0127         in_lambda = itheta;
0128     else
0129         % an element is removed from the primal support
0130         outx_index = find(gamma_xh==itheta);
0131         gamma_xh = gamma_xh([1:outx_index-1 outx_index+1:length(gamma_xh)]);
0132 
0133         outl_index = length(gamma_lh);
0134         new_lambda = gamma_lh(outl_index);
0135         gamma_lh = gamma_lh([1:outl_index-1 outl_index+1:length(gamma_lh)]);
0136         gamma_lh_app = [gamma_lh; new_lambda];
0137 
0138         out_x = itheta;
0139         xp_h(itheta) = 0;
0140 
0141         %%% UPDATE THE MATRIX AND ITS INVERSE HERE
0142         % out_x and new_lambda removed
0143         Q0 = Q_glgx;
0144         R0 = R_glgx;
0145         loc_vector = zeros(length(gamma_lh_old),1);
0146         loc_vector(outl_index) = 1;
0147         rep_row = (A(:,new_lambda)'*A(:,gamma_xh_old))';
0148         [Q1 R1] = qrupdate(Q0, R0, loc_vector,-rep_row);
0149         loc_vector = zeros(length(gamma_xh_old),1);
0150         loc_vector(outx_index) = 1;
0151         rep_vec = A(:,gamma_lh_old)'*A(:,out_x);
0152         rep_vec(outl_index) = 1;
0153         [Q2 R2] = qrupdate(Q1, R1, -rep_vec, loc_vector);
0154         Q_glgx = Q2([1:outl_index-1 outl_index+1:length(gamma_lh_old)],[1:outx_index-1 outx_index+1:length(gamma_xh_old)]);
0155         R_glgx = R2([1:outx_index-1 outx_index+1:length(gamma_xh_old)],[1:outx_index-1 outx_index+1:length(gamma_xh_old)]);
0156     end
0157 
0158     if epsilon >=1
0159         theta_end = (1-e0);
0160         xp_h = xp_old + theta_end*delx_vec;
0161         gamma_xh = gamma_xh_old;
0162         th = cputime-t0;
0163         break;
0164     end
0165 
0166     %%%%%%%%%%%%%%%%%%%
0167     %%% Dual update %%%
0168     %%%%%%%%%%%%%%%%%%%
0169     z_l = sign(pk(new_lambda)+theta*dk(new_lambda));
0170     dl = A(:,gamma_xh)'*A(:,new_lambda)*z_l;
0171 
0172     %del_lambda = -inv(A(:,gamma_xh)'*A(:,gamma_lh))*dl;
0173     del_lambda =  -Q_glgx*(R_glgx'\dl);
0174     del_lam_vec = zeros(N,1);
0175     del_lam_vec(gamma_lh) = del_lambda;
0176     del_lam_vec(new_lambda) = z_l;
0177 
0178     ak = ak_old;
0179     bk = A'*(A*del_lam_vec);
0180     if flag_theta == 0
0181         if sign(ak(out_x))==sign(bk(out_x))
0182             z_l = -z_l;
0183             del_lam_vec = -del_lam_vec;
0184             bk = -bk;
0185         end
0186     end
0187     ak(out_x) = sign(ak(out_x));
0188     
0189     %%% CONTROL THE MACHINE PRECISION ERROR AT EVERY OPERATION: LIKE BELOW.
0190     ak_temp = ak;
0191     gammaX_temp = find(abs(abs(ak)-1)<1e-12);
0192     ak_temp(gammaX_temp) = sign(ak(gammaX_temp));
0193     
0194     % gamma_xc = setdiff([1:N],[gamma_xh]);
0195     temp_gamma = zeros(N,1);
0196     temp_gamma(gamma_xh) = gamma_xh;
0197     gamma_xc = find([1:N]' ~= temp_gamma);
0198 
0199     b_constr1 = (1-ak_temp(gamma_xc))./bk(gamma_xc);
0200     b_constr2 = (1+ak_temp(gamma_xc))./-bk(gamma_xc);
0201     b_constr3 = (-lambda_h(gamma_lh_app)./del_lam_vec(gamma_lh_app));
0202     idelta_1 = find(b_constr1>2*eps);
0203     idelta_2 = find(b_constr2>2*eps);
0204     idelta_3 = find(b_constr3>2*eps);
0205     delta1 = min(b_constr1(idelta_1));
0206     delta2 = min(b_constr2(idelta_2));
0207     delta3 = min(b_constr3(idelta_3));
0208 
0209     if delta1>delta2
0210         delta = delta2;
0211         idelta = gamma_xc(find(b_constr2==delta2));
0212         flag_delta = 1;
0213     else
0214         delta = delta1;
0215         idelta = gamma_xc(find(b_constr1==delta1));
0216         flag_delta = 1;
0217     end
0218     if delta3 < delta
0219         delta = delta3;
0220         idelta = gamma_lh_app(find(b_constr3==delta3));
0221         flag_delta = 0;
0222     end
0223 
0224     lambda_old = lambda_h;
0225     lambda_h = lambda_h + delta*del_lam_vec;
0226     if flag_delta == 1
0227         gamma_lh_old = gamma_lh;
0228         gamma_xh_old = gamma_xh;
0229         gamma_xh = [gamma_xh; idelta];
0230         gamma_lh = gamma_lh_app;
0231         in_x = idelta;
0232 
0233         %%% UPDATE THE MATRIX AND ITS INVERSE HERE
0234         % new_lambda and in_x added
0235         Q0 = Q_glgx;
0236         R0 = R_glgx;
0237         Q0t = [Q0 zeros(length(gamma_xh_old),1); zeros(1,length(gamma_lh_old)) 1];
0238         R0t = [R0 zeros(length(gamma_xh_old),1); zeros(1,length(gamma_lh_old)) 1];
0239         loc_vector = zeros(length(gamma_lh),1);
0240         loc_vector(end) = 1;
0241         rep_row = (A(:,new_lambda)'*A(:,gamma_xh))';
0242         [Q2t R2t] = qrupdate(Q0t, R0t, loc_vector, rep_row);
0243         rep_vec = A(:,gamma_lh)'*A(:,in_x);
0244         rep_vec(end) = -1;
0245         [Q_glgx R_glgx] = qrupdate(Q2t, R2t, rep_vec, loc_vector);
0246     else
0247         % gamma_lh = setdiff(gamma_lh_app,idelta);
0248         outl_index = find(gamma_lh==idelta);
0249         out_lambda = idelta;
0250 
0251         out_lambda = idelta;
0252         lambda_h(idelta) = 0;
0253         if ~isempty(outl_index) % Because otherwise new_lambda is removed
0254             gamma_lh_old = gamma_lh;
0255             gamma_xh_old = gamma_xh;
0256             gamma_lh(outl_index) = new_lambda;
0257 
0258             %%% UPDATE THE MATRIX AND ITS INVERSE HERE
0259             % new_lambda and out_lambda swapped
0260             Q0 = Q_glgx;
0261             R0 = R_glgx;
0262             loc_vector = zeros(length(gamma_xh),1);
0263             loc_vector(outl_index) = 1;
0264             row_replaced =  A(:,new_lambda)'*A(:,gamma_xh)-A(:,out_lambda)'*A(:,gamma_xh);
0265             [Q_glgx, R_glgx] = qrupdate(Q0,R0, loc_vector,row_replaced');
0266         end
0267 
0268     end
0269 
0270     e0 = epsilon;
0271     pk_old = pk+theta*dk;
0272     pk_old([gamma_lh; out_lambda]) = sign(pk_old([gamma_lh; out_lambda]))*tau;
0273     ak_old = ak+delta*bk;
0274     ak_old(gamma_xh) = sign(ak_old(gamma_xh));
0275     
0276     [itheta theta flag_theta idelta delta flag_delta iter epsilon]
0277 
0278 %     figure(10);
0279 %     plot(pk_old,'r.'); shg
0280 %     figure(11);
0281 %     plot(ak_old,'r.'); shg; drawnow;
0282 
0283 end
0284 
0285 th = cputime-t0;