BPDN_homotopy_function_old

 BPDN_homotopy_function.m

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

 Inputs:
 A - m x n measurement matrix
 y - measurement vector
 tau - final value of regularization parameter
 maxiter - maximum number of homotopy iterations

 Outputs:
 x_out - output for BPDN
 gamma_x - support of the solution
 total_iter - number of homotopy iterations taken by the solver
 total_time - time taken by the solver

 Written by: Salman Asif, Georgia Tech
 Email: sasif@ece.gatech.edu

 Modified 
 May 2012: Added qr update feature
 
-------------------------------------------+
 Copyright (c) 2007.  Muhammad Salman Asif
-------------------------------------------+

0001 % BPDN_homotopy_function.m
0002 %
0003 % Solves the following basis pursuit denoising (BPDN) problem
0004 % min_x  \tau ||x||_1 + 1/2*||y-Ax||_2^2
0005 %
0006 % Inputs:
0007 % A - m x n measurement matrix
0008 % y - measurement vector
0009 % tau - final value of regularization parameter
0010 % maxiter - maximum number of homotopy iterations
0011 %
0012 % Outputs:
0013 % x_out - output for BPDN
0014 % gamma_x - support of the solution
0015 % total_iter - number of homotopy iterations taken by the solver
0016 % total_time - time taken by the solver
0017 %
0018 % Written by: Salman Asif, Georgia Tech
0019 % Email: sasif@ece.gatech.edu
0020 %
0021 % Modified
0022 % May 2012: Added qr update feature
0023 %
0024 %-------------------------------------------+
0025 % Copyright (c) 2007.  Muhammad Salman Asif
0026 %-------------------------------------------+
0027 
0028 function out = BPDN_homotopy_function(A, y, in)
0029 
0030 % fprintf('Rewrite this function to add qrupdate \n');
0031 
0032 t0 = cputime;
0033 N = size(A,2);
0034 M = size(A,1);
0035 
0036 tau = in.tau;
0037 Te = inf;
0038 if isfield(in,'Te')
0039     Te = in.Te;
0040 end
0041 maxiter = 4*N;
0042 if isfield(in,'maxiter');
0043     maxiter = in.maxiter;
0044 end
0045 x_orig = in.x_orig;
0046 err_record = in.record;
0047 
0048 % Initialization of primal and dual sign and support
0049 z_x = zeros(N,1);
0050 gamma_x = [];       % Primal support
0051 
0052 % Initial step
0053 pk_old = -A'*y;
0054 [c i] = max(abs(pk_old));
0055 
0056 gamma_xk = i;
0057 
0058 epsilon = c;
0059 xk_1 = zeros(N,1);
0060 
0061 z_x(gamma_xk) = -sign(pk_old(gamma_xk));
0062 pk_old(gamma_xk) = sign(pk_old(gamma_xk))*epsilon;
0063 
0064 z_xk = z_x;
0065 
0066 % loop parameters
0067 done = 0;
0068 iter = 0;
0069 data_precision = eps;   % floating point precision
0070 
0071 old_delta = 0;
0072 out_x = [];
0073 count_delta_stop = 0;
0074 
0075 constraint_plots = 1;
0076 
0077 AtAgx = A(:,gamma_xk)'*A(:,gamma_xk);
0078 iAtAgx = inv(A(:,gamma_xk)'*A(:,gamma_xk));
0079 
0080 G = @(z) A*z;
0081 Gt = @(z) A'*z;
0082 
0083 if err_record
0084     error_table = [epsilon norm(xk_1-x_orig) 1];
0085 end
0086 
0087 while iter < maxiter
0088     iter = iter+1;
0089     % warning('off','MATLAB:divideByZero')
0090     
0091     gamma_x = gamma_xk;
0092     z_x = z_xk;
0093     x_k = xk_1;
0094     
0095     %%%%%%%%%%%%%%%%%%%%%
0096     %%%% update on x %%%%
0097     %%%%%%%%%%%%%%%%%%%%%
0098     
0099     % Update direction
0100     % delx = inv(A(:,gamma_x)'*A(:,gamma_x))*z_x(gamma_x);
0101     %     if iter > 1
0102     %     [delx_cg, res, iter] = cgsolve2(delx_vec(gamma_x), AtAgx, z_x(gamma_x), 1e-8, 100, 1)
0103     %     end
0104     delx = iAtAgx*z_x(gamma_x);
0105     % delx = (A(:,gamma_x)'*A(:,gamma_x))\z_x(gamma_x);
0106     delx_vec = zeros(N,1);
0107     delx_vec(gamma_x) = delx;
0108     %delx_vec = sparse(delx_vec);
0109     
0110     pk = pk_old;
0111     %dk = A'*(A*delx_vec);
0112     % Agdelx = A(:,gamma_x)*delx;
0113     % dk = A'*Agdelx;
0114     % dk = A'*(A*delx_vec);
0115     dk_t = G(delx_vec);
0116     dk = Gt(dk_t);
0117     
0118     %%% CONTROL THE MACHINE PRECISION ERROR AT EVERY OPERATION: LIKE BELOW.
0119     pk_temp = pk_old;
0120     gammaL_temp = find(abs(abs(pk_old)-epsilon)<min(epsilon,2*eps));
0121     %     pk_temp(gammaL_temp) = sign(pk_old(gammaL_temp))*epsilon;
0122     
0123     xk_temp = x_k;
0124     gammaX_temp = find(abs(x_k)<1*eps);
0125     %     xk_temp(gammaX_temp) = 0;
0126     %%%---
0127     
0128     % Compute the step size
0129     [idelta, out_x, delta, chk_x] = update_primal(gamma_x, gamma_x, z_x,  xk_temp, delx_vec, pk_temp, dk, epsilon, out_x);
0130     %     disp(sprintf(['iter = %d, delta = %3.4g, idelta = %d, flag = %d'], iter, delta, idelta, chk_x));
0131     %     pause(1/60)
0132     if old_delta < 4*eps && delta < 4*eps
0133         count_delta_stop = count_delta_stop + 1;
0134     else
0135         count_delta_stop = 0;
0136     end
0137     if count_delta_stop >= 500
0138         disp('stuck somewhere');
0139         break;
0140     end
0141     old_delta = delta;
0142     
0143     xk_1 = x_k+delta*delx_vec;
0144     pk_old = pk+delta*dk;
0145     epsilon_old = epsilon;
0146     epsilon = epsilon-delta;
0147     
0148     %% Check convergence criterion
0149     %     if length(gamma_x) < M/2
0150     %         xhat = zeros(N,1);
0151     %         xhat(gamma_x) = AtAgx\(A(:,gamma_x)'*y);
0152     %         if norm(y-A*xhat) < tau
0153     %             xk_1 = xhat;
0154     %             break;
0155     %         end
0156     %     end
0157     
0158     if epsilon <= tau;
0159         xk_1 = x_k + (epsilon_old-tau)*delx_vec;
0160         total_time= cputime-t0;
0161         break;
0162     end
0163     if length(gamma_x)-chk_x >= Te
0164         total_time = cputime-t0;
0165         break;
0166     end
0167     if chk_x == 1
0168         % If an element is removed from gamma_x
0169         gx_old = gamma_x;
0170         len_gamma = length(gamma_x);
0171         
0172         outx_index = find(gamma_x==out_x);
0173         gamma_x(outx_index) = gamma_x(len_gamma);
0174         gamma_x(len_gamma) = out_x;
0175         gamma_xk = gamma_x(1:len_gamma-1);
0176         
0177         rowi = outx_index; % ith row of A is swapped with last row (out_x)
0178         colj = outx_index; % jth column of A is swapped with last column (out_lambda)
0179         AtAgx_ij = AtAgx;
0180         temp_row = AtAgx_ij(rowi,:);
0181         AtAgx_ij(rowi,:) = AtAgx_ij(len_gamma,:);
0182         AtAgx_ij(len_gamma,:) = temp_row;
0183         temp_col = AtAgx_ij(:,colj);
0184         AtAgx_ij(:,colj) = AtAgx_ij(:,len_gamma);
0185         AtAgx_ij(:,len_gamma) = temp_col;
0186         iAtAgx_ij = iAtAgx;
0187         temp_row = iAtAgx_ij(colj,:);
0188         iAtAgx_ij(colj,:) = iAtAgx_ij(len_gamma,:);
0189         iAtAgx_ij(len_gamma,:) = temp_row;
0190         temp_col = iAtAgx_ij(:,rowi);
0191         iAtAgx_ij(:,rowi) = iAtAgx_ij(:,len_gamma);
0192         iAtAgx_ij(:,len_gamma) = temp_col;
0193         
0194         AtAgx = AtAgx_ij(1:len_gamma-1,1:len_gamma-1);
0195         iAtAgx = update_inverse(AtAgx_ij, iAtAgx_ij,2);
0196         xk_1(out_x) = 0;
0197     else
0198         % If an element is added to gamma_x
0199         gamma_xk = [gamma_x; idelta];
0200         new_x = idelta;
0201         
0202         AtgxAnx = A(:,gamma_x)'*A(:,new_x);
0203         AtAgx_mod = [AtAgx AtgxAnx; AtgxAnx' A(:,new_x)'*A(:,idelta)];
0204         
0205         AtAgx = AtAgx_mod;
0206         iAtAgx = update_inverse(AtAgx, iAtAgx,1);
0207         xk_1(idelta) = 0;
0208         gamma_x = gamma_xk;
0209     end
0210     
0211     z_xk = zeros(N,1);
0212     z_xk(gamma_xk) = -sign(pk_old(gamma_xk));
0213     pk_old([gamma_x]) = sign(pk_old([gamma_x]))*epsilon;
0214     %     figure(1); plot(pk_old)
0215     
0216     if err_record
0217         error_table = [error_table; epsilon norm(xk_1-x_orig) length(gamma_x)];
0218     end
0219 end
0220 if err_record
0221     error_table = [error_table; epsilon norm(xk_1-x_orig) length(gamma_x)];
0222 end
0223 total_iter = iter;
0224 total_time = cputime-t0;
0225 
0226 out = [];
0227 out.x_out = xk_1;
0228 out.gamma = gamma_x;
0229 out.iter = total_iter;
0230 out.time = total_time;
0231 out.error_table = error_table;