DS_homotopy_function

 DS_homotopy_function.m
 
 Solves the following Dantzing selector (DS) problem
 min_x  ||x||_1  subject to  ||A'(Ax-y)||_\infty <= epsilon

 using Primal Dual Pursuit Homotopy method. 

 Usage: 
 [xk_1, lambdak_1, gamma_xk, gamma_lambdak, iter, th] = DS_homotopy_function(A, y, thresh, maxiter);
 
 Inputs: 
 A - MxN measurement matrix.  
 y - Mx1 vector of observations.
 thresh - scalar (desired value of epsilon)
 maxiter - Maximum number of primal-dual iterations.  
           Default = 10*M.
 
 Outputs:
 xk_1 - Final Primal solution
 lambdak_1 - Final Dual solution
 gamma_xk - Primal support at the solution
 gamma_lambdak - Dual support at the solution

0001 % DS_homotopy_function.m
0002 %
0003 % Solves the following Dantzing selector (DS) problem
0004 % min_x  ||x||_1  subject to  ||A'(Ax-y)||_\infty <= epsilon
0005 %
0006 % using Primal Dual Pursuit Homotopy method.
0007 %
0008 % Usage:
0009 % [xk_1, lambdak_1, gamma_xk, gamma_lambdak, iter, th] = DS_homotopy_function(A, y, thresh, maxiter);
0010 %
0011 % Inputs:
0012 % A - MxN measurement matrix.
0013 % y - Mx1 vector of observations.
0014 % thresh - scalar (desired value of epsilon)
0015 % maxiter - Maximum number of primal-dual iterations.
0016 %           Default = 10*M.
0017 %
0018 % Outputs:
0019 % xk_1 - Final Primal solution
0020 % lambdak_1 - Final Dual solution
0021 % gamma_xk - Primal support at the solution
0022 % gamma_lambdak - Dual support at the solution
0023 
0024 % Written by: Salman Asif, Georgia Tech
0025 % Email: sasif@ece.gatech.edu
0026 %
0027 %-------------------------------------------+
0028 % Copyright (c) 2007.  Muhammad Salman Asif
0029 %-------------------------------------------+
0030 
0031 function [xk_1, lambdak_1, gamma_xk, gamma_lambdak, iter, th] = DS_homotopy_function(A, y, thresh, maxiter);
0032 
0033 N = size(A,2);
0034 M = size(A,1);
0035 t0 = cputime;
0036 if (nargin < 4), maxiter = 10*M;  end
0037 
0038 % Initialization of primal and dual sign and support
0039 z_x = zeros(N,1);
0040 z_lambda = zeros(N,1);
0041 gamma_lambda = [];  % Dual support
0042 gamma_x = [];       % Primal support
0043 
0044 % Initial step
0045 Primal_constrk = -A'*y;
0046 [c i] = max(abs(Primal_constrk));
0047 
0048 gamma_lambdak = i;
0049 gamma_xk = gamma_lambdak;
0050 
0051 z_lambda(gamma_lambdak) = sign(Primal_constrk(gamma_lambdak));
0052 epsilon = c;
0053 Primal_constrk(gamma_lambdak) = sign(Primal_constrk(gamma_lambdak))*epsilon;
0054 xk_1 = zeros(N,1);
0055 
0056 lambdak_1 = zeros(N,1);
0057 lambdak_1(gamma_lambdak) = inv(A(:,gamma_lambdak)'*A(:,gamma_lambdak))*z_lambda(gamma_lambdak);
0058 
0059 Dual_constrk = A'*(A*lambdak_1);
0060 
0061 z_x(gamma_xk) = -sign(Dual_constrk(gamma_xk));
0062 Dual_constrk(gamma_xk) = sign(Dual_constrk(gamma_xk));
0063 
0064 z_xk = z_x;
0065 z_lambdak = z_lambda;
0066 
0067 % loop parameters
0068 done = 0;
0069 iter = 0;
0070 data_precision = eps;   % floating point precision
0071 
0072 old_delta = 0;
0073 out_lambda = [];
0074 out_x = [];
0075 count_delta_stop = 0;
0076 
0077 AtglAgx = A(:,gamma_lambdak)'*A(:,gamma_xk);
0078 iAtglAgx = inv(A(:,gamma_lambdak)'*A(:,gamma_xk));
0079 AtgxAgl = AtglAgx';
0080 iAtgxAgl = iAtglAgx';
0081 
0082 while iter < maxiter
0083     iter = iter+1;
0084     % warning('off','MATLAB:divideByZero')
0085 
0086     gamma_x = gamma_xk;
0087     gamma_lambda = gamma_lambdak;
0088     z_lambda = z_lambdak;
0089     z_x = z_xk;
0090     x_k = xk_1;
0091     lambda_k = lambdak_1;
0092 
0093     %%%%%%%%%%%%%%%%%%%%%
0094     %%%% update on x %%%%
0095     %%%%%%%%%%%%%%%%%%%%%
0096     
0097     % Update direction
0098     %del_x = -inv(A(:,gamma_lambda)'*A(:,gamma_x))*z_lambda(gamma_lambda);
0099     del_x = -iAtglAgx*z_lambda(gamma_lambda);
0100     del_x_vec = zeros(N,1);
0101     del_x_vec(gamma_x) = del_x;
0102 
0103     pk = Primal_constrk;
0104     % dk = AtA*del_x_vec;
0105     Agdelx = A(:,gamma_x)*del_x;
0106     dk = A'*Agdelx;
0107     
0108     %%% CONTROL THE MACHINE PRECISION ERROR AT EVERY OPERATION: LIKE BELOW.
0109     pk_temp = Primal_constrk;
0110     gammaL_temp = find(abs(abs(Primal_constrk)-epsilon)<min(epsilon,1e-12));
0111     pk_temp(gammaL_temp) = sign(Primal_constrk(gammaL_temp))*epsilon;
0112     
0113     xk_temp = x_k;
0114     gammaX_temp = find(abs(x_k)<1*eps);
0115     xk_temp(gammaX_temp) = 0;
0116     %%%---
0117     
0118     % Step size computation
0119     [i_delta, out_x, delta, chk_x] = update_primal(gamma_x, gamma_lambda, z_x,  xk_temp, del_x_vec, pk_temp, dk, epsilon, out_lambda);
0120 
0121     if old_delta < 4*eps & delta < 4*eps
0122         count_delta_stop = count_delta_stop + 1;
0123     else
0124         count_delta_stop = 0;
0125     end
0126     if count_delta_stop >= 50
0127         disp('stuck in some corner');
0128         break;
0129     end
0130     old_delta = delta;
0131     
0132     xk_1 = x_k+delta*del_x_vec;
0133     Primal_constrk = pk+delta*dk;
0134     epsilon_old = epsilon;
0135     epsilon = epsilon-delta;
0136     if epsilon < thresh;
0137         xk_1 = x_k + (epsilon_old-thresh)*del_x_vec;
0138         break;
0139     end
0140 
0141     if chk_x == 1
0142         % If an element is removed from gamma_x
0143         out_x = out_x(1);
0144         gl_old = gamma_lambda;
0145         gx_old = gamma_x;
0146 
0147         outx_index = find(gamma_x==out_x);
0148         gamma_x(outx_index) = gamma_x(end);
0149         gamma_x(end) = out_x;
0150         gamma_x = gamma_x(1:end-1);
0151 
0152         %CHECK THIS SINGULARITY CONDITION USING SCHUR COMPLEMENT IDEA !!!
0153         % X = [A B; C D];
0154         % detX = detA detS
0155         % S = D-C A^{-1} B
0156 
0157         % So the idea here is slightly different, we don't want the last
0158         % entry (Q22) in the inverse of modified Gram matrix to be zero, because
0159         % that will make inversion impossible.
0160 
0161         % This is a better thing, pick the index with maximum absolute
0162         % value in the column corresponding to out_x in iAtgxgl. This will
0163         % ensure that new matrix will be invertible.
0164         [max_val outl_index] = max(abs(iAtgxAgl(:,outx_index)));
0165         new_lambda = gamma_lambda(outl_index);
0166         some_lambda = new_lambda;
0167 
0168         outl_index = find(gamma_lambda==some_lambda);
0169         gamma_lambda(outl_index) = gamma_lambda(end);
0170         gamma_lambda(end) = some_lambda;
0171         gamma_lambdak = gamma_lambda;
0172         gamma_lambda = gamma_lambda(1:end-1);
0173 
0174         rowi = outx_index; % ith row of A is swapped with last row (out_x)
0175         colj = outl_index; % jth column of A is swapped with last column (out_lambda)
0176         AtgxAgl_ij = AtgxAgl;
0177         temp_row = AtgxAgl_ij(rowi,:);
0178         AtgxAgl_ij(rowi,:) = AtgxAgl_ij(end,:);
0179         AtgxAgl_ij(end,:) = temp_row;
0180         temp_col = AtgxAgl_ij(:,colj);
0181         AtgxAgl_ij(:,colj) = AtgxAgl_ij(:,end);
0182         AtgxAgl_ij(:,end) = temp_col;
0183         iAtgxAgl_ij = iAtgxAgl;
0184         temp_row = iAtgxAgl_ij(colj,:);
0185         iAtgxAgl_ij(colj,:) = iAtgxAgl_ij(end,:);
0186         iAtgxAgl_ij(end,:) = temp_row;
0187         temp_col = iAtgxAgl_ij(:,rowi);
0188         iAtgxAgl_ij(:,rowi) = iAtgxAgl_ij(:,end);
0189         iAtgxAgl_ij(:,end) = temp_col;
0190 
0191         AtgxAgl = AtgxAgl_ij(1:end-1,1:end-1);
0192         AtglAgx = AtgxAgl';
0193         iAtgxAgl = update_inverse(AtgxAgl_ij, iAtgxAgl_ij,2);
0194         iAtglAgx = iAtgxAgl';
0195         xk_1(out_x) = 0;
0196     else
0197         % If a new element is added to gamma_lambda
0198         gamma_lambdak = [gamma_lambda; i_delta];
0199         new_lambda = i_delta;
0200         lambda_k(new_lambda) = 0;
0201     end
0202 
0203     out_lambda = [];
0204     z_lambdak = zeros(N,1);
0205     z_lambdak(gamma_lambdak) = sign(Primal_constrk(gamma_lambdak));
0206     sgn_new_lambda = sign(Primal_constrk(new_lambda));
0207     Primal_constrk(gamma_lambdak) = sign(Primal_constrk(gamma_lambdak))*epsilon;
0208 
0209 
0210     %%%%%%%%%%%%%%%%%%%%%%%%
0211     %%% update on lambda %%%
0212     %%%%%%%%%%%%%%%%%%%%%%%%
0213     
0214     sgn_new_lambda = sign(Primal_constrk(new_lambda));
0215     
0216     AtgxAnl = A(:,gamma_x)'*A(:,new_lambda);
0217     %AtgxAnl = AtA(gamma_x,new_lambda);
0218 
0219     % Update direction
0220     %del_lambda = -inv(A(:,gamma_x)'*A(:,gamma_lambda))*(A(:,gamma_x)'*A(:,new_lambda));
0221     del_lambda = -iAtgxAgl*AtgxAnl;
0222     del_lambda_p = zeros(N,1);
0223     del_lambda_p(gamma_lambda) = del_lambda*sgn_new_lambda;
0224     del_lambda_p(new_lambda) = 1*sgn_new_lambda;
0225 
0226     ak = Dual_constrk;
0227     % bk = AtA*del_lambda_p;
0228     Agdel_lam = A(:,[gamma_lambda; new_lambda])*del_lambda_p([gamma_lambda; new_lambda]);
0229     bk = A'*Agdel_lam;
0230 
0231     % check if the sign of update direction is correct
0232     if sign(bk(out_x)) == sign(ak(out_x)) & abs(bk(out_x)) >=eps
0233         bk = -bk;
0234         del_lambda_p = -del_lambda_p;
0235     end
0236     
0237     %%% CONTROL THE MACHINE PRECISION ERROR AT EVERY OPERATION: LIKE BELOW.
0238     ak_temp = Dual_constrk;
0239     gammaX_temp = find(abs(abs(Dual_constrk)-1)<2*eps);
0240     ak_temp(gammaX_temp) = sign(Dual_constrk(gammaX_temp));
0241     
0242     lambdak_temp = lambda_k;
0243     gammaL_temp = find(abs(lambda_k)<1*eps);
0244     lambdak_temp(gammaL_temp) = 0;
0245     %%%---
0246     
0247     % Step size computation
0248     [i_theta, out_lambda, theta, chk_lambda] = update_dual(gamma_x, gamma_lambda, z_lambda, lambdak_temp, del_lambda_p, ak_temp, bk, new_lambda, out_x);
0249 
0250     out_x = [];
0251 
0252     lambdak_1 = lambda_k + theta*del_lambda_p;
0253     Dual_constrk = ak+theta*bk;
0254 
0255     if chk_lambda == 1
0256         % If an element is removed from gamma_lambda
0257         outl_index = find(gamma_lambda==out_lambda);
0258         i_theta = gamma_x(end);
0259 
0260         if ~isempty(outl_index) % because otherwise new_lambda is removed.
0261             % Use Sherman-Woodburry-Morrison formula
0262             % (A+BC)^{-1} = A^{-1} - A^{-1}B(I+CA^{-1}B)^{-1}CA^{-1};
0263             iA = iAtgxAgl;
0264             B = AtgxAnl-A(:,gamma_x)'*A(:,out_lambda);
0265             C = zeros(1,length(gamma_lambda));
0266             C(outl_index) = 1;
0267             AtgxAgl(:,outl_index) = AtgxAnl;
0268             AtglAgx = AtgxAgl';
0269 
0270             iAB = iA*B;
0271             CiA = iA(outl_index,:);
0272             iAtgxAgl = iA-iAB*((CiA)/(1+iAB(outl_index)));
0273             iAtglAgx = iAtgxAgl';
0274 
0275             gamma_lambda(outl_index) = new_lambda;
0276         end
0277 
0278         z_xk = zeros(N,1);
0279         z_xk(gamma_x) = -sign(Dual_constrk(gamma_x));
0280         lambdak_1(out_lambda) = 0;
0281     else
0282         % If an element is added to gamma_x
0283 
0284         %%% adding columns
0285         % AtglAgx_mod = zeros(length(gamma_x));
0286         AtglAgx_mod = [AtglAgx A(:,gamma_lambda)'*A(:,i_theta); AtgxAnl' A(:,new_lambda)'*A(:,i_theta)];
0287         %AtglAgx_mod = [AtglAgx AtA(gamma_lambda,i_theta); AtgxAnl' AtA(new_lambda,i_theta)];
0288 
0289         %CHECK THIS SINGULARITY CONDITION USING SCHUR COMPLEMENT IDEA !!!
0290         % X = [A B; C D];
0291         % detX = detA detS
0292         % S = D-C A^{-1} B
0293         A11 = AtglAgx;
0294         A12 = AtglAgx_mod(1:end-1,end);
0295         A21 = AtglAgx_mod(end,1:end-1);
0296         A22 = AtglAgx_mod(end,end);
0297         S = A22 - A21*(iAtglAgx*A12);
0298         if abs(S) < 2*eps
0299             disp('matrix has become singular');
0300             % done = 1;
0301             % break;
0302         end
0303         AtglAgx = AtglAgx_mod;
0304         AtgxAgl = AtglAgx';
0305         iAtglAgx = update_inverse(AtglAgx, iAtglAgx,1);
0306         iAtgxAgl = iAtglAgx';
0307 
0308         out_lambda = [];
0309         gamma_lambda = [gamma_lambda; new_lambda];
0310 
0311         gamma_x = [gamma_x; i_theta];
0312         z_xk = zeros(N,1);
0313         z_xk(gamma_x) = -sign(Dual_constrk(gamma_x));
0314         xk_1(i_theta) = 0;
0315     end
0316     Dual_constrk(gamma_x) = sign(Dual_constrk(gamma_x));
0317     gamma_lambdak = gamma_lambda;
0318     gamma_xk = gamma_x;
0319 end
0320 th = cputime-t0;