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
in - input parameters
tau - final value of regularization parameter
maxiter - maximum number of homotopy iterations
x_orig - original signal for error computation
record - for recording history
delx_mode - rank-one update using matrix inversion lemma (mil) or qr factorization (qr)
Te - maximum number of non-zero elements allowed
Outputs:
out
x_out - output for BPDN
gamma - support of the solution
iter - number of homotopy iterations taken by the solver
time - time taken by the solver
Example usage:
in = [];
in.tau = tau; % regularization parameter
in.maxiter = maxiter; % maximum allowed iterations
in.x_orig = x; % original signal for error computation
in.record = 1; % to record history
in.delx_mode = 'mil'; % rank-one update using matrix inversion lemma (mil) or qr factorization (qr)
in.Te = Te; % maximum number of non-zero elements allowed
out = BPDN_homotopy_function(A, y, in); % solve BPDN using matrix A and measurements y
xh = out.x_out; % output solution
gamma = out.gamma; % support of the solution
iter = out.iter; % number of homotopy iterations used
time = out.time; % computation time
Written by: Salman Asif, Georgia Tech
Email: sasif@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 % in - input parameters 0010 % tau - final value of regularization parameter 0011 % maxiter - maximum number of homotopy iterations 0012 % x_orig - original signal for error computation 0013 % record - for recording history 0014 % delx_mode - rank-one update using matrix inversion lemma (mil) or qr factorization (qr) 0015 % Te - maximum number of non-zero elements allowed 0016 % 0017 % Outputs: 0018 % out 0019 % x_out - output for BPDN 0020 % gamma - support of the solution 0021 % iter - number of homotopy iterations taken by the solver 0022 % time - time taken by the solver 0023 % 0024 % Example usage: 0025 % in = []; 0026 % in.tau = tau; % regularization parameter 0027 % in.maxiter = maxiter; % maximum allowed iterations 0028 % in.x_orig = x; % original signal for error computation 0029 % in.record = 1; % to record history 0030 % in.delx_mode = 'mil'; % rank-one update using matrix inversion lemma (mil) or qr factorization (qr) 0031 % in.Te = Te; % maximum number of non-zero elements allowed 0032 % out = BPDN_homotopy_function(A, y, in); % solve BPDN using matrix A and measurements y 0033 % xh = out.x_out; % output solution 0034 % gamma = out.gamma; % support of the solution 0035 % iter = out.iter; % number of homotopy iterations used 0036 % time = out.time; % computation time 0037 % 0038 % Written by: Salman Asif, Georgia Tech 0039 % Email: sasif@gatech.edu 0040 % 0041 % Modified 0042 % May 2012: Added qr update feature 0043 % 0044 %-------------------------------------------+ 0045 % Copyright (c) 2007. Muhammad Salman Asif 0046 %-------------------------------------------+ 0047 0048 function out = BPDN_homotopy_function(A, y, in) 0049 0050 % fprintf('Rewrite this function to add qrupdate \n'); 0051 0052 t0 = cputime; 0053 N = size(A,2); 0054 M = size(A,1); 0055 0056 tau = in.tau; 0057 Te = inf; 0058 if isfield(in,'Te') 0059 Te = in.Te; 0060 end 0061 maxiter = 4*N; 0062 if isfield(in,'maxiter'); 0063 maxiter = in.maxiter; 0064 end 0065 early_terminate = 0; 0066 if isfield(in,'early_terminate') 0067 early_terminate = in.early_terminate; 0068 end 0069 x_orig = in.x_orig; 0070 err_record = in.record; 0071 0072 % Initialization of sign and support 0073 xk_1 = zeros(N,1); 0074 pk_old = -A'*y; 0075 [c idelta] = max(abs(pk_old)); 0076 0077 gamma_xh = idelta; 0078 temp_gamma = zeros(N,1); 0079 temp_gamma(gamma_xh) = gamma_xh; 0080 gamma_xc = find([1:N]' ~= temp_gamma); 0081 0082 epsilon = c; 0083 z_x = zeros(N,1); 0084 z_x(gamma_xh) = -sign(pk_old(gamma_xh)); 0085 pk_old(gamma_xh) = sign(pk_old(gamma_xh))*epsilon; 0086 0087 % loop parameters 0088 done = 0; 0089 iter = 0; 0090 itr_history = []; 0091 error_table = []; 0092 obj = []; 0093 0094 % initialize delx 0095 in_delx = []; 0096 delx_mode = in.delx_mode; 0097 update_mode = 'init0'; 0098 rhs = z_x; 0099 update_delx; 0100 0101 if err_record 0102 error_table = [epsilon norm(xk_1-x_orig)/norm(x_orig) 1]; 0103 % obj = 1/2*norm(A*xk_1-y)^2+epsilon*norm(xk_1,1); 0104 end 0105 0106 while iter < maxiter 0107 iter = iter+1; 0108 % warning('off','MATLAB:divideByZero') 0109 0110 x_k = xk_1; 0111 0112 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0113 %%%% update support of x %%%% 0114 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0115 % update direction 0116 delx_vec = zeros(N,1); 0117 delx_vec(gamma_xh) = delx; 0118 0119 if sign(delx_vec(idelta)) == sign(pk_old(idelta)) && iter > 1 0120 delta = 0; flag = 0; 0121 else 0122 pk = pk_old; 0123 % dk = A'*(A(:,gamma_xh)*delx); 0124 dk = A'*(A*delx_vec); 0125 0126 %%%--- compute step size 0127 in = []; 0128 0129 % Setting shrinkage_flag to zero shrinks new active constraint towards the 0130 % final value instantly if doing so doesn't disturb the active set 0131 in.delta_flag = 2; 0132 in.pk = pk; in.dk = dk; 0133 in.ak = epsilon; in.bk = -1; 0134 in.gamma = gamma_xh; in.gamma_c = gamma_xc; 0135 in.delx_vec = delx_vec; in.x = xk_1; 0136 out = compute_delta(in); 0137 delta = out.delta; idelta = out.idelta; 0138 flag = out.flag; 0139 end 0140 0141 xk_1 = x_k+delta*delx_vec; 0142 pk_old = pk+delta*dk; 0143 epsilon_old = epsilon; 0144 epsilon = epsilon-delta; 0145 0146 if epsilon <= tau; 0147 xk_1 = x_k + (epsilon_old-tau)*delx_vec; 0148 total_time= cputime-t0; 0149 break; 0150 end 0151 if length(gamma_xh) >= Te 0152 total_time = cputime-t0; 0153 break; 0154 end 0155 0156 % disp(sprintf(['iter = %d, delta = %3.4g, idelta = %d, flag = %d'], iter, delta, idelta, flag)); 0157 itr_history = [itr_history; idelta delta flag]; 0158 0159 %% Check convergence criterion 0160 if early_terminate 0161 if length(gamma_x) < M/2 0162 xhat = zeros(N,1); 0163 rhs = A(:,gamma_xh)'*y; 0164 update_mode = 'recompute'; 0165 delx_update; 0166 xhat(gamma_x) = out; 0167 if norm(y-A*xhat) < tau 0168 xk_1 = xhat; 0169 break; 0170 end 0171 end 0172 end 0173 0174 % update support 0175 update_supp; 0176 0177 pk_old([gamma_xh; idelta]) = sign(pk_old([gamma_xh; idelta]))*epsilon; 0178 0179 % update delx 0180 z_x = -sign(pk_old); 0181 rhs = z_x; 0182 update_mode = 'update'; 0183 update_delx; 0184 0185 0186 constr_violation = nnz((abs(pk_old(gamma_xc))-epsilon)>1e-10); 0187 sign_violation = nnz((sign(pk_old(gamma_xh))+sign(xk_1(gamma_xh)))>1); 0188 if constr_violation 0189 chk = gamma_xc((abs(pk_old(gamma_xc))-epsilon)>1e-10); 0190 stp = 1; 0191 fprintf('problem... with constraint violation -- %s \n', mfilename); 0192 fprintf('Refactorize the matrix... recompute delx \n'); 0193 % some times it comes here due to bad conditioning of AtAgx. 0194 update_mode = 'init0'; 0195 update_delx; 0196 end 0197 if sign_violation>=1 0198 chk = gamma_xh(sign(pk_old(gamma_xh))+sign(xk_1(gamma_xh))>1); 0199 stp = 1; 0200 fprintf('problem... sign mismatch -- %s \n',mfilename); 0201 fprintf('Refactorize the matrix... recompute delx \n'); 0202 update_mode = 'init0'; 0203 update_delx; 0204 end 0205 0206 if err_record 0207 error_table = [error_table; epsilon norm(xk_1-x_orig)/norm(x_orig) length(gamma_xh)]; 0208 % obj = [obj; 1/2*norm(A*xk_1-y)^2+epsilon*norm(xk_1,1)]; 0209 end 0210 end 0211 if err_record 0212 error_table = [error_table; epsilon norm(xk_1-x_orig)/norm(x_orig) length(gamma_xh)]; 0213 % obj = [obj; 1/2*norm(A*xk_1-y)^2+epsilon*norm(xk_1,1)]; 0214 end 0215 total_iter = iter; 0216 total_time = cputime-t0; 0217 0218 out = []; 0219 out.x_out = xk_1; 0220 out.gamma = gamma_xh; 0221 out.iter = total_iter; 0222 out.time = total_time; 0223 out.error_table = error_table; 0224 % out.obj = obj;